# Timing-residual power spectrum of a polarized stochastic gravitational-wave
background in pulsar-timing-array observation

Guo-Chin Liu1 and Kin-Wang Ng2,3 1Department of Physics, Tamkang University,
Tamsui, New Taipei City 25137, Taiwan
2Institute of Physics, Academia Sinica, Taipei 11529, Taiwan
3Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 11529,
Taiwan

###### Abstract

We study the observation of stochastic gravitational-wave background (SGWB)
made by pulsar-timing arrays in the spherical harmonic space. Instead of using
the Shapiro time delay, we keep the Sachs-Wolfe line-of-sight integral for the
timing residual of an observed pulsar. We derive the power spectrum of the
timing residual, from which the overlap reduction functions and the bipolar
spherical harmonics coefficients are constructed for the SGWB intensity and
polarization anisotropies. We have reproduced the previous results, noting
that we have developed a fast algorithm for computing accurate overlap
reduction functions and the bipolar spherical harmonics coefficients for the
linear-polarization anisotropy are worked out for the first time. Our
harmonic-space method is useful for future pulsar-timing-array observation on
a few thousand pulsars and provides optimal estimators for testing the
statistical isotropy of the SGWB.

##  I Introduction

The detection of gravitational waves (GWs) emitted from the coalescence of
binary black holes by the LIGO-Virgo experiment opens up a new era of GW
astronomy and cosmology ligo ; ligo2019 . LIGO/Virgo is a ground-based laser
interferometer that has two perpendicular detector-arms to measure GW strain
amplitudes. Upcoming and future GW interferometry experiments include KAGRA,
GEO600, Einstein Telescope, Cosmic Explorer, as well as space missions such as
LISA, DECIGO, Taiji, and TianQin ligo2050 , aiming to measure GWs at
frequencies ranging from kilohertz to millihertz. Pulsar timing is another
method to detect GWs by monitoring the arrival times of radio pulses from
pulsars with ground-based radio telescopes pta2019 . The line of sight from
the telescope to the precise pulsar clock acts like a detector-arm sensitive
to the passage of nanohertz GWs through the space. Current pulsar timing array
(PTA) experiments, monitoring roughly 100 Galactic millisecond pulsars,
include EPTA EPTA , NANOGrav NANO , and PPTA PPTA . The future SKA project
will observe about 6000 Galactic millisecond pulsars to reach a sensitivity
three to four orders of magnitude better than the current PTAs SKA ; SKA2 .

Stochastic gravitational wave background (SGWB) is a key science goal in GW
experiments. GWs are very weakly interacting. The observation of SGWB enable
us to see directly the physical processes that produce GWs in the early
Universe such as distant compact binary coalescences, early-time phase
transitions, cosmic string or defect networks, second-order primordial scalar
perturbations, and inflationary GWs romano . The SGWB is predicted to be
highly isotropic; however, it has been proposed that it could be anisotropic
and even circularly or linearly polarized alexander ; satoh ; sorbo ;
crowder2013 ; cusin ; bartolo ; pitrou ; liu21 ; cai21 .

In search of SGWB in an interferometry network, the responses of a pair of
detectors to the GW strain amplitude are correlated so as to filter out
detector noises and increase the signal-to-noise ratio romano . In PTA
observation, quadrupolar spatial correlations between pulsar pairs are used to
identify the presence of a SGWB pta2019 . Recently, the NANOGrav Collaboration
nanograv  has found strong evidence of a stochastic common-spectrum process
across 45 ms pulsars, hinting at a SGWB with the spectral energy density of
$\Omega_{\rm GW}\simeq 5.0\times 10^{-9}$ at a reference frequency of
$32\,{\rm nHz}$. A consistent common-spectrum process has also been found in
the second data release of the Parkes Pulsar Timing Array (PPTA) PPTA21 , the
Data Release 2 of the European Pulsar Timing Array (EPTA) covering a timespan
up to 24 years EPTA21 , and the second data release of the International
Pulsar Timing Array (IPTA) synthesizing decadal-length pulsar-timing campaigns
IPTA22 . However, the observations have not found statistically significant
evidence that this process has quadrupolar interpulsar correlations. Lately,
using data from Advanced LIGO’s and Advanced Virgo’s third observing run (O3)
combined with the earlier O1 and O2 runs, upper limits have been derived on an
isotropic SGWB, $\Omega_{\rm GW}<3.9\times 10^{-10}$ at $25\,{\rm Hz}$
ligo2101 , and on anisotropic SGWB, $\Omega_{\rm GW}<(0.56-9.7)\times
10^{-9}{\rm sr}^{-1}$ ligo2103 .

In this paper, we will study the observation of SGWB intensity and
polarization anisotropies in PTAs. There has been a lot of theoretical studies
on the pulsar-timing observation of SGWB pta2019 . Previous works have been
mostly based on spatial correlation functions of the timing residuals from the
Earth-term contribution for the intensity anisotropy anholm ; mingar13 ;
gair14  and for the circular-polarization anisotropy kato16 . The
modifications by the pulsar term to the intensity correlation functions have
been discussed mingar14 . In Ref. qin19 , angular power spectra of the line-
of-sight integral for the timing residual for the SGWB intensity and circular-
polarization anisotropies have been derived in the total-angular-momentum
formalism that uses chiral spherical gravitational waves. Adopting the
technique to expand the polarization basis tensors in terms of spin-weighted
spherical harmonics chu21 , a numerical scheme has been developed to calculate
the correlation functions, being extended to including the linear-polarization
anisotropy chu2107 . Here, following Ref. ng21 , we will formulate the problem
in the spherical harmonic space, directly deriving the power spectrum of the
timing residual of an observed pulsar. We will use the Sachs-Wolfe line-of-
sight integral for the pulsar timing residual, adopted in Refs. ng21 ; chu2107
, which have shown that the integral conveniently incorporates the effects of
the pulsar term.

Many papers have laid out schemes for measuring the SGWB intensity and
polarization anisotropies in pulsar timing data analyses, such as the Bayesian
parameter-estimation pipeline taylor1 ; taylor2 , the spherical-harmonic power
spectrum estimators gair14 ; kam1 ; kam2 , and the Fisher matrix of observed
pulsar pairs smith1 ; smith2 . Some of the limitations of using the spherical-
harmonics approach such as being computationally demanding in the analysis
pipeline and an inhomogeneous sky coverage of pulsars in real PTA observation
have been discussed smith1 ; smith2 . However, the results in this work will
be useful in future PTA observation on a large number of pulsars.

The paper is organized as follows. We will introduce a polarized SGWB and its
Stokes parameters in the next section, followed by a brief account of the
pulsar timing in Sec. III. In Secs. IV and V, the power spectrum of the timing
residual will be derived. We will obtain the overlap reduction functions in
the celestial coordinates in Sec. VI and in the computational frame in Sec.
VII. We will briefly mention the bipolar spherical harmonics coefficients in
Sec. VIII. Sec. IX is our conclusion.

##  II Polarized SGWB

In the Minkowskian spacetime $(t,\vec{x})$, the metric perturbation $h_{ij}$
in the transverse traceless gauge depicts travelling GWs at the speed of light
$c=\omega/k$. It can be expanded by Fourier modes as

$h_{ij}(t,\vec{x})=\sum_{A}\int_{-\infty}^{\infty}\mathbf{d}f\int_{S^{2}}\mathbf{d}\hat{k}\;h_{A}(f,\hat{k})\mathbf{e}^{A}_{ij}(\hat{k})e^{-2\pi
if(t-\hat{k}\cdot\vec{x}/c)}\,,$ (1)

where $A$ stands for the polarization of GWs with basis tensors
$\mathbf{e}^{A}_{ij}(\hat{k})$, which are transverse to the propagation
direction, $\hat{k}$. Here $h_{ij}$ is treated as real, so the Fourier
components with negative frequencies are given by
$h_{A}(-f,\hat{k})=h_{A}^{*}(f,\hat{k})$ for all $f\geq 0$. We define a SGWB
as a collection of GWs satisfying the condition that $h_{ij}$ are random
Gaussian fields with a statistical behavior completely characterized by the
two-point correlation function $\langle
h_{ij}(t,\vec{x}_{1})h_{ij}(t,\vec{x}_{2})\rangle$, where the angle brackets
denote their ensemble averages. The ensemble averages of the Fourier modes
have the following form

$\langle
h_{A}(f,\hat{k})h^{*}_{A^{\prime}}(f^{\prime},\hat{k}^{\prime})\rangle=\delta(f-f^{\prime})\delta(\hat{k}-\hat{k}^{\prime})P_{AA^{\prime}}(f,\hat{k})\,,$
(2)

where the spatial translational invariance dictates the delta function of
their 3-momenta, $\delta(\vec{k}-\vec{k}^{\prime})$. Note that the power
spectra $P_{AA^{\prime}}(f,\hat{k})$ remain to be direction dependent.

For GWs coming from the sky direction $-\hat{k}$ with wave vector $\vec{k}$,
it is customary to write the polarization basis tensors in terms of the basis
vectors in the spherical coordinates,

$\displaystyle\mathbf{e}^{+}(\hat{k})$
$\displaystyle=\hat{\mathbf{e}}_{\theta}\otimes\hat{\mathbf{e}}_{\theta}-\hat{\mathbf{e}}_{\phi}\otimes\hat{\mathbf{e}}_{\phi}\,,$
$\displaystyle\mathbf{e}^{\times}(\hat{k})$
$\displaystyle=\hat{\mathbf{e}}_{\theta}\otimes\hat{\mathbf{e}}_{\phi}+\hat{\mathbf{e}}_{\phi}\otimes\hat{\mathbf{e}}_{\theta}\,,$
(3)

in which $\hat{\mathbf{e}}_{\theta}$, $\hat{\mathbf{e}}_{\phi}$, and $\hat{k}$
form a right-handed orthonormal basis. Also, we can define the complex
circular polarization basis tensors as

$\displaystyle\mathbf{e}_{R}$
$\displaystyle=\frac{(\mathbf{e}_{+}+i\mathbf{e}_{\times})}{\sqrt{2}}\,,$
$\displaystyle\mathbf{e}_{L}$
$\displaystyle=\frac{(\mathbf{e}_{+}-i\mathbf{e}_{\times})}{\sqrt{2}}\,,$ (4)

where $\mathbf{e}_{R}$ stands for the right-handed GW with a positive helicity
while $\mathbf{e}_{L}$ stands for the left-handed GW with a negative helicity.
The corresponding amplitudes in Eq. (1) in the two different bases are related
to each other via

$\displaystyle h_{R}$ $\displaystyle=\frac{(h_{+}-ih_{\times})}{\sqrt{2}}\,,$
$\displaystyle h_{L}$ $\displaystyle=\frac{(h_{+}+ih_{\times})}{\sqrt{2}}\,.$
(5)

Analogous to the case in electromagnetic waves book:BornAndWolf , the
coherency matrix $P_{AA^{\prime}}$ in Eq. (2) is related to the Stokes
parameters, $I$, $Q$, $U$, and $V$ as

