# Conservation of Angular Momentum in the Fast Multipole Method

Oleg Korobkin, Hyun Lim, Irina Sagert, Julien Loiseau, Christopher Mauney, M.
Alexander R. Kaltenborn,
Bing-Jyun Tsao, Wesley P. Even  Los Alamos National Laboratory
Los Alamos, NM 87544 USA

###### Abstract

Smoothed particle hydrodynamics (SPH) is positioned as having ideal
conservation properties. When properly implemented, conservation of total
mass, energy, and both linear and angular momentum is guaranteed exactly, up
to machine precision. This is particularly important for some applications in
computational astrophysics, such as binary dynamics, mergers, and accretion of
compact objects (neutron stars, black holes, and white dwarfs). However, in
astrophysical applications that require the inclusion of gravity, calculating
pairwise particle interactions becomes prohibitively expensive. In the Fast
Multipole Method (FMM), they are, therefore, replaced with symmetric
interactions between distant clusters of particles (contained in the tree
nodes) [1]. Although such an algorithm is linear momentum-conserving, it
introduces spurious torques that violate conservation of angular momentum.

We present a modification of FMM that is free of spurious torques and
conserves angular momentum explicitly. The new method has practically no
computational overhead compared to the standard FMM.

##  I Introduction

To correctly model astrophysical phenomena like the orbital motion of planets,
accretion disks, rotating neutron stars and black holes, compact star mergers,
and galactic disks, numerical approaches have to accurately determine
gravitational interactions and ensure the conservation of energy, linear, and
angular momentum. Many self-gravitating astrophysical systems are well-
described by Newtonian gravity, and there is a variety of numerical methods to
model the corresponding gravitational field. For example, grid-based methods
often determine the gravitational potential by solving the Poisson equation
and using iterative solvers and/or Fourier transforms [2, 3, 4, 5]. Particle-
based methods, like N-body and Smoothed Particle Hydrodynamics (SPH) [6, 7, 8,
9, 10, 11, 12] usually rely on tree codes [13], particle-mesh methods [14], or
the Fast Multipole Method (FMM) [15, 16, 17]. However, when it comes to
angular momentum conservation, some methods, including FMM, struggle. The
latter conserves angular momentum exactly when used in zeroth-order or with a
small error for low values of the so-called multipole acceptance criterion
(MAC) angle. However, this comes with a high computational cost. To reduce
this performance cost, a larger MAC angle is typically chosen, which then
correspondingly increases the error in angular momentum.

In this work, we investigate a propositions within the standard FMM to
conserve the angular momentum _by construction_ , on par with linear momentum.
The spurious torques in standard FMM arise due to the pairwise action-reaction
forces being slightly offset from the line connecting two particles. Here we
propose a new method, in which the misaligned pairwise forces are reprojected
onto the line connecting two particles, such that spurious torques disappear.

The rest of the paper is organized as follows: in Section II, we describe FMM
[1] and the specifics of the problem conserving angular momentum for high-
order FMM. We start by reviewing the symmetric multipole interactions as
introduced in [17, 1]. In Section III, we describe two proposed approach,
which ensures angular momentum conservation. In Section IV, we present the
numerical tests, highlighting conservation of linear and angular momentum; we
conclude with the Section V.

##  II Method

In this section, we briefly outline the standard FMM approach using multipole
expansions in Cartesian coordinates. We employ tensor index notation with
Latin indices ($i,j,k,...$) from the middle of the alphabet, and the unit
metric $g_{ij}\equiv\delta_{ij}$. For such metric, covariant and contravariant
components are the same and there is no need to distinguish lower and upper
indices. The Latin indices from the beginning of the alphabet ($a$, $b$,
$\dots$) will be used to enumerate particles. We also adopt the Einstein
summation convention (repeated indices indicate a sum over that index,
$\sum_{i}a^{i}b_{i}=a^{i}b_{i}$).

Computing gravitational forces for $N$ particles via pairwise interactions
results in an algorithm that scales as $O(N^{2})$ and quickly becomes
computationally unfeasible for growing particle numbers. FMM offers a solution
to this problem by reducing this complexity to almost an $O(N)$ level [16,
18]. This is achieved by replacing pairwise particle interactions with
interactions between distant nodes in a particle tree, where a cluster of
distant particles contained in a node is represented with a single point mass.
With that, pairwise particle interactions can be replaced with long-range
interactions between nodes, provided that the latter are well-separated. A
gravitational potential of a remote node is approximated by a potential of a
single point with gravitational monopole (mass), quadrupole, octopole, and
higher moments up to desired precision.

The remarkable efficiency of FMM is complemented by the property of exact
momentum conservation. A gravitational field of the remote node $B$ acting on
a particle $a$ inside the node $A$ is computed by the Taylor expansion of the
gravitational potential from the center of mass (COM) of $A$ to the position
of $a$ (see Figure 1). For conservation to work, the truncation order of the
Taylor expansion must match the order of multipole moments retained in FMM.

In the following, we follow the presentation by Dehnen[1]. Consider two nodes,
$A$ and $B$, with their respective COMs at ${\bf Z}_{A}$ and ${\bf Z}_{B}$,
and let ${\bf R}:={\bf Z}_{A}-{\bf Z}_{B}$. Let ${\bf X}_{a}$ and ${\bf
Y}_{b}$ be the locations of particles $a$ and $b$ within the nodes $A$ and
$B$, with masses $m_{a}$ and $m_{b}$, respectively. The gravitational
potential $\varphi_{a\leftarrow B}$ at ${\bf X}_{a}$ created by the node $B$
is due to the sum of individual interactions with other particles in $B$:

$\displaystyle\varphi_{a\leftarrow B}=-\sum_{b\in B}G\frac{m_{b}}{|{\bf
X}_{a}-{\bf Y}_{b}|},$ (1)

where $G$ is the Newtonian gravitational constant. This is an exact
expression, however, it scales as $O(N^{2})$.

In FMM, we rewrite the expression ${\bf X}_{a}-{\bf Y}_{b}={\bf R}+({\bf
x}_{a}-{\bf y}_{b})$ with ${\bf x}_{a}:={\bf X}_{a}-{\bf Z}_{A}$ and ${\bf
y}_{b}:={\bf Y}_{b}-{\bf Z}_{B}$ (see Figure 1) and expand the potential at
${\bf X}_{a}$ in Taylor series around the COM of node $A$:

$\displaystyle\varphi_{a\leftarrow b}$ $\displaystyle=\varphi(|{\bf
X}_{a}-{\bf Y}_{b}|)=$ $\displaystyle=\varphi(R)+({\bf x}_{a}-{\bf
y}_{b})\cdot\nabla\varphi(|\bf{r}|)_{\bf{r}={\bf R}}+\dots=$
$\displaystyle=\sum_{p}\frac{1}{p!}\left[[({\bf x}_{a}-{\bf
y}_{b})\cdot\nabla]^{p}\varphi(|\bf{r}|)\right]_{\bf{r}={\bf R}}.$ (2)

Figure 1: Illustration of the FMM method. Depicted are two “well-separated”
nodes $A$ and $B$, each containing a cluster of particles. Gravitational field
of all the particles inside $B$ is approximated by a field of a point mass
$M_{B}$ with optional quadrupole, octopole etc. The field at $\vec{x}_{a}$ is
obtained by using the truncated Taylor expansion of the field from the center
of mass of $A$. The angle $\theta$ shows the MAC angle.

The series converges if $|{\bf x}_{a}-{\bf y}_{b}|$ is small: ${|{\bf
x}_{a}-{\bf y}_{b}|<R}$. This condition is always satisfied if the nodes are
“well-separated”, namely when their size is smaller than the distance between
their COMs:

$\displaystyle\frac{r_{A,{\rm max}}+r_{B,{\rm max}}}{R}\leq\tan\theta,$ (3)

where $r_{A,\textrm{max}}$ and $r_{B,\textrm{max}}$ are the maximum distances
from particles within the nodes to their respective COMs, and $\theta$ is the
MAC angle. The latter controls whether the tree nodes are suitable to use
multipoles. The series converges if $0\leq\tan\theta<1$ (see Fig. 1).

We can truncate the series at a finite $p$ to obtain an approximation for the
potential to arbitrary precision. The symmetry of the potential with respect
to flipping the particles ${a\leftrightarrow b}$ guarantees that the forces
will be antisymmetric (equal in magnitude and opposite in direction), thus
producing exact conservation of linear momentum.

When summed over all particles in node $B$, the potential can be expressed in
terms of the multipoles of $B$:

$\displaystyle\varphi_{a\leftarrow B}$ $\displaystyle=\sum_{b\in
B}\varphi_{a\leftarrow b}=$ $\displaystyle=\sum_{b\in B}\varphi(R)+\sum_{b\in
B}({\bf x}_{a}-{\bf y}_{b})\cdot\nabla\varphi(|\bf{r}|)_{\bf{r}={\bf
R}}+\dots=$
$\displaystyle=-\frac{GM_{B}}{R}-G\Bigg{\\{}\frac{1}{2}Q^{(B)}_{ij}\frac{r^{i}r^{j}}{R^{5}}+\frac{1}{6}H^{(B)}_{ijk}\frac{r^{i}r^{j}r^{k}}{R^{7}}+$
$\displaystyle\qquad\qquad\qquad+\frac{1}{24}X^{(B)}_{ijkl}\frac{r^{i}r^{j}r^{k}r^{l}}{R^{9}}+\dots\Bigg{\\}},$
(4)

where $M_{B}$ is the node mass, $Q^{(B)}_{ij}$, $H^{(B)}_{ijk}$, and
$X^{(B)}_{ijkl}$ are its quadrupole, octopole, and hexadecapole moments,
respectively. Here, the components $\\{r^{i}\\}$ represent vector ${\bf
x}_{a}$, with ${r\equiv|{\bf x}_{a}|}$ being its magnitude.

So, a field created at the location of a particle $a$ by a node $B$ is given
by a relatively simple and fast-to-evaluate expression (4), representing a
truncated Taylor series. We can, therefore, compute gravitational forces
between the nodes $A$ and $B$ in two steps: (i) compute multipole moments of
the nodes; (ii) apply Taylor expansion from their respective COM to the
particle locations. Summing the symmetric expression (2) over the particles in
the node $B$ gives the Taylor expansion from COM of node $A$, and vice versa.
As long as both Taylor expansions are obtained from the same symmetric
expression, the total gravitational forces cancel out exactly.

However, this method also introduces torques generated between pairs of
particles, which in the case of exact $O(N^{2})$ method vanish identically.
Indeed, even though the forces are equal in magnitude and opposite in
direction, they are not necessarily aligned with the line connecting particles
$a$ and $b$, producing small nonzero torques. The existence of these torques
leads to violation of the angular momentum conservation. A possible solution
was found by Marcello [19] by introducing artificial compensating torques to
counteract the spurious FMM torques between pairs of nodes.
Let us illustrate linear momentum conservation versus angular momentum
nonconservation in FMM using the Taylor expansion up to 4th order.

###  II-A Example: expansion up to 4-th order

The monopole $M$ and the quadrupole moments $Q_{ij}$ are given by

$\displaystyle M$ $\displaystyle=\sum_{b}m_{b},$ (5) $\displaystyle Q^{ij}$
$\displaystyle=\sum_{b}m_{b}(3y_{b}^{i}y_{b}^{j}-\delta_{ij}q_{b}^{2}),$ (6)

where $y_{b}^{i}=Y^{i}_{b}-Z_{B}^{i}$ (or, in vector notation, ${\bf
y}_{b}:={\bf Y}_{b}-{\bf Z}_{B}$), and $q_{b}:=|{\bf y}_{b}|$.

The octopole moment $H_{ijk}$ is given by

$\displaystyle
H^{ijk}=\sum_{b}m_{b}\left[15y_{b}^{i}y_{b}^{j}y_{b}^{k}-3q_{b}^{2}(\delta_{ij}y_{b}^{k}+\delta_{jk}y_{b}^{i}+\delta_{ik}y_{b}^{j})\right].$
(7)

