# On the reliability of simulations of collisional stellar systems

Long Wang,1,2, David M. Hernandez 3
1Department of Astronomy, School of Science, The University of Tokyo, 7-3-1
Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan
2RIKEN Center for Computational Science, 7-1-26 Minatojima-minami-machi, Chuo-
ku, Kobe, Hyogo 650-0047, Japan
3Harvard–Smithsonian Center for Astrophysics, 60 Garden St., MS 51, Cambridge,
MA 02138, USA
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

It is well known that numerical errors grow exponentially in $N$-body
simulations of gravitational bound stellar systems, but it is not well
understood how the accuracy parameters of algorithms affect the physical
evolution in simulations. By using the hybrid $N$-body code, petar, we
investigate how escapers and the structure evolution of collisional stellar
systems (e.g., star clusters) depend on the accuracy of long-range and short-
range interactions. We study a group of simulations of ideal low-mass star
clusters in which the accuracy parameters are varied. We find that although
the number of escapers is different in individual simulations, its
distribution from all simulations can be described by Poisson statistics. The
density profile also has a similar feature. By using a self-consistent set-up
of the accuracy parameters for long- and short-range interactions, such that
orbits are resolved well enough, the physical evolution of the models is
identical. But when the short-range accuracy is too low, a nonphysical
dynamical evolution can easily occur; this is not the case for long-range
interactions. This strengthens the need to include a proper algorithm (e.g.
regularization methods) in the realistic modelling of collisional stellar
systems. We also demonstrate that energy-conservation is not a good indicator
to monitor the quality of the simulations. The energy error of the system is
controlled by the hardest binary, and thus, it may not reflect the ensemble
error of the global system.

###### keywords:

methods: numerical – software: simulations – star clusters: general

††pubyear: 2015††pagerange: On the reliability of simulations of collisional
stellar systems–On the reliability of simulations of collisional stellar
systems

##  1 Introduction

The star-by-star numerical $N$-body simulations of collisional stellar
systems, such as star clusters, are difficult due to the existence of compact
few-body systems. “Collisional” means that close encounters between stars
significantly affect the dynamical evolution of the system (e.g., Binney &
Tremaine, 1987). In a smooth potential like the Galactic potential, a star has
a certain orbit around the Galaxy; it will maintain the same orbit unless the
potential evolves. This is a collisionless system. In a collisional system,
close encounters perturb the motions of stars; after a while, the stars will
depart from their original orbits. Such a process is referred as relaxation.
Relaxation drives energy exchange among stars. If we use thermal dynamics to
describe the evolution of gravitational bound systems, we can find that the
heat capacity of the system is negative. As the energy transfers from the core
to the halo of the system, the core contracts while the temperature gradients
along the radial direction of the system increase. This subsequently
accelerates the contracting process and drives the core collapse. However,
when a binary exists in the core, it can encounter stars with separation
comparable to the binary apo-center distance. Then, it can transfer its
binding energy as a heating source. This binary heating prevents an infinite
collapse of the core. Therefore, binaries play a key role to control the long-
term evolution of the system. A realistic numerical simulation of a
collisional stellar system must accurately treat the dynamical evolution of
binaries.

However, as tight binaries shrink after interactions, their periods can become
as short as days, which must be resolved. Since the whole system evolves in a
Gyr timescale, the simulation can become very time-consuming. To solve this
issue, $N$-body methods need to introduce approximations to avoid expensive
computing. For example, the state-of-the-art code, nbody6 (Aarseth, 2003),
ignores the long-range tidal force to binaries. Thus, if a binary has no
neighbor, it is treated as an isolated binary with a pure Kepler orbit. Since
the phase information of the binary is not important for the dynamical
evolution of the system, nbody6 does not evolve the isolated binary until a
neighbor comes. This can significantly reduce the computing cost.

In a recently developed $N$-body code, petar, the slow-down method (Mikkola &
Aarseth, 1996; Wang, Nitadori & Makino, 2020) and the tidal-tensor method
(Wang et al., 2020) are used together. Since the cumulative tidal effect of
the whole system on a binary can be important, petar does not include isolated
binaries like nbody6. Instead, artificial particles are created around each
tight binary to measure the local tidal tensor in a certain time interval.
Then, the tensor is used to evaluate the long-range tidal force to the binary.
Meanwhile, short-range forces from neighbors are directly calculated. Since
there are no isolated binaries, the orbit of each short-period binary needs to
be integrated. The slow-down method is implemented to avoid expensive
computing. Depending on the ratio between the perturbation and the internal
force of a binary, a slow-down factor ($k>1$) is estimated. Then, the internal
motion of the binary is reduced by $k$ times, while the perturbation force is
increased by the same factor. It can be mathematically proved that the slow-
down method ensures a correct secular motion with a cost of losing the phase
information. In the example of a hyperbolic encounter in Wang, Nitadori &
Makino (2020), $k$ can be as large as $3\times 10^{4}$. Thus, the total
integration steps are reduced by a few orders of magnitude compared to that of
a conventional integrator.

These methods cannot solve the problem of the singularity due to the of
inverse square law of the Newtonian force. To follow the motion of an
eccentric binary, a traditional integrator, such as the fourth-order Hermite
method, needs to use an extremely small time step at the peri-center in order
to obtain an enough accuracy. This not only is computationally expensive, but
also leads to a large numerical error after many orbits. Such an error can
completely change the orbit of the binary and even unbind it. To solve this
issue, a group of algorithms, so called “regularization methods”, were
introduced. The famous ones are the KS regularization (Kustaanheimo & Stiefel,
1965), the Burdet-Heggie regularization (Burdet, 1967, 1968) and the time-
transformed symplectic leapfrog integrator (or algorithmic regularization;
Preto & Tremaine, 1999; Mikkola & Tanikawa, 1999). The coordinate
transformations are used in these algorithms to remove the singularity. Thus,
the efficiency and accuracy are significantly improved. However, coupling such
algorithms into $N$-body codes requires introducing a criterion to select the
particles for which the regularization applies. The force from nearby
particles to the regularized group are treated as an external perturbation.
The criterion needs to be well determined so that the perturbation is much
smaller than the internal forces of the binary. This is important for the KS
regularization method, which assumes a weak perturbation. Meanwhile, the start
and termination of the regularization introduces additional numerical errors.
Thus, the criterion also needs to ensure the switching is not frequent.
Determining a suitable criterion to ensure both computing efficiency and
accuracy of integration for all kinds of conditions is extremely challenging.
Neither nbody6 nor petar has a perfect solution.

There are only a few studies that have investigated how the numerical errors
affect the physical properties in the $N$-body simulations of collisional
stellar systems. Heggie (1991) studied how the formation of binaries, the
evolution of half-mass radius and the number of escapers depend on the
accuracy of the integrator. They used a force polynomial combined with
individual time steps (nbody1; Aarseth, 2003). They found that, in an $N$-body
model, when changing the parameter ($\eta$) to determine the time step, the
formation time of a binary, and the evolution of the half-mass radius remain
consistent. However, a noticeable difference in the number of escapers is
seen. Since nbody1 uses a softening potential to avoid singularities, close
encounters were not properly treated. This might affect the number of
escapers. Portegies Zwart & Boekholt (2014) investigates how the accuracy of
integration affects the few-body motions (the decay time) and found that small
errors in individual simulations can finish at completely different orbits due
to chaos, but the statistical mean of many repeating simulations can converge
to the correct result. Hernandez et al. (2020) investigated the long-term
evolution of planetary systems. They found that the statistics of action
variables are accurate as long as all orbits are properly resolved; e.g.,
$~{}16$ steps per effective pericenter period. Furthermore, Hernandez et al.
(2021) finds that ensemble averaging improves measurement of action-like
planetary quantities.

Energy conservation is frequently used to indicate the quality of the $N$-body
simulations. However, it is not a guarantee that good energy conservation
means a correct physical evolution (e.g. Heggie, 1991). In collisional stellar
systems, one tight binary can contain most of the energy in the entire system.
A small energy error from the binary can override any systematic error in the
system. Thus, an energy conservation check can be deceptive.

In this work, by performing a large amount of simulations for a low-mass star
clusters using petar code, we investigate how the statistical properties of
simulations, especially the escapers, depends on the accuracy of long-range
and short-range interactions. The results show that, although individual
models in the groups have significant differences, the distributions of the
models converge. In Section 2, we introduce the numerical algorithms used in
petar. In Section 2.1, we present the initial conditions of the models. Then,
a comparison for the evolution of the half-mass radius, the core radius, the
properties of escapers, and the energy errors are presented in Section 3.
Finally, we make conclude in Section 4.

##  2 Method

In this work, we use petar to carry out numerical simulations (Wang et al.,
2020). petar is a hybrid $N$-body code that combines three integration
methods:

  * •

The Barnes-Hut tree (Barnes & Hut, 1986) is used to calculate long-range
forces between particles, which are integrated with a second-order symplectic
leap-frog integrator (hereafter referred as “PT”).

  * •

The fourth-order Hermite integrator with block time steps (e.g., Aarseth,
2003, hereafter referred as “PP”) is applied to integrate the orbits of stars
and the centers-of-mass of multiple systems with short-range forces.

  * •

The slow-down algorithmic regularization method (SDAR; Wang, Nitadori &
Makino, 2020) is used to integrate the multiple systems, such as hyperbolic
encounters, binaries and hierarchical few-body systems.

When there are no binaries in the stellar system, the long-range forces are
the most expensive computing part in simulations. The direct pair-force
summation for all particles in the PP method requires a computing cost of
$O(N^{2})$, while the cost for the PT method scales as $O(N\log{N})$. However,
using pure PT with a second-order leap-frog integrator to accurately handle
the short-range interactions requires a small shared time step. This can be
even more expensive than the PP method for simulating the collisional stellar
systems. A hybrid algorithm, $\mathrm{P^{3}T}$, combines the PT and the PP
methods to include both their benefits of efficiency and accuracy. This is
done via the Hamiltonian splitting for a system with $N$ particles (e.g.,
Oshino, Funato, & Makino, 2011):

$\displaystyle H_{\mathrm{s}}=$
$\displaystyle\sum_{i=1}^{N}\frac{p_{i}^{2}}{2m_{i}}-\sum_{i<j}^{N}\frac{Gm_{i}m_{j}}{r_{ij}}W(r_{ij})$
(1) $\displaystyle H_{\mathrm{l}}=$
$\displaystyle\sum_{i=1,i<j}^{N}\frac{Gm_{i}m_{j}}{r_{ij}}[1-W(r_{ij})].$

where $H_{\mathrm{s}}$ includes the short-range interactions and kinetic
energies and $H_{\mathrm{l}}$ indicates the long-range interaction; $p_{i}$
and $m_{i}$ are, respectively, the momenta and mass of particle $i$; $r_{ij}$
is the separation between the $i$ and $j$ particles; $G$ is the gravitational
constant; $W(r_{ij})$ is a changeover function to smoothly transfer pieces
between $H_{\mathrm{s}}$ and $H_{\mathrm{l}}$.

