# Angular Resolution of the Search for Anisotropic Stochastic Gravitational-
Wave Background with Terrestrial Gravitational-Wave Detectors

Erik Floden  School of Physics and Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  Vuk Mandic  School of Physics and Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA  Andrew Matas  Max Planck
Institute for Gravitational Physics (Albert Einstein Institute), D-14476
Potsdam, Germany  Leo Tsukada  Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA

###### Abstract

We consider an anisotropic search for the stochastic gravitational-wave (GW)
background by decomposing the gravitational-wave sky into its spherical
harmonics components. Previous analyses have used the diffraction limit to
define the highest-order spherical harmonics components used in this search.
We investigate whether the angular resolution of this search is indeed
diffraction-limited by testing our ability to detect and localize simulated GW
signals. We show that while using low-order spherical harmonics modes is
optimal for initially detecting GW sources, the detected sources can be better
localized with higher-order spherical harmonics than expected based on the
diffraction limit argument. Additionally, we discuss how the ability to
recover simulated GW sources is affected by the number of detectors in the
network, the frequency range over which the search is performed, and the
method by which the covariance matrix of the GW skymap is regularized. While
we primarily consider point-source signals in this study, we briefly apply our
methodology to spatially-extended sources and discuss potential future
modifications of our analysis for such signals.

##  I Introduction

Since the first direct detection of gravitational-waves (GW) from the
collision of two black holes on September 14, 2015 [1], the field of
gravitational-wave astronomy has been an exciting and relatively new way to
probe the far reaches of our universe. The ability to detect gravitational-
waves with detectors like the Laser Interferometer Gravitational-wave
Observatory (LIGO) and Virgo detectors [2, 3] has opened up a wide range of
questions to investigate. Recent developments include the observation of
neutron star-black hole coalescences [4], a binary black hole coalescence with
a total mass of 150 $M_{\odot}$ [5], and a coalescence between a black hole
and a 2.6 $M_{\odot}$ object which is either the lightest black hole or the
heaviest neutron star discovered in a compact binary system [6]. As of the end
of LIGO-Virgo’s third observing run (O3), 90 GW event candidates have been
detected [7], and merger rates have been estimated to be between $17.3$ Gpc-3
yr-1 and $45$ Gpc-3 yr-1 for binary black holes, between $13$ Gpc-3 yr-1 and
$1900$ Gpc-3 yr-1 for binary neutron stars, and between $7.4$ Gpc-3 yr-1 and
$320$ Gpc-3 yr-1 for neutron star-black hole binaries [8].

While these detections are of signals from individual GW events, it is also
possible for a GW background to form as the superposition of many unresolved
GW signals [9, 10]. Such a background may have contributions which are
astrophysical or cosmological in origin. Examples of astrophysical
contributions include binary mergers and supernovae [11, 12, 13, 14, 15, 16,
17] while examples of cosmological contributions include GWs generated during
the inflationary epoch and in phase transitions in the early universe [18, 19,
20] . Moreover, the GW background is likely to be anisotropic. The potential
sources of anisotropy include primordial density fluctuations (e.g. reflected
in the distribution of compact binaries throughout the universe), the local
distribution of GW sources (e.g. pulsars in the Milky Way galactic plane [21,
22]), the velocity of the Solar System, and others [23, 24, 25, 26, 27, 28].

Advanced LIGO and Advanced Virgo data have been used to search for the
anisotropic stochastic gravitational-wave background (SGWB), producing
stringent upper limits on GW energy density across the sky [29, 30, 31].
Traditionally, diffraction limit arguments were used to assess the angular
resolution of these searches. In this paper, we take a closer look at the
intrinsic angular resolution of anisotropic SGWB searches, finding that they
can surpass then resolution expected based on diffraction limit arguments. In
Section II we present the formalism for anisotropic SGWB search. In Section
III we present some intuitive arguments on angular resolution limitiations. In
Sections IV, V, and VI we present the angular resolution in recovery of
simulated individual point sources, simulated multiple point sources, and
simulated extended sources, respectively. We offer concluding remarks in
Section VII.

##  II Spherical Harmonics Decomposition

We consider the anisotropic SGWB search in which we decompose the map of the
gravitational-wave sky into spherical harmonics components. We assume an
unpolarized, Gaussian, and stationary SGWB. The most general quadratic
expectation value of the GW strain $h_{A}(f,\Theta)$ of frequency $f$, sky
direction $\Theta$, and polarization $A$, that satisfies these assumptions is
given by:

$\displaystyle\langle
h^{*}_{A}(f,\Theta)h_{A^{\prime}}(f^{\prime},\Theta^{\prime})\rangle=\frac{1}{4}\mathcal{P}(f,\Theta)\delta_{AA^{\prime}}\delta(f-f^{\prime})\delta(\Theta,\Theta^{\prime})$
(1)

where ${\mathcal{P}}(f,\Theta)$ gives the spectral and angular distribution of
the background [32]. We assume that ${\mathcal{P}}(f,\Theta)$ can be factored
into its separate spectral and angular components

$\displaystyle\mathcal{P}(f,\Theta)=H(f)\,\mathcal{P}(\hat{\Theta}),$ (2)

where $H(f)$ is a dimensionless quantity which we choose to take the form of a
power law both for its simplicity and its ability to approximate most
interesting SGWB models,

$\displaystyle H(f)=\left(\frac{f}{f_{\rm ref}}\right)^{\alpha-3}.$ (3)

Here we use the reference frequency ${f_{\rm ref}}=25$ Hz. The spectral index
$\alpha$ is commonly assumed to take values of 0, 2/3, and 3 corresponding to
a cosmological background, CBC background, and a generic flat strain spectrum,
respectively. We further decompose the angular distribution
$\mathcal{P}(\Theta)$ into a basis of spherical harmonic components $Y_{\ell
m}$:

$\displaystyle\mathcal{P}(\hat{\Theta})=\sum_{\ell=0}^{\ell_{max}}\sum_{m=-\ell}^{\ell}P_{\ell
m}Y_{\ell m}(\hat{\Theta}).$ (4)

The goal of the search is therefore to estimate the values of the $P_{\ell m}$
coefficients. To do this, we must first define the cross-correlation spectrum
between two detectors $I$ and $J$ at time $t$ and frequency $f$:

$\displaystyle C_{IJ}=\frac{2}{\tau}\tilde{s}_{I}^{*}(t;f)\tilde{s}_{J}(t;f).$
(5)

Here $\tau$ is the duration of an observation segment, and we use the finite-
time Fourier transform of each detector time series data, which contains both
the detector noise and the GW signal:
$\tilde{s}_{I}(t;f)=\tilde{n}_{I}(t;f)+\tilde{h}_{I}(t;f)$. If we assume that
the noise between the detectors is uncorrelated, then the expectation value of
the cross-correlation spectrum is:

$\displaystyle\langle
C_{IJ}(t;f)\rangle=\frac{2}{\tau}\langle{h}^{*}_{I}(t;f){h}_{J}(t;f)\rangle.$
(6)

Equations 1 and 2 then lead to [32]:

$\displaystyle\langle
C_{IJ}(t;f)\rangle=H(f)\int_{S^{2}}d\hat{\Theta}\gamma(t;\hat{\Theta},f)\mathcal{P}(\hat{\Theta}).$
(7)

where $\gamma(t;\hat{\Theta},f)$ is a geometric function that takes into
account the response of the detector pair to GW signals given the detectors’
relative locations and orientations [32]:

$\displaystyle\gamma(t;\hat{\Theta},f)=\frac{1}{2}F^{A}_{I}(t;\hat{\Theta})F^{A}_{J}(t;\hat{\Theta})e^{i2\pi
f\hat{\Theta}\cdot(\vec{x}_{I}(t)-\vec{x}_{J}(t))/c}$ (8)

In the above equation we use the detector response functions
$F^{A}_{I}(t;\hat{\Theta})$, detector locations $\vec{x}_{I}(t)$ [32], and
speed of light $c$. We can substitute Eq. 4 into Eq. 7 and integrate over the
two-sphere to obtain

$\displaystyle\langle
C_{IJ}(t;f)\rangle=H(f)\gamma_{lm}(t;f)\mathcal{P}_{lm},$ (9)

where

$\displaystyle\gamma_{lm}(t;f)=\int_{S^{2}}d\hat{\Theta}\gamma(t;\hat{\Theta},f)Y_{lm}(\hat{\Theta})$
(10)