$\displaystyle I$ $\displaystyle=\left[\langle h_{R}h_{R}^{*}\rangle+\langle
h_{L}h_{L}^{*}\rangle\right]/2\,,$ $\displaystyle Q+iU$ $\displaystyle=\langle
h_{L}h_{R}^{*}\rangle\,,$ $\displaystyle Q-iU$ $\displaystyle=\langle
h_{R}h_{L}^{*}\rangle\,,$ $\displaystyle V$ $\displaystyle=\left[\langle
h_{R}h_{R}^{*}\rangle-\langle h_{L}h_{L}^{*}\rangle\right]/2\,,$ (6)

which are functions of the frequency $f$ and the propagation direction
$\hat{k}$. $I$ is the intensity, $Q$ and $U$ represent the linear
polarization, and $V$ is the circular polarization.

##  III Pulsar timing

In the pulsar-timing observation, radio pulses from an array of roughly 100
Galactic millisecond pulsars are being monitored with ground-based radio
telescopes. The redshift fluctuation of a pulsar in the pointing direction
$\hat{e}$ on the sky is given by the Sachs-Wolfe effect sachs ,

$z(\hat{e})=-{1\over
2}\int_{\eta_{e}}^{\eta_{r}}\mathbf{d}\eta\,\hat{e}^{i}\hat{e}^{j}\frac{\partial}{\partial\eta}h_{ij}(\eta,\vec{x})\,,$
(7)

where the lower (upper) limit of integration in the line-of-sight integral
represents the point of emission (reception) of the radio pulse. The physical
distance of the pulsar from the Earth is

$D=c(\eta_{r}-\eta_{e})\,,$ (8)

which is of order $1\,{\rm kpc}$.

The quantity that is actually observed in the pulsar-timing observation is the
timing residual counted as

$r(t)=\int_{0}^{t}\mathbf{d}t^{\prime}z(t^{\prime})\,,$ (9)

where $t^{\prime}$ denotes the laboratory time and $t$ is the duration of the
observation. Using the laboratory time $t^{\prime}$, we rewrite Eq. (7) as

$z(t^{\prime},\hat{e})=-{1\over
2}\int_{t^{\prime}+\eta_{e}}^{t^{\prime}+\eta_{r}}\mathbf{d}\eta\,d^{ij}\frac{\partial}{\partial\eta}h_{ij}(\eta,\vec{x})\,,$
(10)

where the detector tensor is

$d^{ij}=\hat{e}^{i}\hat{e}^{j}\,.$ (11)

##  IV Timing-residual power spectrum

Then, replacing $\vec{x}$ by $c(\eta_{r}-\eta)\hat{e}$ in Eq. (10) and using
the spherical wave expansion (66) for the phane wave (1), Eq. (9) becomes

$r(t,\hat{e})=2\pi\sum_{A}\int_{-\infty}^{\infty}\mathbf{d}f\int_{S^{2}}\mathbf{d}\hat{k}\int_{\eta_{e}}^{\eta_{r}}\mathbf{d}\eta\,(1-e^{-2\pi
ift})h_{A}(f,\hat{k})d^{ij}\mathbf{e}^{A}_{ij}(\hat{k})e^{-2\pi
if\eta}\sum_{LM}i^{L}j_{L}[2\pi
f(\eta_{r}-\eta)]Y_{LM}^{*}(\hat{k})Y_{LM}(\hat{e})\,.$ (12)

We expand

$r(t,\hat{e})=\sum_{\ell m}a_{\ell m}Y_{\ell m}(\hat{e})\,.$ (13)

Defining $x=2\pi f(\eta_{r}-\eta)$, we have

$\displaystyle a_{\ell m}$
$\displaystyle=\int_{-\infty}^{\infty}\frac{\mathbf{d}f}{2\pi f}\,(1-e^{-2\pi
ift})e^{-2\pi
if\eta_{r}}\sum_{A}\int_{S^{2}}\mathbf{d}\hat{k}\,h_{A}(f,\hat{k})J^{A}_{\ell
m}(fD,\hat{k})\,,$ $\displaystyle J^{A}_{\ell m}(fD,\hat{k})$
$\displaystyle\equiv\sum_{LM}2\pi\,i^{L}Y_{LM}^{*}(\hat{k})\int_{0}^{2\pi
fD/c}\mathbf{d}x\,e^{ix}j_{L}(x)\int_{S^{2}}\mathbf{d}\hat{e}\,d^{ij}\mathbf{e}^{A}_{ij}(\hat{k})Y_{LM}(\hat{e})Y^{*}_{\ell
m}(\hat{e})\,.$ (14)

The timing-residual correlation between a pair of Galactic pulsars $a$ and $b$
is constructed as

$\langle
r(t_{a},\hat{e}_{a})r(t_{b},\hat{e}_{b})\rangle=\int_{0}^{t_{a}}\mathbf{d}t^{\prime}\int_{0}^{t_{b}}\mathbf{d}t{{}^{\prime\prime}}\langle
z(t^{\prime},\hat{e}_{a})z(t{{}^{\prime\prime}},\hat{e}_{b})\rangle=\sum_{\ell_{1}m_{1}\ell_{2}m_{2}}\langle
a_{\ell_{1}m_{1}}a^{*}_{\ell_{2}m_{2}}\rangle
Y_{\ell_{1}m_{1}}(\hat{e}_{a})Y^{*}_{\ell_{2}m_{2}}(\hat{e}_{b})\,,$ (15)

where the ensemble average is given by Eq. (2) as

$\langle
a_{\ell_{1}m_{1}}a^{*}_{\ell_{2}m_{2}}\rangle=\int_{-\infty}^{\infty}\frac{\mathbf{d}f}{(2\pi
f)^{2}}\,(1-e^{-2\pi ift_{a}})(1-e^{2\pi
ift_{b}})\sum_{A_{1}A_{2}}\int_{S^{2}}\mathbf{d}\hat{k}\,P_{A_{1}A_{2}}(f,\hat{k})J^{A_{1}}_{\ell_{1}m_{1}}(fD_{a},\hat{k})J^{A_{2}*}_{\ell_{2}m_{2}}(fD_{b},\hat{k})\,.$
(16)

In terms of the Stokes parameters in Eq. (6) and the definitions,

$\displaystyle\mathbb{J}^{I}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})$
$\displaystyle\equiv
J^{R}_{\ell_{1}m_{1}}(fD_{a},\hat{k})J^{R*}_{\ell_{2}m_{2}}(fD_{b},\hat{k})+J^{L}_{\ell_{1}m_{1}}(fD_{a},\hat{k})J^{L*}_{\ell_{2}m_{2}}(fD_{b},\hat{k})\,,$
$\displaystyle\mathbb{J}^{V}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})$
$\displaystyle\equiv
J^{R}_{\ell_{1}m_{1}}(fD_{a},\hat{k})J^{R*}_{\ell_{2}m_{2}}(fD_{b},\hat{k})-J^{L}_{\ell_{1}m_{1}}(fD_{a},\hat{k})J^{L*}_{\ell_{2}m_{2}}(fD_{b},\hat{k})\,,$
$\displaystyle\mathbb{J}^{Q+iU}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})$
$\displaystyle\equiv
J^{L}_{\ell_{1}m_{1}}(fD_{a},\hat{k})J^{R*}_{\ell_{2}m_{2}}(fD_{b},\hat{k})\,,$
$\displaystyle\mathbb{J}^{Q-iU}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})$
$\displaystyle\equiv
J^{R}_{\ell_{1}m_{1}}(fD_{a},\hat{k})J^{L*}_{\ell_{2}m_{2}}(fD_{b},\hat{k})\,,$
(17)

we have

$\langle
a_{\ell_{1}m_{1}}a^{*}_{\ell_{2}m_{2}}\rangle=\int_{-\infty}^{\infty}\frac{\mathbf{d}f}{(2\pi
f)^{2}}\,(1-e^{-2\pi ift_{a}})(1-e^{2\pi ift_{b}})\sum_{X=\\{I,V,Q\pm
iU\\}}\int_{S^{2}}\mathbf{d}\hat{k}\,X(f,\hat{k})\,\mathbb{J}^{X}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})\,.$
(18)

We further expand the Stokes parameters in terms of ordinary and spin-weighted
spherical harmonics as

$\displaystyle I(f,\hat{k})$ $\displaystyle=\sum_{\ell m}I_{\ell
m}(f)\;Y_{\ell m}(\hat{k})\,,$ $\displaystyle V(f,\hat{k})$
$\displaystyle=\sum_{\ell m}V_{\ell m}(f)\;Y_{\ell m}(\hat{k})\,,$
$\displaystyle(Q+iU)(f,\hat{k})$ $\displaystyle=\sum_{\ell m}(Q+iU)_{\ell
m}(f)\;_{+4}Y_{\ell m}(\hat{k})\,,$ $\displaystyle(Q-iU)(f,\hat{k})$
$\displaystyle=\sum_{\ell m}(Q-iU)_{\ell m}(f)\;_{-4}Y_{\ell m}(\hat{k})\,,$
(19)

where the specific combinations, $Q\pm iU$, make them become spin $\pm 4$
objects so that we can expand them nicely by the corresponding spin-weighted
spherical harmonics. A brief introduction to the spin-weighted spherical
harmonics is found in Appendix A.

Hence, we can express the timing-residual correlation in the following form

$\langle
r(t_{a},\hat{e}_{a})r(t_{b},\hat{e}_{b})\rangle=\int_{-\infty}^{\infty}\frac{\mathbf{d}f}{(2\pi
f)^{2}}\,(1-e^{-2\pi ift_{a}})(1-e^{2\pi ift_{b}})\sum_{X=\\{I,V,Q\pm
iU\\}}\sum_{\ell m}X_{\ell m}(f)\gamma_{\ell
m}^{X}(fD_{a},fD_{b};\hat{e}_{a},\hat{e}_{b})\,,$ (20)

where the overlap reduction functions (ORFs) are given by

$\displaystyle\gamma_{\ell m}^{I,V}(fD_{a},fD_{b};\hat{e}_{a},\hat{e}_{b})$
$\displaystyle=\sum_{\ell_{1}m_{1}\ell_{2}m_{2}}Y_{\ell_{1}m_{1}}(\hat{e}_{a})Y^{*}_{\ell_{2}m_{2}}(\hat{e}_{b})\int_{S^{2}}\mathbf{d}\hat{k}\,Y_{\ell
m}(\hat{k})\,\mathbb{J}^{I,V}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})\,,$
(21) $\displaystyle\gamma_{\ell m}^{Q\pm
iU}(fD_{a},fD_{b};\hat{e}_{a},\hat{e}_{b})$
$\displaystyle=\sum_{\ell_{1}m_{1}\ell_{2}m_{2}}Y_{\ell_{1}m_{1}}(\hat{e}_{a})Y^{*}_{\ell_{2}m_{2}}(\hat{e}_{b})\int_{S^{2}}\mathbf{d}\hat{k}\;_{\pm
4}Y_{\ell m}(\hat{k})\,\mathbb{J}^{Q\pm
iU}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})\,.$ (22)