petar uses a mass-dependent eighth-order polynomial changeover function for
each particle. The order of this polynomial can have a significant effect on
the accuracy of chaotic solutions (Hernandez, 2019a, b). In this work, we only
investigate equal-mass star clusters, thus the changeover function is the same
for all particles during one simulation. The simplified changeover function
for the equal-mass condition can be described as

$W(x)=\begin{cases}\beta(1-2x)&(x\leq 0)\\\ \beta(1-2x)-1+f(x)&(0<x<1)\\\
0&(x\geq 1)\\\ \end{cases}$ (2)

where

$\displaystyle f(x)$ $\displaystyle=1+\beta
x^{5}\left(14-28x+20x^{2}-5x^{3}\right)$ (3) $\displaystyle x$
$\displaystyle=\frac{r_{ij}-r_{\mathrm{in}}}{r_{\mathrm{out}}-r_{\mathrm{in}}}$
$\displaystyle\beta$
$\displaystyle=\frac{r_{\mathrm{out}}-r_{\mathrm{in}}}{r_{\mathrm{out}}+r_{\mathrm{in}}}.$

$r_{\mathrm{in}}$ and $r_{\mathrm{out}}$ are the inner and outer boundaries of
the changeover function, respectively. These are determined at the beginning
of the simulation. We fix $r_{\mathrm{out}}/r_{\mathrm{in}}$ to be $10$ in our
simulations, thus we only specify $r_{\mathrm{out}}$.

The PP method uses an individual time step for each particle. It is calculated
by (Aarseth, 2003; Oshino, Funato, & Makino, 2011),

$\Delta
t_{\mathrm{s,i}}=\mathrm{min}\left(\eta\sqrt{\tfrac{\sqrt{\left|\bm{A}^{(0)}_{i}\right|^{2}+A^{2}_{0}}\left|\bm{A}^{(2)}_{i}\right|+\left|\bm{A}^{(1)}_{i}\right|^{2}}{\left|\bm{A}^{(0)}_{i}\right|\left|\bm{A}^{(3)}_{i}\right|+\left|\bm{A}^{(2)}_{i}\right|^{2}}},\Delta
t_{\mathrm{s,max}}\right),\\\ $ (4)

where $\bm{A}^{(0)}_{i}$ is the acceleration of a particle ($i$),
$\bm{A}^{(j)}_{i}$ is its $j$-order time derivative, $A_{0}$ is a smoothing
parameter with the unit of acceleration. It avoids an unnecessarily small step
when all the other terms are close to zero. $\Delta t_{\mathrm{s,max}}$ is the
maximum time step. In our simulations, $\eta=0.1$, a choice used in previous
works (e.g., Aarseth, 2003), and $A_{0}=0.1\langle
m\rangle/r_{\mathrm{out}}^{2}$, where $\langle m\rangle$ is the local average
stellar mass.

To switch on the SDAR method, a mass-dependent radial criterion,
$r_{\mathrm{g}}$, is used to select the members of multiple systems. In this
work $r_{\mathrm{g}}$ is uniform for all particles and
$r_{\mathrm{g}}<r_{\mathrm{out}}$.

To deal with large $N$-body simulations, petar is developed by using the
parallelization framework for developing particle simulation codes (fdps;
Iwasawa et al., 2016; Iwasawa, et al., 2020). fdps provides the parallel-
computing support for the PT part by using MPI and OpenMP software. petar also
supports the use of special accelerators, such as GPUs (using the CUDA
language), SIMD (AVX, AVX2, AVX512) and Fujitsu A64FX (in the Fugaku
supercomputer), to speed up the force calculation.

The hybrid integration method of petar is very suitable for this research
because we can distinguish the impacts of long-range relaxation and short-
range close encounters on the evolution of collisional stellar systems.

###  2.1 Initial conditions

In this work, we use the same initial conditions of a star cluster for all
simulations. The cluster contains $1000$ equal-mass particles. The positions
and the velocities of particles are randomly sampled from the Plummer
distribution (Plummer, 1911). Initially, the system is in virial equilibrium.
We adopt Hénon (1971) units in the simulations (this is frequently referred to
as the $N$-body unit, hereafter we use the abbreviation NB unit). In these
units, the total mass of the cluster ($M$) is unity, the total kinetic energy
is $0.25$, and the total potential energy is $-0.5$. The gravitational
constant is $G=1$. The initial virial radius is one. The corresponding initial
half-mass radius ($R_{\mathrm{h,0}}$) is approximately $0.768$. The initial
crossing time is approximately $0.674$, estimated by using

$T_{\mathrm{cr}}=\sqrt{\frac{R_{\mathrm{h}}^{3}}{GM}}.$ (5)

The initial relaxation time is approximately 31.038, calculated by using the
formula from Spitzer (1987),

$T_{\mathrm{rh}}\approx
0.138\frac{N^{1/2}R_{\mathrm{h}}^{3/2}}{m^{1/2}G^{1/2}\ln\Lambda},$ (6)

where $m=0.001$ and $\Lambda=0.02N$. The factor of $0.02$ follows the
measurement from Giersz & Heggie (1996). For the equal-mass Plummer model,
core collapse happens after approximately $15T_{\mathrm{rh}}$ (Heggie & Hut,
2003), which corresponds to approximately 466 NB time units for our models.
Thus, we stop the simulations at 500 NB time unit to cover the core collapse
process. When the energy of a particle becomes positive and its distance to
the cluster center exceeds 10 NB units, the particle is treated as an escaper.

We carry out 7 groups of models by varying the tree time step ($\Delta
t_{\mathrm{l}}$) of the PT method, $r_{\mathrm{out}}$, and $r_{\mathrm{g}}$.
$\Delta t_{\mathrm{l}}$ represents the accuracy of long-range interactions.
$r_{\mathrm{out}}$ and $r_{\mathrm{g}}$ together represent the accuracy of
short-range interactions. Smaller $r_{\mathrm{out}}$ and $r_{\mathrm{g}}$
indicate low accuracy because less particles are treated with the Hermite or
SDAR method.

Table 1 shows the parameters of each group. For the names of groups, the
prefixes, “L”, “S” and “H”, in the front of “PT” and “PP” indicate the low,
standard, and high accuracy, respectively. The suffix “-C” indicates that the
model follows the automatic determination of $r_{\mathrm{out}}$ as (Equation
41 in Wang et al., 2020):

$r_{\mathrm{out}}=10\Delta t_{\mathrm{l}}\sigma_{\mathrm{1D}},$ (7)

where $\sigma_{\mathrm{1D}}$ is the one-dimensional velocity dispersion of the
system. This relation ensures that $\Delta t_{\mathrm{l}}$ is small enough to
resolve the motion of particles inside the changeover region. Therefore, in
the LPT-HPP-C group, the long-range interactions have a low accuracy (a large
$\Delta t_{\mathrm{l}}$). In balance, the short-range interactions have a high
accuracy ($r_{\mathrm{out}}$ is about $12.8\%$ the initial virial radius).

Meanwhile, for the HPP and SPP groups, $r_{\mathrm{g}}=0.8r_{\mathrm{in}}$.
This value ensures that the SDAR method is used for most multiple systems
inside $r_{\mathrm{in}}$ while the Hermite integrator mainly handles the
interactions inside the changeover region. This choice can best avoid an
abnormal orbital evolution of binaries due to the cumulative error of the
Hermite integrator. In the LPP case, $r_{\mathrm{g}}$ is set to a small value
$10^{-8}$ so that most wide multiple systems and encounters are not integrated
by the SDAR method. We cannot set $r_{\mathrm{g}}$ to zero because the
formation of a tight binary can result in a very small $\Delta
t_{\mathrm{s}}$, and thus, the simulation becomes extremely slow.

Table 1: The accuracy controlling parameters for each group of the $N$-body models. $\Delta t_{\mathrm{l}}$ is tree time step, $r_{\mathrm{out}}$ is the outer boundary of the changeover function. $r_{\mathrm{g}}$ is the radial criterion to switch on the SDAR method. All values are in the NB unit. Groups | STD | LPT-SPP | LPP
---|---|---|---
Model Name | SPT-SPP-C | LPT-HPP-C | HPT-SPP | LPT-SPP | LPT-LPP | SPT-LPP | HPT-LPP
$\Delta t_{\mathrm{l}}$ | $1/128$ | $1/32$ | $1/1024$ | $1/32$ | $1/32$ | $1/128$ | $1/1024$
$r_{\mathrm{out}}$ | 0.03 | 0.128 | 0.03 | 0.03 | 0.0001 | 0.0001 | 0.0001
$r_{\mathrm{g}}$ | 0.0024 | 0.01024 | 0.0024 | 0.0024 | $1\times 10^{-}8$ | $1\times 10^{-}8$ | $1\times 10^{-}8$

The 7 groups can be divided into three classes. The SPT-SPP-C, LPT-HPP-C and
HPT-SPP groups have self-consistent (satisfying the criterion determined by
Equation 7) accuracy parameters for long- and short-range interactions. We
call them as the STD class. The LPT-LPP, SPT-LPP and HPT-LPP groups have a low
accuracy of short-range interactions. We call them as the LPP class. The LPT-
SPP has the same accuracy parameters for short-range interactions as those in
SPT-SPP-C and HPT-SPP groups, but it uses low accuracy parameters for long-
range interactions. Thus we separate it into the third class.

For each group, we carry out $300$ simulations by only varying the values of
$r_{\mathrm{out}}$ via adding a 3 digit number. For example, in the SPT-SPP-C
group, the values of $r_{\mathrm{out}}$ for the first, the second and the last
simulations are $0.03$, $0.03001$, and $0.03299$, respectively. In this way,
we ensure that the initial conditions are exactly the same but only the
accuracy parameter of the integration is slightly perturbed. Since the
numerical error grows exponentially in a crossing time (Heggie, 1991), we can
observe a large divergence of the individual orbits of stars soon after the
simulations start. Thus, we can investigate how sensitively the evolution of
the star cluster depends on the accuracy parameter.

##  3 Results

###  3.1 Structure evolution

Figure 1 shows the evolution of the half-mass radius ($R_{\mathrm{h}}$) and
the core radius ($R_{\mathrm{c}}$) for all groups. The left and right panels
compare the STD/LPT-SPP and the LPP classes. The STD class shows an identical
evolution of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$. The core collapse finishes
around 300 time units. After that, $R_{\mathrm{h}}$ starts to increase. The
HPT-SPP group has a better accuracy of long-range interactions compared to
that of the SPT-SPP-C group; this suggests that when $\Delta t_{\mathrm{l}}$
and $r_{\mathrm{out}}$ are set in a self-consistent way, the statistical
result can converge and does not depend on the accuracy parameters.

