# How Thermal Annealing Process Determines the Inherent Structure Evolution in
Amorphous Silicon: An Investigation from Atomistic Time Scales to Experimental
Time Scales

Yanguang Zhou<EMAIL_ADDRESS>Department of Mechanical and Aerospace
Engineering, The Hong Kong University of Science and Technology, Clear Water
Bay, Kowloon, Hong Kong

(August 28, 2024)

###### Abstract

Abstract

The annealing treatment in the advanced manufacturing process, e.g., laser
assisted manufacturing, determines the final state of glasses which is
critical to its thermal, electrical and mechanical properties. Energy barriers
analysis based on the potential energy surface offers an effective way to
study the microscopic evolution of the inherent structures during the
annealing process in a broadening timescale range, i.e., from atomistic
timescale ($\sim$ ps) to experimental timescale ($\sim$ s). Here, we find the
distribution activation energy barriers in the potential energy surface can be
divided into three regimes: 1, the distribution mainly follows the Rayleigh
distribution when the annealing rate $\dot{R}<5\times{{10}^{15}}\ \text{K/s}$;
2, two different modes, i.e., an exponentially decaying mode and a Rayleigh
distribution mode, are found in the spectra when the annealing rate
$5\times{{10}^{15}}\ \text{K/s}<\dot{R}<instant$; 3, the spectra is almost
following the exponentially decaying mode when the system is under the instant
annealing process. However, the spectra of relaxation energy barriers only
show an exponentially decaying mode with a decreasing decay parameter. A multi
timescale model for any specific annealing rate, which is beyond the limit of
the conventional atomistic simulations, i.e., molecular dynamics simulations,
is then proposed based on the distribution of the energy barriers. Such a
model enables quantitative explanations and predictions of the heat release
during the annealing process of the nanocalorimetry measurements or laser
assisted manufacturing.

Structural defects such as dislocations and grain boundaries in crystal cannot
be well characterized in amorphous material due to its topology of atomic
connectivity varies locally. For such a reason, amorphous material, e.g.,
amorphous silicon, has many promising physical, chemical and mechanical
properties comparing to their crystalline counterparts, e.g., extremely low
thermal conductivity Goldsmid1983 ; Larkin2014 ; Kwon2016 ; He2011 ; Liu2009 ;
Moon2018 ; Saaskilahti2016 ; Zhou2016 ; Zhou2015 ; Zhou20162 , exceptional
charge capacity Ding2015 ; Ding2017  and high strength Patinet2016 ; Fan2013 .
In amorphous material, many properties are strongly related to the degrees of
relaxation in the system which can be quantified via the inherent structural
energy (ISE), i.e., $E_{\text{IS}}$, calculating by removing the kinetic
energy Kallel2010 . Therefore, the properties such as damage evolution
Patinet2016 ; Fan2015  and flow stress Fan20132  of amorphous materials can be
modulated through ISE. One of most popular approach to vary the ISE of
amorphous materials is to adjust the dynamic annealing process Liu2020 ;
Fan2017 ; Ramachandramoorthy2016 . Recent development in ultrafast laser pulse
techniques Shou2019  make the annealing rate as high as 106 K/s which is still
an unresolved challenge to date to the traditional modelling tools, i.e.,
molecular dynamics (MD) techniques Fan2017 ; Deringer2018 ; Fan2014 ;
Chowdhury2016 . The potential energy surface (PES) generated based on the
inherent structure energy has been proved to be capable of characterizing the
complex phenomenology in amorphous materials Fan2015 ; Liu2020 ; Fan2017 ;
Fan2014 ; Yan2016 , and more importantly they allow focusing on the energy
distribution rather than the dynamic process, making it possible to capture
the experimental timescale annealing rates, i.e., $10^{-2}$ $\sim$ $10^{6}$
K/s Shou2019 ; Wan2011 ; Yuan2017 .

To generate the PES, the elementary processes are the hopping between
neighboring local minima. This process consists of two stages: 1, the
activation stage which is from the initial state to the connecting saddle
state with activation barrier ${{E}_{\text{A}}}$; and 2, the relaxation stage,
i.e., from the saddle state to the final state with energy barrier
${{E}_{\text{R}}}$. Then, the time evolution of ${{E}_{\text{IS}}}$ can be
expressed in the form of Liu2020 ; Fan2017 ; Derlet2013

$\displaystyle\frac{d{{{\bar{E}}}_{\text{IS}}}}{dt}=$ $\displaystyle\
\frac{dT}{dt}\frac{d{{{\bar{E}}}_{\text{IS}}}}{dT}\approx\dot{R}\left[\frac{{{{\bar{E}}}_{\text{IS}}}(T+\Delta
T)-{{{\bar{E}}}_{\text{IS}}}(T)}{\Delta T}\right]$ (1) $\displaystyle=$
$\displaystyle
f\cdot\exp\left(-\frac{{{{\bar{E}}}_{\text{A}}}}{{{k}_{b}}T}\right)\left({{{\bar{E}}}_{\text{A}}}-{{{\bar{E}}}_{\text{R}}}\right)$

in which, $f$ is the jump frequency which includes the entropy effects
Goldstein1976 ; Johari1977 ; Berthier2011 , ${{k}_{b}}$ is the Boltzmann
constant, is the temperature of the current inherent structure, $\dot{R}$ is
the annealing rate, and
${{\left[f\cdot\exp\left(-{{{\bar{E}}}_{\text{A}}}/{{k}_{b}}T\right)\right]}^{-1}}$
represents the average residence time in current inherent structure.
${{\bar{E}}_{\text{A}}}$ and ${{\bar{E}}_{\text{R}}}$ are the average
activation and relaxation energy barriers, which consider the effect of
multiple transition pathways in the PES of amorphous materials, and can be
calculated via Fan2017 ; Derlet2013

${{\bar{E}}_{\text{A}}}=-{{k}_{b}}T\cdot\ln\left[\int{P\left({{E}_{\text{A}}}|{{E}_{\text{IS}}}\right){{e}^{-\frac{{{E}_{\text{A}}}}{{{k}_{b}}T}}}d{{E}_{\text{A}}}}\right]$
(2) ${{\bar{E}}_{\text{R}}}=\int{{{E}_{\text{R}}}\cdot
P\left({{E}_{\text{R}}}|{{E}_{\text{IS}}}\right)d{{E}_{\text{R}}}}$ (3)

where $P\left({{E}_{A}}|{{E}_{ISE}}\right)$ and
$P\left({{E}_{R}}|{{E}_{ISE}}\right)$ are the activation and relaxation energy
spectra which are calculated via the activation-relaxation technique (ART)
Cancs2009 ; Barkema1996  – a method of atomistic simulation known to be able
to capturing the saddle point states and providing the standard PES samples
including these high activation energy barrier events that cannot be
accessible by MD simulations in amorphous materials. Here, the crystalline
silicon with 1728 atoms with interaction depicting using Tersoff potential
Tersoff1989  is first equilibrated at a high temperature (4000 K) liquid
state, and then cooled to 10 K with nine annealing rates ranging from 1010 K/s
up to 1016 K/s and the instant quench using MD simulations. The size of the
system is 3.26 nm$\times$3.26 nm$\times$3.26 nm, with periodic boundary
conditions applied to all the three directions. The system with 4096 atoms is
used to validate the size effect in our results can be ignored. Figure 1 shows
the evolution of the ISE at various temperatures in MD simulations at the
annealing rate of 1010 K/s up to 1016 K/s. In the high annealing rate process,
e.g., $\dot{R}={{10}^{16}}\ \text{K/s}$, the system would jump out of the
equilibrium state and eventually turn into the glassy state with an almost
constant ISE (black dots in Figure 1) since the timescale of cooling rate
cannot keep up with the intrinsic timescale in supercooled liquid. With the
decreasing of the annealing rate, e.g., $\dot{R}<{{10}^{16}}\ \text{K/s}$, the
systems at temperature above 1800 K can reach equilibrium state easily and
their corresponding ISEs decrease quickly, and then froze to glassy states
with different ISEs with temperature decreasing below 1800 K. Next, the ISs
generated in the MD runs are going to be used as input for ART simulations to
calculate the ${{E}_{\text{A}}}$ and ${{E}_{\text{R}}}$ spectra. For each of
the ten input structures, 5000 ART searchers with different random
perturbation directions are employed. The magnitude of perturbation
displacement is fixed at 0.5 Å, and the system is relaxed to the saddle point
following the Lanczos algorithm Barkema1996 ; Malek2000  when the curvature of
the PES is less than -0.01 eV$\cdot$Å-2 and the force of the system is less
than 0.05 eV$\cdot$Å-1. After removing the failed and redundant searches,
around 3,000 different searches are left for each IS.

Figure 1: The inherent energy of the systems during the annealing process with
different annealing rates, i.e., 1010 K/s up to 1016 K/s, in molecular
dynamics simulations. The solid dots are results from the models with 1728
atoms, and the open dots are calculated by the system of 4096 atoms. Our
results show the system with 1728 atoms is large enough to ignore the size
effects. Figure 2: The raw histogram data on the activation energy barrier
spectra in different samples generated at various annealing rates (a)
$10^{10}$ K/s, (b) $10^{14}$ K/s, (c) $10^{15}$ K/s, (d) 5$\times$$10^{15}$
K/s, (e) 1016 K/s and (f) instant quench obtained by ART. The spectra can be
fitted into two modes: an Rayleigh distribution (M1 mode, blue lines in the
Figure) with an obvious peak and an exponential decay mode (M2 mode, purple
lines in the Figure), using Eq. (4).

We then focus on the activation energy spectra
$P\left({{E}_{\text{A}}}|{{E}_{\text{IS}}}\right)$ and relaxation energy
spectra $P\left({{E}_{\text{R}}}|{{E}_{\text{IS}}}\right)$. Figure 2 shows the
activation energy spectra obtained for these samples. It is not surprising to
find faster annealing rate yield smaller activation energy range since the
intrinsic timescale is larger than the timescale of cooling rate and therefore
it is difficult for the system to catch up to the high saddle energy states
Fan2017 . More interesting, it is found that
$P\left({{E}_{\text{A}}}|{{E}_{\text{IS}}}\right)$ can be divided into three
regimes: 1, the spectra dominantly follows the Rayleigh distribution with a
clear peak when $\dot{R}<5\times{{10}^{15}}\ \text{K/s}$ (M1 mode in our
manuscript, Figure 2a-c); 2, the spectra can be decomposed into two different
modes: an exponentially decaying mode (M2 mode) and a Rayleigh distribution
mode (M1 mode) when $5\times{{10}^{15}}\ \text{K/s}<\dot{R}<instant$ Figure
2d-e); 3, the spectra is mainly following the M2 mode when the system is under
the instant annealing process (Figure 2f). Based on the probability
distribution of the spectra, the spectra can be represented phenomenologically
using Liu2020 ; Fan2017 :