Eqs. (21) and (22) are the most general ORFs for a pair of Galactic pulsars
$a$ and $b$, respectively, at distances $D_{a}$ and $D_{b}$ from the Earth.

##  V Calculation of $J^{A}_{\ell m}(fD,\hat{k})$

Now we calculate the contribution of a $k$-mode to the redshift fluctuation of
a pulsar, namely $J^{A}_{\ell m}(fD,\hat{k})$ in Eq. (14). For convenience, we
first assume that $\hat{k}$ points to the direction of the polar axis or
$\hat{\mathbf{z}}$-axis. In this case, the basis vectors (II) become

$\displaystyle\mathbf{e}^{+}(\hat{k})$
$\displaystyle=\hat{\mathbf{x}}\otimes\hat{\mathbf{x}}-\hat{\mathbf{y}}\otimes\hat{\mathbf{y}}\,,$
$\displaystyle\mathbf{e}^{\times}(\hat{k})$
$\displaystyle=\hat{\mathbf{x}}\otimes\hat{\mathbf{y}}+\hat{\mathbf{y}}\otimes\hat{\mathbf{x}}\,,$
(23)

and we have

$\displaystyle\hat{e}=\sin\theta\cos\phi\,\hat{\mathbf{x}}+\sin\theta\sin\phi\,\hat{\mathbf{y}}+\cos\theta\,\hat{\mathbf{z}}\,,$
(24) $\displaystyle Y_{LM}(\hat{k})=\sqrt{\frac{2L+1}{4\pi}}\delta_{M0}\,.$
(25)

This gives

$d^{ij}\mathbf{e}^{R,L}_{ij}(\hat{k})=4\sqrt{\frac{\pi}{15}}Y_{2\pm
2}(\hat{e})\,,$ (26)

where the helicity $R$ takes the value of $2$ and $L$ the value of $-2$.
Hence, we obtain

$J^{R,L}_{\ell
m}(fD,\hat{\mathbf{z}})=\sum_{L}4\pi\,i^{L}\sqrt{\frac{2L+1}{15}}\int_{0}^{2\pi
fD/c}\mathbf{d}x\,e^{ix}j_{L}(x)\int_{S^{2}}\mathbf{d}\hat{e}\;Y_{2\pm
2}(\hat{e})\;Y_{L0}(\hat{e})\;Y^{*}_{\ell m}(\hat{e})\,.$ (27)

From Eq. (67), we have

$\int_{S^{2}}\mathbf{d}\hat{e}\;Y_{2\pm
2}(\hat{e})\;Y_{L0}(\hat{e})\;Y^{*}_{\ell
m}(\hat{e})=(-1)^{m}\sqrt{\frac{5(2L+1)(2\ell+1)}{4\pi}}\begin{pmatrix}2&&L&&\ell\\\
0&&0&&0\end{pmatrix}\begin{pmatrix}2&&L&&\ell\\\ \pm 2&&0&&-m\end{pmatrix}\,,$
(28)

which vanishes unless $m=\pm 2$ and $L=\ell-2,\ell,\ell+2$. These nonzero
integral values are given by

$\displaystyle\int_{S^{2}}\mathbf{d}\hat{e}\;Y_{2\pm
2}(\hat{e})\;Y_{\ell-2,0}(\hat{e})\;Y^{*}_{\ell\pm 2}(\hat{e})$
$\displaystyle={1\over
4}\sqrt{\frac{15}{2\pi}}\left[\frac{(\ell-1)\ell(\ell+1)(\ell+2)}{(2\ell-3)(2\ell-1)^{2}(2\ell+1)}\right]^{1\over
2},$ $\displaystyle\int_{S^{2}}\mathbf{d}\hat{e}\;Y_{2\pm 2}(\hat{e})\;Y_{\ell
0}(\hat{e})\;Y^{*}_{\ell\pm 2}(\hat{e})$ $\displaystyle=-{1\over
2}\sqrt{\frac{15}{2\pi}}\left[\frac{(\ell-1)\ell(\ell+1)(\ell+2)}{(2\ell-1)^{2}(2\ell+3)^{2}}\right]^{1\over
2},$ $\displaystyle\int_{S^{2}}\mathbf{d}\hat{e}\;Y_{2\pm
2}(\hat{e})\;Y_{\ell+2,0}(\hat{e})\;Y^{*}_{\ell\pm 2}(\hat{e})$
$\displaystyle={1\over
4}\sqrt{\frac{15}{2\pi}}\left[\frac{(\ell-1)\ell(\ell+1)(\ell+2)}{(2\ell+1)(2\ell+3)^{2}(2\ell+5)}\right]^{1\over
2}.$ (29)

Substituting them in Eq. (27), we obtain

$J^{R,L}_{\ell m}(fD,\hat{\mathbf{z}})=-\delta_{m\pm
2}2\pi\,i^{\ell}\sqrt{\frac{(2\ell+1)}{8\pi}\frac{(\ell+2)!}{(\ell-2)!}}\int_{0}^{2\pi
fD/c}\mathbf{d}x\,e^{ix}\left[\frac{j_{\ell-2}(x)}{(2\ell-1)(2\ell+1)}+\frac{2j_{\ell}(x)}{(2\ell-1)(2\ell+3)}+\frac{j_{\ell+2}(x)}{(2\ell+1)(2\ell+3)}\right].$
(30)

This can be cast into a compact form by using the recursion relation,
$j_{\ell}(x)/x=[j_{\ell-1}(x)+j_{\ell+1}(x)]/(2\ell+1)$, which gives

$\displaystyle J^{R,L}_{\ell m}(fD,\hat{\mathbf{z}})=$
$\displaystyle-\delta_{m\pm
2}2\pi\,i^{\ell}\sqrt{\frac{(2\ell+1)}{8\pi}\frac{(\ell+2)!}{(\ell-2)!}}\;\times$
$\displaystyle\int_{0}^{2\pi
fD/c}\mathbf{d}x\,e^{ix}\frac{j_{\ell}(x)}{x^{2}}.$ (31)

Through a three-dimensional rotation that takes the $\hat{\mathbf{z}}$-axis
into the direction $\hat{k}$, we can relate chu21

$J^{R,L}_{\ell
m}(fD,\hat{k})=\sum_{m^{\prime}}D^{\ell\,*}_{m^{\prime}m}(-\alpha,-\theta,-\phi)J^{R,L}_{\ell
m^{\prime}}(fD,\hat{\mathbf{z}})\,,$ (32)

where $\hat{k}=(\theta,\phi)$. Here, the Wigner-D matrix is given by

$D^{\ell}_{m^{\prime}m}(-\alpha,-\theta,-\phi)=\sqrt{\frac{4\pi}{2\ell+1}}{}_{-m^{\prime}}Y_{\ell
m}(\theta,\phi)e^{im^{\prime}\alpha}\,,$ (33)

where $e^{im^{\prime}\alpha}$ is a redundant phase that reflects a remaining
degree of freedom in the rotation about the $\hat{k}$-direction. Hence, we
have

$\displaystyle J^{R}_{\ell m}(fD,\hat{k})$
$\displaystyle=D^{\ell\,*}_{2m}(-\alpha,-\theta,-\phi)J^{R}_{\ell
2}(fD,\hat{\mathbf{z}})\,,$ $\displaystyle J^{L}_{\ell m}(fD,\hat{k})$
$\displaystyle=D^{\ell\,*}_{-2m}(-\alpha,-\theta,-\phi)J^{L}_{\ell-2}(fD,\hat{\mathbf{z}})\,,$
(34)

noting that $J^{R}_{\ell
2}(fD,\hat{\mathbf{z}})=J^{L}_{\ell-2}(fD,\hat{\mathbf{z}})$ and
$e^{im^{\prime}\alpha}$ will not appear in physical observables.

##  VI Overlap Reduction Functions in the Celestial Coordinates

Inserting the results (34) into Eq. (17), we have

$\displaystyle\mathbb{J}^{I,V}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})$
$\displaystyle=(-1)^{m_{1}}\left[{}_{2}Y_{\ell_{1}-m_{1}}(\hat{k}){}_{-2}Y_{\ell_{2}m_{2}}(\hat{k})\pm{}_{-2}Y_{\ell_{1}-m_{1}}(\hat{k}){}_{2}Y_{\ell_{2}m_{2}}(\hat{k})\right]J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})\,,$
$\displaystyle\mathbb{J}^{Q\pm
iU}_{\ell_{1}m_{1}\ell_{2}m_{2}}(fD_{a},fD_{b},\hat{k})$
$\displaystyle=(-1)^{m_{1}}{}_{\mp 2}Y_{\ell_{1}-m_{1}}(\hat{k}){}_{\mp
2}Y_{\ell_{2}m_{2}}(\hat{k})\,e^{\pm
i4\alpha}J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})\,,$ (35)

where we have defined the function

$J_{\ell}(fD)=\sqrt{2}\,\pi\,i^{\ell}\sqrt{\frac{(\ell+2)!}{(\ell-2)!}}\int_{0}^{2\pi
fD/c}\mathbf{d}x\,e^{ix}\frac{j_{\ell}(x)}{x^{2}},$ (36)

and the phase factors $e^{\pm i4\alpha}$ are resulted from a rotation of angle
$\alpha$ about the $\hat{k}$-direction on the spin-4 objects,
$\mathbb{J}^{Q\pm iU}_{\ell_{1}m_{1}\ell_{2}m_{2}}$, respectively.
Simultaneously, the spin-4 spherical harmonics ${}_{\pm 4}Y_{\ell m}(\hat{k})$
in Eq. (22) are augmented by the same phase factors of opposite signs $e^{\mp
i4\alpha}$ under the rotation, which exactly cancel $e^{\pm i4\alpha}$ from
$\mathbb{J}^{Q\pm iU}_{\ell_{1}m_{1}\ell_{2}m_{2}}$. Using Eq. (67) and the
property (68), the ORFs in Eqs. (21) and (22) become

$\displaystyle\gamma_{\ell
m}^{I,V}(fD_{a},fD_{b};\hat{e}_{a},\hat{e}_{b})=\sum_{\ell_{1}m_{1}\ell_{2}m_{2}}$
$\displaystyle(-1)^{m_{1}}\left[1\pm(-1)^{\ell+\ell_{1}+\ell_{2}}\right]J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})\;Y_{\ell_{1}m_{1}}(\hat{e}_{a})Y^{*}_{\ell_{2}m_{2}}(\hat{e}_{b})\,\times$
$\displaystyle\sqrt{\frac{(2\ell+1)(2\ell_{1}+1)(2\ell_{2}+1)}{4\pi}}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
0&&-2&&2\end{pmatrix}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
m&&-m_{1}&&m_{2}\end{pmatrix}\,,$ (37) $\displaystyle\gamma_{\ell m}^{Q\pm
iU}(fD_{a},fD_{b};\hat{e}_{a},\hat{e}_{b})=\sum_{\ell_{1}m_{1}\ell_{2}m_{2}}$
$\displaystyle(-1)^{m_{1}}J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})\;Y_{\ell_{1}m_{1}}(\hat{e}_{a})Y^{*}_{\ell_{2}m_{2}}(\hat{e}_{b})\,\times$
$\displaystyle\sqrt{\frac{(2\ell+1)(2\ell_{1}+1)(2\ell_{2}+1)}{4\pi}}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
\mp 4&&\pm 2&&\pm 2\end{pmatrix}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
m&&-m_{1}&&m_{2}\end{pmatrix}\,.$ (38)