The LPT-SPP group shows a faster core collapse that finishes at approximately
100 time units. In addition, this group shows a strong decrease of
$R_{\mathrm{h}}$ before the end of core collapse.

In contrast, the LPP class shows a completely different evolution of
$R_{\mathrm{h}}$ and $R_{\mathrm{c}}$. The whole system expands immediately
after the simulations start, and there is no core collapse. The expansion
indicates that the system is not in virial equilibrium. Such behaviour may be
caused by the artificial energy heating coming from the integration errors of
close encounters. This indicates that the star cluster has a nonphysical
evolution. For example, in the STD class, before the tight binary, which acts
as a heating source, forms, virial equilibrium is the natural state of the
star cluster; i.e., $R_{\mathrm{h}}$ does not significantly evolve. If the
cluster is not in virial equilibrium, it will evolve to the virial equilibrium
via phase mixing and violent relaxation in a free-fall timescale (e.g. Binney
& Tremaine, 1987).

This result suggests that without an accurate integration of short-range
interactions, the simulation is completely unreliable. The accuracy of long-
range interactions is less crucial but affects the evolution of core and
$R_{\mathrm{h}}$.

Figure 1: The evolution of the half-mass radius ($R_{\mathrm{h}}$) and the
core radius ($R_{\mathrm{c}}$) for all groups. The left panels include the STD
and the LPT-SPP classes. The right panels include the LPP class. The central
curves are the averaged values of all models in each group. The shaded regions
show the standard deviations.

###  3.2 Number of escapers

When the energy of a star becomes positive after encounters, the star becomes
a potential escaper. Then, it travels through the halo to escape, but during
travel, distant perturbations can still change its orbit and the star may
again become bound. Thus, the number of escapers ($N_{\mathrm{esc}}$) is
sensitive to both long- and short-range interactions. In Figure 2, we show the
probability of $N_{\mathrm{esc}}$ of all models in the SPT-SPP-C group.
Similar to the evolution of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$, a large
scatter of $N_{\mathrm{esc}}$ appears among the models. However, the
probability of $N_{\mathrm{esc}}$ approximately follows a Poisson distribution
with an expected value of $\lambda=123$. This suggests that the behaviour of
escapers has no correlation with the accuracy parameters; the difference of
$N_{\mathrm{esc}}$ is purely due to a random effect. Therefore,
$N_{\mathrm{esc}}$ of individual models is not reliable; to understand the
properties of escapers in a low-mass star cluster, it is necessary to obtain
ensemble statistics. Hernandez et al. (2020) found a similar behaviour for the
phase space structure of planetary dynamics.

Figure 2: The probability of $N_{\mathrm{esc}}$ for all models in the SPT-
SPP-C group. A fit of the Poisson distribution is shown in the brown color.
The expected value is $\lambda=123$.

Subsequently, we compare the distribution of $N_{\mathrm{esc}}$ from different
groups in Figure 3. The three groups in the STD class well agree with each
other. This is consistent with the behaviour of $R_{\mathrm{h}}$ and
$R_{\mathrm{c}}$. The LPT-SPP group has a slightly larger mean
$N_{\mathrm{esc}}$. In contrast, the LPP class with a low-accuracy short-range
interaction and different accuracy long-range interactions (different $\Delta
t_{\mathrm{l}}$) shows a much larger mean $N_{\mathrm{esc}}$. This indicates
that $N_{\mathrm{esc}}$ is much more sensitive to the accuracy of short-range
interactions. Thus, the accurate treatment of close encounters and multiple
systems is crucial for simulating collisional stellar systems.

Figure 3: The distribution of $N_{\mathrm{esc}}$ for the models in individual
groups.

###  3.3 Kinetic energy of escapers

To investigate why the low-accuracy short-range interactions cause a large
number of escapers, we plot the distribution for the kinetic energies of
escapers in Figure 4. The STD class shows the same distribution again. Most
escapers have a low kinetic energy with a peak around $5\times 10^{-5}$. There
is a long and inconspicuous tail of high-energy
($e_{\mathrm{kin,esc}}>10^{-3}$) escapers. The LPT-SPP group has an almost
identical distribution. Thus, the low accuracy of long-range interactions does
not have an obvious effect on the energy distribution of escapers.

In contrast, the LPP class has a very different distribution with a bimodal
shape. There are more low-energy escapers with $e_{\mathrm{kin,esc}}<10^{-4}$
and a large cluster of high-energy escapers with
$10^{-4}<e_{\mathrm{kin,esc}}<10^{-2}$. The peak of the high-energy part
shifts to a higher value when the accuracy of long-range interactions
increases.

Figure 4: The kinetic energy distribution of escapers. The histograms show the
averaged value of all models in each group. The standard deviations are shown
as the error bars. The upper panel includes the STD and the LPT-SPP groups.
The lower panel includes the LPP class.

In order to understand the origin of the bimodal distribution forms, we
analyze how $N_{\mathrm{esc}}$ depends on the time as shown in Figure 5. For
the STD and the LPT-SPP classes, only a few low-energy escapers appear
initially, and then, their number increases. After 300 time units, the high-
energy escapers start to appear. This is due to core collapse, which leads to
the formation of tight binaries. These binaries can generate high-energy
escapers via the decays of unstable multiple systems. Compared to the STD
class, the low- and middle-energy escapers appear slightly earlier in the LPT-
SPP group because of an earlier core collapse.

The LPP class (right panel) shows a completely different behaviour. A large
number of escapers appear immediately after the simulations start. This is the
case for all energy bins in the LPT-LPP group, which has the lowest accuracy
for both the long-range and the short-range interactions. These initial
escapers are corresponding to the right peaks in the distribution of
$e_{\mathrm{kin,esc}}$ (Figure 4). When $\Delta t_{\mathrm{l}}$ decreases, the
numbers of low- and middle-energy escapers decrease. As shown in Figure 1, the
LPP class has immediately left the virial equilibrium at the beginning. The
integration errors of close encounters lead to a boost of escapers. A higher
accuracy of long-range interactions help somewhat.

Meanwhile, if we ignore the initial escapers, the numbers of middle- and high-
energy escapers decrease in the later evolution, in contrast to those in the
STD class. There is no core collapse while the density continues to decrease
in the LPP class. Thus, such behaviour is expected. The low-energy escapers
continue to appear and cause the left peak in the distribution of
$e_{\mathrm{kin,esc}}$.

Figure 5: The number of escapers depends on the time: $N_{\mathrm{esc}}(t)$.
The escapers are collected into to three (low, middle and high) kinetic-energy
bins. The left panels include the STD and the LPT-SPP classes. The right
panels include the LPP class. The step curves are the mean
$N_{\mathrm{esc}}(t)$ of all models in each group. The errorbars are the
standard deviations.

###  3.4 Energy error

The total energy in a Newtonian gravitational $N$-body system is conserved.
Thus, the energy conservation is frequently used to measure the fidelity of
$N$-body models. However, this does not work well for collisional stellar
systems. If a tight binary forms, based on the Heggie-Hills law (Heggie, 1975;
Hills, 1975), the maximum semi-major axis can be described as

$a_{\mathrm{m}}\approx\frac{Gm_{\mathrm{1}}m_{\mathrm{2}}}{\langle
m\rangle\sigma^{2}},$ (8)

where $m_{\mathrm{1}}$, $m_{\mathrm{2}}$, and $\langle m\rangle$ are the
masses of the two components and the locally averaged stellar mass,
respectively; $\sigma$ is the local velocity dispersion. In our models,
$m_{\mathrm{1}}m_{\mathrm{2}}/\langle m\rangle=1.0$ and $\sigma\approx 1$
after core collapse for the STD class. Thus, $a_{\mathrm{m}}\approx 1$ and the
corresponding binding energy of the binary is $0.5$. This is close to the
initial total potential energy of the whole cluster. After many encounters,
the binding energy can even become larger. Then, the total energy of the
system is completely dominated by this binary. The energy error from the
binary can mask all other small energy errors in the system.

Similarly, when there are still no binaries, the error is dominated by close
encounters with a pair separation below $a_{\mathrm{m}}$. If a single
encounter is not accurately followed, it can create a large energy error that
does not affect the dynamical state of the whole system significantly.
However, if every encounter is treated inaccurately, each with small energy
error, the error may not add coherently. The cumulative error may appear small
even when the evolution is nonphysical.

Figure 6 shows the averaged energy error per time units
($\langle\mathcal{E}\rangle$) vs. $N_{\mathrm{esc}}$ for all groups. The three
classes show different properties of $\langle\mathcal{E}\rangle$. The STD
class has a relatively low mean $\langle\mathcal{E}\rangle$, and some models
in the HPT-SPP group have the lowest $\langle\mathcal{E}\rangle$. The LPT-SPP
group shows a narrow distribution of $\langle\mathcal{E}\rangle$ with the mean
value close to the upper boundary of the STD class. The LPP class has the
largest values of $\langle\mathcal{E}\rangle$. For individual models in each
group, there is no obvious correlation between $N_{\mathrm{esc}}$ and
$\langle\mathcal{E}\rangle$. Especially, in the LPP class, the mean
$N_{\mathrm{esc}}$ among the groups vary widely, while the distributions of
$\langle\mathcal{E}\rangle$ are similar; some models in the STD class have a
large $\langle\mathcal{E}\rangle$ similar to those in the LPP class.
Therefore, $\langle\mathcal{E}\rangle$ cannot explain whether a system has a
correct physical evolution.

Figure 6: The average energy error per time unit ($\langle\mathcal{E}\rangle$)
vs. $N_{\mathrm{esc}}$ in each model from all groups. The upper and right
panels show the corresponding histograms.

##  4 Discussion and conclusion

In this work, we carry out 7 groups of $N$-body simulations of an equal-mass
star cluster with 1000 particles by using the petar code. We investigate how
the accuracy parameters of long-range ($\Delta t_{\mathrm{l}}$) and short-
range interactions ($r_{\mathrm{out}}$ and $r_{\mathrm{g}}$) affect the
dynamical evolution of the star cluster. All simulations use the exact same
initial condition. In each group, we perform 300 integrations with slight
modifications to $r_{\mathrm{out}}$ to obtain different starting numerical
errors. We analyze how the structure evolution and the properties of escapers
are related to the accuracy parameters of the integrator.