$\displaystyle P\left({{E}_{\text{A}}}|{{E}_{\text{IS}}}\right)=$
$\displaystyle{{P}_{\text{{M}1}}}\left({{E}_{\text{A}}}|{{E}_{\text{IS}}}\right)+{{P}_{\text{{M}2}}}\left({{E}_{\text{{A}}}}|{{E}_{\text{{IS}}}}\right)$
(4) $\displaystyle=$
$\displaystyle{{W}_{1}}\cdot{{E}_{\text{{A}}}}\cdot\exp\left[-\frac{{{({{E}_{\text{{A}}}}-\mu({{E}_{\text{{IS}}}}))}^{2}}}{2{{\sigma}^{2}}}\right]+$
$\displaystyle{{W}_{2}}\cdot\frac{1}{{{\xi}_{\text{{A}}}}}\cdot\exp\left(-\frac{{{E}_{\text{{A}}}}}{{{\xi}_{\text{{A}}}}}\right)$

where ${{W}_{1}}$ and ${{W}_{2}}$ are the weigh factors related to the
amplitude of M1 and M2 modes, respectively. $\mu$ is the location parameter of
M2, $\sigma$ is the width of the M1 mode, and ${{\xi}_{\text{A}}}$ is the
decay parameter of M2 mode. It is worth noting the integral of
$P\left({{E}_{\text{A}}}|{{E}_{\text{IS}}}\right)$ should be equal to one due
to the normalization condition of the spectra. As seen in Figure 2, when the
inherent structure energy becomes higher, i.e., the higher annealing rate, the
amplitude of M1 becomes smaller and disappear eventually (Figure 2), while the
amplitude of M2 shifts to a larger value and becomes the main mode when the
system is under instant annealing (Figure 2f). In the similar way, we find the
spectra of relaxation energy
$P\left({{E}_{\text{R}}}|{{E}_{\text{IS}}}\right)$ can be fitted in the
phenomenological formula of Kallel2010 ; Liu2020 ; Fan2017 ; Derlet2013

$P\left({{E}_{\text{R}}}|{{E}_{\text{IS}}}\right)=\frac{1}{{{\xi}_{\text{R}}}}\cdot\exp\left(-\frac{{{E}_{\text{R}}}}{{{\xi}_{\text{R}}}}\right)$
(5)

in which, ${{\xi}_{\text{R}}}$ is the decay parameter. Unlike the previous
studies found the relaxation spectra is independent of the annealing rate
Kallel2010 , our results show that higher annealing rate will lead to a broad
range of the relaxation energy barrier (Figure 3). As discussed in Ref.
Fan2015 ; Fan2014 , the higher annealing rate will lead to a faster
exponential decay behavior of relaxation energy barrier distribution
indicating a localized relaxation process which is found in various systems,
and with decreasing of the annealing rate, the probability of the relaxation
energy barrier distribution will become broader due to the combination of
localization and cascade (i.e., self-organized criticality which has been
observed universally in many different systems experimentally Krisponeit2014
and numerically Salerno2012 ) relaxations of ${{E}_{\text{R}}}$ in the PES.

Figure 3: The raw data on the relaxation energy barrier spectra in different
samples generated at various annealing rates (a) 1010 K/s, (b) 1014 K/s, (c)
5$\times$1010 K/s and (d) instant quench obtained by ART. Lines are fitted
using Eq. (5).

To solve Eq. (1) at any given ISE, we need prepare the forms of ${{W}_{1}}$,
$\mu$, $\sigma$, ${{W}_{2}}$, ${{\xi}_{\text{A}}}$ and ${{\xi}_{\text{R}}}$.
In our study, we have nine samples which can provide nine data of each
parameter mentioned above. What is worth noting is that a wider range of ISE,
particularly those ISEs at slower annealing rate, e.g., $\dot{R}<{{10}^{10}}\
\text{K/s}$, can make sure the fitting process of the parameters more
accurately. But it will take prohibiting computing time or even out of the
power of MD simulations. In this study, we firstly fit these parameters using
appropriate functions and then compare their generated average saddle-initial
or final minimal energy differences with MD results directly. Based on the
data generated from the ART results, we can find ${{W}_{1}}=\
-92.4+97.6\cdot\exp({{E}_{\text{IS}}}/81.42)$ (Figure 4a), $\sigma\ =\
0.233+(2.234-0.233)/\left[1+\exp(({{E}_{\text{IS}}}+4.367)/0.019)\right]$
(Figure 4b),
$\mu=2.224-17.25\times{{10}^{6}}\cdot\exp({{E}_{\text{IS}}}/0.25)$ (Figure
4c),
${{W}_{2}}=10\cdot\exp(-48.61-28.30{{E}_{\text{IS}}}-4.19E_{\text{IS}}^{2})$
(Figure 4d),
${{\xi}_{\text{A}}}=1.308+1.370\times{{10}^{-28}}\cdot\exp(-{{E}_{\text{IS}}}/0.066)+1.242\times{{10}^{-28}}\cdot\exp(-{{E}_{\text{IS}}}/0.065)$
(Figure 4e),
${{\xi}_{R}}=3.836-1.075{{E}_{\text{IS}}}-0.985E_{\text{IS}}^{2}-0.125E_{\text{IS}}^{3}$
(Figure 4f). The average saddle-initial or final energy difference is then
calculated for validating our fitting formulas via

$\displaystyle{{\bar{E}}_{{saddle}-initial\ or\ final}}=$ (6)
$\displaystyle\int{{{E}_{\text{A}\ or\ \text{R}}}\cdot P\left({{E}_{\text{A}\
or\ \text{R}}}|{{E}_{\text{IS}}}\right)d{{E}_{\text{A}\ or\ \text{R}}}}$

Figure 4: The magnitude of (a) ${{W}_{1}}$, (b) $\sigma$, (c) $\mu$, (d)
${{W}_{2}}$, (e) ${{\xi}_{\text{A}}}$ in Eq. (4) and (f) ${{\xi}_{\text{R}}}$
in Eq. (5) v.s. the inherent energy $E_{\text{IS}}$. To make Eq. (1) to be
used at any specific annealing rate, we fit the discrete data using empirical
equations: (a) ${{W}_{1}}=-92.4+97.6\cdot\exp({{E}_{\text{IS}}}/81.42)$, (b)
$\sigma=0.233+(2.234-0.233)/\left[1+\exp(({{E}_{\text{IS}}}+4.367)/0.019)\right]$,
(c) $\mu=2.224-17.25\times{{10}^{6}}\cdot\exp({{E}_{\text{IS}}}/0.25)$, (d)
${{W}_{2}}=10\cdot\exp(-48.61-28.30{{E}_{\text{IS}}}-4.19E_{\text{IS}}^{2})$,
(e)
${{\xi}_{\text{A}}}=1.308+1.370\times{{10}^{-28}}\cdot\exp(-{{E}_{\text{IS}}}/0.066)+1.242\times{{10}^{-28}}\cdot\exp(-{{E}_{\text{IS}}}/0.065)$
and (f)
${{\xi}_{R}}=3.836-1.075{{E}_{\text{IS}}}-0.985E_{\text{IS}}^{2}-0.125E_{\text{IS}}^{3}$.

The lines in Figure 5a are calculated using Eq. (6) with the functions of the
parameters fitted above (Figure 4), and the dots in Figure 5a are the results
obtained from discrete ART data. The good agreement between the lines and the
dots in Figure 5a indicates rationality of the fitting forms of ${{W}_{1}}$,
$\mu$, $\sigma$, ${{W}_{2}}$, ${{\xi}_{\text{A}}}$ and ${{\xi}_{\text{R}}}$.
To fully obtain the dynamic process of the ISE during the annealing or heating
process at arbitrary annealing rate, or equivalently solving Eq. (1), the
hypothesis on the jump frequency is still necessary. As firstly proposed by
Goldstein Goldstein1976  and Johari Johari1977 , the jump frequency $f$
depends on the entropy of the system which is related to the temperature and
ISE. Since the ISE is a function of temperature as well, we then assume $f$
only depends on the temperature and can be expressed as
${{f}_{0}}\exp[-1.25\times{{10}^{-4}}(T-300)]/{{(9.6T)}^{2}}$ Fan2017  where
${{f}_{0}}$ is the zero temperature jump frequency and is set as
$1\times{{10}^{9}}$ THz Fan2013 ; Yan2016 ; Fan2012 . Next, the Eq. (1) can be
solved at arbitrary annealing rate $\dot{R}$, and good agreements between our
theoretical predictions with MD simulations are found (Figure 5b).

The model developed above is then used to study the heat release of the system
$Q$ at experimental annealing rates, i.e., ${{10}^{-2}}\sim{{10}^{6}}\
\text{K/s}$ Shou2019 ; Wan2011 ; Yuan2017 , which can be written as:

$\displaystyle Q=$
$\displaystyle\frac{d\bar{E}}{dT}=\frac{d{{{\bar{E}}}_{\text{IS}}}+d{{{\bar{E}}}_{kinetic}}}{dT}$
(7) $\displaystyle\approx$
$\displaystyle\left[\frac{{{{\bar{E}}}_{\text{IS}}}(T+\Delta
T)-{{{\bar{E}}}_{\text{IS}}}(T)}{\Delta T}+1.5{{k}_{b}}\right]$

Combining Eqs. (1) and (7), we can know

$\displaystyle
Q=\frac{f}{{\dot{R}}}\cdot\exp\left(-\frac{{{{\bar{E}}}_{\text{A}}}}{{{k}_{b}}T}\right)\left({{{\bar{E}}}_{\text{A}}}-{{{\bar{E}}}_{\text{R}}}\right)+1.5{{k}_{b}}$
(8)

Figure 5: (a) The average relaxation energy calculated using ART data and
empirical formulas fitted in Figure 4. (b) The inherent energy changes with
temperature at various annealing rates, dots are from MD simulations and lines
are calculated using Eq. (1).