The two lowest moments, $\gamma_{00}^{I}$ and $\gamma_{00}^{V}$, select the
unpolarized and the circularly polarized components of an isotropic SGWB,
respectively. When $\ell=m=0$, the Wigner-3j symbol is proportional to
$\delta_{\ell_{1}\ell_{2}}\delta_{m_{1}m_{2}}$. It immediately gives us

$\gamma_{00}^{V}=0\,,$ (39)

whereas

$\displaystyle\gamma_{00}^{I}$
$\displaystyle=\sum_{l}\frac{2l+1}{4\pi}C_{l}P_{l}(\hat{e}_{a}\cdot\hat{e}_{b})\quad{\rm
with}$ (40) $\displaystyle C_{l}$
$\displaystyle\equiv\frac{1}{\sqrt{\pi}}J_{l}(fD_{a})J^{*}_{l}(fD_{b})\,,$
(41)

which depends solely on the separation angle as expected for an isotropic
SGWB. The power spectrum $C_{l}$ has an analytic form under the limit that
$fD_{a}\gg c$ and $fD_{b}\gg c$ qin19 ; ng21 . Using the integral result

$\int_{0}^{\infty}\mathbf{d}x\,e^{ix}\frac{j_{l}(x)}{x^{2}}=2i^{l-1}\frac{(l-2)!}{(l+2)!}\,,$
(42)

Eq. (36) can be approximated as

$J_{l}(fD)|_{fD/c\rightarrow\infty}=2^{3/2}\pi\,i^{2l-1}\sqrt{\frac{(l-2)!}{(l+2)!}}\,,$
(43)

which gives

$C_{l}=\frac{8\pi^{3/2}}{(l+2)(l+1)l(l-1)}\,.$ (44)

It was shown gair14  that the $\gamma_{00}^{I}$ (40) with this $C_{l}$
reproduces the Hellings and Downs curve for the quadrupolar interpulsar
correlations downs . For finite $fD_{a}$ and $fD_{b}$, one needs to calculate
numerically the power spectrum (41) to get corrections to the Hellings and
Downs curve ng21 .

For higher multipole moments, we can easily set up a numerical scheme, similar
to that in Ref. chu2107 , to compute the ORFs in Eqs. (37) and (38) for any
pair of Galactic pulsars on the sky with known distances and coordinates,
$(D_{a},\theta_{a},\phi_{a})$ and $(D_{b},\theta_{b},\phi_{b})$. The factor
$J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})$ is generally a complex number.
When $fD_{a}\gg c$ and $fD_{b}\gg c$, we can approximate it as

$J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})\simeq
8\pi^{2}(-1)^{\ell_{1}+\ell_{2}}\sqrt{\frac{(\ell_{1}-2)!\,(\ell_{2}-2)!}{(\ell_{1}+2)!\,(\ell_{2}+2)!}}\,.$
(45)

##  VII Overlap Reduction Functions in the Computational Frame

To compare the present method with previous works anholm ; mingar13 ; gair14 ;
kato16 ; mingar14 ; chu2107 , we also compute the ORFs in the so-called
computational frame: pulsar $a$ is placed along the $\hat{\mathbf{z}}$-axis
while pulsar $b$ is in the $\hat{\mathbf{x}}$-$\hat{\mathbf{z}}$ plane. Then,
their polar coordinates are given by

$\hat{e}_{a}=(0,0),\quad\hat{e}_{b}=(\zeta,0)\,,$ (46)

where $\zeta$ is their separation angle, and we have

$\displaystyle Y_{\ell_{1}m_{1}}(\hat{e}_{a})$
$\displaystyle=\sqrt{\frac{2\ell_{1}+1}{4\pi}}\delta_{m_{1}0},$ $\displaystyle
Y^{*}_{\ell_{2}m_{2}}(\hat{e}_{b})$
$\displaystyle=Y^{*}_{\ell_{2}m_{2}}(\zeta,0)=Y_{\ell_{2}m_{2}}(\zeta,0).$
(47)

Hence, Eqs. (37) and (38) simplify to

$\displaystyle\gamma_{\ell
m}^{I,V}(fD_{a},fD_{b},\zeta)=\sum_{\ell_{1}\ell_{2}}$
$\displaystyle(-1)^{m}\frac{2\ell_{1}+1}{4\pi}\left[1\pm(-1)^{\ell+\ell_{1}+\ell_{2}}\right]J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})\;Y_{\ell_{2}m}(\zeta,0)\,\times$
$\displaystyle\sqrt{(2\ell+1)(2\ell_{2}+1)}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
0&&-2&&2\end{pmatrix}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
m&&0&&-m\end{pmatrix}\,,$ (48) $\displaystyle\gamma_{\ell m}^{Q\pm
iU}(fD_{a},fD_{b},\zeta)=\sum_{\ell_{1}\ell_{2}}$
$\displaystyle(-1)^{m}\frac{2\ell_{1}+1}{4\pi}J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})\;Y_{\ell_{2}m}(\zeta,0)\,\times$
$\displaystyle\sqrt{(2\ell+1)(2\ell_{2}+1)}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
\mp 4&&\pm 2&&\pm 2\end{pmatrix}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
m&&0&&-m\end{pmatrix}\,.$ (49)

Using the property (68), it is straightforward to show that the ORFs in the
computational frame have the conjugate relations:

$\displaystyle\gamma_{\ell-m}^{I}$ $\displaystyle=(-1)^{m}\gamma_{\ell
m}^{I}\,,$ (50) $\displaystyle\gamma_{\ell-m}^{V}$
$\displaystyle=(-1)^{m+1}\gamma_{\ell m}^{V}\,,$ (51)
$\displaystyle\gamma_{\ell-m}^{Q\pm iU}$ $\displaystyle=(-1)^{m}\gamma_{\ell
m}^{Q\mp iU}\,.$ (52)

We have computed numerically the ORF multipole moments in Eqs. (48) and (49).
For $fD_{a}\gg c$ and $fD_{b}\gg c$, we reproduce the results for the
intensity and circular-polarization ORFs in Refs. anholm ; mingar13 ; gair14 ;
kato16 . For $fD_{a}=fD_{b}=10c$, we confirm the contribution of the pulsar
term to the ORFs on small angular scales mingar14 ; chu2107  and reproduce the
results for the linear-polarization ORFs chu2107 . We note that the algorithm
in the present work is less complicated and more efficient than that in Ref.
chu2107 .

##  VIII Bipolar spherical harmonics coefficients

We have derived the ORFs using the harmonic-space method. Let us go back to
the power spectrum in the timing-residual correlation function (15), which can
be also expanded in terms of bipolar spherical harmonics (BiPoSHs)
$\left\\{Y_{\ell_{1}}(\hat{e}_{a})\otimes
Y_{\ell_{2}}(\hat{e}_{b})\right\\}_{\ell m}$ as book:Varshalovich

$\langle
r(t_{a},\hat{e}_{a})r(t_{b},\hat{e}_{b})\rangle=\sum_{\ell_{1}\ell_{2}\ell
m}A^{\ell m}_{\ell_{1}\ell_{2}}\left\\{Y_{\ell_{1}}(\hat{e}_{a})\otimes
Y_{\ell_{2}}(\hat{e}_{b})\right\\}_{\ell m}\,,$ (53)

where $A^{\ell m}_{\ell_{1}\ell_{2}}$ are the expansion coefficients and

$\left\\{Y_{\ell_{1}}(\hat{e}_{a})\otimes
Y_{\ell_{2}}(\hat{e}_{b})\right\\}_{\ell
m}=\sum_{m_{1}m_{2}}(-1)^{m_{1}}\sqrt{2\ell+1}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
m&&-m_{1}&&m_{2}\end{pmatrix}Y_{\ell_{1}m_{1}}(\hat{e}_{a})Y^{*}_{\ell_{2}m_{2}}(\hat{e}_{b})\,,$
(54)

which satisfy the orthogonal condition,