are purely geometric factors associated with the $IJ$ detector pair and
repeated indices are summed over.

We proceed to define a likelihood function for the cross-correlation spectrum
given a sky map defined by spherical coefficients $\\{\mathcal{P}_{lm}\\}$
[29]:

$\displaystyle p(C_{IJ}|\\{\mathcal{P}_{lm}\\})\propto\,{\rm
exp}\Big{(}[C_{IJ}(t;f)-H(f)\,\gamma_{lm}(t;f)\,\mathcal{P}_{lm}]^{*}\Big{.}$
$\displaystyle\Big{.}N^{-1}_{ft,f^{\prime}t^{\prime}}[C_{IJ}(t^{\prime};f^{\prime})-H(f^{\prime})\,\gamma_{l^{\prime}m^{\prime}}(t^{\prime};f^{\prime})\,\mathcal{P}_{l^{\prime}m^{\prime}}]\Big{)}$
(11)

where $N_{ft,f^{\prime}t^{\prime}}$ is the covariance matrix of $C_{IJ}(t;f)$
given by [33]:

$\displaystyle
N_{ft,f^{\prime}t^{\prime}}=\delta_{tt^{\prime}}\delta_{ff^{\prime}}P_{I}(t;f)P_{J}(t;f),$
(12)

and $P_{I}(t;f)$ is the one-sided power spectrum of the data from detector
$I$. The spherical harmonic coefficients that maximize this likelihood
function are given by:

$\displaystyle\hat{P}_{lm}=\sum_{l^{\prime}m^{\prime}}\big{(}\Gamma_{R}^{-1}\big{)}_{lm,l^{\prime}m^{\prime}}X_{l^{\prime}m^{\prime}}.$
(13)

where

$\displaystyle
X_{lm}=\sum_{f,t}\gamma_{lm}^{*}(t;f)\frac{H(f)}{P_{I}(t;f)P_{J}(t;f)}C_{IJ}(t;f),$
(14)
$\displaystyle\Gamma_{lm,l^{\prime}m^{\prime}}=\sum_{f,t}\gamma_{lm}^{*}(t;f)\frac{H^{2}(f)}{P_{I}(t;f)P_{J}(t;f)}\gamma_{l^{\prime}m^{\prime}}(t;f).$
(15)

Equation 14 defines the “dirty map” of the SGWB convolved with the response
antenna pattern of the $IJ$ detector pair. The covariance matrix of the dirty
map is given in Eq. 15. Since Eq. 13 represents the deconvolution of the GW
signal from the detector response, we refer to it as the “clean map”. We can
obtain the covariance matrix of this clean map by inverting the covariance
matrix of the dirty map. Because of this property, we refer to the matrix
defined in Eq. 15 as the Fisher information matrix [32].

Inverting the Fisher matrix, however, is nontrivial, as gaps in the detector
network’s sensitivity across the sky cause singular values within the Fisher
matrix. We therefore must regularize the Fisher matrix prior to inversion. In
Eq. 13 the subscript $R$ denotes that the inverted Fisher matrix has been
regularized. We utilize the singular value decomposition (SVD) method for this
regularization, choosing which eigenmodes to discard and which to retain prior
to inversion. In Section IV we consider several options for implementing the
SVD regularization.

We note that the spherical harmonic expansion defined in Eq. 4 is cut off at
some $\ell_{max}$. The choice of $\ell_{max}$ is an assumption of the search.
In past searches, this value was chosen based on a diffraction limit argument
in which a diffraction-limited spot size $\theta$ for a particular baseline is
defined using the distance $d$ between detectors in a baseline pair and the
most sensitive frequency $f_{\alpha}$ of the baseline (which depends on the
assumed spectral model):

$\displaystyle\ell_{\rm max}=\frac{\pi}{\theta},$ (16)

where for the LIGO-Hanford and LIGO-Livingston detectors

$\displaystyle\theta=\frac{c}{2df_{\alpha}}\approx\frac{50\rm{Hz}}{f_{\alpha}}.$
(17)

Whether or not this argument is appropriate or sufficient for making optimal
GW signal detections is a question we address in Sections IV and V.

##  III Source of Angular Resolution Limitation

The diffraction limit relation in Eq. 17 gives a reasonable order-of-magnitude
estimate for the angular resolution and captures the intuition that a larger
baseline should improve the angular resolution [34]. However, there are
important shortcomings of this formula as usually presented. First, the word
“diffraction” gives a misleading picture of the physics, since the angular
resolution limitation actually comes from phase coherence between sources at
different angular separations (as discussed below in more detail). Second,
gravitational-wave sources are typically not monochromatic, and it is not
clear a-priori how to transfer this diffraction limit to broadband sources.
Therefore a more detailed investigation is needed to pin down the angular
resolution of a given GW detector network.

An important subtlety is that for deterministic signal models, such as
transient CBC signals or isolated rotating neutron stars, the angular
resolution of the recovered signal is evidently not limited by the diffraction
limit. The CBC signal model _assumes_ (quite reasonably) that there is one
source explaining the gravitational-wave signal. This allows triangulation of
the arrival times of the signal, along with the antenna patterns, to pin down
the sky location of a source to (in principle) an arbitrary good
localization—currently the best localization is of order 10 square degrees for
the highest SNR events seen to date. On the other hand, the anisotropic SGWB
search is making an unmodeled determination of the GW sky. This generality
comes with a loss in angular resolution, but does not require an assumption of
a single point-like source, which is not appropriate for large scale
structures. In other words, a skymap produced by the CBC search is a
probability distribution for the location of _one_ source, while the
anisotropic SGWB search skymap is an image of the entire population of GW
sources across the sky.

In addition to the phase coherence, there is an additional effect due to blind
spots in the detector network. This problem is a manifestation of the well-
known _bias-variance tradeoff._ As noted above, the Fisher matrix is usually
not invertible, meaning that formally the variance is infinite. Regularizing
the Fisher matrix reduces the variance, but at the cost of introducing a bias
in the recovered signal.

Inspired by this, we briefly describe a toy model that can be used to gain
insight into the origin of angular resolution of the SGWB search. We consider
a pair of detectors rotating on the boundary of a circle of diameter $L$ and
observing scalar waves propagating in the same plane (e.g. two water buoys
measuring the height of passing water waves). The detectors rotate with
frequency $f_{E}$ such that $f_{E}\ll f$. Our detectors will have a response
function $F(\theta)$, where $\theta$ denotes the polar angle in the plane. We
consider two special cases: (a) an isotropic response $F(\theta)\propto 1$ and
(b) a directional “lighthouse” response
$F(\theta)\propto\delta(\theta-\theta_{0})$, where $\theta_{0}$ is the opening
angle of the detector.

In the circular harmonic basis, we can define our overlap reduction function
as (equivalent to Eq. 10):

$\gamma_{n}(t,f)=\frac{1}{2}\int d\theta
e^{in\theta}F_{1}(\theta,t)F_{2}(\theta,t)e^{2\pi
ifL\hat{\Theta}\cdot\hat{n}(t)/c},$ (18)

where $\hat{\theta}$ is the unit vector pointing in the direction of $\theta$
and $\hat{n}$ points from detector 2 to detector 1. In the isotropic response
case, a suitable choice of coordinate axes allows us to rewrite our overlap
reduction function as

$\gamma_{n}(t,f)=\frac{1}{2}e^{in(2\pi f_{E}t)}J_{n}(2\pi fL/c),$ (19)

where $J_{n}$ are the Bessel functions. In the lighthouse response case, this
function reduces to

$\gamma_{n}(t,f)=\frac{1}{2}e^{in2\pi f_{E}t}e^{i2\pi fL\cos(2\pi
f_{E}t)/c}\delta^{2}(\phi_{1},\phi_{2}),$ (20)

where $\phi_{1},\phi_{2}$ are the directions of detectors 1,2 at time $t=0$.
The cross correlation between the two detectors is defined as in Eq. 5, with
its expectation value defined as

$\langle C(t,f)\rangle=\sum_{n}\gamma_{n}(t,f)\mathcal{P}_{n}$ (21)

We are now ready to analyze the angular resolution of this toy model.
Following the anisotropic SGWB search formalism, we define the likelihood
function

$\ln(\mathcal{L})\propto\sum_{t}\frac{|C(t,f_{0})-\gamma_{n}(t,f_{0})\mathcal{P}_{n}|^{2}}{P_{0}^{2}},$
(22)