Figure 6a shows the heat release of the system during the annealing process of
amorphous Si at the rates of ${{10}^{-2}}\ \sim{{10}^{6}}\ \text{K/s}$ using
Eq. (8). At high temperature range ($\geq$ 200 K), the systems can be
equilibrated at very short time scale and therefore the abundant heat will be
released. However, at the low temperature range ($\leq$ 200 K), the annealing
rate is too fast to equilibrate the system, and the system is undergoing a
normal ageing behavior Liu2020 ; Fan2017 . We now try to compare our
predictions with the experimental measurements Karmouch2007 . In the
experiments, the Si ions are implanted into the amorphous samples, and then
the released heat of the system is measured. By comparing the results in
Figure 3 and 4 in Ref. Karmouch2007 , we can assume the number of atoms of the
ion-implanted region is about 1016. By choosing the ISE of the systems and
annealing rate carefully, our theoretical predictions agree very well with the
experimental measurement (Figure 6b). It is worth noting that the higher
density implanted ions will cause higher ISE of the system as we assumed here.
Meanwhile, the annealing rate can be as high as ${{10}^{6}}\ \text{ K/s}$ in
the laser flash experiments Shou2019  and as low as ${{10}^{-2}}\ \text{K/s}$
in the furnace lamp annealing experiments Shrestha2018 , here we use the
medium annealing rate 10 K/s in our model.

Figure 6: (a) the predicted heat release during the annealing process with
different annealing rates using Eq. (8), which we assume the
${{E}_{\text{IS}}}=-4.51\ \text{eV}$. (b) The het release when the annealing
rate is 10 K/s with appropriate inherent energy, experimental measurements are
from Ref. Karmouch2007 , excellent agreement between theoretical predictions
and experimental data is found.

We also note that in Ref. Kallel2010 , the potential energy surface of
activation and relaxation energy barriers in amorphous Si is almost
independent on the inherent structures. The reason for this may be due to the
original structures used in their study are the same, and the inherent energy
is varying via deleting atoms in the system. In our paper, we generate the
inherent structures via mimicing the real experimental annealing process, we
find the distribution of activation energy can be decomposed in two modes
which is found in other systems such as amorphous Cu56Zr44 Fan2017  and the
distribution of relaxation energy is varying with the annealing rate which is
also observed in other systems, e.g., Stillinger-Weber silicon Middleton2001 .
We would also like to emphasize that the accuracy of our theoretical
predictions depends on the accuracy of
$P\left({{E}_{\text{A}}}|{{E}_{\text{IS}}}\right)$ and
$P\left({{E}_{\text{R}}}|{{E}_{\text{IS}}}\right)$ over the entire range of
inherent structure energy. Molecular dynamics simulations can assure the
parameters of $P\left({{E}_{\text{A}}}|{{E}_{\text{IS}}}\right)$ and
$P\left({{E}_{\text{R}}}|{{E}_{\text{IS}}}\right)$ in high annealing rate from
${{10}^{10}}\ \text{ K/s}$ to ${{10}^{16}}\ \text{ K/s}$ and the instant
annealing rate, i.e., inherent surface energy from -4.44 eV/atom to -3.97
eV/atom, in our theoretical model are reasonable, while there is no guarantee
that their extrapolations are still suitable at experimental annealing rates,
i.e., ${{10}^{-1}}\sim{{10}^{4}}\ \text{ K/s}$. However, the underlying
mechanism behind the theoretical model proposed here is still valid since Eq.
(8) describes the physics process of hopping in potential energy surface.

In conclusion, the activation and relaxation energy barriers’distributions are
extracted based on the extensive sampling of the potential energy surface of
a-Si. We find the distribution of the activation energy barriers can be
divided into three regimes: 1, the spectra dominantly follows the Rayleigh
distribution when $\dot{R}<5\times{{10}^{15}}\ {K/s}$; 2, the spectra can be
decomposed into two different modes: an exponentially decaying mode and a
Rayleigh distribution mode when $5\times{{10}^{15}}\ {K/s}<\dot{R}<instant$;
3, the spectra is mainly following the exponentially decaying mode when the
system is under the instant annealing process. The distribution of relaxation
energy barrier is found to follow the exponentially decaying mode with a
decreasing decaying parameter at the whole annealing rate range. Then, a long-
time scale model beyond the limit of the conventional atomistic simulations
used to capture the annealing process observed in experiments of generating
a-Si is proposed based on distributions of the energy barriers. The heat
release during the annealing process is calculated using our proposed long
timescale model, and compared with the nanocalorimetry measurements. Excellent
agreement between our numerical results and experimental measurements, makes
sure the present model not only can unravel the microscopic information of
evolution of the inherent structure but also applicable to many thermodynamics
processes such as laser assisted manufacturing process.

Y.Z. thanks startup fund from Hong Kong University of Science and Technology
(HKUST). Y. Z. gratefully acknowledges Mr. Zhitong Bai (University of
Michigan-Ann Arbor) for valuable discussions.

## References

  * (1) H. J. Goldsmid, M. M. Kaila, and G. L. Paul, Phys. Status Solidi 76, K31 (1983).
  * (2) J. M. Larkin and A. J. H. McGaughey, Phys. Rev. B 89, 1 (2014).
  * (3) S. Kwon, J. Zheng, M. C. Wingert, S. Cui, and R. Chen, ACS Nano 11, 2470 (2016).
  * (4) Y. He, D. Donadio, and G. Galli, Appl. Phys. Lett. 98, 3 (2011).
  * (5) X. Liu, J. L. Feldman, D. G. Cahill, R. S. Crandall, N. Bernstein, D. M. Photiadis, M. J. Mehl, and D. A. Papaconstantopoulos, Phys. Rev. Lett. 102, 1 (2009).
  * (6) J. Moon, B. Latour, and A. J. Minnich, Phys. Rev. B 97, 1 (2018).
  * (7) K. Sääskilahti, J. Oksanen, J. Tulkki, S. Volz, J. and A. J. H. Mcgaughey, AIP Adv. 6, 121904 (2016).
  * (8) Y. Zhou and M. Hu, Nano Lett. 16, 6178 (2016).
  * (9) Y. Zhou and M. Hu, Phys. Rev. B 92, 195205 (2015).
  * (10) Y. Zhou, X. Zhang, and M. Hu, Nanoscale 8, 1994 (2016).
  * (11) B. Ding, X. Li, X. Zhang, H. Wu, Z. Xu, and H. Gao, Nano Energy 18, 89 (2015).
  * (12) B. Ding, H. Wu, Z. Xu, X. Li, and H. Gao, Nano Energy 38, 486 (2017).
  * (13) S. Patinet, D. Vandembroucq, and M. L. Falk, Phys. Rev. Lett. 117, 1 (2016).
  * (14) F. Fan, S. Huang, H. Yang, M. Raju, D. Datta, V. B. Shenoy, A. C. T. Van Duin, S. Zhang, and T. Zhu, Model. Simul. Mater. Sci. Eng. 21, 074002 (2013).
  * (15) H. Kallel, N. Mousseau, and F. Schiettekatte, Phys. Rev. Lett. 105, 045503 (2010).
  * (16) Y. Fan, T. Iwashita, and T. Egami, Phys. Rev. Lett. 115, 1 (2015).
  * (17) Y. Fan, Y. N. Osetskiy, S. Yip, and B. Yildiz, Proc. Natl. Acad. Sci. 110, 17756 (2013).
  * (18) C. Liu, X. Yan, P. Sharma, and Y. Fan, Comput. Mater. Sci. 172, 109347 (2020).
  * (19) Y. Fan, T. Iwashita, and T. Egami, Nat. Commun. 8, 1 (2017).
  * (20) R. Ramachandramoorthy, W. Gao, R. Bernal, and H. Espinosa, Nano Lett. 16, 255 (2016).
  * (21) W. Shou, B. Ludwig, L. Wang, X. Gong, X. Yu, C. P. Grigoropoulos, and H. Pan, ACS Appl. Mater. Interfaces11, 34416 (2019).
  * (22) V. L. Deringer, N. Bernstein, A. P. Bartók, M. J. Cliffe, R. N. Kerber, L. E. Marbella, C. P. Grey, S. R. Elliott, and G. Csányi, J. Phys. Chem. Lett. 9, 2879 (2018).
  * (23) Y. Fan, T. Iwashita, and T. Egami, Nat. Commun. 5, 1 (2014).
  * (24) S. C. Chowdhury, B. Z. Haque, and J. W. Gillespie, J. Mater. Sci. 51, 10139 (2016).
  * (25) X. Yan and P. Sharma, Nano Lett. 16, 3487 (2016).
  * (26) Z. Wan, S. Huang, M. A. Green, and G. Conibeer, Nanoscale Res. Lett. 12, 478 (2011).
  * (27) Z. Yuan, C. Wang, K. Chen, Z. Ni, and Y. Chen, Nanoscale Res. Lett. 12, 1 (2017).
  * (28) P. M. Derlet and R. Maa$\beta$, Philos. Mag. 93, 4232 (2013).
  * (29) M. Goldstein, J. Chem. Phys. 64, 4767 (1976).
  * (30) G. P. Johari, Philos. Mag. 35, 1077 (1977).
  * (31) L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).
  * (32) E. Cancs, F. Legoll, M. C. Marinica, K. Minoukadeh, and F. Willaime, J. Chem. Phys. 130, (2009).
  * (33) G. T. Barkema and N. Mousseau, Phys. Rev. Lett. 77, 4358 (1996).
  * (34) J. Tersoff, Phys. Rev. B 39, 5566 (1989).
  * (35) R. Malek and N. Mousseau, Phys. Rev. E 62, 7723 (2000).
  * (36) J. O. Krisponeit, S. Pitikaris, K. E. Avila, S. Küchemann, A. Krüger, and K. Samwer, Nat. Commun. 5, 3616 (2014).
  * (37) K. M. Salerno, C. E. Maloney, and M. O. Robbins, Phys. Rev. Lett. 109, 1 (2012).
  * (38) Y. Fan, Y. N. Osetsky, S. Yip, and B. Yildiz, Phys. Rev. Lett. 109, 1 (2012).
  * (39) R. Karmouch, Y. Anahory, J. F. Mercure, D. Bouilly, M. Chicoine, G. Bentoumi, R. Leonelli, Y. Q. Wang, and F. Schiettekatte, Phys. Rev. B 75, 075304 (2007).
  * (40) M. Shrestha, K. Wang, B. Zheng, L. Mokrzycki, and Q. H. Fan, Coatings 8, 97 (2018).
  * (41) T. F. Middleton and D. J. Wales, Phys. Rev. B 64, 024205 (2001).