The dynamical evolution, especially the relaxation process, of the clusters
are affected by both short-range interactions (close encounters and multiple
systems) and long-range interactions. If the accuracy parameters are set in a
self-consistent way (the STD class; see Equation 7), the statistical
properties of the evolution are independent of the accuracy parameters. When
the accuracy is not sufficient for long-range interactions (the LPT-SPP group)
compared to that of the STD class, the core collapse occurs earlier and the
number of escapers $N_{\mathrm{esc}}$ is slightly larger. However, when the
short-range interaction is inaccurate (the LPP class), the clusters have
nonphysical evolution: the clusters evolve out of virial equilibrium
immediately after the simulations start, i.e., there is a strong expansion of
$R_{\mathrm{h}}$. Since the phase mixing and the violent relaxation drive a
stellar system into virial equilibrium, a contrary evolution is not correct.
In addition, no core collapse occurs and $N_{\mathrm{esc}}$ is significantly
larger than that of the STD group. Therefore, realistic $N$-body simulations
of collisional stellar systems must accurately treat short-range interactions.

Although individual models in one group have a large difference of
$R_{\mathrm{h}}$, $R_{\mathrm{c}}$ and $N_{\mathrm{esc}}$, the distributions
of them in a group converge. For example, the distribution of
$N_{\mathrm{esc}}$ follows the Poisson distribution function (Figure 2). Thus,
although an initial small difference can grow exponentially so that individual
$N$-body models cannot reproduce the exactly same evolution of a real cluster,
the statistical means are meaningful. A similar result in the orbits of
planets is found in Hernandez et al. (2020).

The energy error has a large scatter among the models, but the physical
properties of the systems do not have a clear dependence on
$\langle\mathcal{E}\rangle$ (Figure 6). This is reasonable because in
collisional stellar systems, the integration errors from tight binaries or
close encounters dominate the total energy error and dominate over the
systematic error generated by the nonphysical evolution. Thus, checking the
energy conservation is not the proper way to measure the quality of $N$-body
simulations for collisional stellar systems. It is much more important to have
a proper algorithm to correctly treat the short-range interactions rather than
achieving a good energy conservation by using an approximation (e.g., a
softening potential).

Hernandez et al. (2020) found that the relevant metric for accuracy was the
number of steps per bound orbit, rather than the energy error. They found that
orbits needed at minimum about $16$ steps per effective period at pericenter,
a number which depends on eccentricity. This number holds for a Wisdom–Holman
map (Wisdom & Holman, 1991), but nonetheless, we assume this approximate
criteria applies to other maps. Using a maximum of $16$ integration steps per
orbit, the minimum semi-major axes of a binary ($a_{\mathrm{min}}$) that the
long-range PT integrator can resolve are approximately $0.002$, $0.009$, and
$0.02$ NB length units for the HPT, SPT and LPT groups, respectively (see
Table 1). In the HPT-LPP group, $r_{\mathrm{out}}=0.0001$, which is $5$ times
less than $a_{\mathrm{min}}$. The situation in the other two groups of the LPP
class are even worse. However, the three groups in the STD class have
$r_{\mathrm{out}}>a_{\mathrm{min}}$. Thus, our result is consistent with what
was found in Hernandez et al. (2020).

In this work, we only investigate the low-$N$ star clusters, which have a
large statistical scatter. In the case of massive stellar systems such as
globular clusters with million stars, the statistical scatter is expected to
be much smaller. Thus, $N_{\mathrm{esc}}$ in one model would be closer to the
statistically expected value.

Our analysis is based on the specific algorithms used in petar. We refer the
PT and PP to the long-range and short-range interaction, respectively. This
does not represent all type of $N$-body algorithms. With different
integrators, the results can change, although the general trend is not
expected to vary significantly.

## Acknowledgements

L.W. thanks the financial support from JSPS International Research Fellow
(School of Science, The university of Tokyo). Numerical computations were in
part carried out on Cray XC50 at Center for Computational Astrophysics,
National Astronomical Observatory of Japan. We thank discussions with Chris
Hamilton.

## Data Availability

The $N$-body simulations underlying this article were generated by using the
petar code on the desktop computer of the corresponding author and the
supercomputer, Cray XC50 at Center for Computational Astrophysics (CfCA),
National Astronomical Observatory of Japan. The data underlying this article
will be shared on reasonable request to the corresponding author.

## References

  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations, Cambridge University Press
  * Barnes & Hut (1986) Barnes J., Hut P., 1986, Natur, 324, 446
  * Binney & Tremaine (1987) Binney J., Tremaine S., 1987, Princeton, NJ, Princeton University Press
  * Burdet (1967) Burdet C. A., 1967, ZaMP, 18, 434
  * Burdet (1968) Burdet C. A., 1968, ZaMP, 19, 345
  * Giersz & Heggie (1996) Giersz M., Heggie D. C., 1996, MNRAS, 279, 1037
  * Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
  * Heggie (1991) Heggie D. C., 1991, ASIB, 272, 47
  * Heggie & Hut (2003) Heggie D., Hut P., 2003, the Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics, Cambridge University Press, 2003
  * Hénon (1971) Hénon M., 1971, Ap&SS, 13, 284
  * Hernandez (2019a) Hernandez David M., 2019, MNRAS, 486, 5231
  * Hernandez (2019b) Hernandez David M., 2019, MNRAS, 490, 4175
  * Hernandez et al. (2020) Hernandez, David M., Hadden, S., Makino, J., 2020, MNRAS, 493, 1913
  * Hernandez et al. (2021) Hernandez D. M., Agol E., Holman M. J., Hadden S., 2021, RNAAS, 5, 77. doi:10.3847/2515-5172/abf4e3
  * Hills (1975) Hills J. G., 1975, AJ, 80, 809
  * Iwasawa et al. (2016) Iwasawa M., Tanikawa A., Hosono N., Nitadori K., Muranushi T., Makino J., 2016, PASJ, 68, 54
  * Iwasawa, et al. (2017) Iwasawa M., Oshino S., Fujii M. S., Hori Y., 2017, PASJ, 69, 81
  * Iwasawa, et al. (2020) Iwasawa M., Namekata D., Nitadori K., Nomura K., Wang L., Tsubouchi M., Makino J., 2020, PASJ, 72, 13
  * Kustaanheimo & Stiefel (1965) Kustaanheimo P., Stiefel E., 1965, J. Reine Angew. Math., 218, 204
  * Mikkola & Aarseth (1996) Mikkola S., Aarseth S. J., 1996, CeMDA, 64, 197
  * Mikkola & Tanikawa (1999) Mikkola S., Tanikawa K., 1999, MNRAS, 310, 745
  * Oshino, Funato, & Makino (2011) Oshino S., Funato Y., Makino J., 2011, PASJ, 63, 881
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Preto & Tremaine (1999) Preto M., Tremaine S., 1999, AJ, 118, 2532
  * Portegies Zwart & Boekholt (2014) Portegies Zwart S., Boekholt T., 2014, ApJ, 785, LL3
  * Spitzer & Hart (1971) Spitzer L., Hart M. H., 1971, ApJ, 164, 399
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters, Princeton University Press
  * Wang, Nitadori & Makino (2020) Wang L., Nitadori K., Makino J., 2020, MNRAS, 493, 3398
  * Wang et al. (2020) Wang L., Iwasawa M., Nitadori K., Makino J., 2020, MNRAS, 497, 536
  * Wisdom & Holman (1991) Wisdom J., Holman M., 1991, AJ, 102, 1528. doi:10.1086/115978


Ali Taani

# Accretion Induced Collapse of White Dwarfs as an Alternative Symbiotic
Channel to Millisecond Pulsars

###### Abstract

Recently, extra motivation has been given to the investigations of an
unresolved problem of millisecond pulsars (MSPs) produced by the recycling
process, as an apparent role of the accretion-induced collapse (AIC) in white
dwarfs (WDs) was suggested to this concern. I have found that the distribution
of the orbital periods of binary MSPs in the Galactic disk ($N_{obs,orb}$)
closely follows an exponential distribution. I have also determined the best-
fit mean value of $N_{obs,orb}$ by fitting our data with an exponential
distribution for the MSP population. As a result, it can be stated that
reaching the Chandrasekhar limit may cause an explosion of a massive WD as a
Type Ia supernova (in the case of a CO WD) or an ignition of a ONeMg WD, and
possibly merging in some CO WDs, all resulting in peculiar MSP systems. A
possible formation scenario, where the system has a circular orbit during this
evolutionary stage, is discussed.

##  1 Introduction

Millisecond pulsars (MSPs) are remarkable objects born via Type II supernovae
(SNe) explosions, they are characterized by short spin periods
(P${}_{spin}\leq$ 20 ms), weak magnetic fields (B $\leq 10^{9}$ G), and
extremely old ages $\sim 10^{9}$ yr, based on the recycling process. The ATNF
catalog counts 500 Galactic MSPs both in the Galactic plane and in the
globular clusters [1, 2]. They are often found in binaries (about 53%) with
white dwarfs (WDs) in circular orbits, having companions with masses of
$\rm\sim 0.15M_{\odot}-0.45M_{\odot}$ (see e.g., [3, 4, 5]).

It has been suggested that MSPs form in low-mass X-ray binary systems. This
argument leads to the recycling process, in which a slowly rotating old
neutron star (NS) may be spun up into an MSP via accretion from a binary
companion [3, 6, 7].

An alternative formation process is an accretion-induced collapse (AIC) of an
ONeMg white dwarf (see, e.g., [8, 9, 10, 11, 12, 13, 14]). In this scenario,
the mass transfer can increase the mass of the WD up to the Chandrasekhar
limit. This causes its collapse and a violent release of gravitational energy,
which might be observable by gravitational-wave observatories such as LIGO,
VIRGO, and GEO [15].