$\int_{S^{2}}\mathbf{d}\hat{e}_{a}\int_{S^{2}}\mathbf{d}\hat{e}_{b}\;\left\\{Y_{\ell_{1}}(\hat{e}_{a})\otimes
Y_{\ell_{2}}(\hat{e}_{b})\right\\}_{\ell
m}\left\\{Y_{\ell^{\prime}_{1}}(\hat{e}_{a})\otimes
Y_{\ell^{\prime}_{2}}(\hat{e}_{b})\right\\}^{*}_{\ell^{\prime}m^{\prime}}=\delta_{\ell_{1}\ell^{\prime}_{1}}\delta_{\ell_{2}\ell^{\prime}_{2}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\,,$
(55)

where we have used the relation (69). The BiPoSH coefficients is thus related
to the power spectrum by

$A^{\ell
m}_{\ell_{1}\ell_{2}}=\sum_{m_{1}m_{2}}(-1)^{m_{1}}\sqrt{2\ell+1}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
m&&-m_{1}&&m_{2}\end{pmatrix}\langle
a_{\ell_{1}m_{1}}a^{*}_{\ell_{2}m_{2}}\rangle\,.$ (56)

From Eqs. (15) and (20), we have

$\langle
a_{\ell_{1}m_{1}}a^{*}_{\ell_{2}m_{2}}\rangle=\int_{-\infty}^{\infty}\frac{\mathbf{d}f}{(2\pi
f)^{2}}\,(1-e^{-2\pi ift_{a}})(1-e^{2\pi ift_{b}})\sum_{X=\\{I,V,Q\pm
iU\\}}\sum_{\ell m}X_{\ell m}(f){\tilde{\gamma}}_{\ell
m,\ell_{1}m_{1}\ell_{2}m_{2}}^{X}\,,$ (57)

where

$\gamma_{\ell
m}^{X}(fD_{a},fD_{b};\hat{e}_{a},\hat{e}_{b})=\sum_{\ell_{1}m_{1}\ell_{2}m_{2}}{\tilde{\gamma}}_{\ell
m,\ell_{1}m_{1}\ell_{2}m_{2}}^{X}Y_{\ell_{1}m_{1}}(\hat{e}_{a})Y^{*}_{\ell_{2}m_{2}}(\hat{e}_{b})\,.$
(58)

Using the ORFs in Eqs. (37) and (38), we can then explicitly write the BiPoSH
coefficients as

$\displaystyle A^{\ell m}_{\ell_{1}\ell_{2}}=$
$\displaystyle\int_{-\infty}^{\infty}\frac{\mathbf{d}f}{(2\pi
f)^{2}}\,(1-e^{-2\pi ift_{a}})(1-e^{2\pi
ift_{b}})\sqrt{\frac{(2\ell_{1}+1)(2\ell_{2}+1)}{4\pi}}J_{\ell_{1}}(fD_{a})J^{*}_{\ell_{2}}(fD_{b})\,\times$
$\displaystyle\left\\{I_{\ell
m}\left[1+(-1)^{\ell+\ell_{1}+\ell_{2}}\right]\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
0&&-2&&2\end{pmatrix}+V_{\ell
m}\left[1-(-1)^{\ell+\ell_{1}+\ell_{2}}\right]\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
0&&-2&&2\end{pmatrix}\right.$ $\displaystyle\left.+(Q+iU)_{\ell
m}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
-4&&2&&2\end{pmatrix}+(Q-iU)_{\ell
m}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
4&&-2&&-2\end{pmatrix}\right\\}\,,$ (59)

which shows that the intensity and circular-polarization anisotropies induce
even-parity and odd-parity BiPoSHs respectively, whose parity is defined by
the sum, $\ell+\ell_{1}+\ell_{2}$. This selection rule has been pointed out in
Refs. kam1 ; kam2 , though the authors considered a limiting case with
$fD_{a}\gg c$ and $fD_{b}\gg c$. Using the property (68), under the
approximation (45) we have $A^{\ell m}_{\ell_{1}\ell_{2}}=A^{\ell
m}_{\ell_{2}\ell_{1}}$. The BiPoSH coefficients induced by the linear-
polarization anisotropy are derived for the first time in this work. These
BiPoSH coefficients can be used to construct optimal estimators for testing
the statistical isotropy of the SGWB intensity and polarization kam1 ; kam2 .

##  IX Conclusion

We have studied the pulsar-timing-array observation of the Stokes parameters
of an anisotropic stochastic gravitational-wave background, based on the
spherical harmonic expansion of the pulsar timing residual. A numerical scheme
to compute the overlap reduction functions (ORFs) and the bipolar spherical
harmonics (BiPoSH) coefficients has been developed. We have used the Sachs-
Wolfe line-of-sight integral for the time residual of a Galactic millisecond
pulsar, which properly accounts for the power contribution of the pulsar term
to the ORF multipoles at small angular separation of the pulsar pair. Using a
spherical harmonic analysis, we can compute an ORF multipole with a desired
angular resolution. Our method also gives ORFs for a Galactic pulsar pair at
different distances from the Earth. Indeed, the Sachs-Wolfe line-of-sight
integral can be generalized to considering extragalactic pulsars at high
redshifts. In future pulsar-timing-array observation on a large number of
pulsars, our method provides a fast algorithm for computing accurate ORF
multipoles and BiPoSH coefficients.

###### Acknowledgements.

This work was supported in part by the Ministry of Science and Technology
(MOST) of Taiwan, Republic of China, under Grants No. MOST 110-2112-M-032 -007
(G.C.L.) and No. MOST 110-2112-M-001-036 (K.W.N.).

##  Appendix A Spin-Weighted Spherical Harmonics

The explicit form of the spin-weighted spherical harmonics that we use is

$\displaystyle{}_{s}Y_{\ell
m}(\theta,\phi)=(-1)^{m}e^{im\phi}\sqrt{\frac{(2\ell+1)}{(4\pi)}\frac{(\ell+m)!(\ell-m)!}{(\ell+s)!(\ell-s)!}}\sin^{2\ell}\\!\left(\frac{\theta}{2}\right)\sum_{r}\binom{\ell-s}{r}\binom{\ell+s}{r+s-m}(-1)^{\ell-
r-s}\cot^{2r+s-m}\\!\left(\frac{\theta}{2}\right)\,.$ (60)

When $s=0$, it reduces to the ordinary spherical harmonics,

$Y_{\ell
m}(\hat{n})=\sqrt{\frac{(2\ell+1)}{(4\pi)}\frac{(\ell-m)!}{(\ell+m)!}}P^{m}_{\ell}(\cos\theta)e^{im\phi}\,.$
(61)

Spin-weighted spherical harmonics satisfy the orthogonal relation,

$\int_{S^{2}}\mathbf{d}{\hat{n}}\;{}_{s}Y^{*}_{\ell
m}(\hat{n}){}_{s}Y_{\ell^{\prime}m^{\prime}}(\hat{n})=\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\,,$
(62)

and the completeness relation,

$\displaystyle\sum_{\ell m}{}_{s}Y^{*}_{\ell m}(\hat{n}){}_{s}Y_{\ell
m}(\hat{n}^{\prime})$ $\displaystyle=\delta(\hat{n}-\hat{n}^{\prime})$
$\displaystyle=\delta(\phi-\phi^{\prime})\delta(\cos\theta-\cos\theta^{\prime})\,.$
(63)

Its complex conjugate is

${}_{s}Y^{*}_{\ell m}(\hat{n})=(-1)^{s+m}{}_{-s}Y_{\ell-m}(\hat{n})\,,$ (64)

and its parity is given by

${}_{s}Y_{\ell m}(-\hat{n})\equiv{}_{s}Y_{\ell
m}(\pi-\theta,\phi+\pi)=(-1)^{\ell}{}_{-s}Y_{\ell m}(\hat{n})\,.$ (65)

Also, we have the spherical wave expansion:

$e^{i\vec{k}\cdot\vec{r}}=4\pi\sum_{l=0}^{\infty}\sum_{m=-\ell}^{\ell}i^{\ell}j_{\ell}(kr)Y_{\ell
m}^{*}(\hat{k})Y_{\ell m}(\hat{r})\,,$ (66)

where $j_{\ell}(x)$ is the spherical Bessel function.

We can calculate the integral of a product of three spin-weighted spherical
harmonics using the formula:

$\int_{S^{2}}\mathbf{d}\hat{e}\;{}_{s_{1}}Y_{l_{1}m_{1}}(\hat{e})\;{}_{s_{2}}\\!Y_{l_{2}m_{2}}(\hat{e})\;{}_{s_{3}}\\!Y_{l_{3}m_{3}}(\hat{e})=\sqrt{\frac{(2l_{1}+1)(2l_{2}+1)(2l_{3}+1)}{4\pi}}\begin{pmatrix}l_{1}&&l_{2}&&l_{3}\\\
-s_{1}&&-s_{2}&&-s_{3}\end{pmatrix}\begin{pmatrix}l_{1}&&l_{2}&&l_{3}\\\
m_{1}&&m_{2}&&m_{3}\end{pmatrix}\,,$ (67)

which involves two Wigner-3j symbols representing the coupling coefficients
between different spherical harmonics book:Varshalovich . The Wigner-3j symbol
is zero unless it satisfies: $l_{1}$, $l_{2}$, and $l_{3}$ have to meet the
triangular condition, i.e. $l_{1}+l_{2}\geq l_{3}\geq|l_{1}-l_{2}|$, while
$m_{1}+m_{2}+m_{3}=0$; when $m_{1}=m_{2}=m_{3}=0$, $l_{1}+l_{2}+l_{3}$ is
even. The Wigner-3j symbols have the reflection property and the summation
relation:

$\displaystyle\begin{pmatrix}l_{1}&&l_{2}&&l_{3}\\\
s_{1}&&s_{2}&&s_{3}\end{pmatrix}=(-1)^{l_{1}+l_{2}+l_{3}}\begin{pmatrix}l_{1}&&l_{2}&&l_{3}\\\
-s_{1}&&-s_{2}&&-s_{3}\end{pmatrix}=(-1)^{l_{1}+l_{2}+l_{3}}\begin{pmatrix}l_{1}&&l_{3}&&l_{2}\\\
s_{1}&&s_{3}&&s_{2}\end{pmatrix}\,,$ (68)
$\displaystyle(2\ell+1)\sum_{m_{1}m_{2}}\begin{pmatrix}\ell&&\ell_{1}&&\ell_{2}\\\
m&&-m_{1}&&m_{2}\end{pmatrix}\begin{pmatrix}\ell^{\prime}&&\ell_{1}&&\ell_{2}\\\
m^{\prime}&&-m_{1}&&m_{2}\end{pmatrix}=\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\,.$
(69)

## References

  * (1) B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), Phys. Rev. Lett. 116, 061102 (2016).
  * (2) B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), Class. Quant. Grav. 37, 055002 (2020).
  * (3) For examples, see M. A. Sedda et al., arXiv:1908.11375; V. Baibhav et al., arXiv:1908.11390; J. Baker et al., arXiv:1908.11410.
  * (4) For a review, see S. Burke-Spolaor et al., Astron. Astrophys. Rev. 27, 5 (2019).
  * (5) L. Lentati et al., Mon. Not. Roy. Astron. Soc. 453, 2576 (2015).
  * (6) Z. Arzoumanian et al. (NANOGrav), Astrophys. J. 821, 13 (2016).
  * (7) R. M. Shannon et al., Science 349, 1522 (2015).
  * (8) G. Janssen et al., Advancing Astrophysics with the Square Kilometre Array, PoS AASKA14 (2015) 037.
  * (9) A. Weltman et al., Publ. Astron. Soc. Austral. 37, e002 (2020).
  * (10) For a review, see J. D. Romano, arXiv:1909.00269.
  * (11) S. H. S. Alexander, M. E. Peskin, and M. M. Sheikh-Jabbari, Phys. Rev. Lett. 96, 081301 (2006).
  * (12) M. Satoh, S. Kanno, and J. Soda, Phys. Rev. D 77, 023526 (2008).
  * (13) L. Sorbo, J. Cosmol. Astropart. Phys. 06 (2011) 003.
  * (14) S. G. Crowder, R. Namba, V. Mandic, S. Mukohyama, and M. Peloso, Phys. Lett. B 726, 66 (2013).
  * (15) G. Cusin, R. Durrer, and P. G. Ferreira, Phys. Rev. D 99, 023534 (2020).
  * (16) N. Bartolo, D. Bertacca, S. Matarrese, M. Peloso, A. Ricciardone, A. Riotto, and G. Tasinato, Phys. Rev. D 102, 023527 (2020).
  * (17) C. Pitrou, G. Cusin, and J.-P. Uzan, Phys. Rev. D 101, 081301(R) (2020).
  * (18) J. Liu, R.-G. Cai, and Z.-K. Guo, Phys. Rev. Lett. 126, 141303 (2021).
  * (19) R.-G. Cai, Z.-K. Guo, and J. Liu, arXiv:2112.10131.
  * (20) NANOGrav Collaboration: Z. Arzoumanian et al., Astrophys. J. Lett. 905, L34 (2020).
  * (21) B. Goncharov et al., Astrophys. J. Lett. 917, L19 (2021).
  * (22) S. Chen et al., Mon. Not. Roy. Astron. Soc. 508, 4970 (2021).
  * (23) J. Antoniadis et al., arXiv:2201.03980.
  * (24) R. Abbott et al. (LIGO Scientific, Virgo, and KAGRA Collaborations), Phys. Rev. D 104, 022004 (2021).
  * (25) R. Abbott et al. (LIGO Scientific, Virgo, and KAGRA Collaborations), Phys. Rev. D 104, 022005 (2021).
  * (26) M. Anholm, S. Ballmer, J. D. E. Creighton, L. R. Price, and X. Siemens, Phys. Rev. D 79, 084030 (2009).
  * (27) C. M. F. Mingarelli, T. Sidery, I. Mandel, and A. Vecchio, Phys. Rev. D 88, 062005 (2013).
  * (28) J. Gair, J. D. Romano, S. Taylor, and C. M. F. Mingarelli, Phys. Rev. D 90, 082001 (2014).
  * (29) R. Kato and J. Soda, Phys. Rev. D 93, 062003 (2016).
  * (30) C. M. F. Mingarelli and T. Sidery, Phys. Rev. D 90, 062011 (2014).
  * (31) W. Qin, K. K. Boddy, M. Kamionkowski, and L. Dai, Phys. Rev. D 99, 063002 (2019).
  * (32) Y.-K. Chu, G.-C. Liu, and K.-W. Ng, Phys. Rev. D 103, 063528 (2021).
  * (33) Y.-K. Chu, G.-C. Liu, and K.-W. Ng, Phys. Rev. D 104, 124018 (2021).
  * (34) K.-W. Ng, arXiv:2106.12843.
  * (35) S. R. Taylor and J. R. Gair, Phys. Rev. D 88, 084001 (2013).
  * (36) S. R. Taylor, R. van Haasteren, and A. Sesana, Phys. Rev. D 102, 084039 (2020).
  * (37) S. C. Hotinli, M. Kamionkowski, and A. H. Jaffe, Open J. Astrophys. 2, 1 (2019).
  * (38) E. Belgacem and M. Kamionkowski, Phys. Rev. D 102, 023004 (2020).
  * (39) Y. Ali-Haïmoud, T. L. Smith, and C. M. F. Mingarelli, Phys. Rev. D 102, 122005 (2020).
  * (40) Y. Ali-Haïmoud, T. L. Smith, and C. M. F. Mingarelli, Phys. Rev. D 103, 042009 (2021).
  * (41) M. Born and E. Wolf, Principles of Optics, 6th ed. (Pergamon Press, New York, 1980).
  * (42) R. K. Sachs and A. M. Wolfe, Astrophys. J. 147, 73 (1967).
  * (43) R. W. Hellings and G. S. Downs, Astrophys. J. 265, L39 (1983).
  * (44) D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii, Quantum Theory of Angular Momentum (World Scientific, Singapore, 1988).