# Two New Black Widow Millisecond Pulsars In M28

Andrew Douglas University of Virginia, 1826 University Avenue,
Charlottesville, VA 22904, USA  Prajwal V. Padmanabh Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Scott M. Ransom
NRAO, 520 Edgemont Road., Charlottesville, VA 22903, USA  Alessandro Ridolfi
INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047
Selargius (CA), Italy  Max-Planck-Institut für Radioastronomie, Auf dem Hügel
69, D-53121 Bonn, Germany  Paulo Freire Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Vivek Venkatraman
Krishnan Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121
Bonn, Germany  Ewan D. Barr Max-Planck-Institut für Radioastronomie, Auf dem
Hügel 69, D-53121 Bonn, Germany  Cristina Pallanca Dipartimento di Fisica e
Astronomia “Augusto Righi”, Università degli Studi di Bologna, Via Gobetti
93/2, 40129 Bologna, Italy  INAF-Osservatorio di Astrofisica e Scienze dello
Spazio di Bologna, Via Gobetti 93/3 I-40129 Bologna, Italy  Mario Cadelano
Dipartimento di Fisica e Astronomia “Augusto Righi”, Università degli Studi di
Bologna, Via Gobetti 93/2, 40129 Bologna, Italy  INAF-Osservatorio di
Astrofisica e Scienze dello Spazio di Bologna, Via Gobetti 93/3 I-40129
Bologna, Italy  Andrea Possenti INAF – Osservatorio Astronomico di Cagliari,
Via della Scienza 5, I-09047 Selargius (CA), Italy  Ingrid Stairs Dept. of
Physics and Astronomy, University of British Columbia, 6224 Agricultural Road,
Vancouver, BC V6T 1Z1, Canada  Jason W. T. Hessels ASTRON, the Netherlands
Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The
Netherlands  Anton Pannekoek Institute for Astronomy, University of Amsterdam,
Postbus 94249, 1090 GE Amsterdam, The Netherlands  Megan E. DeCesar George
Mason University, Fairfax, VA 22030, resident at the U.S. Naval Research
Laboratory, Washington, D.C. 20375, USA  Ryan S. Lynch Green Bank Observatory,
PO Box 2, Green Bank, West Virginia, 24944, USA  Matthew Bailes Centre for
Astrophysics and Supercomputing/OzGrav, Swinburne University of Technology,
Mail 74, PO Box 218, Hawthorn, Vic, 3122, Australia  Marta Burgay INAF –
Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius
(CA), Italy  David J. Champion Max-Planck-Institut für Radioastronomie, Auf
dem Hügel 69, 53121 Bonn, Germany  Ramesh Karuppusamy Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Michael Kramer Max-
Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany
Benjamin Stappers Jodrell Bank Centre for Astrophysics, Department of Physics
and Astronomy, The University of Manchester, Manchester M13 9PL, UK  Laila
Vleeschower Jodrell Bank Centre for Astrophysics, Department of Physics and
Astronomy, The University of Manchester, Manchester M13 9PL, UK

(Received November 10, 2021; Revised December 28, 2021; Accepted December 28,
2021)

###### Abstract

We report the discovery of two Black Widow millisecond pulsars in the globular
cluster M28 with the MeerKAT telescope. PSR J1824$-$2452M (M28M) is a 4.78-ms
pulsar in a 5.82 hour orbit and PSR J1824$-$2452N (M28N) is a 3.35-ms pulsar
in a 4.76 hour orbit. Both pulsars have dispersion measures near 119.30 pc
cm-3 and have low mass companion stars ($\sim$0.01-0.03 $M_{\odot}$), which do
not cause strong radio eclipses or orbital variations. Including these
systems, there are now five known black widow pulsars in M28. The pulsar
searches were conducted as a part of an initial phase of MeerKAT’s globular
cluster census (within the TRAPUM Large Survey Project). These faint
discoveries demonstrate the advantages of MeerKAT’s survey sensitivity over
previous searches and we expect to find additional pulsars in continued
searches of this cluster.

M28 – Pulsar – PRESTO – Jerk Search – MeerKAT – MeerTime – TRAPUM

††journal: Astrophysical Journal††facilities: MeerKAT, GBT††software: PRESTO
(Ransom, 2011b), PSRCHIVE (van Straten et al., 2011), Astropy (Astropy
Collaboration et al., 2013)

##  1 Introduction

Globular clusters (GCs) are rich environments for the creation of millisecond
pulsars (MSPs; e.g. Ransom, 2008). MSPs are created when neutron stars accrete
mass and angular momentum from companion stars and are spun-up to millisecond
periods. Due to the high stellar density in the central regions of globular
clusters between 1,000 and 100,000 $M_{\odot}$ pc-3 111Structural Parameters
of Globular Clusters:
https://people.smp.uq.edu.au/HolgerBaumgardt/globular/parameter.html stellar
interactions produce dynamically-formed binaries that can spin-up neutron
stars (e.g. Verbunt & Freire, 2014). Additionally, these dynamical encounters
can produce MSPs in highly eccentric and exotic orbits. Due to the high
likelihood of finding MSPs within GCs, and especially those in exotic orbits,
GCs have been the focus of many pulsar searches using the largest radio
telescopes and over 232 have been found so far222List of pulsars in globular
clusters: http://www.naic.edu/~pfreire/GCpsr.html. There are 39 pulsars
currently known in Terzan 5 (Lyne et al., 1990, 2000; Ransom et al., 2005;
Hessels et al., 2006; Andersen & Ransom, 2018; Cadelano et al., 2018; Ridolfi
et al., 2021) and 27 in 47 Tucanae (Camilo et al., 2000; Manchester et al.,
1991; Robinson et al., 1995; Pan et al., 2016; Ridolfi et al., 2021). M28 had,
up until this work, 12 known pulsars (Bégin, 2006), ranking it as the third
cluster in number of known pulsars.333Note that NGC 1851 is now tied with M28:
http://www.naic.edu/~pfreire/GCpsr.html

M28 (NGC 6626) is a fairly massive globular cluster in the bulge of the Milky
Way at a distance of $\sim$5.4 kpc (Baumgardt & Vasiliev, 2021). It contains
the first known globular cluster pulsar (PSR B1821$-$24 or M28A; Lyne et al.,
1987), as well as a transitional MSP that behaves alternatively as a radio MSP
and as an active X-ray binary (IGR J18245$–$2452 or M28I; Papitto et al.,
2013). M28 has been the subject of many multi-wavelength searches in the past
few decades. ROSAT observed M28A as a possible x-ray source (Danner et al.,
1994) and the Chandra X-ray Observatory has observed the cluster multiple
times and in 2011 confidently located Pulsars A, C, D, F, H, J and K in
X-rays. (e.g. Bogdanov et al., 2011). The Green Bank Telescope (GBT) also
observed M28 multiple times and uncovered 11 radio MSPs in the early to mid
2000s (Bégin, 2006). One notable point about these MSPs is that they seem to
contain a disproportionately large percentage of Black Widow pulsars, i.e.
systems where the companion star is of very low mass ($\leq$ 0.04
$M_{\odot}$), often causing eclipses of the radio pulsations, and having
orbits of 2$-$10 hrs in duration (Roberts, 2013). M28 also contains two
Redback pulsars, i.e. pulsar binaries with a more massive companion ($\sim$
0.1$-$0.2 $M_{\odot}$) than the Black Widow systems, but also showing eclipses
or timing irregularities, sometime across the entire orbit: Black Widow and
Redback pulsars are collectively dubbed as spiders (Roberts, 2013).

MeerKAT is a newly commissioned 64-dish radio telescope array in South Africa
(Booth & Jonas, 2012). Two pulsar-based large science projects, MeerTime and
TRAPUM, working together as described below, have made observations of several
globular clusters. Initial results from a globular cluster survey targeting
clusters with known pulsars, for the purposes of comparing older detections
with MeerKAT ones, has been presented by Ridolfi et al. (2021). With 15 years
since the last dedicated searches (using the GBT), the high sensitivity of
MeerKAT coupled with better computing power and the availability of higher-
dimension acceleration or “jerk” searches (e.g. Andersen & Ransom, 2018)
motivated new dedicated searches of M28.

In this paper we present the discovery and timing of M28M and M28N. In §2 we
present the radio observations and detection of the MSPs and the subsequent
timing procedures. In §3 we discuss searches for the systems in archival HST
and Chandra data. In §4 we discuss some implications of these new systems and
in §5 we present our conclusions.

##  2 Observations and Data Analysis

###  2.1 Data

The initial MeerKAT data that we used to search for new pulsars in M28 was
part of the series of survey observations using the core MeerKAT antennas
($\sim$42 out of 64 dishes) that targeted globular clusters containing known
pulsars (Ridolfi et al., 2021). These observations were made as part of a
joint effort between the MeerTime project (Bailes et al., 2020), which
conducts long-term pulsar timing, and the TRAPUM project (TRansients And
PUlsars with MeerKAT; Stappers & Kramer, 2016)444TRAPUM: http://trapum.org,
whose goal is to discover new radio pulsars in targeted regions that have
yielded significant discoveries in the past as well as environments that are
conducive for pulsar formation.

Two 9000-second duration MeerTime observations occurred on 2019 July 19 and
2020 February 4, with 642 MHz of nominal bandwidth centered at 1283.58 MHz and
split into 768 frequency channels. Each channel was coherently dedispersed at
a dispersion measure (DM) of 119.892 pc cm-3, which is roughly the average
value of the DMs of the known pulsars in the cluster, and were subsequently
further frequency subbanded, using incoherent dedispersion, into 256 total-
intensity channels with a sampling time of $\sim$76.56 $\mu$s. In addition to
the two MeerTime observations, we had access to the central beam of two
9000-second TRAPUM observations taken on 2020 December 6 and 11, which were
not coherently dedispersed but had 2048 frequency channels over 856 MHz of
bandwidth centered at 1283.86 MHz. These observations used additional antennas
than the earlier MeerTime observations, and were therefore approximately
10$-$20% more sensitive per unit time than the earlier observations. As with
the MeerTime data, these observations were later subbanded, after intra-
channel dispersion removal, into 256 frequency channels and integrated in time
to $\sim$76.56 $\mu$s samples for analysis. Analysis of the non-central TRAPUM
beams will be published elsewhere.