---
Figure 1: Left: a histogram of the orbital periods of binary MSPs in the
Galactic disk (white) and in the globular clusters (black). The solid line
denotes the best fit to the data. Right: the cumulative probability
distribution of the observed data for the Galactic disk (solid line) and the
globular clusters (dashed line). The data were taken from the ATNF catalog
([1] and the Paulo Freire’s GC Pulsar Catalog (www.naic.edu/
pfreire/GCpsr.html). Figure 2:  Final eccentricity $\rm e_{f}$ as a function
of $M_{b}$.

It is worth noting that in contrast to the Galactic field, globular clusters
may provide distinct routes to the formation and evolution of accreting
binaries. As a result, the tidal capture process [16, 17] could be a source of
long orbital period systems. In addition, there is a possibility of the
formation of such long-period binaries through mass exchange or mergers [18,
19, 20, 21]. More observations in globular clusters showed that some fraction
of NSs must have received low kicks at birth. Fig. 1 shows the distribution of
some representative observed data of the orbital periods of binary MSP systems
in the Galactic disk and globular clusters (left) and their cumulative
probability distribution (right). These parameters could play a role in
examining the theoretical scenarios of the formation and evolution of recycled
MSPs. The cumulative probability distribution is best described by a log-
normal law with $\sigma=0.33$, $R^{2}=92\%$, and $\chi^{2}=0.98$. In Fig. 2,
we can see that the eccentricity $e$ is inversely proportional to the binary
system mass $M_{b}$. It turns out that the low eccentricity is significant for
the solutions presented for each system, as the circular orbit approximation
is used. This generally reflects a substantial spread in orbits, as in Figs. 2
& 3.

Figure 3: The relation between masses and orbital periods, showing that
massive WDs recycle to MSPs. The orbital periods show a clear difference (over
two orders of magnitude). It should be noticed that the AIC process leads to
dynamical interactions in the systems during the conversion of a WD into an
MSP. The data for MSPs and WDs are taken from the ATNF catalogue [1] and from
[22] respectively.

##  2 The Evolutionary Path from WDs to MSPs

Fig. 4 is a simple flowchart of the currently favored model that is used to
explain the formation of various types of the systems. The left column
illustrates a binary system in which a more massive star has typically a mass
in the range of $8-11M_{\odot}$ and produces a WD after exhausting its
available thermonuclear fuel.

Later, as the binary system sheds its angular momentum through the wind, the
companion eventually fills its Roche lobe and begins accreting onto the WD.
The accretion causes the mass of the WD to steadily increase until it
approaches the Chandrasekhar limit. Thereafter, the WD becomes unstable and
undergoes an AIC, forming an NS, or detonates completely as a Type Ia SN in a
thermonuclear explosion. Our major assumption is that the WD–MS channel for
producing MSPs is more common than the WD-NS channel, because (1) not all WDs
have the right composition (O, Ne, Mg), and not all of them accumulate a
sufficient mass ($\rm\Delta m=0.1-0.2M_{\odot}$) to undergo the AIC [23]; (2)
in the WD–NS case the progenitor star of the NS is more massive and hence may
bring the binary into the Roche lobe overflow from a larger original
separation.

Figure 4: A diagram illustrating possible binary evolution scenarios leading
to the formation of WD–MSP systems in the Galactic disk.

In the right column of Fig. 4, the companion star is massive enough so that it
explodes directly as a Type II SN, producing an NS. The binary does not stay
intact if more than half the total pre-supernova mass is ejected from the
system during the explosion [24, 25]. After a time of $10^{7}-10^{8}$ yrs, the
old spun-down NS can gain a new lease on life as a pulsar by accreting matter
and angular momentum at the expense of the orbital angular momentum of the
binary system. During this accretion phase, X-rays are emitted by the
frictional heating of the matter infalling to the NS. This makes the system
visible as an X-ray binary. At the end of this spin-up phase, the secondary
sheds its outer layers to become a WD in an orbit around a now rapidly
spinning MSP.

##  3 Conclusions

I have attempted to follow the AIC scenario channel in the considered
evolutionary stage of close binary MSPs since some circular binaries are found
both in the Galactic disk and in globular clusters. I present and analyze a
whole set of terminal evolution tracks from WDs to the MSP stage that are
suitable for the study of their formation and evolution including the mass
exchange phases. With a focus on the combination of stellar parameters such as
mass, Porb, mass accretion rate, and eccentricity, I argue that our data
appear to be more consistent with this scenario. However, knowledge of these
tracks is required to study the cosmological perspective and any effects on
galaxy evolution models. In addition, the existence of MSPs, WDs, and SNe Ia
has expanded the possible evolutionary channels for the formation of near-
Chandrasekhar-mass models and MSPs appearing via the AIC process. In
forthcoming work, I will expand on this research with studies aimed at delving
deeper into these issues.

## References

  * [1] R.N. Manchester, G.B. Hobbs, A. Teoh and M. Hobbs, _The Australia Telescope National Facility Pulsar Catalogue_ , _Astron. J._ 129 (2005) 1993 [astro-ph/0412641].
  * [2] L. Ferrario and D.T. Wickramasinghe, _Binary Millisecond Pulsars: The Accretion Induced Collapse hypothesis revisited_ , in _Astrophysics of Compact Objects_ , Y.-F. Yuan, X.-D. Li and D. Lai, eds., vol. 968 of _American Institute of Physics Conference Series_ , pp. 194–196, Jan., 2008, DOI.
  * [3] M.A. Alpar, A.F. Cheng, M.A. Ruderman and J. Shaham, _A new class of radio pulsars_ , _Nature_ 300 (1982) 728.
  * [4] D. Bhattacharya and E.P.J. van den Heuvel, _Formation and evolution of binary and millisecond radio pulsars_ , _Phys. Rep._ 203 (1991) 1.
  * [5] A. Taani, J.C. Vallejo and M. Abu-Saleem, _Assessing the complexity of orbital parameters after asymmetric kick in binary pulsars_ , _Journal of High Energy Astrophysics_ 35 (2022) 83 [2206.11015].
  * [6] A. Taani, _Systematic comparison of initial velocities for neutron stars in different models_ , _Research in Astronomy and Astrophysics_ 16 (2016) 101.
  * [7] A. Taani, A. Abushattal and M.K. Mardini, _The regular dynamics through the finite-time Lyapunov exponent distributions in 3D Hamiltonian systems_ , _Astronomische Nachrichten_ 340 (2019) 847.
  * [8] K. Nomoto, _Evolution of 8–10 M sun Stars toward Electron Capture Supernovae. II. Collapse of an O + NE + MG Core_, _Astrophys. J._ 322 (1987) 206.
  * [9] A. Taani, C.M. Zhang, M. Al-Wardat and Y.H. Zhao, _Where do the progenitors of millisecond pulsars come from?_ , _Astronomische Nachrichten_ 333 (2012) 53 [1112.1312].
  * [10] A. Taani, C. Zhang, M. Al-Wardat and Y. Zhao, _Investigation of some physical properties of accretion induced collapse in producing millisecond pulsars_ , _Astrophys. Space Sci._ 340 (2012) 147 [1201.3779].
  * [11] A. Taani and A. Khasawneh, _Probing the accretion induced collapse of white dwarfs in millisecond pulsars_ , in _Journal of Physics Conference Series_ , vol. 869 of _Journal of Physics Conference Series_ , p. 012090, June, 2017, DOI [1901.00547].
  * [12] T.M. Tauris, M. Kramer, P.C.C. Freire, N. Wex, H.T. Janka, N. Langer et al., _Formation of Double Neutron Star Systems_ , _Astrophys. J._ 846 (2017) 170 [1706.09438].
  * [13] I. Ablimit, P. Podsiadlowski, R. Hirai and J. Wicker, _Stellar core-merger-induced collapse: new formation pathways for black holes, Thorne-Żytkow objects, magnetars, and superluminous supernovae_ , _Mon. Not. R. Astron. Soc._ 513 (2022) 4802 [2108.08430].
  * [14] M.K. Mardini, A. Frebel, R. Ezzeddine, A. Chiti, Y. Meiron, A.P. Ji et al., _The chemical abundance pattern of the extremely metal-poor thin disk star 2MASS J1808-5104 and its origins_ , _Mon. Not. R. Astron. Soc._ (2022)  [2208.03891].
  * [15] S. Darbha, B.D. Metzger, E. Quataert, D. Kasen, P. Nugent and R. Thomas, _Nickel-rich outflows produced by the accretion-induced collapse of white dwarfs: light curves and spectra_ , _Mon. Not. R. Astron. Soc._ 409 (2010) 846 [1005.1081].
  * [16] L. Fuhrmann, J.A. Zensus, T.P. Krichbaum, E. Angelakis and A.C.S. Readhead, _Simultaneous Radio to (Sub-) mm-Monitoring of Variability and Spectral Shape Evolution of potential GLAST Blazars_ , in _The First GLAST Symposium_ , S. Ritz, P. Michelson and C.A. Meegan, eds., vol. 921 of _American Institute of Physics Conference Series_ , pp. 249–251, July, 2007, DOI [0704.3944].
  * [17] A. Taani, A. Khasawneh, M. Mardini, A. Abushattal and M. Al-Wardat, _Probability Distribution of Magnetic Field Strengths through the Cyclotron Lines in High-Mass X-ray Binaries_ , _arXiv e-prints_ (2020) arXiv:2002.03011 [2002.03011].
  * [18] Z. Dai, P. Szkody, A. Taani, P.M. Garnavich and M. Kennedy, _Quiescent photometric modulations of two low-inclination cataclysmic variables <ASTROBJ>KZ Geminorum</ASTROBJ> and <ASTROBJ>TW Virginis</ASTROBJ>_, _Astron. Astrophys._ 606 (2017) A45 [1708.04948].
  * [19] M.K. Mardini, V.M. Placco, A. Taani, H. Li and G. Zhao, _Metal-poor Stars Observed with the Automated Planet Finder Telescope. II. Chemodynamical Analysis of Six Low-metallicity Stars in the Halo System of the Milky Way_ , _Astrophys. J._ 882 (2019) 27 [1906.08439].
  * [20] M.K. Mardini, H. Li, V.M. Placco, S. Alexeeva, D. Carollo, A. Taani et al., _Metal-poor Stars Observed with the Automated Planet Finder Telescope. I. Discovery of Five Carbon-enhanced Metal-poor Stars from LAMOST_ , _Astrophys. J._ 875 (2019) 89 [1904.09603].
  * [21] M.K. Mardini, A. Frebel, A. Chiti, Y. Meiron, K.V. Brauer and X. Ou, _The Atari Disk, a Metal-poor Stellar Population in the Disk System of the Milky Way_ , _Astrophys. J._ 936 (2022) 78 [2206.08459].
  * [22] H. Ritter and U. Kolb, _VizieR Online Data Catalog: Cataclysmic Binaries, LMXBs, and related objects (Ritter+, 2004)_ , _VizieR Online Data Catalog_ (2011) B/cb.
  * [23] Y.-C. Wei, A. Taani, Y.-Y. Pan, J. Wang, Y. Cai, G.-C. Liu et al., _Neutron Star Motion in the Disk Galaxy_ , _Chinese Physics Letters_ 27 (2010) 119801.
  * [24] D.R. Lorimer, _Radio Pulsar Statistics_ , in _Astrophysics and Space Science Library_ , W. Becker, ed., vol. 357 of _Astrophysics and Space Science Library_ , p. 1, Jan., 2009, DOI.
  * [25] M.K. Mardini, V.M. Placco, Y. Meiron, M. Ishchenko, B. Avramov, M. Mazzarini et al., _Cosmological Insights into the Early Accretion of r-process-enhanced Stars. I. A Comprehensive Chemodynamical Analysis of LAMOST J1109+0754_ , _Astrophys. J._ 903 (2020) 88 [2009.12142].


# Razmecheno: Named Entity Recognition from Digital Archive of Diaries
“Prozhito”

###### ?abstractname?

The vast majority of existing datasets for Named Entity Recognition (NER) are
built primarily on news, research papers and Wikipedia with a few exceptions,
created from historical and literary texts. What is more, English is the main
source for data for further labelling. This paper aims to fill in multiple
gaps by creating a novel dataset “Razmecheno”, gathered from the diary texts
of the project “Prozhito” in Russian. Our dataset is of interest for multiple
research lines: literary studies of diary texts, transfer learning from other
domains, low-resource or cross-lingual named entity recognition.
Razmecheno comprises $1331$ sentences and $14119$ tokens, sampled from
diaries, written during the Perestroika. The annotation schema consists of
five commonly used entity tags: person, characteristics, location,
organisation, and facility. The labelling is carried out on the crowdsourcing
platfrom Yandex.Toloka in two stages. First, workers selected sentences, which
contain an entity of particular type. Second, they marked up entity spans. As
a result $1113$ entities were obtained. Empirical evaluation of Razmecheno is
carried out with off-the-shelf NER tools and by fine-tuning pre-trained
contextualized encoders. We release the annotated dataset for open access.

Keywords: named entity recognition, text annotation, datasets

languageresourceLanguage Resources

Razmecheno: Named Entity Recognition from Digital Archive of Diaries
“Prozhito”

Timofey Atnashev♡, Veronika Ganeeva♡, Roman Kazakov♡
---
Daria Matyash${}^{\heartsuit\text{\textdollar}}$, Michael Sonkin♡, Ekaterina
Voloshina♡
Oleg Serikov${}^{\heartsuit\Diamond\text{\textdaggerdbl}\sharp}$, Ekaterina
Artemova♡†♠
♡ HSE University ◇ DeepPavlov lab, MIPT ${}^{\text{\textdaggerdbl}}$ AIRI ♯
The Institute of Linguistics RAS
† Huawei Noah’s Ark Lab ♠ Lomonosov Moscow State University
${}^{\text{\textdollar}}$ Sber AI Centre
{taatnashev, vaganeeva, rmkazakov, dsmatyash, mvsonkin,
<EMAIL_ADDRESS>
{oserikov<EMAIL_ADDRESS>
Moscow, Russia

Abstract content

## 1\. Introduction

Modern Named Entity Recognition (NER) systems are typically evaluated on
datasets such as ACE, OntoNotes and CoNLL 2003, collected from news or
Wikipedia. Other common setups to test NER systems include cross-lingual
evaluation [Liang et al., 2020] and evaluation in domains, other than general,
such as biomedical domain [Weber et al., 2020, Wang et al., 2019].

Additionally, the vast majority of NER dataset are in English. A few large-
scale datasets for other languages are NoSta-D benikova2014nosta (German),
NorNE jorgensen2020norne (Norwegian), AQMAR mohit2012recall (Arabic),
OntoNotes hovy2006ontonotes (Arabic, Chinese), FactRuEval
starostin2016factrueval (Russian).

We present in this work a new annotated dataset for named entity recognition
from diaries, written in Russian, – “Razmecheno”111“Got annotated”. The short
form of the past participle neuter singular of the verb размечать (“to
annotate”). https://github.com/hse-cl-masterskaya-prozhito/main. The texts are
provided by the project “Prozhito”222“Got lived”. The short form of the past
participle neuter singular of the verb прожить (“to live”).
https://prozhito.org/, which digitizes and publishes personal diaries. Diaries
exhibit different surface and style features, such as a complex narrative
structure, an author-centricity, mostly expressed in simple sentences with
predominance of verbs and noun phrases.

Design choices, made for the corpus construction, are the following. We follow
the standard guidelines of named entity annotation and adopt four commonly-
used types Person (per), Location (loc), Organization (org), Facility (fac).
We add one more type, char, which is used for personal characteristic (e.g.
nationality, social group, occupation). Texts, used in the corpus, are sampled
from the diaries, written in the late 1980’s, the time period addressed as
Perestroika. We utilized crowd-sourcing to label texts.

Our dataset enables assessing performance of the NER models in a new domain or
in a cross-domain transferring. We make the following contributions:

  1. 1.

We present a new dataset for Named Entity Recognition of 14119 tokens from 124
diaries from Prozhito. Entity types, used in the dataset, follow standard
guidelines. The dataset will be freely available for download under a Creative
Commons ShareAlike 4.0 license at https://github.com/hse-cl-masterskaya-
prozhito/main;

  2. 2.

We assess the performance of the off-the-shelf NER taggers and fine-tuned
BERT-based model on this data.

## 2\. Related work

Most of the standard datasets for named entity recognition, as ACE
walker2005ace and CoNLL sang2003introduction, consist of general domain news
texts in English. For our study there are two related research lines: NER for
the Russian language and NER in Digital Humanities domain.

### 2.1. NER for Russian language

The largest dataset for Russian was introduced by loukachevitch-
etal-2021-nerel. In NEREL, entities of types per, org, loc, fac, gpe
(Geopolitocal entity), and family were annotated, and the total number of
entities accounts to 56K.

starostin2016factrueval presented FactRuEval for NER competition. The dataset
included news and analytical texts, and the annotation was made manually for
the following types: per, org and loc. As of now, it is one of the largest
datasets for NER in Russian as it includes 4907 sentences and 7630 entities.

Several other datasets for Russian NER, such as Named Entities 5, WikiNER, are
included into project Corus333https://github.com/natasha/corus. Its annotation
schema consists of 4 types: PER, LOC, Geolit (geopolitical entity), and MEDIA
(source of information). Another golden dataset for Russian was collected by
gareev2013introducing. The dataset of 250 sentences was annotated for per and
org. For the BSNLP-2019 shared task, a manually annotated dataset of 450
sentences was introduced piskorski2019second. The annotation includes per,
org, loc, pro (products), and evt (events).

Several silver datasets exist for Russian NER. WikiNEuRal tedeschi-
etal-2021-wikineural-combined uses multilingual knowledge base and transfomer-
based models to create an automatic annotation for per, loc, prg, and misc. It
includes 123,000 sentences and 2,39 million tokens. In Natasha project, a
silver annotation corpus for Russian Nerus444https://github.com/natasha/nerus
was introduced. The corpus contains news articles and is annotated with three
tags: per, loc, and org. For Corus project, an automatical corpus WikiNER was
created, based on Russian Wikipedia and methodology of WiNER ghaddar2017winer.

### 2.2. NER applications to Digital Humanities

bamman2019annotated introduced LitBank, a dataset built on literary texts. The
annotation was based on ACE types of named entities and it includes the
following types: per, org, fac, loc, gpe (geo-political entity) and veh
(Vehicle). The annotation was made by two of the authors for 100 texts. The
experiments with models trained on ACE and on LitBank showed that NER models
trained on the news-based datasets decrease significantly in the quality on
literary texts. brooke2016bootstrapped trained unsupervised system for named
entity recognition on literary texts, which bootstraps a model from term
clusters. For evaluation they annotated 1000 examples from the corpus.
Compared to NER systems, the model shows better results on the literary corpus
data.

Apart from English LitBank, a dataset for Chinese literary texts was created
and described by xu2017discourse. The dataset for Chinese literature texts had
both rule-based annotation and machine auxiliary tagging, hence, only examples
where gold labels and predicted labels differ were annotated manually. The
corpus of 726 articles were annotated by five people. Besides standard tags,
as per, loc, and org, the authors used tags thing, time, metric, and abstract.

Another approach to annotation was presented by wohlgenannt2016extracting. The
authors’ purpose was to extract social networks of book characters from
literary texts. To prepare an evaluation dataset, the authors used paid micro-
task crowd-sourcing. The crowd-sourcing showed high quality results and
appeared to be a suitable method for digital humanities tasks.

## 3\. Dataset collection

### 3.1. Annotation schema

Our tag set consists of five types of entities. This tag set was designed
empirically for texts of diaries from common tags used in related works
walker2005ace,bamman2019annotated.

  * •

per: names/surnames of people, famous people and characters are included (see
Example 3.1);

  * •

char: characteristics of people, such as titles, ranks, professions,
nationalities, belonging to the social group (see Example 3.1);

  * •

loc: locations/places, this tag includes geographical and geopolitical objects
such as countries, cities, states, districts, rivers, seas, mountains,
islands, roads etc. (see Example 3.1);

  * •

org: official organizations, companies, associations, etc. (see Example 3.1);

  * •

fac: facilities that were built by people, such as schools, museums, airports,
etc. (see Example 3.1);

  * •

misc: other miscellaneous named entities.

These five entity types can be clearly divided into two groups: the first one,
per-char, is related to people and the second one, org-loc-fac, is related to
places and institutions.

We annotated flat entities, so that the overlap between two entities is not
possible. The main principle of the annotation is to choose as long as
possible text span for each entity, not to divide them when not required,
because our schema does not assume multi-level annotation, when one entity can
include another ones. For example, a name and a surname coalesce in one PER
entity, rather than being two different ones (see Example 3.1).

А ведь Леон просил меня отозваться лишь о Жаке Ланге// And really
$\underbrace{\textbf{Leon}}_{\text{{per}}}$ asked me to.talk only about
$\underbrace{\textbf{Jack Lang}}_{\text{{per}}}$// ‘And Leon asked me to talk
only about Jack Lang’.//

Орёл самый литературный город в России //
$\underbrace{\textbf{Orel}}_{\text{{loc}}}$ the.most literary city in
$\underbrace{\textbf{Russia}}_{\text{LOC}}$// ‘Orel is the most literary city
in Russia’.//

Позвонил в “Урал”: надо все-таки дать им знать о моем прилете. // called in
$\underbrace{\textbf{``Ural''}}_{\text{{org}}}$ need after.all give them know
about my arrival // ‘I called the “Ural”: after all I have to let them know
about my arrival’.//

Солдаты живут в вагоне на этой станции. //
$\underbrace{\textbf{soldiers}}_{\text{{char}}}$ live in car on this
$\underbrace{\textbf{station}}_{\text{{fac}}}$// ‘Soldiers live in a car at
this station’.//

?figurename? 1: Annotation of a phrase given in Yandex.Toloka: Ira brought
socks as small presents for Sasha.
Available annotations are: Person (blue), Characteristics (green), Misc
(grey), no entities present (checkbox)

In ambiguous cases entity tags were identified based on the context, so the
same entity in different sentences could be tagged as two different types, for
instance, university could be annotated as org or fac . If an entity was used
in a metaphorical sense, it would not be annotated with any tag.

Будет и на нашей улице праздник // will and on our street a.festival// ‘Every
dog has its own day’.//

### 3.2. Preliminary markup

We performed preliminary analysis of the random subsets of the “Prozhito”
corpus. The analysis revealed that most of the sentences contain no entities
at all. To avoid costly looping over all sentences, we developed a two-stage
annotation pipeline. The first stage aims at selecting sentence candidates,
which may include entities of interest. This helps to reduce the amount of
sentences, sent to crowd workers and exclude sentences with no entities at
all. During the second stage entity spans are labeled in the pre-selected
candidates from the first stage.

Two classifiers were trained on a small manually annotated training set — for
per-char and org-loc-fac groups, respectively. The task of these classifiers
is to predict, whether an entity from a group is present in a sentence, or
not. These classifiers do not aim at entity recognition, but rather at binary
entity detection.

We leverage upon four possible base models as classifiers: ruBERT-
tiny555https://huggingface.co/cointegrated/rubert-tiny,
ruBERT666https://huggingface.co/DeepPavlov/rubert-base-cased,
ruRoBERTa777https://huggingface.co/sberbank-ai/ruRoberta-large, XLM-
RoBERTa888https://huggingface.co/xlm-roberta-base. Table 1 presents with the
classification scores. A small number of marked up sentences (198) were taken
as test sample.

Models | Precision | Recall | Micro f1-score
---|---|---|---
ruBERT-tiny | 0.81 | 0.88 | 0.84
ruBERT | 0.89 | 0.91 | 0.90
ruRoBERTa | 0.90 | 0.88 | 0.89
XLM-RoBERTa | 0.80 | 0.99 | 0.89

?tablename? 1: Transformer-based binary classifiers scores

As a result, ruRoBERTa was chosen as the base model. In this task, the
precision is more important than the recall, since we markup only part of the
corpus and, therefore, we still miss some of the information, but at the same
time we want to have any entities in the selected sentences with a high
probability. It is worth mentioning that this approach allowed to save big
amount of money, because the most optimal solution for maximizing the recall
would be to give all the data to assessors for markup.

To train both classifiers, a random sample of size 1500 was taken from diaries
belonging to Perestroika period. Texts were independently marked up by
assessors for the presence of org-loc-fac and per-char. Due to the fact that
it was important to achieve a balance of classes in the training sample, and
there were more texts with per-char then org-loc-fac , the training samples
for org-loc-fac and per-char turned out to be different - 829 and 1465 records
accordingly (see Table 2 for the validation set scores).

All available sentences were marked up by binary classifier and after that
were chosen sentences with following conditions:

  1. 1.

In sentence there are entities from per-char and org-loc-fac groups,
respectively;

  2. 2.

Classifier was the most confident on these sentences.

Entity Type | Precision | Recall | F1-score
---|---|---|---
org-loc-fac | 0.94 | 0.92 | 0.94
per-char | 0.89 | 0.81 | 0.82
?tablename? 2: ruRoBERTa scores in the binary classification task

Most confidence here means the average probabilities of the each entities
groups. Finally, the sentences selected this way were given to the assessors
for further marking.

### 3.3. Crowd-sourcing annotation

#### Annotation setup

For annotation we used Russian crowd-sourcing platform Yandex.Toloka
999https://toloka.yandex.ru/. We prepared two tasks for assessors:
determination of per-char and of org-loc-fac in “Prozhito” texts. The task
required native knowledge of the Russian language for satisfactory solution.
Before annotation it is necessary to get through the learning pool with hints
(20 sentences) and an exam (10 sentences) that show whether assessors
understand the meaning of the given NE tags. The sentences were tokenized with
Razdel tokenizer101010https://github.com/natasha/razdel.

The tasks for learning, exam and control were initially annotated by the co-
authors with help of annotation tool BRAT111111https://brat.nlplab.org/.

If an annotator succeeded in learning and exam (mark $\geq$ 50% for learning
and $\geq$ 80% for exam), he/she got a special skill which allowed to start
assessment of sentences in the main pool. Our main pools in both tasks consist
of approximately $1500$ tasks and $400$ control sentences. Tasks were given to
annotators on pages, Figure 1 depicts the task interface. Each page consisted
of $4$ normal tasks and $1$ control task. A fee for one page was 0.05$. The
average time of completion of page was about one minute. Overall, the fee per
hour exceeded minimum wage in Russia. The overlap for each sentence given in
Toloka is 3 in order to choose the most popular variant of markup as a correct
one. Control tasks are necessary for monitoring of an annotation quality. We
banned users if they skipped more than $7$ task suites in a row or if they had
less than $30\%$ correct control responses.

#### Annotators agreement analysis

While in most of the cases annotators had no dispute, voting mechanism has
been involved in nearly one third of cases provided in the corpus ($38\%$ in
the org-loc-fac task, $36\%$ in per-char tasks, respectively).

In both tasks typical annotators disagreement pattern was two competing
annotation hypotheses. In the org-loc-fac task, that was mostly caused by
different labels plausible for certain rare events. The ability to correctly
disambiguate such terms relied onto rather rare factual knowledge, thus
provoquing annotation errors (as in Сижу в гостинице “Одесса”. (‘Staying in
the hotel “Odessa”’.), the challenging choice is ‘hotel “Odessa”’ is a fac or
an org entity). While the same group of annonators disagreements was found in
the per-char task, there also emerged two more disagreements patterns: (i)
identifying the proper span for the characteristics (annotating the whole
полковник в отставке (‘the retired colonel’) or only полковник (‘colonel’) )
and (ii) inaccurate boundaries detection for persons initials, which mostly
emerged when the annontators missed to higlight the dot in the name
shortenings (as with М . С . in М . С . его очень ценил поначалу. (‘M.S.
valued him a lot in the beginning’) ) .

Type | # Entities | % Entities | # Mentions | % Mentions
---|---|---|---|---
char | 282 | 25.0% | 290 | 19.7%
fac | 71 | 6.4% | 106 | 7.2%
loc | 186 | 16.7% | 221 | 15.0%
org | 73 | 6.6% | 137 | 9.3%
per | 490 | 44.0% | 708 | 48.0%
misc | 11 | 1.0% | 12 | 0.8%
Total | 1113 | 100.0% | 1474 | 100.0%

?tablename? 3: Dataset entities statistics

Rare cases with more than two competing annotations were mostly of random
nature (as with birds being annotated as per), or caused by the appearance of
rare words (as with calzones being annotated as Person).

### 3.4. Dataset statistics

Entity Type |  Top-10 mentions
---|---
char |  ребёнок (‘a child’), женщина (‘a woman’), президент (‘a president’), друг (‘a friend’), поэт (‘a poet’), папа (‘a dad’), писатель (‘a writer’), жена (‘a wife’), отец (‘a father’), военный (‘a military’)
fac |  театр (‘a theatre’), аэропорт (‘an airport’), дом (‘a house’), школа (‘a school’), музей (‘a museum’), кафе (‘a cafe’), станция (‘a station’), библиотека (‘a library’), посольство (‘an embassy’), тюрьма (‘a prison’)
loc |  город (‘a city’), Москва (‘Moscow’), Россия (‘Russia’), улица (‘a street’), Ленинград (‘Leningrad’), проспект (‘an avenue’), Кандагар (‘Kandagar’), озеро (‘a lake’), страна (‘a country’), запад (‘west’)
org |  ЦК (‘Central Committee’), совет (‘a council’), парламент (‘a parliament’), Политбюро (‘Politburo’), Правда (‘Pravda’), КПСС (‘the Communist Party of the Soviet Union’), издательство (‘a publishing house’), верховный (‘supreme’), Мосфильм (‘Mosfilm’), союз (‘a union’)
per |  Горбачев (‘Gorbachev’), Борис (‘Boris’), Ельцин (‘Yeltsin’), Володя (‘Volodya’), Таня (‘ Tanya’), Витя (‘Vitya’), Рыжков (‘Ryzhkov’), Яковлев (‘Yakovlev’), Сергей (‘Sergey’), Иван (‘Ivan’)

?tablename? 4: Top-10 mentions for each entity type ?figurename? 2:
Distribution of tags

The total number of sentences in the dataset is $1331$ and the total number of
tokens is $14119$. The average sentence length is $10.61$ tokens (see Figure
2). $1113$ entities were identified at all ($1474$ mentions). The average
length of entity in tokens is $1.32$ token.

?figurename? 3: Distribution of entity types

Table 3 and Figure 3 describe dataset statistics. per is the most frequent
tag, a little less than a half of all entities are of this type. Persons are
often provided via a few tokens. The rest of types do not represent the same
variance between mentions and entities. MISC entities are only 1% of all
entities.

121212Notably, PER class contains lots of multiple words sentences

Popular mentions of entities actually represent concepts and personalities of
Perestroika period (see Table 4). As we can see, there are main politic
figures in the list (e.g. Boris Yeltsin, Mikhail Gorbachev, Nikolai Ryzhkov)
as well as old soviet political authorities (e.g. Central Committee, the
Communist Party of the Soviet Union, Politburo). Some words that were new at
that time, such as ‘a president’ (since Gorbachev became the first president
of USSR in 1990) are among the most frequent words. Another important trend is
the discussion of the Soviet-Afghan war, as Kardagan was one of the centres of
soviet troops’ dislocation.

## 4\. Evaluation

### 4.1. Off-the-shelf tools

We use a selection of of publicly available, NER systems: Natasha-Slovnet,
Stanza, SpaCy.

Slovnet is a neural network based tool for NLP tasks, including NER
annotation. Slovnet is a part of Natasha project.
131313https://github.com/natasha/slovnet#ner Slovnet’s annotation includes
per, loc and org.

Stanza is a Stanford state-of-art model
141414https://stanfordnlp.github.io/stanza/. Stanza is based on Bi-LSTM model
and CRF-decoder. Stanza for Russian is a 4-entity system, which includes per,
loc, org and misc.

NER system developed by SpaCy is a transition-based named entity recognition
component. We use natasha-SpaCy 151515https://github.com/natasha/natasha-spacy
model trained on two resources - Nerus 161616https://github.com/natasha/nerus
and Navec 171717https://github.com/natasha/navec. Natasha-SpaCy model can
detect per, loc and org entities in our dataset.

We have compared results of these models on our dataset.

Models | per | loc | org | Overall
---|---|---|---|---
SpaCy | 0.64 | 0.54 | 0.16 | 0.95
Stanza | 0.69 | 0.4 | 0.11 | 0.94
Natasha | 0.77 | 0.54 | 0.14 | 0.96
?tablename? 5: The performance of off-the-shelf tools

As seen from the table 5, Natasha-Slovnet showed the best perfomance on our
dataset for per and loc, while SpaCy was the best on loc and org detection.
However, the results of all models are significantly worse than the results on
other datasets (Appendix A). Such results prove our hypothesis that off-the-
shelf tools do not recognise entities on a diary-based dataset, for they were
trained on news data.

?figurename? 4: Confusion matrix for off-the-shelf tools per token in relative
weights

A closer look at the models’ performance (figure 4) reveals what caused entity
recognition issues. Most of the models often detect false loc and per
entities. In this case SpaCy shows the best results. Natasha-Slovnet has the
greatest recall especially on loc and per. All models often annotated org as a
non-entity. As our texts come from diaries written in 1990s, some
organisations could not exist anymore, and models do not recognise them.

fac and char were not on the entity lists of the models, therefore, the models
did not recognise these tags. However, we would expect the models to mark char
as per and fac as loc or org because those tags are related. Indeed this
happens for fac but not for char. This happens as most of named entities are
proper nouns and start with capitalized letters, unlike char. All models
annotated fac more often as org than as non entity.

Another problem is caused by false detection of named entities’ span
boundaries. To account for this, we introduced the following approach. We
counted all cases when models did not found entities at all, detected false
entities or used a wrong tag (combined as ‘false detected’) or models included
more or less words from one or both sides. Natasha showed the best results,
for it detects right boundaries for the most of the spans. The most common
error though for all models was not finding an entity. Other mistakes include
a shift of boundaries to the left and including more or less words on the left
side especially for per recognition. It could be possibly explained that
textscchar entity proceeds per entity (for instance, профессор Иванов
(‘professor Ivanov’) where ‘professor’ is char). Off-the-shelf models do not
include textscchar entity and could annotate them as per. Problems of narrower
boundaries could be caused by excluding quoting markers in automatic
annotation.

### 4.2. Fine-tuned models

We fine-tuned multiple Transformer models for NER: ruBERT, ruBERT-tiny,
ruRoBERTa, XLM-RoBERTa. The performance was evaluated according to F1-scores
per named entity and overall micro F1-score.

We used weighted cross-entropy as a loss function. Inverse tag frequency was
taken as weights for cross-entropy, which helped us gain better better results
on unbalanced data. We also sorted the dataset by the length of tokens and
than split it to batches, which slightly improved models’ performance. Models
were trained in an unfrozen manner. The detailed hyperparameters values used
to train the models are provided in the Appendix C. The performance was
evaluated according to per-class and overall micro-averaged F1-score.

### 4.3. Results

Natasha had the best F1-score among all off-the-shelf tools. Nevertheless,
results achieved for our corpus are below Natasha’s results on news-based
datasets.

Fine-tuned transformers showed better results than off-the-shelf tools.
Predictions made by XLM-RoBERTa had the highest overall F1-score, ruBERT
model’s performance had the best F1-scores for most tags (fac, loc, org) and
top-3 best results for char and per tags. According to Table 6, we can
consider ruBERT the best model for our datasets as it successfully predicts
major and minor classes.

The number of epochs was chosen according to the following criteria: the model
does not overfit on the train data and shows high results on the development
data. To this end, we used early-stopping. For ruBERT-tiny even 50 epochs were
not sufficient for reaching results comparable to other models’ performances.

Models | char | fac | loc | org | per | Overall
---|---|---|---|---|---|---
ruBERT-tiny | 0.712 | 0.8 | 0.748 | 0.4 | 0.738 | 0.731
ruBERT | 0.757 | 1.0 | 0.793 | 0.4 | 0.854 | 0.813
ruRoBERTa | 0.703 | 0.333 | 0.729 | 0.166 | 0.795 | 0.739
XLM-RoBERTa | 0.817 | 0.363 | 0.742 | 0.333 | 0.825 | 0.8

?tablename? 6: Transformer architectures F1-scores

According to Figure 4, char and per entities were mostly wrongly detected as o
by natasha, SpaCy and Stanza annotators. org tags were also erroneously
detected by these parsers, which was quite similar to the results of
transformer models’ results. loc tags almost in all cases were detected
correctly both by pre-trained parsers transformer models, while fac tags were
significantly better found by the former ones.

According to Figure 6, XLM-RoBERTa’s performance could be considered quite
successful: char tags, as well as per and loc, were almost infallibly
predicted. More exactly, per entity was never predicted as another entity on
test data. fac entity was mixed with ORG tag in XLM-RoBERTa’s predictions
while org tag itself is nearly in all cases is considered as o tag by the
model.

Figure 6 also presented ruBERT-tiny’s performance: char and org entities were
erroneously predicted as o more often, if compared to XLM-RoBERTa.
Nebertheless, in most cases the model predicts correctly. ruBERT-tiny
extracted all fac and almost all per tags without major errors.

As for ruBERT’s results, o tags were rarely misclassified as char, while all
other tags were predicted entirely correctly or with inconsequential mistakes.

ruRoBERTa’s performance was far from being perfect, as o-entities were heavily
confused with other tags, but most predictions of other entities were correct.

As for major tendencies in models’ predictions, we can notice that org entity
in most cases was detected as o tag which although was not desired, but still
can encourage us to re-analyze org entities and collect substantially more
examples of org tag occurrence. fac entities were either (in most cases)
correctly predicted, or mispredicted as org. o tags were sometimes detected as
per entity.

## 5\. Conclusion

This paper introduces Razmecheno, a novel dataset for Named Entity
Recognition. The texts in the dataset are sampled from the project “Prozhito”,
which comprises personal diaries, written in Russian, from the 17th century up
to the end of the 20th century. In particular, texts, marked up in Razmecheno
belong to the mid-1980 years, the period in Russia, commonly known as
Perestroika. Razmecheno is a middle-scale dataset so that it contains enough
data to carry out literal and historical studies.

The annotation schema, used in Razmecheno, is simplistic. It consists of five
named entity types, of which four are commonly used in NER datasets, namely,
persons, locations, organization, and facilities. An only named entity type,
introduced in this project, is characteristics of the different groups of
people. The annotations are flat; overlapped, or nested entities are not
allowed at the moment.

As our annotation schema matches a commonly used inventory of named entity
types, it is possible to leverage upon pre-trained models and transfer
learning techniques. The experimental evaluation of Razmecheno is two-fold.
First, we carry out an extensive analysis of how available off-the-shelf NER
tools cope with the task. The results reveal, that Natasha outperforms other
tools under consideration by a small margin. However, of five named entity
types, the off-the-shelf tools used to support only three. Next, we experiment
with four state-of-the-art pre-trained Transformers. A monolingual model,
ruBERT significantly outperforms other Transformers, followed by a
multilingual model XLM-RoBERTa.

There are a few directions for Razhmecheno development. We plan to annotate
the collected sentences for other information extraction tasks, including co-
reference resolution, relation extraction, and entity linking. This way,
Razhmecheno could serve as a test-bed for end-to-end information extraction
models. Experiments in domain adaptation and cross-lingual transfer from other
languages are another research line. Finally, we have set up the whole
environment to annotate texts from “Prozhito”, so that diaries from other
periods can be marked up with a little effort.

## 6\. Bibliographical References

##  ?refname?

  * Liang et al., 2020 Liang, Y., Duan, N., Gong, Y., Wu, N., Guo, F., Qi, W., Gong, M., Shou, L., Jiang, D., Cao, G., et al.  (2020).  Xglue: A new benchmark datasetfor cross-lingual pre-training, understanding and generation.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 6008–6018.
  * Wang et al., 2019 Wang, X., Zhang, Y., Ren, X., Zhang, Y., Zitnik, M., Shang, J., Langlotz, C., and Han, J.  (2019).  Cross-type biomedical named entity recognition with deep multi-task learning.  Bioinformatics, 35(10):1745–1752.
  * Weber et al., 2020 Weber, L., Sänger, M., Münchmeyer, J., Habibi, M., Leser, U., and Akbik, A.  (2020).  Hunflair: an easy-to-use tool for state-of-the-art biomedical named entity recognition.  arXiv preprint arXiv:2008.07347.

## 7\. Language Resource References

lrec2022-bib custom

##  ?appendixname? A Models performance on different datasets

Models | factru | ne5 | bsnlp | razmecheno
---|---|---|---|---
per | loc | org | per | loc | org | per | loc | org | per | loc | org
SpaCy | 0.901 | 0.886 | 0.765 | 0.967 | 0.928 | 0.918 | 0.919 | 0.823 | 0.693 | 0.64 | 0.54 | 0.16
Stanza | 0.943 | 0.865 | 0.687 | 0.923 | 0.753 | 0.734 | 0.938 | 0.838 | 0.724 | 0.69 | 0.4 | 0.11
Natasha | 0.959 | 0.915 | 0.825 | 0.984 | 0.973 | 0.951 | 0.944 | 0.834 | 0.718 | 0.77 | 0.54 | 0.14
ruBERT-tiny | 0.619 | 0.395 | 0.558 | 0.619 | 0.414 | 0.564 | 0.318 | 0.333 | 0.180 | 0.738 | 0.748 | 0.4
ruBERT | 0.548 | 0.358 | 0.461 | 0.883 | 0.777 | 0.856 | 0.483 | 0.451 | 0.423 | 0.854 | 0.793 | 0.4
ruRoBERTa | 0.468 | 0.261 | 0.406 | 0.768 | 0.593 | 0.687 | 0.192 | 0 | 0 | 0.795 | 0.729 | 0.166
XLM-RoBERTa | 0.879 | 0.763 | 0.78 | 0.963 | 0.936 | 0.944 | 0.762 | 0.899 | 0.726 | 0.825 | 0.742 | 0.333

?tablename? 7: See Section 2.1 for the review of these corpora in the Nerus
suite. The data on the performance for off-the-shelf were taken from Natasha
project 191919https://github.com/natasha/slovnet#ner

##  ?appendixname? B Off-the-Shelf models’ span recognition

To evaluate how precise off-the-shelf models are in span recognition, we
divide all cases of recognition in 11 groups:

  * •

left more: the right border of a span was detected correctly but on the left
border a model included more words than in our annotation;

  * •

right more: more words were included into a span on the right side;

  * •

left less: the right border was correctly detected but on the left side one or
more words were missing;

  * •

right less: the left border was detected but on the right side less words were
included;

  * •

more: on both sides a model annotated more words than in the data;

  * •

less: on the both sides a model detected a smaller span;

  * •

equal: a model detected a span correctly;

  * •

left right: the borders of a span were shifted from left to right, i.e. on the
left side less words were included and on the right side a model detected some
extra words;

  * •

right left: the borders of a span were shifted from right to left;

  * •

not found: models did not find a span or annotated it with a wrong tag;

  * •

false detected: models found spans that were not in the manual annotation.

Figure 5 shows the absolute number of cases of each type described above.

?figurename? 5: Off-the-shelf tools’ mistakes in span recognition for each
entity

##  ?appendixname? C Transformers hyper-parameters

Models | Number of epochs | Learning rate | Weight decay
---|---|---|---
ruBERT-tiny | 50 | 1e-5 | 3e-5
ruBERT | 10 | 1e-4 | 2e-5
ruRoBERTa | 5 | 1e-5 | 2e-5
XLM-RoBERTa | 10 | 3e-5 | 1e-4
?tablename? 8: Transformer architectures’ hyperparameters

##  ?appendixname? D Transformers confusion matrix

?figurename? 6: Confusion matrix of ruBERT, ruBERT-tiny, ruRoBERTa and XLM-
RoBERTa models’ results on the test dataset