# On Searching for Minimal Integer Representation of Undirected Graphs

Victor Parque
Department of Modern Mechanical Engineering
Waseda University
3-4-1 Okubo, Shinjuku, Tokyo, Japan 169-8555
<EMAIL_ADDRESS>
& Tomoyuki Miyashita
Department of Modern Mechanical Engineering
Waseda University
3-4-1 Okubo, Shinjuku, Tokyo, Japan 169-8555
Corresponding author.

###### Abstract

Minimal and efficient graph representations are key to store, communicate, and
sample the search space of graphs and networks while meeting user-defined
criteria. In this paper, we investigate the feasibility of gradient-free
optimization heuristics based on Differential Evolution to search for minimal
integer representations of undirected graphs. The class of Differential
Evolution algorithms are population-based gradient-free optimization
heuristics having found a relevant attention in the nonconvex and nonlinear
optimization communities. Our computational experiments using eight classes of
Differential Evolution schemes and graph instances with varying degrees of
sparsity have shown the merit of attaining minimal numbers for graph
encoding/representation rendered by exploration-oriented strategies within few
function evaluations. Our results have the potential to elucidate new number-
based encoding and sample-based algorithms for graph representation, network
design and optimization.

_Keywords_ undirected graphs $\cdot$ graph representation $\cdot$ enumerative
representation $\cdot$ differential evolution $\cdot$ optimization.

##  1 Introduction

Graphs allow to compute with dependencies in several fields. Devising succinct
graph representations is essential to realize efficient mechanisms for
storing, transmission and sampling graphical structures ubiquitously. Although
matrices and lists are useful data structures allowing the simple and
straightforward representation of graphs, they are unable to meet the
information-theoretic tight bounds. In a seminal paper, Turan outlined the
concept of succinctness when dealing with graph representations, and proposed
the succinct representations of both labeled and unlabeled planar graphs using
size up to $12n$ bits ($n$ is the number of nodes in the graph)[Turán(1984)].
And Farzan and Munro proposed the representation of arbitrary graphs with the
size being a multiplicative factor $1+\epsilon$ above the information-
theoretic minimum, for small positive $\epsilon$[Farzan and Munro(2008)].
Finding efficient representation of arbitrary graphs implies finding
regularities in graphs that can be polynomially encoded/decoded. As such, the
community has used regular features in the graph topology such as
planarity[Turán(1984), Jacobson(1989), Barbay et al.(2012)Barbay, Castelli
Aleardi, He, and Munro, Aleardi et al.(2006)Aleardi, Devillers, and Schaeffer,
He et al.(1999)He, Kao, and Lu], triangularity[Barbay et al.(2012)Barbay,
Castelli Aleardi, He, and Munro, Aleardi et al.(2006)Aleardi, Devillers, and
Schaeffer, He et al.(2000)He, Kao, and Lu], separability[He et al.(2000)He,
Kao, and Lu, Blelloch and Farzan(2010), Blandford et al.(2003)Blandford,
Blelloch, and Kash] and symmetry [Katebi et al.(2012)Katebi, Sakallah, and
Markov]. When using regularity in graph representation, it is often useful
find label encodings that allow tests for node adjacency[Erdös and
Evans(1989), Gallian(2018)]. Generally speaking, when regularity in the graph
topology is known, the search space of tailored graph representations is
reduced; thus it becomes feasible to investigate tight bounds on the
representation mechanism.

The community has also explored the feasibility of compressing real-world
networks[Claude and Navarro(2010), Boldi et al.(2011)Boldi, Rosa, Santini, and
Vigna, Buehrer and Chellapilla(2008), Brisaboa et al.(2014)Brisaboa, Ladra,
and Navarro], in which compression efficiency is often reported in bits per
edge. It is also possible to realize the number-based representation of graphs
by using combinatorial ideas[Kreher and Stinson(1999), Parque et
al.(2014)Parque, Kobayashi, and Higashi, Parque and Miyashita(2017a), Parque
and Miyashita(2018a)], in which a graph coding scheme computes an integer
number representing the graph, and a decoding scheme generates the graph
topology from the integer number. Representing graphs with numbers is
advantageous when a-priori knowledge of graph regularity is unknown. Also, the
number representing the graph meets the information-theoretic tight bound over
all possible graph representations in the same class. On the other hand, when
the problem is to find an optimal graph topology, the number based
representation enables the one-dimensional (integer) search space and the
implicit parallelization of search.

One of the key questions in number-based graph representation is whether
exploring the landscape of enumerative encodings would render representations
with fewer digits, and thus smallest possible integers. Although approaches
for succinct graph encodings have relied on apriori knowledge of graph
invariants such as planarity, separability, symmetry, and sparsity, the study
of representing arbitrary graphs with the smallest possible integers has
received little attention in the literature. In this paper, we study the
feasibility of using gradient-free population-based optimization heuristics
inspired by Differential Evolution for smallest/minimal integer representation
of graphs, and present an analysis of the convergence performance. In
particular, our contributions are as follows: (1) We formulate the problem of
searching for minimal integer representation of graphs, portray examples of
the search space, and evaluate the performance of eight gradient-free
optimization algorithms with varying forms of exploration-exploitation based
on Differential Evolution schemes. (2) Our computational experiments using
graph instances with varying degrees of sparsity has shown the merit of
exploration strategies to attain better convergence with few function
evaluations.

The rest of this paper is organized as follows. Section 2 provides the
preliminaries on enumerative and minimal graph representation, section 3
describes the computational experiments, and section 4 concludes our paper.

##  2 Graphs with Minimal Integers

###  2.1 Enumerative Representation of Graphs

Let $\mathscr{G}_{n,m}$ be the class of undirected graphs with $n$ nodes and
$m$ edges. When graphs in the class $\mathscr{G}_{n,m}$ are represented by
integer numbers $g$, the following condition holds:
$g\in[0,|\mathscr{G}_{n,m}|-1]$, in which,

$|\mathscr{G}_{n,m}|=\binom{\binom{n}{2}}{m}.$ (1)

The above implies that graphs in the class $\mathscr{G}_{n,m}$ are
combinatorial structures that involves the combinations of $\binom{n}{2}$
edges taken $m$ at the time. Let the graph instance
$G=(V,E)\in\mathscr{G}_{n,m}$, with $n=|V|$ and $m=|E|$; a numerical
representation (or numerical encoding) of $G$ is the tuple
$(g,~{}CODE,~{}DECODE)$ as shown by Fig. 1:

$CODE:G\rightarrow g,$ (2)

$DECODE(CODE(g))\rightarrow G.$ (3)

Figure 1: Basic concept of representing graphs as numbers. Given an graph, the
_CODE_ function generates the number $g$ representing the graph and,
viceversa, the _DECODE_ function renders the graph from the integer number
$g$.

For a defined labeling of nodes with the tuple $(1,~{}2,~{}3,...,~{}n)$, it is
possible to define the _CODE_ function as the following algebraic
relation[Kreher and Stinson(1999), Parque et al.(2014)Parque, Kobayashi, and
Higashi]:

$g(x)=\sum_{c=1}^{m}(-1)^{m-c}\Bigg{[}\binom{i_{c}}{c}-1\Bigg{]},$ (4)

where $i_{c}$ is the numerical label (integer number) of the $c$-th edge
$(u,v)\in E$, computed as follows[Parque et al.(2014)Parque, Kobayashi, and
Higashi]:

$i_{c}=\binom{u-1}{2}+v,~{}\text{for }u>v,$ (5)

where $u\in[1,n]$ and $v\in[1,n]$ are labels of nodes from the set $V$.

The reader may note that the number representation of a graph is contingent
upon the values of $n$ and $m$, and the above-mentioned definitions consider a
defined labeling order of the set of nodes $V$. It is possible to define the
_DECODE_ function with a decoding algorithm based on the revolving door
order[Kreher and Stinson(1999), Parque et al.(2014)Parque, Kobayashi, and
Higashi].

###  2.2 Minimal Integer Representation

Graphs can be represented with minimal integers by solving the following:

$\displaystyle\underset{x}{\text{Minimize}}$ $\displaystyle L(x)$ (6) subject
to $\displaystyle x\in\mathbf{F}(n)$

where the objective function is

$L(x)=\log(g(x)),$ (7)

where $x$ denotes the representation of the ordering of labels in the graph,
and $\mathbf{F}(n)$ is the search space related to the factorial numbering
system of $n$. The problem described above is NP-hard, and without further
knowledge of the graph topology or regularities in the structure, computing
the gradient of $L$ is unfeasible; thus, the use of a gradient-free and
stochastic optimization heuristic is a desirable choice to sample the search
space of $F(x)$. The logarithmic function in Eq. (7) facilitates computing the
number of digits in number-based representations of graphs, thus it provides
an intuitive idea on the space needed to encode undirected graph. Furthermore,
the use of the factorial (or _factoradic_) space $\mathbf{F}(n)$ in Eq. (6) is
due to the permutation of node labels in terms of $n!$.