###  2.2 Searches

We used PRESTO (Ransom, 2011a) to perform Fourier domain acceleration and jerk
searches on the M28 data using its routine accelsearch (Ransom, 2001). The
jerk search functionality is a relatively recent addition to PRESTO (Andersen
& Ransom, 2018). Jerk searching aims to find pulsars whose accelerations are
changing linearly during an observation. The jerk term improves search
sensitivities to binary pulsars when the integration times (or duration of
searched segments) are between 5$-$15% of the orbital periods of the pulsars
(Bagchi et al., 2013; Andersen & Ransom, 2018).

We searched a range of dispersion measures (DMs) from 115.0 pc cm-3 to 125.0
pc cm-3, with a step size of 0.1 pc cm-3. We chose this range as it is roughly
centered on the average DM of the pulsars in the cluster, and it would be
unlikely that we would find any MSPs significantly outside this range. We
conducted differing levels of maximal acceleration and jerk searching (i.e. up
to -zmax 200 -wmax 600), where $z$ represents the number of Fourier bins a
signal can drifted during an observation and $w$ is the derivative of that. We
also conducted acceleration-only searches on shorter segments of the MeerTime
data with 29 and 58 million samples, corresponding to around one third (50
minutes) and two thirds (100 minutes) of the 117 million sample data file,
respectively. This segmented method of searching allows one to detect more
compact orbits that would be missed with full-duration acceleration searches.
With the TRAPUM data, we conducted segmented searches of 45 and 90 million
samples, very roughly one quarter (60 minutes) and one half (120 minutes) of
the 187 million sample data file, respectively.

We first searched the full-duration observations and segments with
acceleration only (-zmax 200), hoping to find a pulsar through the increased
sensitivity of MeerKAT, but without success. We then switched to jerk
searching, with -zmax 200 and -wmax 600. In the full search we detected the
new M28M at a DM of 119.30 pc cm-3 in the 2019 July 19 MeerTime observation at
$z$=19.50 and $w$=10. We were unable to confirm the pulsar in the 2020
February 4 MeerTime observation (until much later, via folding with an orbital
ephemeris) and thus searched the more sensitive TRAPUM observations from 2020
December 6th and 11th, and confirmed M28M in those data. We continued
searching TRAPUM and MeerTime data at -zmax 200 and -wmax 600, and soon
detected a second candidate that we confirmed to be M28N at $z$=3 and $w$=60
in the 2020 February 4 data. Since both pulsars were detected in the central
TRAPUM beam (0.246′ radius for the December 6 observation and 0.210′ radius
for the December 11 observation), we knew the MSPs must be very close to the
cluster center.

###  2.3 Timing

With multiple detections of both new pulsars over the four MeerKAT
observations, we were able to determine circular orbits and average
barycentric spin periods for each pulsar using the fit_circular_orbit.py
routine from PRESTO. These initial ephemerides allowed us to iteratively fold,
detect, and determine pulse arrival times from the large number of archival
observations of M28 from the GBT, starting in 2005 and continuing to the
present day. We note that the spin period of M28M might actually be 9.56 ms,
twice the period that we report. The pulsar was initially discovered at that
period, but folding of the discovery, subsequent, and archival observations at
the longer period show two gaussian-like peaks with similar amplitudes and
widths, and separated by $\sim$180∘. We therefore timed the pulsar assuming it
was a 4.78 ms rotator. This ambiguity will be resolved by future MeerTime
observations taken with higher time resolution and polarization information.

The GBT data from 2005 to 2008 was a mixture of L-band (centered at 1.5 GHz)
and S-band (centered at 2 GHz) observations initially using the Pulsar Spigot
(Kaplan et al., 2005) and having effectively $\sim$600 MHz of available
bandwidth with 1536 frequency channels sampled every 81.92 $\mu$s (see Bégin,
2006, for more details). Starting in 2010 and continuing through mid-2020,
observations were made with the Green Bank Ultimate Pulsar Processing
Instrument (GUPPI; DuPlain et al., 2008), using the same bands, but with 512
coherently dedispersed channels, subsequently partially integrated to 128
channels with 40.96 $\mu$s sampling. Finally, for the past year, we used VEGAS
(Bussa & VEGAS Development Team, 2012), in an identical mode to the GUPPI
observations.

Both the new pulsars were detected, albeit often with very low signal-to-
noise, in the vast majority of the archival GBT observations. The total GBT
integration time, sampled roughly quarterly, is 49.6 hours at L-band and 73.8
hours at S-band.

We determined times of arrival (TOAs) for each pulsar using the get_TOAs.py
routine in PRESTO after fitting gaussians with PRESTO’s pygaussfit.py to the
highest signal-to-noise TRAPUM observation (see Fig. 1) to use as noiseless
pulse templates. We were typically able to get roughly one TOA per hour of
telescope time with the GBT and one per half hour with MeerKAT. These TOAs
enabled us to phase connect each pulsar back to 2005 using TEMPO
555http://tempo.sourceforge.net/, although the slightly fainter M28M required
us to use the DRACULA666https://github.com/pfreire163/Dracula software (Freire
& Ridolfi, 2018) to establish phase connection. The reduced $\chi^{2}$ values
of the fits for both pulsars are near one (1.20 for M28M and 1.01 for M28N),
although for M28M, we needed to increase the nominal TOA errors (i.e. using
the EFAC parameter) by 50% for Spigot data and 30% for GUPPI/VEGAS/MeerKAT
data to achieve that value. The final timing solutions, as well as the plots
of the associated timing residuals, are provided in Table 2 and Fig. 2.

Figure 1: The best detections of M28M and M28N, respectively (left and right),
from the MeerKAT TRAPUM observation on 2020 December 11. The pulses, each with
64 bins in the profiles, show two full rotations of the pulsar and the
greyscale shows the intensity of the emission as a function of time. The error
bar under the summed profiles indicates the approximate time resolution of the
data, including dispersive effects, $\sim$350 $\mu$s.  Figure 2: Timing
residuals for M28M and M28N. The top two panels show timing residuals as a
function of time. The bottom two panels show the timing residuals as a
function of orbital phase. Eclipses of the pulsed radio emission would
typically occur at superior conjunction, which happens at orbital phase 0.25.

##  3 Optical and X-ray observations

###  3.1 Data set and data reduction

To search for the optical counterparts of the two newly discovered binaries
M28M and M28N, we use the same archival observations used to discover the
optical counterparts to M28H and M28I (Pallanca et al., 2010, 2013). The data
set is composed of deep and high-resolution optical images acquired with the
Wide Field Camera 3 (WFC3) and the Advanced Camera for Surveys (ACS) mounted
on the Hubble Space Telescope (HST). We refer to the aforementioned papers for
more details about the adopted data set. The data reduction is performed using
DAOPHOT (Stetson, 1987) following the prescriptions by Cadelano et al. (2019,
2020a). Instrumental magnitudes are calibrated to the VEGAMAG photometric
system using the stars in common with the catalogs presented in Pallanca et
al. (2010, 2013). The instrumental positions, after correction for geometric
distortions, are transformed in the absolute coordinate systems using the
stars in common with the Gaia Early Data Release 3 (Gaia Collaboration et al.,
2020). The resulting astrometric uncertainty is about $0.01\arcsec$. We verify
the presence of a possible offset between the radio and optical frames using
M28H as reference, whose position is very well measured in both frames. The
position of M28H in the Gaia optical frame is RA=$18^{h}24^{m}31.61112^{s}$
and Dec=$-24^{\circ}52\arcmin 17.39856\arcsec$ and thus differs of about
$-0.01\arcsec$ and $0.18\arcsec$ in RA and Dec, respectively, from the radio
position quoted in Pallanca et al. (2010). We therefore apply this offset to
the optical frame. We cannot include M28I in the offset calculation due to its
large uncertainty in the radio timing position.

###  3.2 Searching for optical counterparts

We investigate all of the sources in the area surrounding the positions of the
two binaries. In both cases, we find a stellar object in a position
compatible, within the positional uncertainty, with that derived from radio
data. The finding charts are presented in Figure 3. The candidate optical
counterparts to M28M and M28N are displaced from the radio pulsar positions by
only $\sim 0.05\arcsec$ and $\sim 0.07\arcsec$, respectively.

In the color-magnitude diagrams (CMDs), the candidate companion to M28M is
located along the cluster main sequence, $\sim 2.9$ mags below the turn-off
(see Figure 4). Interestingly, it shows a hint of red excess in the
F390W(U-band)$-$F814W(I-band) filter combination, while a small blue excess is
detected in the F606W(V-band)$-$F814W combination. The observations in the
$H_{\alpha}$ narrow filter (i.e. F658N) show no sign of $H_{\alpha}$ excess.

On the other hand, the candidate companion to M28N is located at $\sim 5.4$
mag below the turn-off and shows a blue excess in all the available filter
combinations. Such an excess is commonly observed in exotic binaries like
cataclysmic variables and MSPs (e.g. Rivera Sandoval et al., 2018; Zhao et
al., 2019, 2020), and it is usually related to the presence of an accretion
disk and/or heating of the secondary star by the primary. Also in this case,
no evidence of $H_{\alpha}$ excess is detected.

Figure 3: Left-hand panel: $2\arcsec\times 2\arcsec$ finding chart of the
region around the radio position of M28M. The image is a combination of the
F814W images. The red cross is centered on the MSP position, and the red
circle has a radius equal to the combined radio and optical positional
uncertainties. Right-hand panel: same as in the left-hand panel, but for the
case of M28N.

Figure 4: Left-hand panel: CMD of M28 obtained in the F390W$-$F814W filter
combination. Magnitudes have been corrected for the effects of differential
reddening following the prescriptions by Cadelano et al. (2020b). The
positions of the candidate counterparts to M28M and M28N are highlighted with
a red point and a blue square, respectively. Right-hand panel: same as in the
left-hand panel, but in the F606W$-$F814W filter combination.