where for this model we will only be considering a single frequency bin and
$P_{0}$ is the noise variance in this bin. Our maximum likelihood estimator is
then given as

$\hat{\mathcal{P}}_{n}=(\Gamma^{-1})_{nm}X_{m},$ (23)

with

$X_{n}=\frac{1}{P_{0}^{2}}\sum_{t}\gamma^{*}_{n}(t,f_{0})C(t,f_{0})$ (24)
$\Gamma_{nm}=\frac{1}{P_{0}^{2}}\sum_{t}\gamma^{*}_{n}(t,f_{0})\gamma_{m}(t,f_{0}).$
(25)

We first consider the case of an isotropic response. The Fisher matrix
$\Gamma_{mn}$, after summing over $t$, reduces to

$\Gamma_{mn}=\frac{1}{4P_{0}^{2}}J_{n}^{2}(2\pi f_{0}L/c)N_{S}\delta_{mn},$
(26)

where $N_{S}$ is the number of time segments. Since the Fisher matrix is
diagonal, we see that all eigenvalues are nonzero for general $2\pi f_{0}L/c$.
However, $J_{n}(2\pi f_{0}L/c)$ peaks around $2\pi f_{0}L/c=n_{0}$, so wave
numbers $n>2\pi f_{0}L/c$ are suppressed. The effect of this eigenvalue
suppression is quantitatively similar to the diffraction limit argument (c.f.
Eq. 17), and it applies even for arbitrarily long observation times (and, by
extension, high SNRs) since all eigenvalues scale equally with these
quantities. We therefore see that for an isotropic response, a detector pair’s
resolution is limited to angular scales similar to those given by the
diffraction limit argument, albeit for a very different physical reason.

Next, we consider the lighthouse detector response function and apply the same
analysis. We again obtain a diagonal Fisher matrix that reduces to

$\Gamma_{mn}=\frac{1}{4P_{0}^{2}}\delta_{mn}\delta^{2}(\phi_{1},\phi_{2}),$
(27)

for which all eigenvalues are equal implying that there is no limit to the
angular resolution, i.e. the resolution can be arbitrarily improved by having
a larger SNR signal. We conclude therefore that detector networks with
perfectly isotropic responses have angular resolutions similar to those set by
the diffraction limit argument, while perfectly directional detectors have
unlimited angular resolution. In practice, however, the LIGO-Virgo detectors
have anisotropic responses but not perfectly so. That is, they have responses
somewhere between the two extreme cases considered above. Consequently, we
should expect the anisotropic SGWB search angular resolution to be limited,
but also better than what is predicted by the diffraction limit argument.

Unfortunately, the above analysis cannot be simply extended to GW detectors
operating in three dimensions and across a wide frequency band with a colored
noise power spectrum. Consequently, to assess the angular resolution of the
anisotropic SGWB search, we resort to a series of simulations described in the
following Sections.

##  IV Single Point Source

As noted above, the diffraction-limit argument in Eq. 17 may give a reasonable
order-of-magnitude estimate of the angular resolution for the anisotropic SGWB
search. However, this argument is inaccurate for multiple reasons. First, the
angular resolution of this search is not limited by a diffraction process of
GWs interacting with the detectors but rather by the phase coherence between
detectors. As discussed in the previous Section, the directionality of the
detector response plays an important role in the angular resolution, which is
not accounted for in the diffraction limit argument. Second, when considering
a network of multiple baselines, each with different separation distances and
peak sensitivities, Eq. 17 cannot be directly applied. Third, the question of
regularization of the Fisher matrix (which is also related to the number of
detector baselines in the network and to their relative sensitivities) also
directly impacts the angular resolution of the search.

To assess the angular resolution limitations, we resort to a series of
simulations and we cast the problem in terms of finding the search parameters
that give optimal results when performing a spherical harmonics decomposition
(SHD) SGWB search. These parameters include the frequency range over which the
search is performed, the highest-order spherical harmonic modes that are
considered, and the method/threshold used for Fisher matrix regularization.
These parameters will differ depending on the signal model, and we will
address what constitutes optimal results.

We begin by considering a broadband point-source GW simulation. While the SHD
analysis is particularly well-suited for the recovery of spatially-extended
sources, a point-source recovery is still possible with the SHD analysis, and
the angular resolution of the recovery can be quantified in terms of
$\ell_{max}$. In this analysis, we use the same data processing pipeline used
to perform the O3 anisotropic search [29] and outlined in Section II. The
point-source simulation assumes a power-law spectral shape as defined in Eq.
28:

$AH(f)=A\left(\frac{f}{f_{\rm ref}}\right)^{\alpha-3},$ (28)

where we choose $\alpha=2/3$ and vary the amplitude $A$. Fig. 1 illustrates
the strain power spectrum of this simulation for $A=5\times 10^{-49}$,
alongside the power spectral density (PSD) of Advanced LIGO’s design
sensitivity. The single point-source was simulated at $0^{\degree}$
declination and $12^{h}$ right ascension. Each simulation was 24 hours long,
and for each detector the SGWB was added to the detector noise corresponding
to the Advanced LIGO or Advanced Virgo design sensitivities (which are not
taken to be the same). Recoveries were performed using a frequency range of
20-500 Hz and spectral index of $\alpha=2/3$.

Figure 1: Power spectral density (PSD) of Advanced LIGO design sensitivity
(blue, upper) is shown in comparison with the strain power spectrum of the
simulated point source with $\alpha=2/3$ (orange, lower).

To quantify the localization of the recovered point source, we count the
number of pixels in the recovered SNR skymap whose SNR is larger than 75% of
the map’s peak SNR value. The spot size defined in this way is shown in Fig. 2
(left) as a function of $\ell_{max}$. As $\ell_{max}$ is increased, the spot
size is reduced and the resolution of the recovered map is improved. For
comparison, for $\alpha=2/3$ the most sensitive frequency is about 65 Hz,
resulting in the diffraction limit in Eq. 17 of $\theta\approx 45^{\circ}$ and
the diffraction limit spot size of about 4% of the sky. We can also consider
the fact that we make statistically significant recoveries of the simulated
source at an $\ell_{max}$ value of 20, and the diffraction-limited spot size
corresponding to this $\ell_{max}$ value is $\theta\approx 9^{\circ}$, which
represents about 0.2% of the sky. It is evident from Fig. 2 that the
anisotropic SGWB search can substantially surpass the diffraction limit, as
the smallest spot size measurements we obtain are about 0.2% of the sky map.

Figure 2 (left) also compares different Fisher matrix regularization methods.
The first regularization threshold we apply is to keep two thirds of all
eigenvalues of the Fisher matrix ($f_{keep}=2/3$). This is the method that has
been in use for all previous SHD searches for anisotropic SGWB [29, 30, 31].
The second regularization threshold we apply removes the eigenvalues that are
smaller that $10^{-3}$ times the largest eigenvalue 111This is often stated as
the condition number of an eigenvalue being smaller than $10^{-3}$. A
comparison between these two methods is shown in Fig. 2 (left)—while the two
regularization methods yield similar results, the $10^{-3}$ threshold performs
slightly better. For the $10^{-3}$ threshold, a high $\ell_{max}$ makes little
difference to the spot size since most new eigenvalues introduced by raising
$\ell_{max}$ are small and therefore removed by the regularization threshold.
This behavior can be seen in the relatively constant spot size measurements
for high $\ell_{max}$ values in Fig. 2 (left). We also include spot size
measurements that result from using a more aggressive regularization scheme in
which we apply a $10^{-2}$ threshold. This threshold results in more
eigenmodes being discarded during the Fisher matrix regularization, giving
larger spot size measurements.

Figure 2: Size of recovered point-source signal as a function of $\ell_{max}$
for multiple regularization thresholds for the Fisher matrix. The left plot
shows results for the Hanford-Livingston baseline, and the right plot shows
results from the HLV network. For the HLV case, we also consider the
unregularized Fisher matrix (keep all). Data points are excluded if the point-
source was not recovered.

We also consider a modified version of this search in which we include the
Virgo detector, resulting in a network of three detector baselines. We denote
this network as HLV. A single baseline will not be sensitive to all sky
directions, and these gaps in sensitivity manifest as singular Fisher matrix
eigenvalues. When more detector baselines are included in a network, it is
possible for the Fisher matrix to become naturally regularized as the network
gains sensitivity to more sky directions. A different Fisher matrix
regularization strategy may therefore be better suited for multi-baseline
detector networks. The size of the recovered point-source with this multi-
baseline network is shown in Fig. 2 (right). In addition to the above two
regularization thresholds, we also consider using the entire (unregularized)
Fisher matrix, i.e. keeping all of its eigenvalues.