The hexadecapole moment $X_{ijkl}$ is given by

$\displaystyle X^{ijkl}$
$\displaystyle=\sum_{b}m_{b}\big{[}105y_{b}^{i}y_{b}^{j}y_{b}^{k}y_{b}^{l}-15q_{b}^{2}(\delta_{ij}y_{b}^{k}y_{b}^{l}+\delta_{il}y_{b}^{j}y_{b}^{k}$
$\displaystyle\qquad\qquad+\delta_{ik}y_{b}^{j}y_{b}^{l}+\delta_{jl}y_{b}^{i}y_{b}^{k}+\delta_{jk}y_{b}^{i}y_{b}^{l}+\delta_{lk}y_{b}^{i}y_{b}^{j})$
$\displaystyle\qquad\qquad+3q_{b}^{4}(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\big{]}.$
(8)

The gradient and higher partial derivatives of the gravitational potential
$\varphi$, $\nabla^{p}\varphi(R)=-\nabla^{p}\frac{1}{R}$, are [19]:

$\displaystyle D$ $\displaystyle=-\frac{1}{R}\;,$ (9) $\displaystyle D_{i}$
$\displaystyle=\frac{R^{i}}{R^{3}}\;,$ (10) $\displaystyle D_{ij}$
$\displaystyle=\frac{-3R^{i}R^{j}+\delta_{ij}R^{2}}{R^{5}}\;,$ (11)
$\displaystyle D_{ijk}$
$\displaystyle=\frac{15R^{i}R^{j}R^{k}-3R^{2}(\delta_{ij}R^{k}+\delta_{jk}R^{i}+\delta_{ik}R^{j})}{R^{7}}\;,$
(12) $\displaystyle D_{ijkl}$
$\displaystyle=\frac{1}{R^{9}}\big{[}-105R^{i}R^{j}R^{k}R^{l}+$
$\displaystyle\qquad+15R^{2}(\delta_{ij}R^{k}R^{l}+\delta_{il}R^{j}R^{k}+\delta_{ik}R^{j}R^{l}+$
$\displaystyle\qquad\qquad+\delta_{jl}R^{i}R^{k}+\delta_{jk}R^{i}R^{l}+\delta_{lk}R^{i}R^{j})$
$\displaystyle\qquad-3R^{4}(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\big{]}.$
(13)

The 4-th order truncated Taylor expansion of the gravitational potential at
${\bf X}_{a}$ generated by the node $B$ gives:

$\displaystyle\varphi_{a\leftarrow B}({\bf X}_{a})\approx
M_{B}D+D_{i}M_{B}x^{i}+\frac{1}{2}D_{ij}(Q^{ij}_{B}+M_{B}x^{i}x^{j})$
$\displaystyle+\frac{1}{6}D_{ijk}(H^{ijk}_{B}+3Q^{ij}_{B}x^{k}+M_{B}x^{i}x^{j}x^{k})$
$\displaystyle+\frac{1}{24}D_{ijkl}(X^{ijkl}_{B}+4H^{ijk}_{B}x^{l}+6Q^{ij}_{B}x^{k}x^{l}+M_{B}x^{i}x^{j}x^{k}x^{l}).$
(14)

Using Eqn. 14, the gravitational acceleration (specific force) exerted by the
node $B$ at a point ${\bf X}_{a}$ in the node $A$ can be computed as follows
(using $\vec{a}=-\nabla\varphi$):

$\displaystyle a^{i}_{a\leftarrow B}({\bf X}_{a})$ $\displaystyle\approx-
M_{B}D_{i}-D_{ij}M_{B}x^{j}$
$\displaystyle-\frac{1}{2}D_{ijk}(Q^{jk}_{B}+M_{B}x^{j}x^{k})$
$\displaystyle-\frac{1}{6}D_{ijkl}(H^{ijk}_{B}+3Q^{jk}_{B}x^{l}+M_{B}x^{j}x^{k}x^{l}).$
(15)

It is straightforward to check that the net force between the nodes $A$ and
$B$ vanishes identically:

$\displaystyle{\bf F}_{B\leftarrow A}+{\bf F}_{A\leftarrow B}=\sum_{a\in
A}m_{a}{\bf a}_{a\leftarrow B}+\sum_{b\in B}m_{b}{\bf a}_{b\leftarrow A}=0.$
(16)

Zero net force implies exact conservation for linear momentum. For the angular
momentum, on the other hand, this result does not hold. The net torque
$\mathbf{\tau}_{AB}$ taken at the COM of the node $B$ due to the mutual
gravitation of $A$ and $B$ is:

$\displaystyle\tau_{AB}^{i}=\sum_{a\in
A}m_{a}\epsilon_{ijk}(R^{j}+x_{a}^{j})a^{k}_{a\leftarrow B}+\sum_{b\in
B}m_{b}\epsilon_{ijk}y_{b}^{j}a^{k}_{b\leftarrow A},$ (17)

where $\epsilon_{ijk}$ is the Levi-Civita symbol. Using the previous
expressions for acceleration, we can obtain the net nonzero torque up to 4-th
order (check also [19] for a simplified result):

$\displaystyle\tau^{i}_{AB}$
$\displaystyle\approx\epsilon_{ijk}D_{kl}(M_{A}Q_{B}^{jl}-M_{B}Q_{A}^{jl})$
$\displaystyle+\frac{1}{2}\epsilon_{ijk}D_{klm}(M_{A}H_{B}^{jlm}-M_{B}H_{A}^{jlm})$
$\displaystyle+\frac{1}{6}\epsilon_{ijk}D_{klmp}\Big{[}3(Q_{A}^{mp}Q_{B}^{jl}-Q_{B}^{mp}Q_{A}^{jl})$
$\displaystyle\quad+M_{A}X_{B}^{jlmp}-M_{B}X_{A}^{jlmp}\Big{]}$
$\displaystyle-\frac{1}{2}\epsilon_{ijk}D_{klm}(M_{A}Q_{B}^{lm}+M_{B}Q_{A}^{lm})R^{j}$
$\displaystyle-\frac{1}{6}\epsilon_{ijk}D_{klmp}(M_{A}H_{B}^{lmp}+M_{B}H_{A}^{lmp})R^{j}.$
(18)

##  III Improving the Method

We consider two different approaches to conserve angular momentum in FMM.

###  III-A Zeroth Order FMM

One of the simplest solutions to conserve angular momentum in FMM is to
refrain from using FMM orders higher than zeroth. Indeed, when only the
constant term is left in the expression for acceleration (15), there is no net
torque. Every particle experiences the same acceleration from the other node,
both in magnitude and in direction. As a result, the net spin torque on both
nodes is zero. The torque between COMs of the nodes is also zero, because the
action and reaction are exactly aligned along the line connecting the COMs.
This can also be inferred from the expression (II-A) where there is no
monopole contribution in the torques.

The problem with low-order expansion is that it lacks accuracy. The accuracy
can be improved by decreasing MAC, but this requires compromising the
efficiency. On the other hand, simplicity of implementation makes this
approach ideal for a class of problems where angular momentum conservation may
be more important than the accuracy of interactions.

###  III-B Realigning pairwise forces

In the exact Newtonian gravity of $N$ interacting particles, conservation of
angular momentum comes not only from the balance of action and reaction in
pairwise interactions but also from the fact that they are collinear with the
line connecting the particles. For the same reason, in FMM non-conservation of
angular momentum ultimately comes from the misalignment of the pairwise forces
between particles on different nodes, as approximated by the truncated Taylor
series (2). One way to restore it, therefore, would be to re-align these
forces, for instance by projecting them onto the unit vector ${\bf n}_{ab}$ in
the direction from $a$ to $b$:

$\displaystyle{\bf a}^{\|}_{a\leftarrow b}=({\bf a}_{a\leftarrow b}\cdot{\bf
n}_{ab}){\bf n}_{ab}.$ (19)

Another possibility is to simply multiply the unit vector ${\bf n}_{ab}$ by
the acceleration magnitude:

$\displaystyle{\bf a}^{(0),\|}_{a\leftarrow b}=-|{\bf a}^{(0)}_{a\leftarrow
b}|{\bf n}_{ab}.$ (20)

The unit vector ${\bf n}_{ab}$ in (19) and (20) can be expressed in terms of
${{\bf n}_{R}:={\bf R}/R}$, the unit vector along the line connecting COMs of
the nodes:

$\displaystyle{\bf n}_{ab}=\frac{{\bf X}_{a}-{\bf Y}_{b}}{|{\bf X}_{a}-{\bf
Y}_{b}|}=\frac{{\bf n}_{R}+{\boldsymbol{\delta}}_{ab}}{|{\bf
n}_{R}+{\boldsymbol{\delta}}_{ab}|},$ (21)

where ${\boldsymbol{\delta}}_{ab}:=({\bf x}_{a}-{\bf y}_{b})/R$ must be small
such that ${\bf n}_{ab}$ can be expanded around ${\bf n}_{R}$ in powers of
${\boldsymbol{\delta}}_{ab}\propto{\bf x}_{a}-{\bf y}_{b}$, in the spirit of
FMM:

$\displaystyle{\bf n}_{ab}\approx({\bf
n}_{R}+{\boldsymbol{\delta}}_{ab})(1-{\bf
n}_{R}\cdot{\boldsymbol{\delta}}_{ab}+\dots)=$ $\displaystyle\qquad={\bf
n}_{R}+{\boldsymbol{\delta}}_{ab}-({\bf
n}_{R}\cdot{\boldsymbol{\delta}}_{ab}){\bf n}_{R}-({\bf
n}_{R}\cdot{\boldsymbol{\delta}}_{ab}){\boldsymbol{\delta}}_{ab},$ (22)
$\displaystyle({\bf a}_{a\leftarrow b}\cdot{\bf n}_{ab}){\bf
n}_{ab}=\frac{({\bf a}_{a\leftarrow b},{\bf
n}_{R}+{\boldsymbol{\delta}}_{ab})}{|{\bf
n}_{R}+{\boldsymbol{\delta}}_{ab}|^{2}}({\bf
n}_{R}+{\boldsymbol{\delta}}_{ab})\approx$ $\displaystyle\approx\bigg{[}({\bf
a}_{a\leftarrow b}\cdot{\bf n}_{R})+({\bf a}_{a\leftarrow
b}\cdot{\boldsymbol{\delta}}_{ab})$ $\displaystyle\qquad-2({\bf
a}_{a\leftarrow b}\cdot{\bf n}_{R})({\boldsymbol{\delta}}_{ab}\cdot{\bf
n}_{R})+\dots\bigg{]}({\bf n}_{R}+{\boldsymbol{\delta}}_{ab}).$ (23)

Both expressions (22) and (23) have the form of a scalar quantity times $({\bf
n}_{R}+{\boldsymbol{\delta}}_{ab})$, which is collinear with ${\bf n}_{ab}$.
In these scalar quantities, we only keep terms up to the first order, but
higher terms can be added for better accuracy. The overall expressions (22)
and (23) are quadratic in ${\boldsymbol{\delta}}_{ab}$, so, when summed over
all particles $b\in B$, they produce terms up to a quadrupole.

At zeroth order, the net acceleration on a particle $a$ from the node $B$
according to equations (20) and (22) is:

$\displaystyle{\bf a}^{(0),\|}_{a\leftarrow B}$ $\displaystyle=\sum_{b\in
B}{\bf a}^{(0),\|}_{a\leftarrow b}=-|\nabla\varphi^{(0)}(R)|\sum_{b\in
B}m_{b}{\bf n}_{ab}$ $\displaystyle\approx-\frac{G}{R^{2}}\sum_{b\in
B}m_{b}\bigg{[}{\bf n}_{R}+{\boldsymbol{\delta}}_{ab}-$ $\displaystyle-({\bf
n}_{R}\cdot{\boldsymbol{\delta}}_{ab}){\bf n}_{R}-({\bf
n}_{R}\cdot{\boldsymbol{\delta}}_{ab}){\boldsymbol{\delta}}_{ab}\bigg{]}=$
$\displaystyle=-\frac{GM_{B}}{R^{2}}\bigg{[}{\bf n}_{R}+\frac{1}{R}{\bf
x}_{a}^{\perp}-$ $\displaystyle\quad-\frac{1}{R^{2}}({\bf x}_{a}\cdot{\bf
n}_{R}){\bf x}_{a}-\frac{1}{M_{B}R^{2}}{\bf
n}_{R}\cdot\bar{\bar{Q}}_{B}\bigg{]},$ (24)

where ${\bf x}_{a}^{\perp}:={\bf x}_{a}-{\bf n}_{R}({\bf n}_{R}\cdot{\bf
x}_{a})$ is the component of ${\bf x}_{a}$ orthogonal to ${\bf n}_{R}$. Since
we made an approximation (22), the length of the vector on which we are
projecting is incorrect. However, because its direction is along the line
connecting the two particles, the artificial torque automatically vanishes,
and angular momentum conservation is achieved. Moreover, because the projected
pairwise force is antisymmetric in the particles $a$ and $b$, the linear
momentum conservation is also respected.

The first three terms inside the square brackets in (24) represent a force
reoriented towards the COM of $B$, and the last term is a quadrupole
interaction:

$\displaystyle{\bf a}^{(0),\|}_{a\leftarrow B}\approx-\frac{GM_{B}}{R^{2}}{\bf
n}_{aB}+\frac{G}{R^{4}}\bar{\bar{Q}}_{B}\cdot{\bf n}_{R}.$ (25)

|
---|---
Figure 2: Left panel: setup for a simple test. Two shapes are placed at a
varying distance $d$ between their centers of mass. The shapes consist of
$4\times 10^{3}$ particles each. Right panel: relative error between the exact
Newtonian force and approximate force, as a function of distance, computed
using the 0th and 1st order FMM, with and without realignment. The error is
measured for every particle and averaged over the shapes using the $L_{2}$
norm.

However, as discussed above, zeroth-order FMM already conserves angular
momentum. Non-conservation issues appear at the first order and higher. With
the notation ${\boldsymbol{\xi}_{a}}:={\bf x}_{a}/R$,
${\boldsymbol{\eta}_{b}}:={\bf y}_{b}/R$, and the matrix
$\bar{\bar{D}}:=\|D_{ij}\|$ from Eq.(11), the first-order correction for the
acceleration of a particle $a$ due to particle $b$ is:

$\displaystyle{\bf a}^{(1)}_{a\leftarrow
b}=-m_{b}R\bar{\bar{D}}\cdot{\boldsymbol{\xi}_{a}}.$ (26)

In the subsequent derivation, it is convenient to introduce a “scalar product”
based on the matrix $\bar{\bar{D}}$:

$\displaystyle\langle\boldsymbol{\xi},\boldsymbol{\eta}\rangle_{D}:=RD_{ij}\xi^{i}\eta^{j}.$
(27)

To get the net realigned first-order correction, we can use (19) and (23):

$\displaystyle{\bf a}^{(0),\|}_{a\leftarrow B}+{\bf a}^{(1),\|}_{a\leftarrow
B}=\sum_{b\in B}({\bf a}^{(0)}_{a\leftarrow b}+{\bf a}^{(1)}_{a\leftarrow
b}\cdot{\bf n}_{ab}){\bf n}_{ab}\approx$ $\displaystyle\approx-
M_{b}\big{[}D^{2}(1-{\boldsymbol{\xi}_{a}}\cdot{\bf
n}_{R})+\langle{\boldsymbol{\xi}_{a}},{\bf n}_{R}\rangle_{D}\big{]}({\bf
n}_{R}+{\boldsymbol{\xi}_{a}})$ $\displaystyle\quad-\langle{\bf
n}_{R},\bar{\bar{Q}}\rangle_{D}D^{2}+({\bf
n}_{R}\cdot\bar{\bar{Q}}_{B})D^{4}.$ (28)

The equation (28) gives an expression for the net acceleration, accurate up to
$O(\delta^{2})$ and conserving both linear and angular momenta.

##  IV Results

A simple test of the new method is presented in Figure 2. In the test, we
explored the error in computing the Newtonian force acting on particles
arranged in two irregular shapes. Figure 2 shows the error as a function of
distance between the shapes. The distance is shown in the units of
characteristic shape size (3 cm). Red and black curves display errors for
zeroth and first-order FMM, with (dashed) and without (solid) realignment. It
demonstrates that the accuracy of the FMM method is not affected. This is
expected, since the expansions (24) and (28) in our realignment method are
carried out to second order in $\delta$, $O(\delta^{2})$, higher than the FMM
orders to which it is applied. Most importantly, both the net force and the
net torque on the particles vanish at machine precision.

The zeroth order FMM method has been implemented in FleCSPH [20].
FleCSPH111https://github.com/laristra/flecsph is a smoothed particle
hydrodynamics simulation tool based on FleCSI [21].
FleCSI222https://github.com/laristra/flecsi is a compile-time configurable
framework providing distributed and parallel topologies, such as structured
and unstructured mesh, Narray, and Ntree. FleCSPH is based on the MPI
implementation of the Ntree topology implementing distributed binary, quad,
and octree in 1, 2, and 3 dimensions respectively.

To demonstrate the implementation, we test evolution of a stable isolated star
in equilibrium. This test checks consistency and conservation properties for
the coupled hydrodynamics and gravity.

Figure 3: Evolution of angular momentum for a self-gravitating model of a
white dwarf with two different MAC values: $\tan\theta_{\textrm{MAC}}=0.3$,
and $0.5$.

For the initial data, we solve Lane-Emden equation [22] for polytrope with
$\Gamma=5/3$, $K=10^{12}$ (in CGS units), and central density
$\rho_{c}=5.2\times 10^{6}\,{\rm g}\,{\rm cm}^{-3}$. This results in a
polytrope resembling a white dwarf, with mass $0.2\,M_{\odot}$ and radius 4790
km.

Fig. 3 shows the evolution of angular momentum for a self-gravitating model of
a white dwarf, with zeroth-order FMM method. Two different MAC angles,
$\tan\theta_{\textrm{MAC}}=0.3$ and $0.5$, were chosen and compared with exact
$N$-body scheme. Both methods conserve angular momentum up to machine
precision, as expected.

##  V Conclusion and Discussion

We present a new modification to the FMM that allows to conserve angular
momentum by construction. The main idea of the method is to realign the
approximate pairwise forces in such a way that they are parallel to the line
connecting two interacting particles. The unit vector along this direction can
be expanded in terms of the difference between the local particle position
vectors relative to the COMs of their nodes,
${{\boldsymbol{\delta}}\propto({\bf x}_{a}-{\bf y}_{b})}$, and then resummed
to arbitrary multipole order, while retaining the momentum conservation
property in the spirit of the standard FMM. The method is demonstrated for
zeroth and first orders with a simple setup of gravitational interaction
between two irregular shapes. We also show an excellent conservation of
angular momentum in the evolution of a single star with zeroth order FMM (see
Figure 3).

The method can be extended to higher multipole moments. The underlying algebra
of multipole operators is that of symmetric tensors, which can be compactly
represented. Software that manipulates and translates symbolic expressions,
such as SymPy [23], can greatly simplify the generation of code to evaluate
multipole operators. We will therefore seek opportunities to automate our
method before advancing to higher multipoles case in the future work.

## Acknowledgment

The authors would like to thank Christopher Fryer and Stephan Rosswog for
valuable feedback that helped improve the paper. This work is supported by the
LANL ASC Program and LDRD grants 20200145ER and 20190021DR. This work used
resources provided by the LANL Institutional Computing Program. LANL is
operated by Triad National Security, LLC, for the National Nuclear Security
Administration of the U.S.DOE (Contract No. 89233218CNA000001). This work is
authorized for unlimited release under LA-UR-21-24198.

## References

  * [1] W. Dehnen, “A Very Fast and Momentum-conserving Tree Code,” _ApJL_ , vol. 536, pp. L39–L42, Jun. 2000.
  * [2] B. Fryxell, K. Olson, P. Ricker, F. X. Timmes, M. Zingale, D. Q. Lamb, P. MacNeice, R. Rosner, J. W. Truran, and H. Tufo, “FLASH: An adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes,” _The Astrophysical Journal Supplement Series_ , vol. 131, no. 1, pp. 273–334, nov 2000. [Online]. Available: https://doi.org/10.1086/317361
  * [3] M. C. R. D’Souza, P. M. Motl, J. E. Tohline, and J. Frank, “Numerical simulations of the onset and stability of dynamical mass transfer in binaries,” _The Astrophysical Journal_ , vol. 643, no. 1, pp. 381–401, may 2006. [Online]. Available: https://doi.org/10.1086/500384
  * [4] J. M. Stone, T. A. Gardiner, P. Teuben, J. F. Hawley, and J. B. Simon, “Athena: A new code for astrophysical MHD,” _The Astrophysical Journal Supplement Series_ , vol. 178, no. 1, pp. 137–177, sep 2008. [Online]. Available: https://doi.org/10.1086/588755
  * [5] A. S. Almgren, V. E. Beckner, J. B. Bell, M. S. Day, L. H. Howell, C. C. Joggerst, M. J. Lijewski, A. Nonaka, M. Singer, and M. Zingale, “CASTRO: A NEW COMPRESSIBLE ASTROPHYSICAL SOLVER. i. HYDRODYNAMICS AND SELF-GRAVITY,” _The Astrophysical Journal_ , vol. 715, no. 2, pp. 1221–1238, may 2010. [Online]. Available: https://doi.org/10.1088/0004-637x/715/2/1221
  * [6] L. Hernquist and N. Katz, “TREESPH: A Unification of SPH with the Hierarchical Tree Method,” _Astrophysical Journal Supplement_ , vol. 70, p. 419, Jun. 1989.
  * [7] V. Springel, “The cosmological simulation code gadget-2,” _Monthly Notices of the Royal Astronomical Society_ , vol. 364, no. 4, pp. 1105–1134, 12 2005. [Online]. Available: https://doi.org/10.1111/j.1365-2966.2005.09655.x
  * [8] S. Vanaverbeke, R. Keppens, S. Poedts, and H. Boffin, “Gradsph: A parallel smoothed particle hydrodynamics code for self-gravitating astrophysical fluid dynamics,” _Computer Physics Communications_ , vol. 180, no. 7, pp. 1164–1182, 2009. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0010465509000071
  * [9] R. Yokota and L. Barba, “Hierarchical n-body simulations with autotuning for heterogeneous systems,” _Computing in Science Engineering_ , vol. 14, no. 3, pp. 30–39, 2012.
  * [10] P. Lorén-Aguilar, J. Isern, and E. García-Berro, “Smoothed particle hydrodynamics simulations of white dwarf collisions and close encounters,” _Monthly Notices of the Royal Astronomical Society_ , vol. 406, no. 4, pp. 2749–2763, 08 2010. [Online]. Available: https://doi.org/10.1111/j.1365-2966.2010.16878.x
  * [11] D. J. Price and et al., “Phantom: A Smoothed Particle Hydrodynamics and Magnetohydrodynamics Code for Astrophysics,” _Publications of the Astronomical Society of Australia_ , vol. 35, p. e031, Sep. 2018.
  * [12] S. Rosswog, “The Lagrangian hydrodynamics code magma2,” _Monthly Notices of the Royal Astronomical Society_ , vol. 498, no. 3, pp. 4230–4255, 08 2020. [Online]. Available: https://doi.org/10.1093/mnras/staa2591
  * [13] J. Barnes and P. Hut, “A hierarchical o(n log n) force-calculation algorithm,” _Nature_ , vol. 324, no. 6096, pp. 446–449, Dec 1986. [Online]. Available: https://doi.org/10.1038/324446a0
  * [14] J. S. Bagla, “Treepm: A code for cosmological n-body simulations,” _Journal of Astrophysics and Astronomy_ , vol. 23, no. 3, pp. 185–196, Dec 2002. [Online]. Available: https://doi.org/10.1007/BF02702282
  * [15] L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” _Journal of Computational Physics_ , vol. 135, no. 2, pp. 280–292, 1997. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999197957065
  * [16] H. Cheng, L. Greengard, and V. Rokhlin, “A fast adaptive multipole algorithm in three dimensions,” _Journal of Computational Physics_ , vol. 155, no. 2, pp. 468–498, 1999. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999199963556
  * [17] M. S. Warren and J. K. Salmon, “A portable parallel particle program,” _Computer Physics Communications_ , vol. 87, no. 1, pp. 266–290, 1995, particle Simulation Methods. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0010465594001774
  * [18] R. Capuzzo-Dolcetta and P. Miocchi, “A comparison between the fast multipole algorithm and the tree-code to evaluate gravitational forces in 3-d,” _Journal of Computational Physics_ , vol. 143, no. 1, pp. 29–48, 1998. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999198959496
  * [19] D. C. Marcello, “A Very Fast and Angular Momentum Conserving Tree Code,” _Astrophysical Journal_ , vol. 154, no. 3, p. 92, 2017.
  * [20] J. Loiseau and et al., “FleCSPH: The next generation FleCSIble parallel computational infrastructure for smoothed particle hydrodynamics,” _SoftwareX_ , vol. 12, p. 100602, 2020. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S2352711020303150
  * [21] B. K. Bergen and M. R. J. Charest, “Flexible computational science infrastructure (FleCSI): Overview & application progress,” Los Alamos National Lab.(LANL), Los Alamos, NM (United States), Tech. Rep., 2017.
  * [22] S. Chandrasekhar, _An Introduction to the Study of Stellar Structure_ , ser. Astrophysical monographs. Dover Publications, 1957. [Online]. Available: https://books.google.com/books?id=VwPLAgAAQBAJ
  * [23] A. Meurer, C. P. Smith, M. Paprocki, O. Čertík, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh _et al._ , “Sympy: symbolic computing in python,” _PeerJ Computer Science_ , vol. 3, p. e103, 2017.


# Searching for Mass Concentrations with Precision Pulsar Timing

John M. LoSecco Physics and Astronomy Department, University of Notre Dame,
Notre Dame, IN 46556-5670, USA John M. LoSecco<EMAIL_ADDRESS>
###### Abstract

This papers searches for evidence of mass concentrations along the path of
radio pulses in the PPTA2 survey data release. Radio pulse travel times are
influenced via gravitational fields along the path from the source to the
observer. Transient time delays in transit are a useful measure of the matter
distribution along the path. Many pulsars have very well understood timing
solutions with predicable arrival times and can be used to sample the mass
variation. Changes in the source, observer and mass concentration positions
produce changes in arrival times which can be significant for precision pulsar
times. Twelve candidates are reported from this search.

dark matter pulsars: general proper motions time gravitation

revtex4-1.clsrevtex4-2.cls

##  1 Introduction

Gravitational interactions have been suggested as a source of transient
distortions in pulsar arrival time. Larchenkova (1995). Recently efforts have
been made to understand and correct for a broad class of timing distortions in
pulsar arrival times with results approaching microsecond resolution. This
suggests that searching these data may help strengthen models of galactic mass
distributions.

In this note we search for the relativistic delay in pulsar pulses due to
encounters with the gravitational fields of masses distributed close the path
of source to observer, as illustrated in figure 1. A static mass distribution
would add an unmeasurable constant time delay to the signal, so just as in
other mass concentration searches Wyrzykowski (2022) one looks for a time
delay caused by a transient encounter.

Figure 1: An illustration of the source to observer geometry. The impact
parameter $D$ and the gravitating mass $\mathcal{M}$ are shown.

The accumulated time delay is given by Shapiro (1964); Larchenkova (2007):

$\Delta T=-\frac{2G\mathcal{M}}{c^{3}}\ln(1-\hat{R}\cdot\hat{s})$ (1)

where $\frac{2G\mathcal{M}}{c^{3}}$ is the Schwarzschild radius divided by the
speed of light and $\hat{R}\cdot\hat{s}$ is the cosine of the angle, as viewed
from the observer, between the source of the pulse and the source of the
gravitational field created by mass $\mathcal{M}$.

Pulsar timing arrays Kerr (2021) have improved arrival time fits for pulses
and modeled the various known sources of time deviations in the transit time.
Orbital motion of the source and the observer can be fit and removed.
Observing at multiple frequencies permits removal of frequency dependent
plasma delays. Even the relativistic delay caused by a massive companion of
the pulsar can be fit and removed.

##  2 Gravitational Time Delay

Figure 2: The time dependent geometry of the pulsar - mass concentration
cosine as viewed by the observer. The Z axis goes from the observer to the
pulsar (dot on the left). The point of closest approach is a distance $D$
along the X axis. The Y coordinate is upward along the projected speed of the
mass concentration. The displacement along Y is taken as $Vt$. The mass
concentration is in the plane at (D,Vt).

In this paper we fit the pulsar timing residuals of the PPTA2 survey sources
Kerr (2021) to a 4 parameter model of an encounter of the beam with a mass
source moving with constant speed in the plane of the observation.

Larchenkova (1995) suggest the 4 parameter fit of equation 2.

$\Delta T=-M~{}\ln{(1-\frac{1}{\sqrt{1+(v(t-t_{0}))^{2}+d^{2}}})}$ (2)

with $M=\frac{2G\mathcal{M}}{c^{3}}$ as the mass of the source of the field
(in units of microseconds), $t_{0}$ is the time of the closest encounter, $v$
is the (assumed) constant projected speed of the encounter and $d$ is the
distance of closest approach, at $t=t_{0}$. For reference $M_{sun}$ is 9.85
microseconds. To be more precise $d$ is the tangent of the angle of closest
approach and $v(t-t_{0})$ is the tangent of the angle from the point of
closest approach to the location at time $t$

Equation 2 is found by explicitly calculating the value of
$\hat{R}\cdot\hat{s}$ in the geometry shown in figure 2. Take the $\hat{z}$
axis along the direction from the observer to the pulsar. The distance from
the pulsar beam path to the mass concentration candidate point of closest
approach is $D$. The direction from the pulsar beam path to the mass
concentration candidate point of closest approach defines the $\hat{x}$ axis.
The vertical axis in the figure is the $\hat{y}$ axis and is the direction of
motion of the mass concentration candidate. The Z distance is arbitrary but we
can take it as the distance along the $\hat{z}$ axis to the point of closest
approach of the beam to the mass concentration candidate. Figure 2 is drawn in
the plane perpendicular to the $\hat{z}$ axis at distance Z. In this
coordinate system $\vec{R}=(0,0,Z)$ and $\vec{s}=(D,Vt,Z)$.

$\begin{array}[]{c}\displaystyle\hat{R}\cdot\hat{s}=\frac{Z^{2}}{Z\times\sqrt{Z^{2}+(Vt)^{2}+D^{2}}}\\\
\displaystyle=\frac{Z^{2}}{Z^{2}\sqrt{1+(vt)^{2}+d^{2}}}\\\
\displaystyle=\frac{1}{\sqrt{1+(vt)^{2}+d^{2}}}\end{array}$ (3)

which is equation 2 with $d=D/Z$ and $v=V/Z$.

##  3 The Data

The PPTA2 release Antoniadis (2022) includes 65 millisecond pulsars and their
fits. This data was processed with Tempo2 Edwards (2006) to provide residuals,
the deviation of the pulse arrival time from the time estimated by the model.
These residuals hold information about potential encounters with massive
bodies on route to the detector.

The IPTA data release 2 combines data from prior releases, EPTA 1, NANOGrav
9-year, PPTA and their extensions. The observations have been made by seven
different radio telescopes (Effelsberg, Lovell, Nançay, Westerbork, Green
Bank, Arecibo, Parkes) The pulsar data samples span varying time intervals
from 231 to 10753 days and have from 116 to 17487 fitted pulses. A total of
210148 pulse residuals over a total time period of 255327 MJD are provided in
the release.

##  4 Method

The method used has two parts. The first part searches for signal delay
candidates over a limited portion of the recorded data. Once found a local fit
is done on that candidate to recover the fit parameters.

Based on gravitational lensing results the encounters are expected to occur
over a period of weeks or months. The data sample span periods from less than
a year to almost 30 years. Reliable fits also require many more data points to
be fit compared to the four parameters to be measured. The method adopted is
to fit overlapping regions of at about 120 consecutive residuals, less if
there is not enough data to provide 120 samples. A modest event candidate is
refit with the candidate time at the center of the interval. When possible the
time width of the fit is chosen to be the point where the peak has dropped to
10% of its maximum height. The peak position and width are based on the
parameters extracted from the search step. The fit is sensitive to the tails
of the fit function so having two tails to fit helps reject noise. If the
candidate time is near either end of the data sample that tail is truncated.
Using a sample size large compared to the encounter time makes the fit more
sensitive to fluctuations in the underlying tempo2 fit since the fit curve
returns to zero on both sides of an encounter. An example of a fit is in
figure 3

Figure 3: A fit to pulsar J0034-0534 showing a significant delay at 52371.4
MJD.

To proceed from the search fit to the centered fit the search result must have
a log likelihood ratio test at the 95% level. The fit itself should have an
Anderson-Darling test for gaussian residuals at the 95% level and the absolute
value of the fitted mass, $|M|$ must be at least 3 time its $\sigma$.
Accepting $M<-3\sigma$ is included to study potential background fluctuations
indicated by negative time delays with the same time structure expected from a
point mass.

The same selection plus the requirement that $\chi^{2}/DOF$ probability is
greater than 10% is applied to the centered fit before it can be accepted as a
Shapiro delay candidate.

##  5 Technical Details

From equation 1 it is clear that the fit can diverge as
$\hat{R}\cdot\hat{s}\rightarrow 1$. Since the residual data is finite it is
unlikely that the solution will fall on this point but the fitter may still
explore this region to approach the best fit. Due to finite precision of
computer representations of numbers $1-\cos{\theta}$ may approach 0 more
rapidly due to loss of precision. The fitter uses the “log1p” function of
$-\cos{\theta}$ to delay the round off error in the logarithm function
evaluation. At even smaller values of $\theta$ the approximation of
$1-\cos{\theta}=\theta^{2}/2$ is used as the argument to the logarithm.
Eventually that is replaced with $\ln(\theta^{2}/2)=2\ln(\theta)-\ln(2)$. This
provides access to a maximum amplification of about 1490 with standard double
precision calculations.

The transition from search to fit uses the preliminary fitted velocity from
the search to determine the temporal range of the fit. Ideally one wants a
time range in which the fit is dominated by the signal and not the asymptotic
tails which should go to zero. The first attempt is to construct a time
symmetric region about the candidate time with the amplitude dropping to 10%
of the peak height. In the case of a very local fluctuation the search fitted
velocity is very large and one ends up with very few timing residuals in the
narrow time range. To get more statistics the refit region is expanded to get
at least 100 points.

Alternatively a systematic shift in the residuals raises the baseline. There
is no baseline parameter in the fit so this is turned into a very small
velocity. The small velocity results in a very large time range for the refit,
sometimes including the whole data sample for the pulsar. The refitting for
nearby search regions then have a large overlap, reducing their independence.
In this case the refit is shortened to under 500 residuals.

In both cases the refit sample is centered on the candidate to the extent
possible and the initial values for the refit parameters are those taken from
the search. “to the extent possible” is needed since the data itself is not
uniformly sampled in time and is a merger of data from multiple instruments.
The search time is put in the middle of the refit sample. Gaps in the data
sample and different sampling rates may cause the centered fit to not have the
same number of residuals before and after the event time. Lower or upper time
limits of the whole data sample also restrict the range for the centered fit.

In the case of low velocities this adjustment procedure may produce a small
bias against a true signal with the direction of proper motion close to the
line of sight. On the other hand the search samples may have a non-zero
offset, perhaps coming from very long period signals. A symptom of such an
offset is a very small value for $v$. If $v\Delta t<10$ mrad the offset is
removed by subtracting the sample mean from the sample. $\Delta t$ is the
actual time length of the fit sample. This baseline adjustment was not need
for any of the candidate encounters in this paper.

An MCMC based search was used to explore the parameter space near the solution
to improve the error estimates. A fit was accepted as a candidate if
$\mathcal{M}/\sigma_{\mathcal{M}}>3$ and the $\chi^{2}$/DOF probability exceed
0.1.

##  6 Results

Pulsar | Sample | Duration | Count | RMS | Prob | $\mathcal{M}\pm\sigma_{\mathcal{M}}$ | $\mathcal{M}/\sigma_{\mathcal{M}}$
---|---|---|---|---|---|---|---
J0613-0200 | 19 | 236.8 | 457 | 3.551 | 0.9995 | 0.01644 $\pm$ 0.004938 | 3.329
J1024-0719 | 29 | 111.7 | 414 | 10.54 | 0.9139 | -0.4137 $\pm$ 0.07403 | -5.588
J0034-0534 | 0 | 2011.4 | 109 | 9.095 | 0.8481 | 1.656 $\pm$ 0.3839 | 4.314
J1024-0719 | 75 | 94.7 | 446 | 7.891 | 0.6924 | 0.3192 $\pm$ 0.1049 | 3.042
J1713+0747 | 211 | 53.3 | 452 | 2.111 | 0.6369 | 0.03642 $\pm$ 0.009973 | 3.652
J1640+2224 | 3 | 84.1 | 153 | 4.339 | 0.6173 | 0.06589 $\pm$ 0.01322 | 4.984
J1857+0943 | 46 | 127.4 | 403 | 2.747 | 0.4474 | 0.05748 $\pm$ 0.009663 | 5.948
J1909-3744 | 175 | 54.2 | 476 | 1.273 | 0.2458 | 0.01408 $\pm$ 0.003733 | 3.772
J1643-1224 | 17 | 230.7 | 399 | 3.563 | 0.1902 | 0.4636 $\pm$ 0.09352 | 4.957
J1024-0719 | 39 | 123.7 | 442 | 10.72 | 0.1343 | -0.2983 $\pm$ 0.07129 | -4.184
J0613-0200 | 72 | 70 | 436 | 2.661 | 0.1122 | 0.03849 $\pm$ 0.01255 | 3.068
J1643-1224 | 1 | 1765.6 | 103 | 10.78 | 0.1056 | -2.206 $\pm$ 0.4336 | -5.087
Table 1:  There are 12 candidates after duplicates have been removed. Eight
pulsars are represented. The first 4 columns give details of the fit sample
including the time duration of the sample and the number of points fit. The
fifth column gives the RMS value of the residuals used in the fit as a measure
of the noise level in the sample. The last 3 columns give the result of the
centered 10% fit including the $\chi^{2}$ confidence of the fit, the mass
fitted and the significance of the mass fit in $\sigma$’s. The RMS and the
mass are given in microseconds.

2084 candidates pass the search criteria to be recentered and refit. After the
refit 108 of these have $\mathcal{M}/\sigma_{\mathcal{M}}>3$. 16 of these have
a fit with a $\chi^{2}/DOF$ probability $>0.1$, the criteria in the centered
fit to become a final event candidate. Since the search uses overlapping
regions and the centered fit can extend the fit into regions already fit there
are duplicate candidates in some cases. Since the data samples do not overlap
completely the duplicates are not identical but the date is very close and the
other parameters agree to fractions of a $\sigma$. The 4 duplicates have been
removed from the table.

Table 1 lists the results of the search. The table includes details about the
sample and results of the fit including the fitted mass and the $\chi^{2}/DoF$
probability. The number of residuals fit, the duration of the sample (in MJD)
and the RMS of the sample in microseconds are also listed

Table 2 includes the full list of the fitted parameters including $d$, $v$ and
the time of closest approach $t_{0}$. There are no selection constraints on
the values or significance of these parameters. These data are also printed on
the plots.

##  7 Figures

Figures showing all the fits are displayed in this section. The figures are
ordered from high to low by $\chi^{2}/DOF$ probability. Each figure has an
upper panel showing the data and the fitted curve. The lower panel shows the
data minus the fitted value of the delay curve, the new residuals. The dashed
line on the lower curve is zero and the large point marks the position of the
peak in the upper plot. The vertical axes on each plot are individually scaled
and range from microseconds to tens of microseconds. The auto-scaling
sometimes makes the signal peak hard to see due to large error bars on some of
the points. The values for “vel”=$v$ and “ang”=$d$ printed on the plots are
scaled by 10,000 to reduce the number of zeros to the right of the decimal.

If the fitted curve is small compared to the vertical scale of the plot a
dashed red curve, a rescaled copy of the signal fit is also plotted. The
legend gives the rescaling.

Pulsar | Sample | $\mathcal{M}\pm\sigma_{\mathcal{M}}$ | $d\pm\sigma_{d}$ | $v\pm\sigma_{v}$ | $t_{0}\pm\sigma_{t_{0}}$
---|---|---|---|---|---
J0613-0200 | 19 | 0.01644$\pm$0.004938 | -1.314$\times 10^{-15}$$\pm$4.172$\times 10^{-13}$ | -5.321$\times 10^{-10}$$\pm$3.357$\times 10^{9}$ | 53859$\pm$0.0062922
J1024-0719 | 29 | -0.4137$\pm$0.07403 | 1.186$\times 10^{-9}$$\pm$0.0538 | 2.584$\times 10^{-6}$$\pm$4.426$\times 10^{-6}$ | 55561$\pm$3.185
J0034-0534 | 0 | 1.656$\pm$0.3839 | -1.466$\times 10^{-8}$$\pm$654.6 | 1.236$\times 10^{-3}$$\pm$0.0004958 | 52371$\pm$23.047
J1024-0719 | 75 | 0.3192$\pm$0.1049 | -2.109$\times 10^{-9}$$\pm$0.00826 | 2.841$\times 10^{-4}$$\pm$0.0004075 | 56396$\pm$6.8566
J1713+0747 | 211 | 0.03642$\pm$0.009973 | -2.757$\times 10^{-8}$$\pm$0.00000118 | 2.035$\times 10^{-2}$$\pm$0.01307 | 56271$\pm$0.0002075
J1640+2224 | 3 | 0.06589$\pm$0.01322 | -7.233$\times 10^{-9}$$\pm$0.00000815 | 4.909$\times 10^{-4}$$\pm$0.0008107 | 53355$\pm$0.0074092
J1857+0943 | 46 | 0.05748$\pm$0.009663 | -9.726$\times 10^{-12}$$\pm$0.0001744 | 7.528$\times 10^{-9}$$\pm$2.009$\times 10^{-8}$ | 56120$\pm$12.117
J1909-3744 | 175 | 0.01408$\pm$0.003733 | 1.534$\times 10^{7}$$\pm$0.0004855 | -2.204$\times 10^{-2}$$\pm$0.014 | 56522$\pm$0.022584
J1643-1224 | 17 | 0.4636$\pm$0.09352 | -1.05$\times 10^{6}$$\pm$7.695 | 5.348$\times 10^{-3}$$\pm$0.002258 | 53994$\pm$1.4459
J1024-0719 | 39 | -0.2983$\pm$0.07129 | -2.229$\times 10^{-12}$$\pm$7.177 | 1.422$\times 10^{-6}$$\pm$3.241$\times 10^{-6}$ | 55692$\pm$2.2011
J0613-0200 | 72 | 0.03849$\pm$0.01255 | 1.936$\times 10^{-11}$$\pm$2.007$\times 10^{-10}$ | -1.55$\times 10^{-5}$$\pm$0.00005461 | 55578$\pm$0.00043247
J1643-1224 | 1 | -2.206$\pm$0.4336 | 1.619$\times 10^{-6}$$\pm$15.02 | 6.51$\times 10^{-4}$$\pm$0.0002622 | 51024$\pm$17.367
Table 2: Summary of the fit parameters for the 12 events Pulsar | Sample Length (MJD) | Sample Size | Candidates
---|---|---|---
J0034-0534 | 4935.41 | 276 | 1
J0613-0200 | 5863.88 | 9322 | 2
J1024-0719 | 6647.34 | 5865 | 3
J1640+2224 | 6301.78 | 3097 | 1
J1643-1224 | 7356.24 | 8136 | 2
J1713+0747 | 8201.93 | 17487 | 1
J1857+0943 | 10344.48 | 5004 | 1
J1909-3744 | 3951.20 | 11483 | 1
Table 3: Summary of the data sample of the 8 pulsars

##  8 Conclusions

This article presents some evidence of detection of massive objects via their
gravitational fields. There has been no attempt made to determine if the
objects detected via timing delay are also visible via some other detection
method. The measured masses may suggest a visible object. Since this paper has
measured transients one may not be able to find a visible match in the same
direction with data not collected at the same time. The fit yields an angular
velocity on the sky but the direction in the plane is not known. This method
does not provide a distance estimate for these objects except that it is
limited by the distance to the pulsar source.

Three of the events in table 1 have $M<0$ which is not compatible with the
local mass hypothesis on which the fit function is based. These could be due
to statistical fluctuations, noise, a systematic problem or something new. All
three of these negative fits are found for pulsars, J1024-0719 and J1643-1224,
that also have positive events in the table.

The cumulative distribution function of the $chi^{2}$ probability distribution
is consistent with flat except for the region 0.1 to 0.2 which has 4 of the 12
events. This suggests that the sample has relatively little background except
for the region 0.1 to 0.2 which shows some evidence for noise content.

One might expect noise events to be distributed uniformly in time or sample
size. Table 3 lists exposure time and sample size for pulsars showing signals.
The total data set spans 255327 MJD with 210148 timing residuals. The pulsar
samples in table 3 represent 19% of the exposure time and 29 % of the timing
residuals.

Some objects, such as rotating, Kerr-Newman, solutions have the potential to
shorten travel times Pei (2021). But it is not clear if they would have the
time profile fit for a relativistically induced time delay.

##  9 Future Possibilities

Classic gravitational lensing Refsdal (1964) provides three observables, the
time delay used in this analysis, amplification of the signal used in
microlensing and multiple images. It may be possible to add these to the
pulsar observables. The recorded energy may show changes correlated with the
time delays or one may be able to discern changes in the apparent source size
correlated with the delay. Even modest additional input may help reject noise
events and better characterize true lensing sources.

## Acknowledgments

This work would not have been possible without the efforts of the maintainers
of the PPTA second data release and the Tempo2 code used to fit the residuals.
Jim Rich has been an invaluable aid in understanding the microlensing method
of observing dark matter that has many concepts in common with those used
here.

## References

  * Larchenkova (1995) T.I Larchenkova and O.V. Doroshenko, “Pulsars as a Tool for Detection of Dark Matter in the Galaxy” Astron. Astrophys. 297, 607-609 (1995)
  * Wyrzykowski (2022) L. Wyrzykowski et al., “Gaia Data Release 3 Microlensing Events from All Over the Sky” Astron. Astrophys. abc, 1-45 (2022). (arXiv:2206.06121)
  * Shapiro (1964) Irwin I. Shapiro (1964), “Fourth Test of General Relativity”, Phys. Rev, Lett. 13 (26): 789–791 (1964).
S. Weinberg, “Gravitation and Cosmology”, Wiley & Sons, New York, (1972).

  * Larchenkova (2007) T.I Larchenkova and A.A. Lutovinov (2007), “On the Possibility of Observing the Shapiro Effect for Pulsars in Globular Clusters”, Astro. Lett. 33 (7): 455-467 (2007)
S. Desai and E.O. Kahya (2016), “Galactic one-way Shapiro delay to PSR
B1937+21”, Mod. Phys. Lett. A 31 (13): 1650083-1 - 1650083-14 (2016).

  * Kerr (2021) M. Kerr et al., “The Parkes Pulsar Timing Array project: second data release”, Publications of the Astronomical Society of Australia, 37, e020 (2020)
D.J. Reardon et al., “The Parkes pulsar timing array second data release:
Timing analysis”. MNRAS, 507, 2137-2153 (2021).

  * Antoniadis (2022) J. Antoniadis et al., “The International Pulsar Timing Array second data release: Search for an isotropic Gravitational Wave Background”, MNRAS, 510 4873-4887 (2022).
  * Edwards (2006) R.T.Edwards, G.B.Hobbs, R.N.Manchester, “tempo2, a new pulsar timing package – II. The timing model and precision estimates”, MNRAS, 372, 1549-1574 (2006).
  * Shapiro (1968) Irwin I. Shapiro et al., “Fourth Test of General Relativity: Preliminary Results”, Phys. Rev, Lett. 20, 1265–1269 (1968).
  * Pei (2021) Ting-Hang Pei, “The Superluminal Phenomenon of Light Near the Kerr–Newman Black Hole or Super-Gravitational Source”, Front. Phys., Volume 9 - 2021 | https://doi.org/10.3389/fphy.2021.701619 (2021).
Tien Hsieh, Da-Shin Lee and Chi-Yong Lin. “Gravitational time delay effects by
Kerr and Kerr-Newman black holes in strong field limits”, Phys. Rev D 104.
104013 (2021).

  * Refsdal (1964) Sjur Refsdal, “The Gravitational Lens Effect”, MNRAS, 128 295-306 (1964).


11institutetext: Computer Science Department, Ariel University, Israel
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
22institutetext: Computer Science Department, Bar Ilan University, Israel
22email<EMAIL_ADDRESS>
# Weighted Burrows-Wheeler Compression

Aharon Fruchtman 11 Yoav Gross 11
Shmuel T. Klein 22 0000-0002-9478-3303 Dana Shapira 11 0000-0002-2320-9064

###### Abstract

A weight based dynamic compression method has recently been proposed, which is
especially suitable for the encoding of files with locally skewed
distributions. Its main idea is to assign larger weights to closer to be
encoded symbols by means of an increasing weight function, rather than
considering each position in the text evenly. A well known transformation that
tends to convert input files into files with a more skewed distribution is the
Burrows-Wheeler Transform. This paper employs the weighted approach on
Burrows-Wheeler transformed files and provides empirical evidence of the
efficiency of this combination.

##  1 Introduction

The Burrows-Wheeler Transform (BWT) [1] is the basis of the popular
compression method bzip2, yielding, on many types of possible input files,
better compression than gzip and other competitors. As a matter of fact, BWT
itself is not a compression method: its output is a permutation of its input,
which has obviously the same size. The usefulness of the transformation is
that it has a tendency to reorganize the data into what seems to be a more
coherent form, grouping many, though not all, identical characters together.
The output is therefore usually more compressible, by applying as simple
methods as run-length coding and move-to-front.

The combination of different methods to be applied one after the other, an
action known as cascading, is not new to data compression. It can, for
example, be found in gzip, which first parses the input using LZ77 [10], and
then applies Huffman coding [7] to the parsed elements. It therefore seems
natural to try to apply BWT in a pre-processing stage, and then compress the
transformed string by means of more sophisticated compression schemes. This
strategy fails, however, for static or dynamic Huffman or arithmetic coding.
The contribution of this paper is to give empirical evidence that cascading
BWT with the recently developed weighted compression schemes implies
significant savings.

A family of dynamic compression algorithms, named weighted coding, has
recently been proposed, which is especially suitable for the encoding of files
with locally skewed distributions. The main idea of the weighted approach is
to assign larger weights to closer to be encoded symbols by means of an
increasing weight function, rather than considering each position in the text
evenly.

Traditional dynamic algorithms use the distribution of the symbols in the
already processed portion of the input file as an estimate for the
distribution of the elements still to come later in the text. However, the
assumption that the past is a good approximation for the future, is not
necessarily true. A Forward-looking dynamic algorithm, using the true
distribution of the remaining portion of the file, was suggested in [6]. In
this method the frequencies of the symbols in the entire file are prepended to
the compressed file. In the encoding process, these frequencies are gradually
updated to reflect the number of occurrences in the remaining part of the file
by decrementing the frequency of the character that is currently being
processed. A hybrid method, combining both traditional and forward-looking
approaches, is proposed in [5]: in this method the frequencies for all
characters are not transmitted at the beginning of the file but rather
progressively, each time a new character is encountered.

Pushing this approach even further, a family of Forward weighted coding
schemes is proposed in [3], in which the encoding is based on the distribution
derived from index-dependent weights. That is, this weighted method assigns
higher priorities to positions that are closer to the currently processed one
in the encoding process, rather than treating all positions in the input file
equally. The weights assigned to the positions are generated by a non
increasing function $f$, and the weight for each symbol $\sigma$ is the sum of
the values of $f$ over all the positions where $\sigma$ occurs in the portion
of the input file that is still to be encoded.

Recently, a specific variant named the Backward Weighted coding has been
studied [4], which suggests a heuristic based on a weighted distribution,
calculated only over positions that have already been processed. The core
advantage of such a heuristic approach is a negligible header, relatively to a
costly header used in Forward Looking and all Forward Weighted variants.
Empirical tests have shown that backward weighted techniques can improve
beyond the lower bound given by the entropy for static encoding.

The paper is constructed as follows. Section 2 reviews the notation and
formulation needed for the weighted encodings using a running example. Section
3 proves some properties involving the combination of general dynamic coding
techniques and BWT. Section 4 presents empirical outcomes supporting the
compression efficiency of the proposed method even in practice.

##  2 Notation and Formulation

For the completeness of the paper, we include the definitions given in [3],
which formalize entropy based compression methods.

Let $T=T[1,n]$ be an input file of size $n$ over an alphabet $\Sigma$ of size
$m$. A weight $W$ is defined based on the following parameters,

  * -

A non negative function $g$, $g:[1,n]\longrightarrow{\rm I\\!R}^{+}$, which
assigns a positive real number to integers, seen as an assignment of a weight
to each position $i\in[1,n]$ within $T$;

  * -

A symbol of the alphabet, $\sigma\in\Sigma$;

  * -

An interval $[\ell,u]$, $1\leq\ell\leq u\leq n$ for restricting the domain of
the function $g$.

The value of $W(g,\sigma,\ell,u)$ is defined for each symbol $\sigma$, as the
sum of the values of the function $g$ for all positions $j$ in the range
$[\ell,u]$ at which $\sigma$ occurs. Formally,

$W(g,\sigma,\ell,u)\;=\sum_{\\{\ell\leq j\leq u\ \mid\ T[j]=\sigma\\}}{g(j)}.$

As a special case of weighted coding, Backward Weighted considers all the
positions that have already been processed, that is, the interval is of the
form $[\ell,u]=[1,i-1]$, and

$W(g,\sigma,1,i-1)\;=\sum_{\\{1\leq j\leq i-1\ \mid\ T[j]=\sigma\\}}{g(j)}.$

Traditional encoding methods can be reformulated as special instances of $W$
for which $g=\mathds{1}\equiv g(i)=1$ for all $i$. For example, static
compression refers to weights for which
$W(g,\sigma,\ell,u)=W(\mathds{1},\sigma,1,n)$ is constant for all indices.

| $i$ |  | 1 | 2 | 3 | 4 | $\cdots$ | 13 | 14 | 15 | 16 | $\cdots$ | 35 | 36 | 37 | 38 | $\cdots$ | 49 | 50 | avg
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| $T$ |  | a | t | a | t | $\cdots$ | a | t | c | g | $\cdots$ | c | g | a | t | $\cdots$ | a | t |
|  | $W$ | 14 | 14 | 14 | 14 | … | 14 | 14 | 11 | 11 | $\cdots$ | 11 | 11 | 14 | 14 | $\cdots$ | 14 | 14 |
static |  | TotIndx | 50 | 50 | 50 | 50 | … | 50 | 50 | 50 | 50 | $\cdots$ | 50 | 50 | 50 | 50 | $\cdots$ | 50 | 50 | 1.990
|  | IC | 1.84 | 1.84 | 1.84 | 1.84 | $\cdots$ | 1.84 | 1.84 | 2.18 | 2.18 | $\cdots$ | 2.18 | 2.18 | 1.84 | 1.84 | $\cdots$ | 1.84 | 1.84 | (+0.291)
|  | Indx$W$ | 1 | 1 | 1 | 1 | $\cdots$ | 1 | 1 | 1 | 1 | $\cdots$ | 1 | 1 | 1 | 1 | $\cdots$ | 1 | 1 |
|  | $W$ | 1 | 1 | 2 | 2 | $\cdots$ | 7 | 7 | 1 | 1 | $\cdots$ | 11 | 11 | 8 | 8 | $\cdots$ | 14 | 14 |
b-adp |  | TotIndx | 4 | 5 | 6 | 7 | $\cdots$ | 16 | 17 | 18 | 19 | $\cdots$ | 38 | 39 | 40 | 41 | $\cdots$ | 52 | 53 | 2.111
|  | IC | 2.00 | 2.32 | 1.58 | 1.81 | $\cdots$ | 1.19 | 1.28 | 4.17 | 4.25 | $\cdots$ | 1.79 | 1.83 | 2.32 | 2.36 | $\cdots$ | 1.89 | 1.92 |
|  | Indx$W$ | 1 | 1 | 1 | 1 | $\cdots$ | 1 | 1 | 1 | 1 | $\cdots$ | 1 | 1 | 1 | 1 | $\cdots$ | 1 | 1 |
|  | $W$ | 14 | 14 | 13 | 13 | $\cdots$ | 8 | 8 | 11 | 11 | $\cdots$ | 1 | 1 | 7 | 7 | $\cdots$ | 1 | 1 |
f-adp |  | TotIndx | 50 | 49 | 48 | 47 | $\cdots$ | 38 | 37 | 36 | 35 | $\cdots$ | 16 | 15 | 14 | 13 | $\cdots$ | 2 | 1 | 1.820
|  | IC | 1.82 | 1.84 | 1.81 | 1.88 | $\cdots$ | 2.25 | 2.21 | 1.71 | 1.67 | $\cdots$ | 4.00 | 3.91 | 1.00 | 0.89 | $\cdots$ | 1.00 | 0.00 | (+0.291)
|  | Indx$W$ | 1 | 1 | 1 | 1 | $\cdots$ | 4 | 4 | 4 | 8 | $\cdots$ | 64 | 128 | 128 | 128 | $\cdots$ | 512 | 512 |
b-2 |  | $W$ | 1 | 1 | 2 | 2 | $\cdots$ | 12 | 13 | 1 | 1 | $\cdots$ | 261 | 281 | 16 | 17 | $\cdots$ | 1552 | 1809 | 1.981
| TotIndx | 4 | 5 | 6 | 7 | $\cdots$ | 27 | 31 | 35 | 39 | $\cdots$ | 575 | 639 | 767 | 895 | $\cdots$ | 4095 | 4607
|  | IC | 2.00 | 2.32 | 1.58 | 1.81 | $\cdots$ | 1.17 | 1.25 | 5.13 | 5.29 | $\cdots$ | 1.14 | 1.19 | 5.58 | 5.72 | $\cdots$ | 1.40 | 1.35 |
|  | Indx$W$ | 1.00 | 1.15 | 1.32 | 1.52 | $\cdots$ | 5.28 | 6.06 | 6.96 | 8.00 | $\cdots$ | 111.43 | 128.00 | 147.03 | 168.90 | $\cdots$ | 776.05 | 891.44 |
b-weight |  | $W$ | 1.00 | 1.00 | 2.00 | 2.15 | $\cdots$ | 14.39 | 16.38 | 1.00 | 1.00 | $\cdots$ | 327.96 | 376.58 | 19.67 | 22.44 | $\cdots$ | 1988.36 | 2283.88 | 1.989
| TotIndx | 4.00 | 5.00 | 6.15 | 7.47 | $\cdots$ | 32.77 | 38.05 | 44.11 | 51.08 | $\cdots$ | 746.65 | 858.08 | 986.08 | 1113.11 | $\cdots$ | 5216.21 | 5992.26
|  | IC | 2.00 | 2.32 | 1.62 | 1.80 | $\cdots$ | 1.19 | 1.22 | 5.46 | 5.67 | $\cdots$ | 1.19 | 1.19 | 5.65 | 5.66 | $\cdots$ | 1.39 | 1.39 |
Table 1: Coding example for $T={\tt{(at)}^{7}{(cg)}^{11}{(at)}^{7}}$.

As a short illustration for the weighted approach, Table 1 brings a
comparative chart for the encoding of a small example of 50 characters:

$T=x_{1}\cdots x_{50}={\tt(at)}^{7}{\tt(cg)}^{11}{\tt(at)}^{7},$

shown in the second row of the table for several representative portions of
$T$, just underneath the indices. The static compression for $T$ considers the
probabilities $\frac{14}{50}$ for a and t and $\frac{11}{50}$ for c and g.
These probabilities can be calculated from the two first rows corresponding to
the static method in Table 1, the first entitled $W$ and representing the
weight of the specific character, and the second entitled TotIndx representing
the cumulative weights of all the characters. The ratio of these weights can
be considered as a probability $p_{i}$, and the corresponding Information
content, $-\log p_{i}$ for each position $i$, is shown in the line entitled
IC. For static, the IC values are 1.84 for a and t and 2.18 for c and g. The
last column of the table, headed avg, gives the weighted average of these IC
values, which is in fact the entropy.

The classic adaptive coding, b-adp [9], is a specific backward weight method
in which $g(i)=1$ for all $i$, and the backward weights refer to all positions
$i$ with $1\leq i<n$ by:

$W(\mathds{1},\sigma,1,i-1)=\sum_{\\{1\leq j\leq i-1\ \mid\
T[j]=\sigma\\}}{1},$

which is simply the number of occurrences of the current character $\sigma$ up
to that point, i.e., in $T[1,i-1]$. The details appear in the second part of
Table 1, headed b-adp. The line entitled $W$ refers now, at position $i$, to
specific weight of the character $\sigma=T[i]$ up to the given column $i$ of
the table, that is, the sum of the index weights Indx$W$ for those indices
$j<i$ at which the character $\sigma$ occurs, including the initial 1 values.
For b-adp, as well as for the following method f-adp, the values of Indx$W$
are just 1 for every $i$. The cumulative $W$ values for all the characters
$\sigma\in\Sigma$ are given in the row entitled TotIndx.

The symmetric counterpart of b-adp is the forward looking method, f-adp,
which, at each location $i$, considers the positions yet to come $[i,n]$
rather than those already processed $[1,i-1]$ as the backward looking b-adp.
That is,

$W(\mathds{1},\sigma,i,n)=\sum_{\\{i\leq j\leq n\ \mid\ T[j]=\sigma\\}}{1}.$

For our running example, f-adp initializes the weights of the characters a, t
to 14 and c, g to 11. The counts are then gradually decremented reflecting the
remaining number of occurrences for each $\sigma$ from position $i$ to the end
of $T$. The header for f-adp describes the exact frequencies of the involved
characters, and its size is 0.291 for our example, as for static.

A simple adaptive weighted coding, denoted by b-2, has been proposed in [4].
It is inspired by Nelson and Gailly [8], who rescale the frequencies in order
to cope with hardware constraints like the representation of the number of
occurrences as 16-bit integers. b-2 divides all the frequencies at the end of
every block of $k$ characters, for a given parameter $k$, regardless of
computer hardware restrictions. That is, the occurrences of characters at the
beginning of the input file contribute to $W$ less than those closer to the
current position. Furthermore, all positions within the same block contribute
equally to $W$, and their contribution weight is twice as large as the weight
assigned to the indices in the preceding block. Formally, for each pair of
indices $i$ and $i+k$, the function $g$ for b-2, denoted by
$g_{\mbox{\sf\scriptsize b-2}}$, maintains the equation
$g_{\mbox{\sf\scriptsize b-2}}(i+k)=2g_{\mbox{\sf\scriptsize b-2}}(i)$. The
first line of the block headed b-2 shows the index weight, Indx$W$, chosen
here with parameter $k=5$. Starting with 1, the value doubles after each block
of 5 positions. The other lines, entitled $W$, TotIndx and IC, are then
defined as above.

A refinement of b-2, named b-weight, is another special weighted coding [4]
inspired by the division by 2, but based on the function
$g_{\mbox{\sf\scriptsize b-weight}}(i)={(\sqrt[k]{2})}^{i-1}$ for $i\geq 1$,
for a given parameter $k$. The function $g_{\mbox{\sf\scriptsize b-weight}}$
still maintains a fixed ratio of 2 between blocks but with a smooth hierarchy
between all indices, rather than sharp differences at block boundaries. The
ratio of 2 for indices that are $k$ apart can be seen by:

$g_{\mbox{\sf\scriptsize
b-weight}}(i+k)={(\sqrt[k]{2})}^{i+k-1}={{(\sqrt[k]{2})}^{k}\cdot(\sqrt[k]{2})}^{i-1}=2\cdot
g_{\mbox{\sf\scriptsize b-weight}}(i).$

The index weight Indx$W$ for the b-weight method consists of real numbers
rather than integers as above. The shown values correspond to $k=5$, so that
$g_{\mbox{\sf\scriptsize b-weight}}(i)={(\sqrt[5]{2})}^{i-1}=1.149^{i-1}$.
This yields an average codeword length of 1.989 bits per symbol.

The weighted approach should be applied only on a text that has skewed
probability distributions in different portions of the file: there is a price
for adjusting the model in the transition between regions of different
distributions, and this overhead gets negligible only when the text becomes
long enough, or if the difference between the distributions is sufficiently
sharp as in this short example.

One of the features of the BWT is that it has a tendency to reorganize the
text $T$ such that BWT$(T)$ contains several runs of repeated characters. In
particular, for our running example, $T$, the BWT is

${\sf BWT}(T)={\tt t}^{7}{\tt g}{\tt t}^{6}{\tt a}^{14}{\tt g}^{10}{\tt t}{\tt
c}^{11},$

even though in the original $T$, there is not a single pair of identical
adjacent characters. This tendency implies that the local distributions seem
more skewed, which is advantageous to the weighted approach.

| $T$ | ${\sf BWT}(T)$
---|---|---
|  $k$ | header | bps |  $k$ | header | bps
static | – | 0.291 | 2.281 | – | 0.404 | 2.394
b-adp | – | – | 2.111 | – | 0.113 | 2.224
f-adp | – | 0.291 | 2.111 | – | 0.404 | 2.224
b-weight | 5 | – | 1.989 | 3 | 0.113 | 1.567
b-2 | 5 | – | 1.981 | 3 | 0.113 | 1.562
Table 2: Storage requirements of the encoding methods on ${\sf BWT}(T)={\tt
t}^{7}{\tt g}{\tt t}^{6}{\tt a}^{14}{\tt g}^{10}{\tt t}{\tt c}^{11}$.

We applied the same 5 methods on ${\sf BWT}(T)$. A small amendment is
necessary because the BWT on its own is not reversible—it needs in addition a
pointer to a starting point, actually the index of the last character of $T$
within BWT$(T)$, which requires $\log n$ bits, or $\frac{\log n}{n}$ per
symbol. We thus included this additional overhead of $\frac{\log
50}{50}=0.113$ bits in the header. Note that this overhead is constant for all
methods. On the other hand, the static and the forward compression methods
require information about the distribution probabilities, which should be
prepended to the compressed file independently from applying BWT or not. We
use the lower bound approximation given by the IC to express this information
in this example. Yet, for our experiments, we encoded the meta data
information by means of a universal Elias [2] $C_{\delta}$ code.

The comparison of the final compression results, before and after applying the
BWT, is presented in Table 2 and includes the parameter $k$, the header size
and the corresponding total storage costs including the header in bits per
symbol (bps). The first three columns refer to the encoding variants on the
original file $T$, and the last three columns are the results on BWT$(T)$. As
can be seen, while the improvement of the weighted methods with BWT is about
21%, the net encoding, excluding the header size, is the same for static,
b-adp and f-adp, regardless whether BWT has been applied or not. The total bps
results, including the header size, are identical for b-adp and f-adp, because
the gain in the net compression by f-adp is exactly the same as the loss
incurred by the additional overhead due to the exact frequencies. In fact,
these are not a coincidences, and in the following section we show that the
compression performance is preserved for static, b-adp and f-adp under any
permutations of the symbols of $T$, hence in particular for BWT. Moreover, we
show that the size of the compressed file including the header for b-adp and
f-adp is the same.

$0$$10$$20$$30$$40$$50$$0$$2$$4$$6$$8$$10$$12$$14$$i$Information
ContentstaticBWT staticf-adpBWT f-adpb-adpBWT b-adpb-2BWT b-2b-weightBWT-b-
weight Figure 1: Information content per index before and after applying BWT
on $T={\tt{(at)}^{7}{(cg)}^{11}{(at)}^{7}}$.

Figure 1 is a visualization of the differences between the methods, plotting
for each the information content as a function of the position $i$ for our
running example. The net encoding results before and after applying BWT are
depicted in cold and warm colors, respectively. We see that the fluctuations
are much more accentuated after applying BWT. Figure 2 displays the cumulative
values of the same data. As can be seen from both figures, the backward
weighted methods are more sensitive to fluctuations in the distribution, but
adjust faster to changes. The final points of the accumulated values for the
traditional methods differ only by the size of the header, and the advantage
of BWT for b-2 and b-weight can be seen by the fact that the corresponding
plots are below their counterparts without BWT. The differences between b-2
and b-weight, with and without BWT, are so small that their plots seem to be
overlapping.

$0$$10$$20$$30$$40$$50$$0$$20$$40$$60$$80$$100$$120$$i$Information
ContentstaticBWT staticf-adpBWT f-adpb-adpBWT b-adpb-2BWT b-2b-weightBWT-b-
weight Figure 2: Accumulated information content as a function of the size of
the processed prefix, before and after applying BWT on
$T={\tt{(at)}^{7}{(cg)}^{11}{(at)}^{7}}$.

##  3 Properties

We show in this section that for the three first mentioned methods, static,
b-adp and f-adp, applying BWT, or any other permutation, does not have any
influence of the compression by arithmetic coding. The next section then
brings empirical evidence that on the other hand, for the weighted methods, a
pre-processing stage by BWT does significantly improve the compression
performance.

Arithmetic coding starts with an interval [0,1), and repeatedly narrows it as
the text $T$ is being processed. The narrowing procedure is a proportional
refinement of the present interval into sub-portions according to the
probability distribution of the symbols of $\Sigma$. The encoding is a real
number that can be selected randomly within the final interval. Static
arithmetic coding uses a fixed probability distribution throughout the process
for the interval partition, while the (backward) adaptive method updates the
proportions for the corresponding partitions according to what has already
been seen. It is well known that the static variant achieves the entropy of
order 0. A straightforward corollary of this fact is that permuting $T$ does
not change the size of the output, though the output file itself will of
course be altered. This is stated in the following lemma.

###### Lemma 1

The size of the compressed file, after having applied static arithmetic
coding, is invariant under permutations of the original input.

###### Proof

Suppose we are encoding the text $T[1,n]$ using the static arithmetic method.
The notation of the weight introduced above for $\sigma=T[i]$,
$W(\mathds{1},T[i],1,n)$, refers to the number of occurrences of $T[i]$ within
$T[1,n]$. For simplicity, we shall use ${\sf occ}(\sigma)$ to denote
$W(\mathds{1},\sigma,1,n)$.

Each processed letter $T[i]$ narrows the current sub-interval of $[0,1)$ by a
factor equal to the probability of $T[i]$ in $P$, that is, by
$\frac{1}{n}W(\mathds{1},T[i],1,n)$. The range, $r$, of the final interval
after processing $T=T[1,n]$ is the product of the sizes of these intervals:

$r=\prod^{n}_{i=1}{\frac{1}{n}W(\mathds{1},T[i],1,n)}=\frac{1}{n^{n}}\displaystyle\prod_{\sigma\in\Sigma}\prod_{i=1}^{{\sf
occ}(\sigma)}{{\sf
occ}(\sigma)}=\frac{1}{n^{n}}\displaystyle\prod_{\sigma\in\Sigma}{\sf
occ}(\sigma)^{{\sf occ}(\sigma)},$

where the middle equality is obtained by reordering the multiplication factors
by character. The size of the compressed file is the information content of
choosing an interval of size $r$ given by $-\log_{2}r$, and it is independent
of the order in which the letters appear.

A similar property can be proven for traditional adaptive arithmetic coding,
as follows.

###### Lemma 2

The size of the compressed file, after having applied adaptive arithmetic
coding, is invariant under permutations of the original input.

###### Proof

To avoid the zero-frequency problem for encoding the first occurrence of a
letter $\sigma$ in $T$, the number of occurrences for each $\sigma$ is
initialized by 1. For adaptive arithmetic coding, b-adp, each processed symbol
$T[i]$, $1<i\leq n$, narrows the current sub-interval of $[0,1)$, representing
the processed prefix $T[1,i-1]$ of $T$ of size $i-1$, by a factor equal to the
probability of $T[i]$ in $T[1,i-1]$. This probability is equal to

$\frac{W(\mathds{1},T[i],1,i-1)+1}{m+i-1},$

where $m=|\Sigma|$ taking the initial 1-values of all the characters into
account. Multiplying the sizes of these intervals contributed by the
occurrences of $T[i]$, $1\leq i\leq n$, we get

$\displaystyle r$
$\displaystyle=\prod^{n}_{i=1}{\frac{W(\mathds{1},T[i],1,i-1)+1}{m+i-1}}=\Big{(}\prod^{n}_{i=1}{\frac{1}{m+i-1}}\Big{)}\prod_{\sigma\in\Sigma}\prod_{i=1}^{{\sf
occ}(\sigma)}{i}$
$\displaystyle=\frac{(m-1)!}{(m+n-1)!}\prod_{\sigma\in\Sigma}{\sf
occ}(\sigma)!$

The size of the compressed file is accordingly $-\log_{2}r$, so that permuting
the input text will not change the size of the output.

A similar proposition can be proven for f-adp, stated as follows.

###### Lemma 3

The size of the compressed file, after having applied forward arithmetic
coding, f-adp, is invariant under permutations of the original input.

###### Proof

By a similar argument to Lemma 2, each processed letter $T[i]$, $1<i\leq n$
narrows the current sub-interval, representing the processed suffix $T[i,n]$
of $T$ of size $n-i+1$, by a factor equal to the probability of $T[i]$ in
$T[i,n]$, which is equal to

$\frac{W(\mathds{1},T[i],i,n)}{n-i+1}.$

Multiplying the sizes of these intervals contributed by the occurrences of
$T[i]$, $1\leq i\leq n$,

$r=\prod^{n}_{i=1}{\frac{W(\mathds{1},T[i],i,n)}{n-i+1}}=\Big{(}\prod^{n}_{i=1}{\frac{1}{n-i+1}}\Big{)}\prod_{\sigma\in\Sigma}\prod_{i=1}^{{\sf
occ}(\sigma)}{i}=\frac{1}{n!}\prod_{\sigma\in\Sigma}{\sf occ}(\sigma)!$

As before, the size of the compressed file is $-\log_{2}r$, which is not
affected by permutation of the input file.

Besides proving that the order of the characters does not matter, only their
quantity, we also get an exact evaluation of the difference in size between
the encoded files:

###### Corollary 1

The forward looking encoding, f-adp, is better than the backward looking
encoding, b-adp by $\log{m+n-1\choose{n}}$ bits.

###### Proof

We use $r_{b}$ and $r_{\\!f}$ to denote the size of the final ranges of b-adp
and f-adp, given by Lemma 3 and Lemma 4, respectively. The difference between
the sizes of the compressed files of b-adp and f-adp is as follows.

$\displaystyle-\log r_{b}+\log r_{\\!f}$
$\displaystyle=\log{\Big{(}\frac{1}{n!}\prod_{\sigma\in\Sigma}{\sf
occ}(\sigma)!\Big{)}}-\log{\Big{(}\frac{(m-1)!}{(m+n-1)!}\prod_{\sigma\in\Sigma}{\sf
occ}(\sigma)!\Big{)}}$
$\displaystyle=\log{\Big{(}\frac{(m+n-1)!}{n!\cdot(m-1)!}}\Big{)}=\log{m+n-1\choose{n}}$

Interestingly, the difference between b-adp and f-adp only depends on the size
of the alphabet and the size of the file and is blind to the content itself.
In fact, on our datasets with $n=4M$ and $m=257$, including the EOF sign, the
difference between the compressed files, discarding the prelude was the
constant 494 bytes, as expected. However, as mentioned above, the prelude for
f-adp is more costly than that for b-adp, as it should include the exact
frequencies of the characters.

The number of sequences $\Big{(}{\sf occ(\sigma_{1})},{\sf
occ(\sigma_{2})},\ldots,{\sf occ(\sigma_{m})}\Big{)}$ such that
$\sum_{i=1}^{m}{{\sf occ}(\sigma_{i})}=n$ is equal to the number of ways to
select $n$ numbers out of $n+m-1$, i.e., $n+m-1\choose{n}$.

The information content of selecting one of these choices,
$\log{n+m-1\choose{n}}$, coincides therefore with the difference between f-adp
and b-adp, which shows that basically there is no much difference between
their compression efficiency.

##  4 Experimental results

In order to study the performance of the weighted methods on BWT transformed
texts on real-life, rather than artificial data, we considered the datasets
from the Pizza & Chili corpus111http://pizzachili.dcc.uchile.cl/, a collection
of files of different nature and alphabets. We compared all methods before and
after BWT has been applied.

###  4.1 Compression Performance

| static | adp | b-2 | b-weight
---|---|---|---|---
|  |  | $T$ | (k) | ${\sf BWT}$ | (k) | $T$ | (k) | ${\sf BWT}$ | (k)
sources | 69.48 | 69.48 | 64.20 | (444) | 26.44 | (24) | 64.17 | (438) | 26.30 | (22)
xml | 65.54 | 65.54 | 65.27 | (9355) | 13.11 | (24) | 65.27 | (8799) | 13.05 | (24)
dna | 24.97 | 24.98 | 24.64 | (62) | 23.02 | (34) | 24.64 | (60) | 23.01 | (34)
english | 57.07 | 57.07 | 56.55 | (4352) | 29.41 | (36) | 56.54 | (3230) | 29.33 | (36)
pitches | 69.32 | 69.32 | 56.17 | (85) | 53.34 | (119) | 56.09 | (83) | 53.26 | (120)
proteins | 52.64 | 52.65 | 51.69 | (329) | 46.71 | (39) | 51.68 | (317) | 46.62 | (38)
Table 3: Compression performance on BWT transformed files (%) of the different
methods.

For our first experiments we used a 4MB prefix of each of our dataset files.
Table 3 shows the compression performance on the original and transformed
files, defined as the size of the compressed file divided by the size of the
original file, in percent. The column headed adp refers to both b-adp and
f-adp, which give identical results, as expected. As static and adp are not
affected by BWT by Lemma 1 and Lemma 2, their results are reported only once
in the first and second column, respectively. The following four columns refer
to b-2 and the last four columns to b-weight. The values of $k$ used to
achieve the results of b-2 and b-weight are displayed in parentheses.

As can be seen, b-2 and b-weight outperform the traditional variants on all
datasets even on the original file $T$. While they are better by up to 5% when
applied on the original file, they become better by more than 40% in some of
the cases, when applied on the BWT transformed file. More precisely, b-2 is
better than the traditional variants by 24.50% on the average, while b-weight
is better by 24.58% on the average. In fact, b-weight outperforms b-2 by 0.11%
on average on all files.

###  4.2 Iterative BWT

The BWT tends to rearrange the symbols of the text so that longer runs are
formed than in the original file. Since applying BWT turned out to be
beneficial for compressing afterwards according to methods based on the
weighted paradigms, this naturally raises the question whether re-applying BWT
iteratively will arrange the symbols even better in this regard. In other
words, given a text $T$, we shall try our methods on BWT${}^{i}(T)$, for
$i\geq 0$, in which $i=0$ refers to the original file $T$ itself.

00.51.01.52.02.53.0$0$$0.2$$0.4$$0.6$$0.8$$1$012301230123012301230123Compressed
SizeNNRDNAENGLISHPITCHESPROTEINSSOURCESXML Figure 3: File size in million
bytes and NNR.

To quantify this intuitive notion of being more ordered, we suggest the
following measure. A run is a sequence of maximal length of the same symbol.
For example, the string aaabaacccc consists of 4 runs, of lengths 3, 1, 2 and
4, respectively. The number of runs in the text $T$, divided by its length
$n$, may serve as a measure of its orderedness; we call it the normalized
number of runs, NNR, which is between 0 and 1. The rationale is that the
larger the number of runs in the input file, the shorter they are,
consequently causing a longer encoding, and vice versa.

The following graph displays, for each of the input files, the compression
performance on the $x$-axis and the NNR on the $y$-axis of b-2 when applied on
BWT${}^{i}(T)$, for $i=0,1,2,3$. For better visibility, the points
corresponding to a given input file appear in the same color and are connected
and numbered by $i$, the number of times BWT has been applied. As can be seen,
the compression efficiency improves for $i=1$ and 2, but becomes worse when
BWT is applied a third time on all of the input files. The smallest NNR is
always achieved for $i=1$, and it increases for higher $i$.

###  4.3 Weighted BWT Block Variant

Some drawbacks of the BWT are its streamless nature and its processing time
cost. These weaknesses stem from the fact that the entire text $T$ should be
known for the permuting process. To improve execution time, we suggest a block
based variant and examine the block-size vs. compression performance trade-
off. Figure 4 displays the compression efficiency as a function of the block
size for b-2. The results for b-weight produce practically the same plot. As
expected, larger blocks yield better compression for textual input, whereas
for special data like dna or pitches, the block size has only a minor
influence.

SOURCESXMLDNAENGLISHPITCHESPROTEINS$0.2$$0.4$Compressed
Ratio8K16K32K64K128K256K512K1M2M4M Figure 4: Compression efficiency of b-2 on
BWT as a function of block size.

##  5 Conclusion

This paper studies the compression performance of the weighted coding on
Burrows-Wheeler transformed files. We have shown that statistical methods
which treat all positions in the files evenly are indifferent to permutations
in the input file, and to BWT in particular. On the other hand, the weighted
approach, being more suitable to skewed files, has been shown empirically to
gain additional savings when applied to the BWT form of the text.

## References

  * [1] Michael Burrows and David J. Wheeler.  A block-sorting lossless data compression algorithm.  Technical Report 124, Digital Equipment Corporation, 1994.
  * [2] Peter Elias.  Universal codeword sets and representations of the integers.  IEEE Trans. Information Theory, 21(2):194–203, 1975.
  * [3] Aharon Fruchtman, Yoav Gross, Shmuel T. Klein, and Dana Shapira.  Weighted adaptive coding.  CoRR, abs/2005.08232, 2020.
  * [4] Aharon Fruchtman, Yoav Gross, Shmuel T. Klein, and Dana Shapira.  Backward weighted coding.  In To appear in Data Compression Conference, DCC 2021. IEEE, 2021.
  * [5] Aharon Fruchtman, Shmuel T. Klein, and Dana Shapira.  Bidirectional adaptive compression.  In Proceedings of the Prague Stringology Conference 2019, pages 92–101, 2019.
  * [6] Shmuel T. Klein, Shoham Saadia, and Dana Shapira.  Forward looking Huffman coding.  Theory of Computing Systems, pages 1–20, 2020.
  * [7] Alistair Moffat.  Huffman coding.  ACM Comput. Surv., 52(4):85:1–85:35, 2019.
  * [8] Mark Nelson and Jean-Loup Gailly.  The Data Compression Book.  M & T Books, 1996.
  * [9] Jeffrey S. Vitter.  Design and analysis of dynamic Huffman codes.  JACM, 34(4):825–845, 1987.
  * [10] Jacob Ziv and Abraham Lempel.  A universal algorithm for sequential data compression.  IEEE Trans. Inf. Theory, 23(3):337–343, 1977.