Despite the good positional agreement and the peculiar CMD positions, both the
objects show no sign of photometric variability in any of the available
filters. However, companions to black widows identified both in globular
clusters and in the Galactic field are often characterized by strong magnitude
variations over their short orbital periods. The phenomenon is commonly due to
the heating of the companion side exposed to the pulsar wind (e.g. Pallanca et
al., 2014; Cadelano et al., 2015; Kaplan et al., 2018; Zharikov et al., 2019).
The absence of any magnitude modulation in the proposed candidate companions
might suggest that either the two objects are not associated with the pulsars
or the heating phenomenon is significantly less efficient than observed in
typical black widow systems (possibly also due to geometrical effects).

Given all this, we cannot safely conclude that the two proposed counterparts
are associated with the MSPs. In fact, the offset positions with respect to
the cluster main-sequence might also indicate that these two objects are field
interlopers rather than cluster members. Such a possibility could be confirmed
in the future through their proper motion analysis. Moreover, the only two
black widow optical counterparts identified in a globular cluster so far (i.e.
PSR J1518$+$0204C and PSR J1953$+$1846A; Pallanca et al. 2014; Cadelano et al.
2015) are located 6 to 9 magnitudes below their cluster turn-off (depending on
the orbital phase) in the F606W filter. Therefore, both the objects are
fainter than the two optical counterparts here presented. This strengthens the
likelihood that these two objects are physically unassociated with the pulsars
but are unrelated stars belonging either to the M28 population or to the
foreground.

###  3.3 Searching for x-ray counterparts

While we did not perform an independent analysis of the multiple Chandra
observations of M28, we did query the online Chandra Source Catalog (CSC) v2.0
to see if any point sources matched the positions of the two pulsars. Since
the pulsars are only separated by 3.7″, the two closest sources to each pulsar
turn out to be the same two CSC sources, 2CXO J182433.0$-$245210, which was
source 29 in Becker et al. (2003), and 2CXO J182433.2$-$245209, although
neither x-ray source is within 1.7″(more than twice the radius of the formal
x-ray positional errors) of the timing position of either pulsar. The limiting
sensitivity of the combined observations at the pulsar positions is
$\sim$1.2$\times$10-15 erg s-1cm-2 for M28M and $\sim$1.0$\times$10-15 erg
s-1cm-2 for M28N, in the ACIS 0.5$-$7 keV band, corresponding to luminosity
limits of $\lesssim$4.0$\times$1030 erg s-1 and $\lesssim$3.4$\times$1030 erg
s-1 for M28M and M28N, respectively, at the 5.37 kpc distance of M28
(Baumgardt & Vasiliev, 2021). These limits are similar to the measured ACIS-
band luminosities of the many detected millisecond pulsars, including several
black widows, in the globular cluster 47 Tucanae (e.g. Bogdanov et al., 2006;
Bhattacharya et al., 2017), and are therefore not exceptionally constraining.

##  4 Discussion

There are many interesting implications from the timing of both M28M and M28N.
In this section, we compare the exotic binary population in M28 to that in
other globular clusters, we examine the implications of orbital period
derivative measurements (or limits), and we compare the timing-based proper
motion measurements with the Gaia-based proper motion of M28 as a whole.

###  4.1 The “Spider” population

A peculiarity of M28 is its significantly larger fraction of spiders (see
Section 1) than 47 Tucanae (hereafter 47Tuc) and Terzan 5 (hereafter Ter5),
the two globular clusters with the richest content in pulsars to date. In
particular, introducing the ratio $\mathcal{R}$ between the population of
spiders and the total known pulsar population in a globular cluster, we have
$\mathcal{R}_{\rm M28}=50\%$ (14 known pulsars, among which are five black
widows and two redbacks (Bégin, 2006)), $\mathcal{R}_{\rm 47Tuc}=33\%$ (27
known pulsars with six black widows and three redbacks), and $\mathcal{R}_{\rm
Ter5}=13\%$ (two black widows and three redbacks out of a total of 39 known
pulsars). Summing over the data for all the other 33 globular clusters
containing known pulsars, the ratio is $\mathcal{R}_{\rm 33gc}=14\%$, close to
the value seen in Ter5777All the numbers are taken from
https://www.naic.edu/~pfreire/GCpsr.html. However, we focus only on the three
aforementioned most populated clusters, as we believe that the sample size for
them is large enough that the fraction of spiders can be seen as a trend, at
least to first approximation.

Using our knowledge of the numbers of spiders in M28, 47Tuc and Ter5, we can
make a guess as to how many of these types of systems might be potentially
observable in the three clusters. Bagchi et al. (2011) utilized Monte Carlo
simulations to predict the number of observable pulsars using six different
models. These models predicted M28 to have between 42 and 952 observable
pulsars, with an average between the models of 320 $\pm$ 60\. From the same
analysis, Ter5 should have a model averaged pulsar content of 380 $\pm$ 40
pulsars and 47Tuc of 160 $\pm$ 20 pulsars. It must be noted that there is a
large discrepancy in the absolute number of predicted pulsars across the range
of these models, but the relative ratios of predicted pulsars in the 3
globular clusters are roughly compatible between models, with Ter5 predicted
to have a total pulsar population 15%-40% larger than M28 and 120%-150% larger
than 47Tuc. Combining these percentage intervals with the values of
$\mathcal{R}$, we find that M28 should be by far the leading globular cluster
in terms of absolute number of spiders, containing about three times the
number of spiders than either of the other two globular clusters, the spider
content of which should be similar to each other. In absolute terms, according
to the entire range of models explored by Bagchi et al. (2011), M28 might host
from $\sim 20$ observable spiders up to $\sim 400$, while assuming a
conservative log-normal distribution in the pulsar luminosity, the simulations
of Bagchi et al. (2011) indicate a number of potentially observable spiders in
M28 ranging from $\sim 20$ to $\sim 80.$

These results naturally raise the question about the reasons for the
overabundance of spiders in M28 with respect to other globular clusters, and
in particular 47Tuc and Ter5. It is well known that the formation of spiders
in a globular cluster can occur via two channels: ordinary binary evolution
(as in the Galactic field), or through the occurrence of 3- or 4-body
dynamical interactions, which can support the creation of new binaries with
suitable parameters to later become spiders. The former pathway certainly
depends on the total mass $M_{\rm GC}$ of the globular cluster, while the
latter process is in turn favored by a crowded and massive stellar
environment, suggesting an exploration of both the role of the rate
$\Gamma_{c}$ of stellar interactions in the core (which depends on the central
density, $\rho_{c}$ and the core radius of the globular cluster, $r_{c}$,
$\Gamma_{c}=\rho_{c}^{3/2}r_{c}^{2}$, see e.g. Verbunt & Freire 2014), and of
the mass density within the half-mass radius $\rho_{hm}$. Using the data
provided in the website by Holger and Baumgardt888Structural Parameters of
Globular Clusters:
https://people.smp.uq.edu.au/HolgerBaumgardt/globular/parameter.html, as far
as the values of the total mass, Ter5, 47Tuc, and M28 (sorted in this way)
appear in the top 23% of the listed globular clusters. They are also in the
top 22% in terms of interaction rate in the core (with Ter5 leading over M28
and 47Tuc). As for $\rho_{hm},$ M28 (sixth position) and Ter5 (ninth position)
fell in the top 10% of the globular clusters, with 47Tuc in the top 27%. The
high ranking in all the mentioned lists is in agreement with the observed
occurrence of a high number of millisecond pulsars in the three clusters.
However, M28 prevails among the three clusters only in terms of mass density
within the half-mass radius, with a limited edge on Ter5, i.e.
$\log\rho_{hm,\rm M28}=3.49,$ and $\log\rho_{hm,\rm Ter5}=3.34.$ On the other
hand, despite a similar predicted number of spiders (see above), Ter5 and
47Tuc show significantly different values of the latter parameter
($\log\rho_{hm,\rm 47Tuc}=2.65$), as well as of $\Gamma_{c}$ (three times
larger in Ter5 than in 47Tuc), having only similar values for the total mass.

In summary, besides the suggested impact of $\rho_{hm}$, which could be
explored with detailed simulations of its effects on the production of the
spiders, some other still unassessed and/or as yet unmeasured parameters are
at play in determining the population of spiders in a cluster, e.g. one can
think about the fraction of binaries, or the past dynamical history of the
cluster. An interesting indication could arise from the distinct rarity of
eclipses in the sample of the M28 spiders. In fact, of the five black widow
systems, only one shows obvious eclipses, M28G (Bégin, 2006). The timing
analyses of these systems also do not display strong orbital variability,
which is very commonly detected in similar systems. In principle, these
peculiarities could be due to a geometrical effect, i.e. we are looking onto
the systems nearly face on, or at an inclination of near $0^{\circ}$. However,
this is statistically unlikely to occur for many binary systems at once. An
alternate option is that the relative absence of eclipses and the orbital
stability are indications that the companion star is stable and there is
little stellar wind being ablated. This often implies that the systems have
had time to stabilize and hence are older than the typical population of
spiders in globular clusters. If this hypothesis is supported by additional
observations and further discoveries in M28, it would suggest that the high
number of spiders in this globular cluster is sustained by an accumulation
mechanism, where the spiders were produced far in the past yet have managed to
survive intact in the cluster environment.

Table 1: Comparative Projected Pulsar Populations in M28, 47Tuc, and Ter5

Average Potential Observable MSPs | 320 $\pm$ 60 | 160 $\pm$ 20 | 380 $\pm$ 40
---|---|---|---
Ratio of ”Spider” Systems To Known MSPs ($\mathcal{R}$) | 0.50 | 0.33 | 0.13
Mass Density Within Half-mass Radius ($\log\rho_{hm}$) | 3.49 | 2.65 | 3.34

Note. — These are the relevant values taken from section 4.1

###  4.2 Corrections for Pulse Period Derivative

Pulsars in Globular clusters are often found to have seemingly nonphysical
values of $\dot{P}$ compared to pulsars in the Galaxy. This contamination is a
result of primarily the acceleration of the pulsar in the gravitational field
of the cluster and leaves no vestige of the intrinsic period derivative in the
observed period derivative (Phinney, 1992). In order to solve for the
intrinsic period derivative $\dot{P}_{\rm int}$, we can (e.g. Freire et al.,
2017) use the equation

$\frac{\dot{P}_{\rm obs}}{P}=\frac{\dot{P}_{\rm
int}}{P}+\frac{\mu^{2}d}{c}+\frac{a_{l,{\rm GC}}}{c}+\frac{{a}_{l,{\rm
gal}}}{c},$ (1)