Figure 3: Fisher matrix conditioning plots for $\ell_{max}=4,8,$ and $12$
comparing the HL and HLV detector networks using simulated data. Vertical
lines mark the largest 2/3 of eigenvalues. Horizontal lines mark the $10^{-3}$
threshold.

For $\ell_{max}<12$, using the $10^{-3}$ threshold is equivalent to keeping
all eigenvalues. The $10^{-3}$ threshold generally gives the most localized
recoveries of the three regularization methods used. For $\ell_{max}\geq 12$,
the $10^{-3}$ threshold performs similarly to the $f_{keep}=2/3$ method. For
$\ell_{max}>13$, the recovered peak SNR values are not statistically
significant when we keep all eigenvalues, so we conclude that the source is
not properly recovered in this case, and we do not report spot size values. We
can conclude therefore that some regularization is still necessary for this
multi-baseline network for $\ell_{max}>13$. If we decrease the amplitude of
our simulated source from $A=5\times 10^{-49}$ to $A=1\times 10^{-49}$, we
only make statistically significant recoveries of the source for
$\ell_{max}<6$, and when we do recover this lower-amplitude source, its spot
size measurements are larger than those of the higher-amplitude source as seen
in Fig. 2 (right). This is consistent with our expectation that larger SNR
allows recovery with better angular resolution.

Figure 3 directly compares Fisher matrix conditioning in the HL and HLV cases.
We see that for a low $\ell_{max}$ value like $\ell_{max}=4$, there is little
difference in the Fisher matrix eigenvalues between the HL and HLV networks if
we only consider eigenvalues kept using the $f_{keep}=2/3$ threshold. If we
consider all the eigenvalues, the smallest eigenvalue in the HLV Fisher matrix
is an order of magnitude larger than that of the HL Fisher matrix. Using
higher values of $\ell_{max}$ shows a widening in the differences between the
HL and HLV Fisher matrix conditioning, with the HLV network having larger
eigenvalues even among those that are accepted by the $f_{keep}=2/3$
regularization method. This improvement in conditioning is likely responsible
for the slightly smaller spot size measurements for HLV (as compared to HL)
shown in Fig. 2.

We next examine the impact of the chosen frequency band on the conditioning of
the Fisher matrix. For the spectral index $\alpha=2/3$, the most sensitive
frequency is $\sim 65$ Hz [36]. Further, 99% of the HL baseline’s sensitivity
lies within the 20-120 Hz band for the search for an isotropic SGWB with
$\alpha=2/3$ [13]. We therefore start with the 20-100 Hz frequency band that
captures most of the sensitivity of the SHD search. We compare the
conditioning of the Fisher matrix in this band with the conditioning in wider
bands (20-200 Hz and 20-500 Hz)—these wider bands include higher frequencies
and therefore may be sensitive to smaller angular scales on the sky.

Figure 4: Fisher matrix conditioning is shown for $\ell_{max}=4,8,$ and $12$,
comparing several different frequency bands for the HL baseline. Vertical
lines mark the lowest 1/3 of eigenvalues and horizontal lines mark the
$10^{-3}$ threshold.

Figure 4 shows that as the $\ell_{max}$ value increases, the regularization of
the Fisher matrix worsens drastically for the 20-100 Hz band. Extending the
frequency band to 20-200 Hz gives a marked improvement, but the conditioning
is the best when using the widest frequency band of 20-500 Hz. This difference
in conditioning implies that including higher frequencies allows the detector
network to resolve smaller-scale structures in the skymap, thereby increasing
the angular resolution of the recovery.

Finally, we examine our ability to detect a point source rather than our
ability to localize it. We quantify the significance of a GW source recovery
in a manner similar to that used in previous directional searches [29]. We
draw from a multivariate Gaussian distribution with zero mean and covariance
given by the inverted Fisher matrix, resulting in a set of $P_{\ell m}$ values
corresponding to a noise skymap. The procedure is repeated 5000 times and an
SNR map is made for each such random realization of noise. In each SNR map we
define a region centered around the location of the simulated source,
extending $20^{\degree}$ above and below the simulated point source and
$2^{h}$ on either side. We then find the mean pixel value among the top 5% of
the pixels within this region and bin this value into a histogram. We make the
same measurement using a skymap containing the simulated point source signal,
like the one shown in Fig. 5 for $\ell_{max}=8$. We attribute a (significance)
p-value to the peak SNR in this map by comparing it to the noise peak-SNR
histogram. An example is shown in Fig. 6, where the vertical line denotes the
peak SNR measurement of the map in Fig. 5.

Figure 5: An example SNR map is shown of a simulated point-source with
amplitude $A=10^{-49}$ and location $0^{\degree}$ declination, $12^{h}$ right
ascension, recovered using the HL baseline, $f_{keep}=2/3$, and
$\ell_{max}=8$. Figure 6: Peak SNR from the map in Fig. 5 is marked by a
vertical line in comparison with the histogram of the peak SNRs obtained from
purely noise simulated skymaps.

We find that increasing $\ell_{max}$ causes an increase in the p-value
associated with its peak SNR and therefore a decrease in the ability to detect
the simulated source. These results are shown in Fig. 7. For each additional
order of spherical harmonics modes we include in our analysis by increasing
$\ell_{max}$, we increase the number of model parameters used to fit to the
data, which ultimately results in a decrease in the SNR of the recovered
source. We therefore find that low $\ell_{max}$ values lend themselves to
optimal point source signal detection, while high $\ell_{max}$ values may be
better suited for optimal localization of the point source signal, depending
on the signal strength. In other words, search parameters that are optimal for
detecting a GW source may differ from those that are optimal for localizing
it.

Figure 7: We show the significance (p-value) of the recovered simulated peak
SNR vs $\ell_{max}$ using HL baseline, $f_{keep}=2/3$, and point-source
amplitude $A=10^{-49}$. p-values equal to zero are denoted with downward-
pointing arrows at the bottom of the plot. We also show results for the HLV
network, for two different regularization methods. We see that using higher
$\ell_{max}$ gives recovered peak SNRs with lower significance (larger
p-values).

##  V Multiple Point Sources

Another way to characterize the angular resolution of the SHD search is by its
ability to distinguish between multiple point sources. We begin by simulating
two point-sources of equal magnitude with declinations of $+30^{\degree}$ and
$-30^{\degree}$, and we measure the spot size of each of the two sources by
the same process as in Fig. 2, but applied separately to the upper and lower
halves of the map. We note that the choice to look at the two halves of each
map separately requires knowledge of the sources’ actual locations, which
would not be available in real GW searches. However, the goal of this test it
to investigate the ability to resolve separated point-sources, and not to
provide the optimal methodology for identifying multiple point sources in the
skymap. Amplitudes of simulated point sources in this test were louder than
those used previously in this paper. This choice was necessary so as to
investigate whether neighboring GW signals with sufficiently high SNRs are
better recovered when we use sufficiently high-order spherical harmonics
modes.

Figure 8: Sizes of two recovered point-source signals as a function of
$l_{max}$ are shown for multiple source amplitudes. The upper and lower point
sources are denoted with upward- and downward-pointing arrows, respectively.
Here we use the HL baseline. Error bars are given for one set of points to
reduce clutter. Other data points have comparable errors.  Figure 9: SNR maps
of two recovered point-source signals using the HL baseline are shown for
point source amplitudes of $10^{-48}$ and for $l_{max}=4,8,12$ (top to
bottom). In the bottom plot we also show the ratio of the angular separation
to the total spot size of the upper and lower sources as we increase
$l_{max}$.

Figure 8 shows each source’s spot size as a function of $\ell_{max}$. We
observe behavior similar to those of the single point-source recoveries,
namely that the spot sizes of the two sources do in fact decrease as we use
higher-order spherical harmonics modes. Also, we see that the upper and lower
spot size measurements are very similar in most cases. However, this does not
necessarily mean that we are properly resolving the two signals.

Figure 9 shows the SNR maps corresponding to dual point-source simulations of
amplitude $A=10^{-48}$. We see that for low $\ell_{max}$, the two sources are
not resolved. Only at higher $\ell_{max}$ do we see distinct recoveries. The
bottom plot of Fig. 9 shows that the resolution of the two point sources
improves with increasing $\ell_{max}$.