To exemplify the landscape of the above-mentioned objective function, we
consider the following examples of graphs:

(a) $G_{1}$, $n=5,m=5$ (b) $G_{2}$, $n=5,m=8$

Figure 2: Example of undirected graph topologies.

  * •

$G_{1}$: an undirected graph with $n=5$ nodes and $m=5$ edges, and

  * •

$G_{2}$: an undirected graph with $n=5$ nodes and $m=8$ edges.

Fig. 2 shows the topologies of graphs $G_{1}$ and $G_{2}$. For both graphs,
$n=5$, yet the reader may note that a node with label 5 in $G_{1}$ is not
connected to any other node. Then, one can compute the objective function $L$
over the factorial space $\mathbf{F}(n)$, for each permutation of the node
labels. For $n=5$, there exists $5!=120$ possible orders of the labels, thus
120 possible values of the metric $L$.

(a) Landscape of $L(x)$ for $G_{1}$, $n=5,m=5$ (b) Landscape of $L(x)$ for
$G_{2}$, $n=5,m=8$

Figure 3: Example of the landscapes of the fitness functions for representing
of graphs by smallest numbers. Green marks imply the lowest smallest integers.

In order to show the landscape of the objective function, Fig. 3 shows the
rendering of the metric $L$ for both graphs $G_{1}$ and $G_{2}$ when
evaluating all possible 120 permutations of labels. Here, the x-axis shows the
factorial number system corresponding to the ordering of the label (for $n=5$,
the x-axis shows numbers in the range 1-120), and the y-axis shows the metric
$L$. By observing Fig. 3, one can note that the landscape is nonsmooth, and it
is possible to find regions with smallest values of $L$. Such regions (labeled
with green color) are of interest to find node and edge labels that render
smallest/minimal integer representation. Examples of such node labels are
provided in Fig. 3-(a) and Fig. 3-(b), for each graph $G_{1}$ and $G_{2}$,
respectively. Sampling through the factorial number system allows to explore
possible configurations of graphs whose representation may lead to the use of
fewer number of digits.

(a) Landscape of $L(x)$ for arbitrary graphs with $n=5,~{}m=2,...,8$ (b)
Landscape of $L(x)$ for arbitrary graphs with $n=6,~{}m=2,...,14$

Figure 4: Examples of landscapes of the function $L$ for arbitrary graphs. (a)
The case of $n=5,m\in[2,8]$. (b) The case of $n=6$ and $m\in[2,14]$.

In order to show the overall landscape for distinct and arbitrary number of
edge configurations, Fig. 4 shows the landscapes of the metric $L$ for
arbitrary graphs when considering varying number of edges. Here, we show the
the case of undirected graphs with $n=5,~{}m\in[2,8]$, and the case of
undirected graphs with $n=6$ and number of edges $m\in[2,14]$. Compared to
Fig. 3, the landscape of in Fig. 4 shows the larger number of instances in
which it is possible to attain smallest integer representation. This is due to
the graph configurations arising from larger number of edges and nodes. Also,
sparse undirected graphs with fewer edges compared to number of nodes will
allow the larger number of combinatorial permutations with smallest integer
representation. Furthermore, although the minimization of the objective
function in Eq. (7) implies searching for smallest numbers, it is possible to
explore other objective functions to search for graphs with tailored
properties (e.g., by using factorization and the fundamental theorem of
arithmetic).

Figure 5: Graph instances used in this study.

##  3 Computational Experiments

To evaluate the performance of the proposed approach, we conducted
computational experiments using diverse graph topologies, configuration of
nodes and edges, with varying degrees of sparsity. Fig. 5 shows the set of
graph instances. For ease of reference, plots correspond to arbitrary
ordering/labeling of nodes (in the range of 10 to 20), considering the
modeling of labeled trees as well. The key motivation for using the settings
in Fig. 5 is due to our interest in finding compact representations for
multiagent communication networks, where the number of agents is expected to
be distributed in indoor environments (e.g., robots, drones, sensors
distributed in the map). As such, exploring the compact representation of
graphs is desirable to communicate and share the topology of the graph among
members of a robot swarm. The reader may also note that some graphs in Fig. 5
are not connected to the entire network. This observation (and situation) is
especially relevant when agents (nodes) are independent of the group (swarm);
thus, communication and information gathering is not essential at all times.

To evaluate the feasibility in attaining quick convergence towards graph
representations with smaller integers, we use a representative class of
optimization algorithms derived from Differential Evolution, which is a
population-based gradient-free optimization algorithm extending the difference
of vectors:

  1. 1.

DESPS: Differential Evolution with Successful Parent Selection[Guo et
al.(2015)Guo, Yang, Hsu, and Tsai]

  2. 2.

OBDE: Differential Evolution with Opposition-based Learning[Rahnamayan et
al.(2008)Rahnamayan, Tizhoosh, and Salama]

  3. 3.

JADE: Adaptive Differential Evolution External Archive[Zhang and
Sanderson(2009)]

  4. 4.

DCMAEA: Differential Covariance Matrix Adaptation[Ghosh et al.(2012)Ghosh,
Das, Roy, Minhazul Islam, and Suganthan]

  5. 5.

DERAND: DE/rand/1/bin Strategy[Storn and Price(1997)],

  6. 6.

DEBEST: DE/best/1/bin Strategy[Storn and Price(1997)],

  7. 7.

RBDE: Rank-based Differential Evolution[Sutton et al.(2007)Sutton, Lunacek,
and Whitley],

  8. 8.

DESIM: Differential Evolution with Similarity Based Mutation[Segredo et
al.(2018)Segredo, Lalla-Ruiz, and Hart].

The key motivation of using the above set is due to our interest in evaluating
distinct forms of selection pressure, parameter adaptation, exploration and
explotation mechanisms during search. As for dimensionality of the problem,
since $x\in\mathbf{F}(n)$, the dimensionality is related to $n$ due to node
labelings in the factorial search space. For simplicity and without loss of
generality, we used the factoradic representation to encode solutions in the
factorial search space. As for algorithm parameters, we used the probability
of crossover $CR=0.5$, scaling factor $F=0.7$, population sizes $POP=10$, the
bias term in RBDE $\rho=3$, and the termination criterion is $\Omega=1000$
function evaluations. In DESIM, the similarity is based on the Euclidean
distance, the learning parameters $F$ and $CR$ are adapted using JADE[Zhang
and Sanderson(2009)], and initialization uses Opposition Based
Learning[Rahnamayan et al.(2008)Rahnamayan, Tizhoosh, and Salama].

Also, due to the stochastic nature of the above mentioned metaheuristics, 10
independent runs were evaluated for each algorithm and each configuration.
Other parameters followed the suggested values of the above-mentioned
references. The fine-tuning of optimization parameters is out of the scope of
this paper. The key motivations for using the above parameters are as follows:
Crossover probability with $CR=0.5$ implements the equal importance and
consideration to historical search directions up to the current number of
iterations $t$. Small population size $NP=10$ and number of evaluations up to
1000 allow evaluating the (frontier) performance of the gradient-free
algorithms under tight computational budgets.

Computational experiments were conducted in Matlab 2018b. Since computing the
number $g$ that represents a graph involves the summation of binomial
coefficients, we rather use the logarithm of $g$ to efficiently compute the
reduction of large binomial coefficients. In this paper, the computation of
logarithmic function of the number $g$ is realized through the parallel
reductions in a GPU (NVIDIA GeForce GTX TITAN X). The reason of using a GPU is
due to the efficiency in the reduction scheme; however, since we are
constrained by the GPU hardware, we can only explore the possibility of
computing the logarithms within the bound of `realmax` in Matlab ($1.8\times
10^{308}$). As such, exploring the compact number-based representations for
large graphs is out of the scope of this paper (due to hardware limitations)
and is left for future work in our agenda (potentially through a parallel
distributed computing environment).

###  3.1 Results and Discussion

Fig. 6 shows the mean convergence performance through distinct graph instances
and independent runs. Also, the x-axis of each plot in Fig. 6 shows the number
of function evaluations, whereas the y-axis shows the value of the objective
function $L$. By observing Fig. 6, we note the following facts: (1)
convergence and stagnation in flat regions often occur after 200 function
evaluations with graphs with ten nodes and 20 edges, whereas convergence
occurs after 400 functions in other graph instances; (2) the minimal objective
values corresponding to each graph instance differ from one to another. It is
possible to observe the pattern that a large (small) number of edges implies
higher (lower) objective values. This is due to the fact that graphs encode
dependencies; as such, more information will be necessary to encode graphs
with a larger number of edges, and the integer number representing the graph
will be larger due to Eq. 4 depending on $m$ on the summation component.

Instance 1 Instance 2 Instance 3 Instance 4
---
Instance 5 Instance 6 Instance 7 Instance 8
Instance 9 Instance 10 Instance 11 Instance 12
Instance 13 Instance 14 Instance 15 Instance 16
Instance 17 Instance 18 Instance 19 Instance 20
Figure 6: Mean convergence of the objective function $L$ across independent
runs.

Among all studied graph instances, the objective value of graph instance 10-11
is lower compared to the objective values of other graph instances. Yet, for
all cases, it is possible to attain lower objective values for all cases.
Overall cases, the algorithms relying on exploration, such as DERAND, and the
algorithm relying on heavy exploration during the early stages of the search,
such as DESIM show better convergence performance. This observation is due to
the noisy landscape of the problem, such as those depicted in Fig. 3 and Fig.
4. Despite the large factorial search space, the convergence figures show that
it is possible to find potential areas of low objective values using few
function evaluations.

Figure 7: Number of cases of outperformance, equal performance and
underperformance derived from pair-wise Wilcoxon tests at 5% significance
level.

The observations above show it is possible to find convergence regions quickly
and that exploration strategies are relevant for better performance. To show
the comparison of the performance of the algorithms across graph instances, we
performed statistical test at 5% significance level (Wilcoxon). Here, our goal
is to compute instances in which an algorithm is significantly better(+),
equal(=), or worse(-) compared to another algorithm. As such, Fig. 7 shows the
overall performance comparison among algorithms based on the statistical
significance. Here, the bars in the plot show the number of times an algorithm
is significantly better (outperformance), equal (equal performance), or worse
(underperformance) compared to other algorithms across all graph instances
overall independent runs. For clarity of exposition, the counts (numbers
depicted over each bar) are computed by accumulating the number of events in
which [+], [=], and [-] occur from the Wilcoxon test. By observing Fig. 7, we
note that algorithms with an explicit exploration strategy outperform other
algorithms. For instance, DESPS, DERAND, OBDE and DESIM outperform other
algorithms in 12-14 problem instances. Also, algorithms with exploitation
strategies such as DEBEST and RBDE perform equally in most problem scenarios,
yet underperform in 13-26 instances.