where $\dot{P}_{\rm obs}$ is the observed spin-period derivative. The second
term on the right represents the Shklovskii effect (Shklovskii, 1970) due to
the motion of the pulsar in the plane of the sky; $\mu$ is the proper motion
of the cluster, $d$ is the distance to the cluster, and $c$ is the speed of
light. The third term, $a_{l,{\rm GC}}$ represents the line-of-sight
acceleration of the pulsar in the globular cluster. The fourth term represents
the difference in acceleration, $a_{l,{\rm gal}}$, between the solar system
and the cluster in the field of the galaxy.

The Shklovskii effect may be calculated as $3.12\times 10^{-10}$ m s-2 using
the total proper motion of the cluster (8.9 mas/year) and its 5.37 kpc
distance (Vasiliev & Baumgardt, 2021; Baumgardt & Vasiliev, 2021)). The
galactic acceleration may be computed as $3.08\times 10^{-10}$ m s-2 using the
galactic coordinates of M28 and its distance (e.g. Freire et al., 2017). The
line-of-sight acceleration from the cluster may then be solved for using

$a_{l,{\rm GC}}=\frac{\dot{P}_{b,{\rm obs}}}{P_{b}}c-\mu^{2}d-a_{l,{\rm gal}}$
(2)

where $P_{b}$ and $\dot{P}_{b,{\rm obs}}$ are the orbital period and its
observed period derivative. For M28N, which has a more significant measurement
of $\dot{P}_{b,{\rm obs}}$, the cluster acceleration is $1.02\times 10^{-8}$ m
s-2. With that value, we can then solve for $\dot{P}_{\rm int}$ with Eqn. 1,
giving $\sim 3.8\times 10^{-20}$ for M28N, and allowing estimates of the spin-
down based physical parameters in Table 2.

With a less precise measurement of $\dot{P}_{b,{\rm obs}}$ for M28M, we can
assume a 95% lower limit of the true $\dot{P}_{b,{\rm obs}}$ of $4.4\times
10^{-13}$ given the estimated error on $\dot{P}_{b,{\rm obs}}$ and then
proceed in a similar manner for M28N to get an upper limit of $\dot{P_{\rm
int}}<2.2\times 10^{-20}$ with corresponding upper limits (magnetic field and
$\dot{E}$) or lower limits (characteristic age) on the resulting physical
spin-down parameters in Table 2 (denoted by †).

###  4.3 Proper Motions

The long timing baseline provided by the GBT data allows low-precision
estimates of the individual pulsar proper motions as well as the orbital
period derivatives mentioned in the previous section. The Gaia mission has
measured the proper motion of M28 as $\mu_{\alpha}=-0.28(3)$ mas yr-1 and
$\mu_{\delta}=-8.92(3)$ mas yr-1, for a total proper motion of 8.92(3) mas
yr-1 (Vasiliev & Baumgardt, 2021)). The timing-based proper motions (or limits
on measurements) are consistent with the Gaia values within one sigma,
supporting the association of both binary systems with M28 and their being
bound to the cluster.

##  5 Summary and ramifications

We used sensitive new TRAPUM and MeerTime data to detect two new black widow
pulsars in the globular cluster M28, and exploited archival GBT observations
of the cluster to determine long-term coherent timing solutions for both
binary pulsars. Their orbital parameters reflect the typical values seen in
black widow systems in other globular clusters. However, their discovery
supports the growing evidence that the structural parameters and the internal
dynamics of M28 make it a very effective factory of exotic pulsar binaries,
here defined as the joint class of black widow and redback systems. Half of
the known pulsars in M28 belong to these categories, and projecting this high
fraction on preceding pulsar population studies, M28 could harbor several tens
or even some hundreds of these binaries. Some of them could be unveiled during
the ongoing search of the outer beams of the TRAPUM observations. We also
searched for the companions (or overall systems) of M28M and M28N in the
optical and X-ray bands, but were unable to conclusively identify either
system in either band.