Since we used loud simulations, the sources are still recoverable at high
$\ell_{max}$ despite the decreasing SNR behavior discussed in the previous
Section. If we were to use sources with lower SNRs, they would only be
statistically distinct from noise if we use low spherical harmonics modes,
which in turn reduces the angular resolution of the map, preventing the
sources from being resolved. Therefore we find that in order to recover
multiple GW point sources, they must be sufficiently loud so as to support
analysis with sufficiently high-order spherical harmonics components. This
result also reinforces the notion that low $\ell_{max}$ should be considered
when initially performing a search in order to detect the presence of a GW
source, before attempting analysis with higher order modes to optimize
localization.

We note that the recovered spot sizes shown in Fig. 8 are even smaller (by
about $2\times$ in area) than those shown for single point-source simulations
in Fig. 2. This is a consequence of the higher amplitude, $A=10^{-48}$, used
in the simulations shown in Fig. 8. In fact this amplitude has already been
excluded by past SGWB searches with Advanced LIGO and Advanced Virgo data.

##  VI Extended Sources

We extend the methodology presented in previous Sections to recovery of
spatially extended sources. Extended sources come in many varieties, and we do
not attempt an exhaustive study here. Instead, we chose to study an extended
source in the shape of a two-dimensional Gaussian distribution centered at
$0^{\degree}$ declination and $12^{h}$ right ascension with amplitude $5\times
10^{-49}$ and standard deviations of $30^{\degree}$ in each direction. We
first create a clean map of this signal in the HEALPix basis [37]. We then
perform the SHD search on noise-only data to produce a Fisher matrix and dirty
map corresponding to only noise. Multiplying this noise-only Fisher matrix by
the clean signal map gives the dirty map corresponding to the extended source
(c.f. Eq. 13). This dirty map is added to the noise-only dirty map to give a
final dirty map containing both noise and the simulated source. Fig. 10 shows
the clean signal map, the dirty signal map, and the final combined dirty map.

Figure 10: From top to bottom, clean map of 2-D Gaussian simulation, its
corresponding dirty map, and the sum of this dirty map with the noise-only
dirty map.

The final combined dirty map and the Fisher matrix are then used to perform
the same analysis described in Sec. IV. Namely, we obtain spot-size and
significance measurements using the same methodology as in the point-source
case, using the HL baseline and $f_{keep}=2/3$. Figure 11 shows the
significance of the recovered signal as a function of $\ell_{max}$. As in the
point-source case, we see the significance is the largest (p-values are the
smallest) for low $\ell_{max}$ values. For this particular simulation,
$\ell_{max}=8$ seems to be the cutoff beyond which our recoveries are not
statistically significant.

Figure 11: Significance of the recovered extended source as a function of
$\ell_{max}$, using the HL baseline and Gaussian source of amplitude of
$5\times 10^{-49}$.

In Fig. 12 we show the spot size as a function of $\ell_{max}$. As in the
point-source case, we observe the spot size decreasing as $\ell_{max}$
increases. The smallest size is reached for $\ell_{max}=7$, at which point the
recovered spot makes up about 1% of the map, roughly consistent with the size
of the simulated extended source. Attempting recovery for $\ell_{max}>7$ is
not successful since the SNR of the simulation is not large enough to support
analyses with larger numbers of free parameters. This is illustrated in Fig.
13, where e.g. for $\ell_{max}=10$ we see no discernible recovery of the
extended source.

Figure 12: Size of the recovered extended source as a function of
$\ell_{max}$, using the HL baseline and Gaussian source of amplitude of
$5\times 10^{-49}$. Figure 13: SNR maps of the recovered extended source using
the HL baseline and Gaussian source of amplitude of $5\times 10^{-49}$. From
top to bottom, $\ell_{max}=4,6,8,10$.

Overall, we observe that recoveries of the simple spatially extended source
follow the same patterns as recoveries of point sources. Specifically, low
values of $\ell_{max}$ are appropriate for the first detection of an
anisotropic SGWB, as they offer the largest statistical significance of the
recovery. Larger values of $\ell_{max}$ should then be used in follow-up
analyses to better localize the source, to the extent supported by the
amplitude (SNR) of the source. Future work is needed for a more exhaustive
study of spatially extended sources of different shapes and morphologies.

##  VII Conclusions

Observations of GWs from individual compact binary coalescences have raised
the possibility of observing the SGWB in the upcoming observation runs of
terrestrial GW detectors. Observation of SGWB anisotropy promises to be an
important source of information about the population of systems that give rise
to the SGWB. Consequently, understanding the angular resolution of the
anisotropic SGWB search and how it changes as a function of search parameters
will be increasingly important in future analyses of GW data.

In this paper we have investigated the impact of several search parameters on
the angular resolution of the anistropic SGWB search. We have found that the
use of a wide frequency band, such as 20-500 Hz, and of a detector network
containing multiple baselines both help increase the number of directions (or
spherical harmonic modes) to which our network is sensitive. Consequently,
this improves the conditioning of the Fisher matrix, simplifies its inversion,
and overall improves the angular resolution of the search.

We have also demonstrated that for detecting the presence of an anisotropic
SGWB, one should search using low-order spherical harmonics modes. Once the
detection is made, follow-up analyses using higher-order spherical harmonic
modes should be pursued to improve the angular resolution of the analysis and
better localize the source. Since searching with increasingly higher-order
modes reduces the significance of the recovery, the maximum number of modes
(i.e. $\ell_{max}$) will be ultimately set by the amplitude (or SNR) of the GW
source. In other words, the optimal localization of a GW source is achieved by
using the highest number of spherical harmonics modes that still give
statistically significant recoveries. Further, for a sufficiently strong SGWB
it is possible to surpass the traditional “diffraction limit” on the angular
resolution of the anisotropic SGWB search.

Similarly, we have shown that multiple point sources can be recovered and
resolved by the anisotropic SGWB search. If the sources’ amplitudes are
sufficiently large, the search could be conducted with a sufficiently large
$\ell_{max}$ so as to resolve the point sources, while still providing
statistically significant recovery. Spatially extended sources follow a
similar pattern, as we have demonstrated on a simple case of a Gaussian-
distributed extended source.

Future studies should consider more complete sampling of shapes and
morphologies of spatially extended sources. We also note that the results
presented in this paper could be altered when a different spectral model is
used (e.g. for spectral index $\alpha$ different from $2/3$), as well as when
using larger GW detector networks that include more than the three detectors
considered here.

##  VIII Acknowledgements