The above results show that using gradient-free heuristics can enable the
search for compact integer representations of graphs using few evaluations,
which may find use in modeling and sampling graphs and networks with utmost
efficiency. The future work will aim at studying the tailored optimization
algorithms extending adaptation and explorative strategies, and the modular
schemes for large graphs[Parque and Miyashita(2018b)]. It also remains to be
answered whether different labeling strategies than the revolving door[Parque
et al.(2014)Parque, Kobayashi, and Higashi] are able to achieve smaller values
of the objective function (smaller numbers representing graphs). Furthermore,
it remains to be answered how to achieve minimal number representation for
weighted graphs, directed graphs[Parque and Miyashita(2017a)], loopy
graphs[Parque and Miyashita(2017b)], and other relevant graph architectures.
On the application side, we aim at using the insights from this study to
extrapolate towards large-scale minimal networks in the plane[Parque(2021),
Parque and Miyashita(2018c)], and the further works in network design and
optimization.

##  4 Conclusion

We have studied the feasibility of gradient-free population-based optimization
heuristics based on Differential Evolution variants for minimal integer
representations of graphs. Our computational experiments using graph instances
with varying degrees of sparsity and eight optimization algorithms with
distinct forms of adaptation and exploration-exploitation have shown the merit
of exploration-oriented strategies to attain better convergence with fewer
function evaluations. Investigating the role of adaptation, exploration and
modular strategies in large graphs, and evaluating other forms of graph
encoding strategies is left for future work. Our results have the potential to
elucidate new number-based approaches for graph representation, network design
and optimization.

## Acknowledgment

This work was supported by JSPS KAKENHI Grant Number 20K11998.

## References

  * [Turán(1984)] György Turán.  On the succinct representation of graphs.  _Discrete Applied Mathematics_ , 8(3):289–294, 1984.
  * [Farzan and Munro(2008)] Arash Farzan and J Ian Munro.  Succinct representations of arbitrary graphs.  In _Esa_ , volume 5193, pages 393–404. Springer, 2008.
  * [Jacobson(1989)] Guy Jacobson.  Space-efficient static trees and graphs.  In _30th annual symposium on foundations of computer science_ , pages 549–554. IEEE Computer Society, 1989.
  * [Barbay et al.(2012)Barbay, Castelli Aleardi, He, and Munro] Jérémy Barbay, Luca Castelli Aleardi, Meng He, and J Ian Munro.  Succinct representation of labeled graphs.  _Algorithmica_ , 62:224–257, 2012.
  * [Aleardi et al.(2006)Aleardi, Devillers, and Schaeffer] Luca Castelli Aleardi, Olivier Devillers, and Gilles Schaeffer.  Optimal succinct representations of planar maps.  In _Proceedings of the twenty-second annual symposium on Computational geometry_ , pages 309–318, 2006.
  * [He et al.(1999)He, Kao, and Lu] Xin He, Ming-Yang Kao, and Hsueh-I Lu.  Linear-time succinct encodings of planar graphs via canonical orderings.  _SIAM Journal on Discrete Mathematics_ , 12(3):317–325, 1999.
  * [He et al.(2000)He, Kao, and Lu] Xin He, Ming-Yang Kao, and Hsueh-I Lu.  A fast general methodology for information-theoretically optimal encodings of graphs.  _SIAM Journal on Computing_ , 30(3):838–846, 2000\.
  * [Blelloch and Farzan(2010)] Guy E Blelloch and Arash Farzan.  Succinct representations of separable graphs.  In _Combinatorial Pattern Matching: 21st Annual Symposium, CPM 2010, New York, NY, USA, June 21-23, 2010. Proceedings 21_ , pages 138–150. Springer, 2010.
  * [Blandford et al.(2003)Blandford, Blelloch, and Kash] Daniel K. Blandford, Guy E. Blelloch, and Ian A. Kash.  Compact representations of separable graphs.  In _Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms_ , SODA ’03, page 679–688, USA, 2003. Society for Industrial and Applied Mathematics.  ISBN 0898715385.
  * [Katebi et al.(2012)Katebi, Sakallah, and Markov] Hadi Katebi, Karem A. Sakallah, and Igor L. Markov.  Graph symmetry detection and canonical labeling: Differences and synergies.  In Andrei Voronkov, editor, _Turing-100 - The Alan Turing Centenary, Manchester, UK, June 22-25, 2012_ , volume 10 of _EPiC Series in Computing_ , pages 181–195. EasyChair, 2012.  doi: 10.29007/gzc1.
  * [Erdös and Evans(1989)] Paul Erdös and Anthony B Evans.  Representations of graphs and orthogonal latin square graphs.  _Journal of Graph Theory_ , 13(5):593–595, 1989\.
  * [Gallian(2018)] Joseph A Gallian.  A dynamic survey of graph labeling.  _Electronic Journal of combinatorics_ , 1(DynamicSurveys):DS6, 2018.
  * [Claude and Navarro(2010)] Francisco Claude and Gonzalo Navarro.  Fast and compact web graph representations.  _ACM Transactions on the Web (TWEB)_ , 4(4):1–31, 2010.
  * [Boldi et al.(2011)Boldi, Rosa, Santini, and Vigna] Paolo Boldi, Marco Rosa, Massimo Santini, and Sebastiano Vigna.  Layered label propagation: A multiresolution coordinate-free ordering for compressing social networks.  In _Proceedings of the 20th international conference on World Wide Web_ , pages 587–596, 2011.
  * [Buehrer and Chellapilla(2008)] Gregory Buehrer and Kumar Chellapilla.  A scalable pattern mining approach to web graph compression with communities.  In _Proceedings of the 2008 international conference on web search and data mining_ , pages 95–106, 2008.
  * [Brisaboa et al.(2014)Brisaboa, Ladra, and Navarro] Nieves R Brisaboa, Susana Ladra, and Gonzalo Navarro.  Compact representation of web graphs with extended functionality.  _Information Systems_ , 39:152–174, 2014.
  * [Kreher and Stinson(1999)] Donald L Kreher and Douglas R Stinson.  Combinatorial algorithms: generation, enumeration, and search.  _ACM SIGACT News_ , 30(1):33–35, 1999.
  * [Parque et al.(2014)Parque, Kobayashi, and Higashi] Victor Parque, Masakazu Kobayashi, and Masatake Higashi.  Bijections for the numeric representation of labeled graphs.  In _2014 IEEE International Conference on Systems, Man, and Cybernetics, SMC 2014, San Diego, CA, USA, October 5-8, 2014_ , pages 447–452. IEEE, 2014.
  * [Parque and Miyashita(2017a)] Victor Parque and Tomoyuki Miyashita.  On succinct representation of directed graphs.  In _2017 IEEE International Conference on Big Data and Smart Computing, BigComp 2017, Jeju Island, South Korea, February 13-16, 2017_ , pages 199–205. IEEE, 2017a.
  * [Parque and Miyashita(2018a)] Victor Parque and Tomoyuki Miyashita.  On graph representation with smallest numerical encoding.  In _2018 IEEE 42nd Annual Computer Software and Applications Conference, COMPSAC 2018, Tokyo, Japan, 23-27 July 2018, Volume 1_ , pages 817–818. IEEE Computer Society, 2018a.  doi: 10.1109/COMPSAC.2018.00135.
  * [Guo et al.(2015)Guo, Yang, Hsu, and Tsai] Shu-Mei Guo, Chin-Chang Yang, Pang-Han Hsu, and Jason S.-H. Tsai.  Improving differential evolution with a successful-parent-selecting framework.  _IEEE Transactions on Evolutionary Computation_ , 19(5):717–730, 2015.  doi: 10.1109/TEVC.2014.2375933.
  * [Rahnamayan et al.(2008)Rahnamayan, Tizhoosh, and Salama] Shahryar Rahnamayan, Hamid R Tizhoosh, and Magdy MA Salama.  Opposition-based differential evolution.  _IEEE Transactions on Evolutionary computation_ , 12(1):64–79, 2008.
  * [Zhang and Sanderson(2009)] Jingqiao Zhang and Arthur C Sanderson.  Jade: adaptive differential evolution with optional external archive.  _IEEE Transactions on evolutionary computation_ , 13(5):945–958, 2009.
  * [Ghosh et al.(2012)Ghosh, Das, Roy, Minhazul Islam, and Suganthan] Saurav Ghosh, Swagatam Das, Subhrajit Roy, S.K. Minhazul Islam, and P.N. Suganthan.  A differential covariance matrix adaptation evolutionary algorithm for real parameter optimization.  _Information Sciences_ , 182(1):199–219, 2012\.  ISSN 0020-0255.  Nature-Inspired Collective Intelligence in Theory and Practice.
  * [Storn and Price(1997)] Rainer Storn and Kenneth Price.  Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces.  _Journal of global optimization_ , 11(4):341, 1997\.
  * [Sutton et al.(2007)Sutton, Lunacek, and Whitley] Andrew M Sutton, Monte Lunacek, and L Darrell Whitley.  Differential evolution and non-separability: using selective pressure to focus search.  In _Proceedings of the 9th annual conference on Genetic and evolutionary computation_ , pages 1428–1435, 2007.
  * [Segredo et al.(2018)Segredo, Lalla-Ruiz, and Hart] Eduardo Segredo, Eduardo Lalla-Ruiz, and Emma Hart.  A novel similarity-based mutant vector generation strategy for differential evolution.  In _Proceedings of the Genetic and Evolutionary Computation Conference_ , pages 881–888, 2018.
  * [Parque and Miyashita(2018b)] Victor Parque and Tomoyuki Miyashita.  Numerical representation of modular graphs.  In _2018 IEEE 42nd Annual Computer Software and Applications Conference, COMPSAC 2018, Tokyo, Japan, 23-27 July 2018, Volume 1_ , pages 819–820. IEEE Computer Society, 2018b.  doi: 10.1109/COMPSAC.2018.00136.
  * [Parque and Miyashita(2017b)] Victor Parque and Tomoyuki Miyashita.  On the numerical representation of labeled graphs with self-loops.  In _29th IEEE International Conference on Tools with Artificial Intelligence, ICTAI 2017, Boston, MA, USA, November 6-8, 2017_ , pages 342–349. IEEE Computer Society, 2017b.  doi: 10.1109/ICTAI.2017.00061.
  * [Parque(2021)] Victor Parque.  On hybrid heuristics for steiner trees on the plane with obstacles.  In _Evolutionary Computation in Combinatorial Optimization - 21st European Conference, EvoCOP 2021, Held as Part of EvoStar 2021, Virtual Event, April 7-9, 2021, Proceedings_ , volume 12692 of _Lecture Notes in Computer Science_ , pages 120–135. Springer, 2021.  doi: 10.1007/978-3-030-72904-2\\_8.
  * [Parque and Miyashita(2018c)] Victor Parque and Tomoyuki Miyashita.  Obstacle-avoiding euclidean steiner trees by n-star bundles.  In _IEEE 30th International Conference on Tools with Artificial Intelligence, ICTAI 2018, 5-7 November 2018, Volos, Greece_ , pages 315–319. IEEE, 2018c.  doi: 10.1109/ICTAI.2018.00057.