Besides confirming the sensitivity and the potential of the MeerKAT telescope
to find new pulsars in globular clusters (with these discoveries the total
count for MeerKAT reaches 36 new
pulsars999http://trapum.org/discoveries.html), the case of M28M and M28N also
highlights the value of the long-term monitoring and archiving of globular
cluster pulsar data. As in the case of Terzan5 (Ridolfi et al., 2021), the use
of archival GBT data allowed us to establish for M28M and M28N a $\sim$
16-year-long timing solution soon after the discovery of these pulsars. This
is especially important for the data from large single-dish telescopes, like
the GBT, which have wider beams and very likely include the positions of newly
discovered pulsars.

In this context, it is worth noting that M28M and M28N were not discovered in
the GBT data because they are, on most occasions, too faint to be detected by
blind search algorithms. However, after their orbits had been characterized in
the MeerKAT observations, the number of trials necessary for a detection
decreased dramatically. That implied that the pulsed signals from these
objects could be recovered in the archival GBT data with a signal-to-noise
high enough for determining a suitable set of pulse times of arrival and then
deriving a phase-connected timing solution with highly precise timing
parameters. For these reasons, it is clear that the current archival data of
many globular clusters will continue to be of great use for fast
characterization of many MeerKAT discoveries. The same also applies, of
course, for the archival data from HST and from high-energy telescopes, which
are also important to immediately investigate the new pulsars at those
wavelengths.

The MeerKAT telescope is operated by the South African Radio Astronomy
Observatory, which is a facility of the National Research Foundation, an
agency of the Department of Science and Innovation. SARAO acknowledges the
ongoing advice and calibration of GPS systems by the National Metrology
Institute of South Africa (NMISA) and the time space reference systems
department of the Paris Observatory. MeerTime data is housed on the OzSTAR
supercomputer at Swinburne University of Technology. The OzSTAR program
receives funding in part from the Astronomy National Collaborative Research
Infrastructure Strategy (NCRIS) allocation provided by the Australian
Government. The Green Bank Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Associated Universities,
Inc. SMR is a CIFAR Fellow and is supported by the NSF Physics Frontiers
Center awards 1430284 and 2020265. AR, AP, and MB acknowledge the support from
the Ministero degli Affari Esteri e della Cooperazione Internazionale -
Direzione Generale per la Promozione del Sistema Paese - Progetto di Grande
Rilevanza ZA18GR02. Part of this work has been funded using resources from the
research grant “iPeska” (P.I. Andrea Possenti) funded under the INAF national
call Prin-SKA/CTA approved with the Presidential Decree 70/2016. Pulsar
research at UBC is supported by an NSERC Discover Grant and by the Canadian
Institute for Advanced Research. AR acknowledges continuing valuable support
from the Max-Planck Society. J.W.T.H. acknowledges funding from an NWO Vici
grant (”AstroFlash”; VI.C.192.045). M.D. acknowledges support from the
National Science Foundation (NSF) Physics Frontier Center award 1430284, and
from several NASA grants at the Naval Research Laboratory. TRAPUM observations
usedthe FBFUSE and APSUSE computing clusters for data acquisition,storage and
analysis. These clusters were funded and installed bythe Max-Planck-Institut
für Radioastronomie and the Max-Planck-Gesellschaft. L.V. acknowledges
financial support from the Dean’s Doctoral Scholar Award from the University
of Manchester.  Table 2: New MeerKAT Discovered M28 Pulsars

Pulsar Name J1824$-$2452M J1824$-$2452N Right Ascension (RA, J2000) $18^{\rm
h}\;24^{\rm m}\;33\fs 1835(5)$ $18^{\rm h}\;24^{\rm m}\;33\fs 1418(12)$
Declination (DEC, J2000) $-24^{\circ}\;52^{\prime}\;08\farcs 2(23)$
$-24^{\circ}\;52^{\prime}\;11\farcs 89(3)$ Proper Motion Right Ascension
(mas/yr) $0.0(1.4)$ $0.5(0.3)$ Proper Motion Declination (mas/yr) $-23(28)$
$-15(6)$ Pulsar Period (ms) 4.7842842982785(7) 3.35287243686099(10) Pulsar
Frequency (Hz) 209.01767906222(3) 298.251728579395(9) Frequency Derivative (Hz
s-1) $-5.36572(12)\times 10^{-15}$ $-1.416227(3)\times 10^{-14}$ Frequency
Second Derivative (Hz s-2) $-1.4(4)\times 10^{-26}$ $-9.4(1.2)\times 10^{-27}$
Reference Epoch (MJD) 56450 56450 Dispersion Measure (pc cm-3) 119.352(11)
119.321(3) Orbital Period (days) 0.2425192190(14) 0.19849331510(11) Orbital
Period Derivative $3.1(1.6)\times 10^{-12}$ $6.2(1.4)\times 10^{-13}$
Projected Semi-Major Axis (lt-s) 0.032463(8) 0.0497303(16) Time of Ascending
Node 56451.272704(15) 56451.2896713(15) Properties of Fit Span of Timing Data
(MJD) 53629$-$59274 53629$-$59274 Number of TOAs 161 230 RMS TOA Residual
($\mu$s) 71.11 16.45 Derived Parameters Mass Function ($M_{\odot}$)
$6.245(4)\times 10^{-7}$ $3.3516(3)\times 10^{-6}$ Min Companion Mass
($M_{\odot}$) $\geq$ 0.011 $\geq$ 0.019 Pulsar Magnetic Field (G) $\dagger$
$<4\times 10^{8}$ $3.6\times 10^{8}$ Pulsar Age (years) $\dagger$ $>3\times
10^{9}$ $1.4\times 10^{9}$ Pulsar $\dot{E}$ (erg s-1) $\dagger$ $<8\times
10^{33}$ $4.0\times 10^{34}$ Flux Density at 2 GHz ($\mu$Jy) 6(3) 10(5)
Angular Offset of Cluster Center (mas) 5.092 3.464 Total Proper Motion
(mas/yr) ${\ddagger}$ 10(22) 14(7)

Note. — Numbers in parentheses represent 1-$\sigma$ uncertainties in the last
digit. The timing solution was determined using TEMPO with the DE440 (Park et
al., 2021) Solar System Ephemeris and the ELL1 binary model. The time system
used is Barycentric Dynamical Time (TDB), referenced to TT via BIPM. The
minimum companion mass was calculated assuming a pulsar mass of 1.4
$M_{\odot}$. $\dagger$ Here we assumed a 95% lower limit to $\dot{P}_{b}$
which allowed us to solve for an upper limit on the intrinsic $\dot{P}$ and
therefore a lower limit on the age and upper limits on the magnetic field and
$\dot{E}$ (see §4.2). ${\ddagger}$ The relatively poorly constrained proper
motions are nevertheless consistent with overall cluster proper motion
measured by Gaia (see §4.3).

## References

  * Andersen & Ransom (2018) Andersen, B. C., & Ransom, S. M. 2018, ApJ, 863, L13, doi: 10.3847/2041-8213/aad59f
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Bagchi et al. (2011) Bagchi, M., Lorimer, D. R., & Chennamangalam, J. 2011, Monthly Notices of the Royal Astronomical Society, 418, 477–489, doi: 10.1111/j.1365-2966.2011.19498.x
  * Bagchi et al. (2013) Bagchi, M., Lorimer, D. R., & Wolfe, S. 2013, Monthly Notices of the Royal Astronomical Society, 432, 1303–1314, doi: 10.1093/mnras/stt559
  * Bailes et al. (2020) Bailes, M., Jameson, A., Abbate, F., et al. 2020, PASA, 37, e028, doi: 10.1017/pasa.2020.19
  * Baumgardt & Vasiliev (2021) Baumgardt, H., & Vasiliev, E. 2021, MNRAS, 505, 5957, doi: 10.1093/mnras/stab1474
  * Becker et al. (2003) Becker, W., Swartz, D. A., Pavlov, G. G., et al. 2003, ApJ, 594, 798, doi: 10.1086/376967
  * Bégin (2006) Bégin, S. 2006, Master’s thesis, University of British Columbia, doi: 10.14288/1.0085181
  * Bhattacharya et al. (2017) Bhattacharya, S., Heinke, C. O., Chugunov, A. I., et al. 2017, MNRAS, 472, 3706, doi: 10.1093/mnras/stx2241
  * Bogdanov et al. (2006) Bogdanov, S., Grindlay, J. E., Heinke, C. O., et al. 2006, ApJ, 646, 1104, doi: 10.1086/505133
  * Bogdanov et al. (2011) Bogdanov, S., van den Berg, M., Servillat, M., et al. 2011, ApJ, 730, 81, doi: 10.1088/0004-637X/730/2/81
  * Booth & Jonas (2012) Booth, R. S., & Jonas, J. L. 2012, African Skies, 16, 101
  * Bussa & VEGAS Development Team (2012) Bussa, S., & VEGAS Development Team. 2012, in American Astronomical Society Meeting Abstracts, Vol. 219, American Astronomical Society Meeting Abstracts #219, 446.10
  * Cadelano et al. (2020a) Cadelano, M., Chen, J., Pallanca, C., et al. 2020a, ApJ, 905, 63, doi: 10.3847/1538-4357/abc345
  * Cadelano et al. (2019) Cadelano, M., Ferraro, F. R., Istrate, A. G., et al. 2019, ApJ, 875, 25, doi: 10.3847/1538-4357/ab0e6b
  * Cadelano et al. (2018) Cadelano, M., Ransom, S. M., Freire, P. C. C., et al. 2018, ApJ, 855, 125, doi: 10.3847/1538-4357/aaac2a
  * Cadelano et al. (2020b) Cadelano, M., Saracino, S., Dalessandro, E., et al. 2020b, ApJ, 895, 54, doi: 10.3847/1538-4357/ab88b3
  * Cadelano et al. (2015) Cadelano, M., Pallanca, C., Ferraro, F. R., et al. 2015, ApJ, 807, 91, doi: 10.1088/0004-637X/807/1/91
  * Camilo et al. (2000) Camilo, F., Lorimer, D. R., Freire, P., Lyne, A. G., & Manchester, R. N. 2000, ApJ, 535, 975, doi: 10.1086/308859
  * Danner et al. (1994) Danner, R., Kulkarni, S. R., & Thorsett, S. E. 1994, ApJ, 436, L153, doi: 10.1086/187655
  * DuPlain et al. (2008) DuPlain, R., Ransom, S., Demorest, P., et al. 2008, in Advanced Software and Control for Astronomy II, ed. A. Bridger & N. M. Radziwill, Vol. 7019, International Society for Optics and Photonics (SPIE), 496 – 505, doi: 10.1117/12.790003
  * Freire & Ridolfi (2018) Freire, P. C. C., & Ridolfi, A. 2018, MNRAS, 476, 4794, doi: 10.1093/mnras/sty524
  * Freire et al. (2017) Freire, P. C. C., Ridolfi, A., Kramer, M., et al. 2017, MNRAS, 471, 857, doi: 10.1093/mnras/stx1533
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv e-prints, arXiv:2012.01533.  https://arxiv.org/abs/2012.01533
  * Hessels et al. (2006) Hessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Science, 311, 1901, doi: 10.1126/science.1123430
  * Kaplan et al. (2018) Kaplan, D. L., Stovall, K., van Kerkwijk, M. H., Fremling, C., & Istrate, A. G. 2018, ApJ, 864, 15, doi: 10.3847/1538-4357/aad54c
  * Kaplan et al. (2005) Kaplan, D. L., Escoffier, R. P., Lacasse, R. J., et al. 2005, Publications of the Astronomical Society of the Pacific, 117, 643.  http://www.jstor.org/stable/10.1086/430368
  * Lyne et al. (1987) Lyne, A. G., Brinklow, A., Middleditch, J., Kulkarni, S. R., & Backer, D. C. 1987, Nature, 328, 399, doi: 10.1038/328399a0
  * Lyne et al. (2000) Lyne, A. G., Mankelow, S. H., Bell, J. F., & Manchester, R. N. 2000, MNRAS, 316, 491, doi: 10.1046/j.1365-8711.2000.03517.x
  * Lyne et al. (1990) Lyne, A. G., Manchester, R. N., D’Amico, N., et al. 1990, Nature, 347, 650, doi: 10.1038/347650a0
  * Manchester et al. (1991) Manchester, R. N., Lyne, A. G., Robinson, C., et al. 1991, Nature, 352, 219, doi: 10.1038/352219a0
  * Pallanca et al. (2013) Pallanca, C., Dalessandro, E., Ferraro, F. R., Lanzoni, B., & Beccari, G. 2013, ApJ, 773, 122, doi: 10.1088/0004-637X/773/2/122
  * Pallanca et al. (2014) Pallanca, C., Ransom, S. M., Ferraro, F. R., et al. 2014, ApJ, 795, 29, doi: 10.1088/0004-637X/795/1/29
  * Pallanca et al. (2010) Pallanca, C., Dalessandro, E., Ferraro, F. R., et al. 2010, ApJ, 725, 1165, doi: 10.1088/0004-637X/725/1/1165
  * Pan et al. (2016) Pan, Z., Hobbs, G., Li, D., et al. 2016, MNRAS, 459, L26, doi: 10.1093/mnrasl/slw037
  * Papitto et al. (2013) Papitto, A., Ferrigno, C., Bozzo, E., et al. 2013, Nature, 501, 517, doi: 10.1038/nature12470
  * Park et al. (2021) Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. 2021, AJ, 161, 105, doi: 10.3847/1538-3881/abd414
  * Phinney (1992) Phinney, E. S. 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 39, doi: 10.1098/rsta.1992.0084
  * Ransom (2011a) Ransom, S. 2011a, PRESTO: PulsaR Exploration and Search TOolkit.  http://ascl.net/1107.017
  * Ransom (2011b) —. 2011b, PRESTO: PulsaR Exploration and Search TOolkit.  http://ascl.net/1107.017
  * Ransom (2001) Ransom, S. M. 2001, PhD thesis, Harvard University
  * Ransom (2008) Ransom, S. M. 2008, in Dynamical Evolution of Dense Stellar Systems, ed. E. Vesperini, M. Giersz, & A. Sills, Vol. 246, 291–300, doi: 10.1017/S1743921308015810
  * Ransom et al. (2005) Ransom, S. M., Hessels, J. W. T., Stairs, I. H., et al. 2005, Science, 307, 892, doi: 10.1126/science.1108632
  * Ridolfi et al. (2021) Ridolfi, A., Gautam, T., Freire, P. C. C., et al. 2021, MNRAS, doi: 10.1093/mnras/stab790
  * Rivera Sandoval et al. (2018) Rivera Sandoval, L. E., van den Berg, M., Heinke, C. O., et al. 2018, MNRAS, 475, 4841, doi: 10.1093/mnras/sty058
  * Roberts (2013) Roberts, M. S. E. 2013, in Neutron Stars and Pulsars: Challenges and Opportunities after 80 years, ed. J. van Leeuwen, Vol. 291, 127–132, doi: 10.1017/S174392131202337X
  * Robinson et al. (1995) Robinson, C., Lyne, A. G., Manchester, R. N., et al. 1995, MNRAS, 274, 547, doi: 10.1093/mnras/274.2.547
  * Shklovskii (1970) Shklovskii, I. S. 1970, Soviet Ast., 13, 562
  * Stappers & Kramer (2016) Stappers, B., & Kramer, M. 2016, in MeerKAT Science: On the Pathway to the SKA, 9
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191, doi: 10.1086/131977
  * van Straten et al. (2011) van Straten, W., Demorest, P., Khoo, J., et al. 2011, PSRCHIVE: Development Library for the Analysis of Pulsar Astronomical Data.  http://ascl.net/1105.014
  * Vasiliev & Baumgardt (2021) Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978, doi: 10.1093/mnras/stab1475
  * Verbunt & Freire (2014) Verbunt, F., & Freire, P. C. C. 2014, A&A, 561, A11, doi: 10.1051/0004-6361/201321177
  * Zhao et al. (2019) Zhao, Y., Heinke, C. O., Cohn, H. N., Lugger, P. M., & Cool, A. M. 2019, MNRAS, 483, 4560, doi: 10.1093/mnras/sty3384
  * Zhao et al. (2020) Zhao, Y., Heinke, C. O., Cohn, H. N., et al. 2020, MNRAS, 499, 3338, doi: 10.1093/mnras/staa2927
  * Zharikov et al. (2019) Zharikov, S., Kirichenko, A., Zyuzin, D., Shibanov, Y., & Deneva, J. S. 2019, MNRAS, 489, 5547, doi: 10.1093/mnras/stz2475