The work of EF and VM was in part supported by NSF grants PHY-2110238 and
PHY-1806630. The authors are grateful for computational resources provided by
the LIGO Laboratory and supported by NSF Grants PHY-0757058 and PHY-0823459.
This material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation.

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. . 116, 061102 (2016).
  * Aasi  _et al._ [2015] J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, and et al., CQGra 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, and et al., CQGra 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (KAGRA Collabration, LIGO Scientific Collaboration, and Virgo Collaboration), (2021a), arXiv:2106.15163 [astro-ph.HE] .
  * Abbott  _et al._ [2020a] R. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), (2020a), arXiv:2009.01075 [gr-qc] .
  * Abbott  _et al._ [2020b] R. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), (2020b), arXiv:2006.12611 [astro-ph.HE] .
  * Collaboration  _et al._ [2021a] T. L. S. Collaboration, the Virgo Collaboration, and the KAGRA Collaboration, (2021a), arXiv:2111.03606 [gr-qc] .
  * Collaboration  _et al._ [2021b] T. L. S. Collaboration, the Virgo Collaboration, and the KAGRA Collaboration, (2021b), arXiv:2111.03634 [astro-ph.HE] .
  * Maggiore [2000] M. Maggiore, Phys. Rep. 331, 283 (2000).
  * Christensen [2019] N. Christensen, Rept. Prog. Phys. 82, 016903 (2019), arXiv:1811.08797 [gr-qc] .
  * Regimbau and Mandic [2008] T. Regimbau and V. Mandic, Class. Quantum Gravity 25, 184018 (2008).
  * Marassi  _et al._ [2009] S. Marassi, R. Schneider, and V. Ferrari, Mon. Not. R. Ast. Soc. 398, 293 (2009).
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (KAGRA Collabration, LIGO Scientific Collaboration, and Virgo Collaboration), (2021b), arXiv:2101.12130 [gr-qc] .
  * Abbott  _et al._ [2019a] B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, V. Adya, C. Affeldt, M. Agathos, and et al., Physical Review D 100 (2019a), 10.1103/physrevd.100.061101.
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 131102 (2016).
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Physical Review Letters 120, 091101 (2018), arXiv:1710.05837 [gr-qc] .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Physical Review Letters 118, 121101 (2017a), arXiv:1612.02029 [gr-qc] .
  * Romero  _et al._ [2021] A. Romero, K. Martinovic, T. A. Callister, H.-K. Guo, M. Martínez, M. Sakellariadou, F.-W. Yang, and Y. Zhao, Phys. Rev. Lett. 126, 151301 (2021).
  * Geller  _et al._ [2018] M. Geller, A. Hook, R. Sundrum, and Y. Tsai, Phys. Rev. Lett. 121, 201303 (2018), arXiv:1803.10780 [hep-ph] .
  * Von Harling  _et al._ [2020] B. Von Harling, A. Pomarol, O. Pujolàs, and F. Rompineve, JHEP 04, 195 (2020), arXiv:1912.07587 [hep-ph] .
  * Talukder  _et al._ [2014] D. Talukder, E. Thrane, S. Bose, and T. Regimbau, Physical Review D 89, 123008 (2014).
  * Mazumder  _et al._ [2014] N. Mazumder, S. Mitra, and S. Dhurandhar, Physical Review D 89, 084076 (2014).
  * Jenkins and Sakellariadou [2018] A. C. Jenkins and M. Sakellariadou, Phys. Rev. D98, 063509 (2018), arXiv:1802.06046 [astro-ph.CO] .
  * Jenkins  _et al._ [2018] A. C. Jenkins, M. Sakellariadou, T. Regimbau, and E. Slezak, Phys. Rev. D98, 063501 (2018), arXiv:1806.01718 [astro-ph.CO] .
  * Jenkins  _et al._ [2019] A. C. Jenkins, O. Richard, M. Sakellariadou, and D. Wysocki, Physical Review Letters 122, 111101 (2019).
  * Cusin  _et al._ [2018a] G. Cusin, C. Pitrou, and J.-P. Uzan, Phys. Rev. D97, 123527 (2018a), arXiv:1711.11345 [astro-ph.CO] .
  * Cusin  _et al._ [2018b] G. Cusin, I. Dvorkin, C. Pitrou, and J.-P. Uzan, Phys. Rev. Lett. 120, 231101 (2018b), arXiv:1803.03236 [astro-ph.CO] .
  * Cusin  _et al._ [2019] G. Cusin, I. Dvorkin, C. Pitrou, and J.-P. Uzan, Phys. Rev. D 100, 063004 (2019), arXiv:1904.07797 [astro-ph.CO] .
  * Abbott  _et al._ [2021c] R. Abbott _et al._ , Phys. Rev. D 104, 022005 (2021c).
  * Abbott  _et al._ [2019b] R. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), (2019b), arXiv:1903.08844 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Physical Review Letters 118, 121102 (2017b), arXiv:1612.02030 [gr-qc] .
  * Thrane  _et al._ [2009] E. Thrane, S. Ballmer, J. D. Romano, S. Mitra, D. Talukder, S. Bose, and V. Mandic, Phys. Rev. D 80, 122002 (2009).
  * Romano and Cornish [2017] J. D. Romano and N. J. Cornish, Living reviews in relativity 20, 2 (2017).
  * Mukherjee and Silk [2019] S. Mukherjee and J. Silk, Monthly Notices of the Royal Astronomical Society (2019), 10.1093/mnras/stz3226.
  * Note [1] This is often stated as the condition number of an eigenvalue being smaller than $10^{-3}$.
  * Abbott  _et al._ [2017] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, _et al._ , Physical review letters 118, 121102 (2017).
  * Gorski  _et al._ [2005] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelman, Astrophys. J. 622, 759 (2005), arXiv:astro-ph/0409513 [astro-ph] .


# Staircase palindromic polynomials

Rabi K.C.111Mathematics, Ohio University, Athens, OH USA 222Corresponding
author<EMAIL_ADDRESS>, Abdalnaser Algoud111Mathematics, Ohio University,
Athens, OH USA

###### Abstract

We study a class of monic-palindromic polynomials that we call staircase
palindromic polynomials. Specifically, suppose $S(x,n,h)$ is a polynomial of
degree $n$ with the special form:

$\begin{split}S(x,n,h)&=x^{n}+2x^{n-1}+3x^{n-2}+\ldots+(h-1)x^{n-h+2}\\\
&\hskip 28.45274pt+hx^{n-h+1}+\ldots+hx^{h-1}\\\ &\hskip
56.9055pt+(h-1)x^{h-2}+\ldots+2x+1.\end{split}$

Then $S(x,n,h)$ can be factored as a product of cyclotomic polynomials.
Moreover, for any given n, there are $\lceil{\frac{n+1}{2}\rceil}$ staircase
polynomials, all of whose factors can be derived using two parameter $n$ and
$h$ with the help of cyclotomic polynomials. After that we explore some
classes of polynomials that can be converted to staircase polynomials.

##  1 Introduction

In this paper, we will study a class of monic-palindromic polynomials that we
call “staircase palindromic polynomials”. For simplicity, we will refer to a
staircase palindromic polynomial as a SP polynomial. Finding the roots of any
given polynomial with higher degree is generally a difficult task. However, in
this paper we will present a formula that gives all the irreducible factors of
a given SP polynomial. Moreover, zeros of these factors are roots of unity.
Because of the symmetric property of coefficients of SP polynomials, we
observe that these polynomials can be written as the product of cyclotomic
polynomials and hence all the roots of the SP polynomial are zeros of
cyclotomic polynomials. In addition to SP polynomials, we will present some
special types of polynomials that behave like SP polynomials.

Interestingly, SP polynomials arise in the study of a certain class of
periodic solutions found in a population model of the cell cycle [1] (see [2]
for background information). They appear as part of the characteristic
polynomial of the linearization about the periodic solutions. A polynomial is
said to be Schur stable if all its roots lie in the open unit disk. A periodic
orbit is asymptotically stable if the corresponding characteristic polynomial
is Schur stable. So, the analysis of SP polynomials plays an important role in
determining the stability of the periodic solution of the model.

##  2 Basic results and notations

First we will go over some notation, few standard definitions and results.

###  2.1 Roots of unity

A $n^{th}$ root of unity is a number $x$ satisfying the equation $x^{n}-1=0$.
For given $n$, the distinct roots of unity for $x^{n}-1$ can be written as

$1,\zeta,\zeta^{2},\ldots,\zeta^{n-1},$

where $\zeta=\cos(2\pi/n)+i\sin(2\pi/n)=e^{\frac{2\pi\cdot 1}{n}}$. These
roots are equally spaced on the unit circle and are the vertices of a regular
polygon, called the unit n-gon.

Next we consider the cyclotomic polynomials, which are closely associated with
roots of unity.

######  Definition 2.1.

The $n^{th}$ cyclotomic polynomial $\Psi_{n}(x)$ for any positive integer $n$,
is the unique irreducible polynomial with integer coefficients that is a
divisor of $x^{n}-1$ and is not a divisor of $x^{k}-1$ for any $k<n$.

Note: For $n>2$, the cyclotomic polynomials are monic-polynomials with integer
coefficients that are irreducible over the field of the rational numbers.

######  Lemma 2.2.

If $k$ is prime, then the cyclotomic polynomial of $x^{k}-1$ is

$\Psi_{k}(x)=x^{k-1}+x^{k-2}+\ldots+x+1.$

######  Theorem 2.3.

Let $n$ be a positive integer, then

$x^{n}-1=\displaystyle\prod_{d|n}\Psi_{d}(x).$

###  2.2 Some elementary definitions

Consider a polynomial

$P(x)=a_{n}x^{n}+a_{n-1}x^{n-1}+\cdots+a_{1}x+a_{0}$

of degree $n$, where $n\geq 1~{}\textnormal{and}~{}a_{i}\in\mathbb{R}$.

######  Definition 2.4.

The polynomial $P(x)$ is called monic-polynomial if $a_{n}=1.$

######  Definition 2.5.

The polynomial $P(x)$ is called all all-one polynomial if $a_{i}=1$ for all
$i$.

######  Definition 2.6.

The polynomial $P(x)$ is called palindromic if $a_{i}=a_{n-i}$ for every
$i=~{}0,1,\ldots,n.$

##  3 Staircase palindromic polynomials

######  Definition 3.1.

For given $n$, consider monic-palindrome polynomials

$P(x)=a_{n}x^{n}+a_{n-1}x^{n-1}+\cdots+a_{1}x+a_{0}$

of degree $n$ and let $h$ be a positive integer satisfying $1\leq
h\leq\lceil{\frac{n+1}{2}\rceil}$. For fixed $h$, if the coefficients of
$P(x)$ satisfy the following conditions:

$\displaystyle a_{\ell}$
$\displaystyle=a_{n-\ell}=\ell+1~{}\textnormal{for}~{}\ell=0,1,\ldots
h-1~{}\textnormal{and}$ $\displaystyle a_{h-1}$
$\displaystyle=a_{h}=\ldots=a_{n-h+1}=h.$

Then the class of such polynomials $P(x)$of degree $n$ will be called a
staircase palindromic polynomials.

For fixed height $h$ an staircase palindromic polynomial of degree $n$ will be
called SP polynomial and is denoted by $S(x,n,h)$.

Note: For any positive integer $n$, there are $\lceil\frac{n+1}{2}\rceil$ SP
polynomials of degree n.

######  Example 1.

For $n=7$, the possible values of the height $h$ are $1,2,3$ and $4$. The
following polynomials are all SP polynomials of degree $7$:

$\displaystyle S(x,1,7)$
$\displaystyle=x^{7}+x^{6}+x^{5}+x^{4}+x^{3}+x^{2}+x+1.$ $\displaystyle
S(x,2,7)$ $\displaystyle=x^{7}+2x^{6}+2x^{5}+2x^{4}+2x^{3}+2x^{2}+2x+1.$
$\displaystyle S(x,3,7)$
$\displaystyle=x^{7}+2x^{6}+3x^{5}+3x^{4}+3x^{3}+3x^{2}+2x+1.$ $\displaystyle
S(x,4,7)$ $\displaystyle=x^{7}+2x^{6}+3x^{5}+4x^{4}+4x^{3}+3x^{2}+2x+1.$

That is, the staircase palindromic polynomials of degree $7$ is given by
$\\{S(x,1,7),S(x,2,7),S(x,3,7),S(x,4,7)\\}$.

######  Lemma 3.2.

Consider an SP polynomial

$S(x,n,h)=\displaystyle\sum_{i=1}^{h-1}ix^{i-1}+h\sum_{j=h-1}^{n+1-h}x^{j}+\sum_{i=1}^{h}(h-i)x^{i+n}.$

Then there exist $\ell,m$ such that $S(x,n,h)$ is equal to the product of two
all-one polynomials, i.e.,

$\displaystyle S(x,n,h)=\sum_{i=1}^{l}x^{i}\cdot\sum_{j=1}^{m}x^{j}.$

###### Proof.

Consider the product

$\displaystyle\displaystyle\sum_{i=1}^{\ell}x^{i}\cdot\sum_{j=1}^{m}x^{j}$
$\displaystyle=(1+x+x^{2}+\ldots+x^{\ell})\cdot(1+x+x^{2}+\ldots+x^{m})$

$=\begin{array}[]{c@{}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}ccl}&\makebox[0.0pt]{\raisebox{-6.0pt}[0.0pt][0.0pt]{}}1+&x+&x^{2}+&\ldots+&x^{\ell-1}+&x^{\ell}&{}\hfil&{}\hfil&{}\hfil&{}\hfil&{}\hfil&\\\
&{}\hfil&x+&x^{2}+&\ldots+&x^{\ell-1}+&x^{\ell}+&x^{\ell+1}&{}\hfil&{}\hfil&{}\hfil&\\\
&{}\hfil&{}\hfil&x^{2}+&\ldots+&x^{\ell-1}+&x^{\ell}+&x^{\ell+1}+&x^{\ell+2}&{}\hfil&{}\hfil&\\\
&&&&&&\vdots&&&&&\\\
&{}\hfil&{}\hfil&{}\hfil&{}\hfil&x^{m-1}+&x^{m}+&x^{m+1}+&\ldots&+x^{\ell+m-2}&+x^{\ell+m-1}&\\\
&{}\hfil&{}\hfil&{}\hfil&{}\hfil&{}\hfil&x^{m}+&x^{m+1}+&\ldots&+x^{\ell+m-2}&+x^{\ell+m-1}&+x^{\ell+m}&\\\
\cline{1-12}\cr\\\\[-13.07996pt] \end{array}$
$\displaystyle=\displaystyle\sum_{i=1}^{\ell}ix^{i-1}+(\ell+1)\sum_{j=\ell}^{m}x^{j}+\sum_{i=1}^{\ell}(\ell+1-i)x^{m+i}$
$\displaystyle=\displaystyle\sum_{i=1}^{h-1}ix^{i-1}+h\sum_{j=h-1}^{n+1-h}x^{j}+\sum_{i=1}^{h-1}(h-i)x^{(n+1-h)+i},$

where $\ell=h-1$ and $m+\ell=n$. ∎

######  Theorem 3.3.

Consider an SP polynomial $S(x,n,h)$, with $h>1$,

$S(x,n,h)=\prod_{\begin{subarray}{c}{\delta|h}\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(x)\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(x),$ (3.1)

where $\delta$ and $\tau$ range over the set of divisors of $h$ and $n+2-h$,
respectively. And

$S(x,n,1)=\prod_{\begin{subarray}{c}\tau_{1}|n+1\\\ \tau_{1}\neq
1\end{subarray}}\Psi_{\tau_{1}}(x),$ (3.2)

where $\tau_{1}$ ranges over the set of divisors of $n+1$. The symbol
$\Psi_{i}(z)$ denotes the $i^{th}$ cyclotomic polynomial.

###### Proof.

If $S(x,n,h)$ is an SP polynomial, we have

$a_{n}=a_{0}=1,a_{n-1}=a_{1}=2,\ldots,a_{n-h+1}=a_{h-1}=h,$

that is,

$S(x,n,h)=1+2x+\cdots+hx^{h-1}+\cdots+hx^{n+1-h}+\cdots+2x^{n-1}+x^{n}.$

Then, using Lemma 3.2 we have

$S(x,n,h)=(1+\cdots+x^{h-1})\cdot(1+\cdots+x^{n-h+1}).$ (3.3)

Since all the roots of all-one polynomials are roots of unity other than unity
itself, so

$S(x,n,h)=\frac{x^{h}-1}{x-1}\cdot\frac{x^{n+2-h}-1}{x-1}.$

Now from Theorem 2.3,

$\displaystyle=\frac{\displaystyle\prod_{\delta|h}\Psi_{\delta}(x)}{x-1}\cdot\frac{\displaystyle\prod_{\tau|n+2-h}\Psi_{\tau}(x)}{x-1}$
$\displaystyle=\displaystyle\prod_{\begin{subarray}{c}{\delta}|h\\\
{\delta}\neq
1\end{subarray}}\Psi_{\delta}(x)\cdot\prod_{\begin{subarray}{c}{\tau}|n+2-h\\\
{\tau}\neq 1\end{subarray}}\Psi_{\tau}(x),$

where $\delta$ and $\tau$ are divisors greater than $1$ of $h$ and $n+2-h$
respectively.

Next for $h=1$, the first all one polynomial on the right hand side of
Equation 3.3 is equal to $1$ and $n+2-h=n+1$. So we have

$S(x,n,1)=\prod_{\begin{subarray}{c}\tau_{1}|n+1\\\ \tau_{1}\neq
1\end{subarray}}\Psi_{\tau_{1}}(x),$

where $\tau_{1}$ is an element from the set of divisors of $n+1$. ∎

######  Example 2.

Consider an SP polynomial

$S(x,4,3)=x^{4}+2x^{3}+3x^{2}+2x+1$

of degree $4$. Then the divisor of both $h=3$ and $n+2-h=4+2-3=3$ is $3$, so

$\displaystyle S(x,4,3)$ $\displaystyle=\Psi_{3}(x)\cdot\Psi_{3}(x)$
$\displaystyle=(x^{2}+x+1)^{2}.$

######  Remark 3.4.

If $\zeta$ is a root of the polynomial $S(x,n,h)$, then the conjugate
$\overline{\zeta}$ is also a root of $S(x,n,h)$. Note that the only possible
real root of an SP polynomials is $-1$. Also, if $\zeta=-1$ is the root of an
SP polynomial $S(x,n,h)$, then either $n$ is odd or $n$ is even and $h=2$.

##  4 Extensions to related polynomials

In this section we will derive some extensions of SP polynomial which are
motivated by following properties of cyclotomic polynomials:

###  4.1 Some properties of cyclotomic polynomials

  1. 1.

For odd number $m>1$ , $\Psi_{m}(-x)=\Psi_{2m}(x).$

  2. 2.

If $p$ is a prime and $m$ is a positive integer. If $p$ divides $m$, then
$\Psi_{mp}(x)=\Psi_{m}(x^{p}).$

  3. 3.

If $p$ is prime and $m$ is a positive integer. If $p$ does not divide $m$,
then

$\Psi_{m}(x^{p})=\Psi_{pm}(x)\Psi_{m}(x)$

###  4.2 Alternating sign SP polynomials (ASP)

######  Definition 4.1.

Consider a polynomial

$P(x)=(-1)^{n}x^{n}+(-1)^{n-1}a_{n-1}x^{n-1}+...+(-1)a_{1}x+a_{0}$

with degree $n$. If the coefficients $a_{i}$ satisfy the conditions of the
coefficients of a SP polynomial then polynomial is called alternating sign SP
polynomial (ASP) of height $h$ and will be denoted by $T(x,n,h)$.

Using the properties of cyclotomic polynomials and Theorem 3.3, we can derive
a formula for ASP analogous to SP polynomials.

######  Lemma 4.2.

Let $T(x,n,h)$ be an $ASP$, then

$T(x,n,h)=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(-x)\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(-x),$

where $\delta$ and $\tau$ are divisors of $h$ and $n+2-h$ respectively. And

$T(x,n,1)=\prod_{\begin{subarray}{c}\tau_{1}|n+1\\\ \tau_{1}\neq
1\end{subarray}}\Psi_{\tau_{1}}(-x),$

where $\tau_{1}$ ranges over the set of divisors of $n+1$.

###### Proof.

Let $y=-x.$ Then $T(x,n,h)=S(-y,n,h)$ is a SP polynomial. Now using Theorem
3.3,

$\displaystyle T(x,n,h)$ $\displaystyle=\prod_{\begin{subarray}{c}\delta|h\\\
\delta\neq
1\end{subarray}}\Psi_{\delta}(y)\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(y)$
$\displaystyle=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(-x)\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(-x).$

Therefore,

$T(x,n,h)=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(-x)\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(-x),$

where $\Psi_{i}(-x)$ are cyclotomic polynomials obtained by plugging $-x$ in
the $i^{th}-$cyclotomic polynomials. Similarly, for $h=1$ using Theorem 3.3 we
get

$T(x,n,1)=\prod_{\begin{subarray}{c}\tau_{1}|n+1\\\ \tau_{1}\neq
1\end{subarray}}\Psi_{\tau_{1}}(-x).$

∎

###  4.3 Missing terms SP polynomials (MSP)

######  Definition 4.3.

Consider a polynomial

$P(x)=x^{n}+a_{(n-d)}x^{n-d}+...+a_{d}x^{d}+1$

of degree $n$. If $d$ is a divisor of $n$ and the coefficients of $P(x)$
satisfy the conditions of an SP polynomial of degree $n$ and height
$h=a_{\lceil\frac{n+1}{2}\rceil}$ then the polynomial will be called missing
terms SP polynomial (MSP) and is denoted by $M(x,n,h).$

######  Lemma 4.4.

Let $M(x,n,h)$ be an MSP polynomial, then

$M(x,n,h)=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(x^{d})\cdot\prod_{\begin{subarray}{c}\tau|\frac{n}{d}+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(x^{d}).$

where $\delta$ and $\tau$ are divisors of $h$ and $n+2-h$ respectively. And

$M(x,n,1)=\prod_{\begin{subarray}{c}\tau_{1}|n+1\\\ \tau_{1}\neq
1\end{subarray}}\Psi_{\tau_{1}}(x^{d}),$

where $\tau_{1}$ ranges over the set of divisors of $n+1$.

###### Proof.

Let $y=x^{d}.$ Since $d|n$ so we know $M(x,n,h)=S(y,n/d,h)$ is an SP
polynomial. Now using Theorem 3.3,

$\displaystyle M(x,n,h)$ $\displaystyle=S(y,n/d,h)$
$\displaystyle=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(y)\cdot\prod_{\begin{subarray}{c}\tau|\frac{n}{d}+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(y)$
$\displaystyle=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(x^{d})\cdot\prod_{\begin{subarray}{c}\tau|\frac{n}{d}+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(x^{d}).$

Therefore,

$M(x,n,h)=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(x^{d})\cdot\prod_{\begin{subarray}{c}\tau|\frac{n}{d}+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(x^{d}),$

where $\Psi_{i}(x^{d})$ are polynomials obtained by plugging $x^{d}$ in the
$i^{th}-$cyclotomic polynomials. Similarly, for $h=1$ using Theorem 3.3 we get

$M(x,n,1)=\prod_{\begin{subarray}{c}\tau_{1}|n+1\\\ \tau_{1}\neq
1\end{subarray}}\Psi_{\tau_{1}}(x^{d}).$

∎

###  4.4 Geometric SP polynomials (GSP)

Next, we will look at polynomials whose coefficients are not palindromic but
can be made palindromic by a particular substitution, such polynomials will be
called a geometric SP polynomial.

######  Definition 4.5.

Consider a polynomial

$G(x)=b_{n}x^{n}+b_{n-1}x^{n-1}+...+b_{1}x+1$

of degree $n$. If $a_{i}=\frac{b_{i}}{\alpha^{i}}$ are coefficients of an SP
polynomial, where $\alpha={b_{n}}^{\frac{1}{n}}$, then the polynomial $G(x)$
will be called a geometric SP polynomial (GSP).

######  Remark 4.6.

Note, if $\alpha=1$ then GSP is an SP polynomial and if $\alpha=-1$ then GSP
is an ASP polynomial.

######  Lemma 4.7.

Let $G(x)$ be a GSP polynomial of degree $n$, then

$G(x)=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(\alpha{x})\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(\alpha{x}).$

where $h=a_{\lceil\frac{n+1}{2}\rceil}>1.$ And if $h=1$

$G(x)=\prod_{\begin{subarray}{c}\tau_{1}|n+1\\\ \tau_{1}\neq
1\end{subarray}}\Psi_{\tau_{1}}(\alpha x),$

where $\tau_{1}$ ranges over the set of divisors of $n+1$.

###### Proof.

Let $y=\alpha{x}.$ Then

$\displaystyle G(x)$ $\displaystyle=b_{n}x^{n}+b_{n-1}x^{n-1}+...+b_{1}x+1.$
$\displaystyle=a_{n}(\alpha x)^{n}+a_{n-1}(\alpha x)^{n-1}+...+a_{1}(\alpha
x)+1.$ $\displaystyle=S(y,n,h).$

Now using Theorem 3.3, we have

$\displaystyle G(x)$ $\displaystyle=S(y,n,h)$
$\displaystyle=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(y)\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(y)$
$\displaystyle=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(\alpha{x})\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(\alpha{x}).$

Therefore,

$G(x)=\prod_{\begin{subarray}{c}\delta|h\\\ \delta\neq
1\end{subarray}}\Psi_{\delta}(\alpha{x})\cdot\prod_{\begin{subarray}{c}\tau|n+2-h\\\
\tau\neq 1\end{subarray}}\Psi_{\tau}(\alpha{x}),$

where $\Psi_{i}(\alpha{x})$ are polynomials obtained by plugging $\alpha{x}$
in the $i^{th}$cyclotomic polynomials. Similarly, for $h=1$ using Theorem 3.3
we get

$G(x)=\prod_{\begin{subarray}{c}\tau_{1}|n+1\\\ \tau_{1}\neq
1\end{subarray}}\Psi_{\tau_{1}}(\alpha x).$

∎

##  5 Acknowledgement

## References

  * [1] R. K.C., A. Algoud, & T.R. Young. “(In)stability of Cyclic Solutions Under Negative Feedback, To be submitted 2020\.
  * [2] T.R. Young, K. Prathom, and J. Rombouts, “Temporal clustering in cell cycle dynamics”, Dynamical Systems Magazine, 2019. Online at: https://dsweb.siam.org/
  * [3] S. Lang, “Algebra;”, Springer, pp. 276-277, 2002.
  * [4] J. Rotman, A. Cuoco, “Learning Modern Algebra”, Mathematical Association of America, (2013).