# MeerKAT discovery of 13 new pulsars in Omega Centauri

W. Chen,1 P. C. C. Freire,1 A. Ridolfi,3,1 E. D. Barr,1 B. Stappers,2 M.
Kramer,1,2 A. Possenti,3 E-mail<EMAIL_ADDRESS>S. M. Ransom,4 L.
Levin,2 R. P. Breton,2 M. Burgay,3 F. Camilo,5 S. Buchner,5 D. J. Champion,1
F. Abbate,1 V. Venkatraman Krishnan,1 P. V. Padmanabh,1,8,9 T. Gautam,1 L.
Vleeschower,2 M. Geyer,5 J-M. Grießmeier,6,7 Y. P. Men,1 V. Balakrishnan,1 M.
C. Bezuidenhout2
1Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn,
Germany
2Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
The University of Manchester, Manchester M13 9PL, UK
3INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047
Selargius (CA), Italy
4National Radio Astronomy Observatory (NRAO), 520 Edgemont Rd.,
Charlottesville, VA 22903 USA
5South African Radio Astronomy Observatory (SARAO), 2 Fir Street, Black River
Park, Observatory, Cape Town 7925, South Africa
6LPC2E - Université d’Orléans / CNRS, 45071 Orléans cedex 2, France
7Observatoire Radioastronomique de Nançay (ORN), Observatoire de Paris,
Université PSL, Univ Orléans, CNRS, 18330 Nancay, France
8 Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
D-30167 Hannover, Germany
9 Leibniz Universität Hannover, D-30167 Hannover, Germany

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The most massive globular cluster in our Galaxy, Omega Centauri, is an
interesting target for pulsar searches, because of its multiple stellar
populations and the intriguing possibility that it was once the nucleus of a
galaxy that was absorbed into the Milky Way. The recent discoveries of pulsars
in this globular cluster and their association with known X-ray sources was a
hint that, given the large number of known X-ray sources, there is a much
larger undiscovered pulsar population. We used the superior sensitivity of the
MeerKAT radio telescope to search for pulsars in Omega Centauri. In this
paper, we present some of the first results of this survey, including the
discovery of 13 new pulsars; the total number of known pulsars in this cluster
currently stands at 18. At least half of them are in binary systems and
preliminary orbital constraints suggest that most of the binaries have light
companions. We also discuss the ratio between isolated and binaries pulsars
and how they were formed in this cluster.

###### keywords:

Pulsar – Globular cluster – Binary

††pubyear: 2015††pagerange: MeerKAT discovery of 13 new pulsars in Omega
Centauri–LABEL:lastpage

##  1 Introduction

Pulsar surveys conducted in Globular clusters (GCs) have yielded fruitful
rewards in recent decades, with the discovery of a total of 261
pulsars111https://www3.mpifr-bonn.mpg.de/staff/pfreire/GCpsr.html. Per unit of
stellar mass, GC have three orders of magnitude more pulsars than the Galactic
disk. The reason for this is their large stellar densities: these prompt
stellar interactions (Verbunt & Hut, 1987) where old, dead Neutron stars (NSs)
acquire new main sequence (MS) companions. These then evolve, forming X-ray
binaries (also exceptionally numerous in GCs, Clark 1975), where the NS is
being spun up by accretion of matter from the MS star. When the accretion
stops, these systems become millisecond pulsar (MSP) binaries, which have
nearly circular orbits and low-mass companions (Bhattacharya & van den Heuvel,
1991). Indeed, the pulsar population in GCs is dominated by such binaries.

However, in some cases, additional exchange encounters can originate different
end products. If they happen during the X-ray binary phase, they can disrupt
the binary and produce many single and/or partially recycled pulsars, which
are not only slower, but appear to be much younger than the GC population
(Verbunt & Freire, 2014). Such encounters can also replace a pulsar’s
companion by a much more massive degenerate object, resulting in massive,
eccentric binary MSPs unlike any seen in the Galactic disk (e.g., Freire et
al. 2004; Lynch et al. 2012; DeCesar et al. 2015; Ridolfi et al. 2021, 2022).
If the massive companions happen to be stellar-mass black holes, these systems
could be superb test-beds for fundamental physics (Liu et al., 2014).

Such exotic pulsar binaries are generally observed in GCs with very dense
cores, especially core-collapsed clusters; these are the types of environments
where each particular binary is likely to go through more than one disruptive
stellar encounter (Verbunt & Freire, 2014). Thus, the pulsar population in a
GC, and the types of binaries the pulsars find themselves in, reflects not
only its current dynamical status, but also the cluster’s previous evolution
(Benacquista & Downing, 2013).

###  1.1 The Omega Centauri globular cluster

Omega Centauri ($\omega$-Cen, also known as NGC 5139), the largest GC in our
Galaxy, is a natural target to search for pulsars. It is located in the
constellation of Centaurus and is 5.2 kpc away from the Sun, with an age of
11.52 Gyr (Forbes & Bridges, 2010). Besides its size and large number of
stars, it differs from other GCs because of its intricate composition of
different populations of stars (Bedin et al., 2004). This could indicate that
$\omega$-Cen is the merger of several clusters, like the Sagittarius dwarf
galaxy (Ibata et al., 1994). Moreover, it was found to be rich in calcium and
heavy metals (Lee et al., 2009), which is a tracer of supernovae explosions.
However these materials ejected by the explosion could not be sustained by the
current gravitational potential of the cluster. This along with the multiple
stellar populations supports the long established idea that $\omega$-Cen is
the relic of a former disrupted dwarf galaxy (Hilker & Richtler, 2000; Ibata
et al., 2019).

Unassociated high-energy emission in GCs is thought to originate from MSPs
(Venter et al., 2009) as supported by observations, (e.g. Abdo et al. 2009).
$\gamma$-ray emission has been detected in several GCs by the Fermi Large Area
Telescope (Abdo et al., 2010), including those that, at the time, had no
previously known pulsars, such as $\omega$-Cen.

Additionally, $\sim$30 unassociated X-ray sources have been found within the
core of $\omega$-Cen to have luminosities similar to the emission of pulsars
in other GCs (Henleywillis et al., 2018). However, previous searches for
pulsars in this cluster turned out to be unsuccessful (Edwards et al., 2001;
Possenti et al., 2005; Camilo et al., 2015).

Notwithstanding, in 2019, Dai et al. (2020) carried out a search for pulsars
using the new Ultra-wide Bandwidth receiver (UWL, Hobbs et al., 2020) of the
64-m "Murriyang" radio telescope at Parkes, NSW, Australia. This finally
allowed the discovery of the first 5 MSPs in $\omega$-Cen. Among these, PSR
J1326$-$4728B is in an eclipsing binary system with a light companion, and it
is associated with an X-ray source (Henleywillis et al., 2018). The authors
suggested that the non-detection of pulsars from the previous searches on this
cluster was caused by the lack of sensitivity of previous surveys. Given the
large distance of $\omega$-Cen (and other GCs in general) it is clear that we
are only detecting the very brightest pulsars in them. This means that more
sensitive telescopes would in principle detect many more pulsars in this and
other globular clusters.

###  1.2 The MeerKAT survey

The MeerKAT 64-antenna array, located in the Karoo desert in South Africa,
(Jonas & MeerKAT Team, 2016; Camilo, 2018) started scientific observations in
late 2019, becoming, by far, the most radio sensitive radio telescope in the
Southern Hemisphere. Numerous targeted observations and surveys of pulsars
have been carried out since and are producing significant results. For
example, the TRAnsients and PUlsars with MeerKAT (TRAPUM, Stappers & Kramer,
2016) Large Survey Project aims to increase the total known population of
pulsars, and discover peculiar binary pulsars that might be suitable for
studies of fundamental physics and stellar evolution. This project has so far
discovered 156 pulsars222http://www.trapum.org/discoveries/. About one-third
of these (e.g. Ridolfi et al. 2021; Douglas et al. 2022; Ridolfi et al. 2022;
Vleeschower et al. 2022; Abbate et al. 2022), are found in various GCs,
including 47 Tucanae (Chen et al., in prep.), a cluster that has been searched
for over 20 years.

Given the small beams produced by the phased interferometer, pulsars surveys
need a beamformer in order to cover as much sky simultaneously as possible.
The beam former developed by Barr (2018) and Chen et al. (2021) can generate
up to 1000 coherent beams that cover a significant part of the primary field
of view of an individual antenna, with high time and frequency resolution.
Such a large number of beams is essential for blind surveys (which cover
mostly the Galactic plane), but they are also important for GCs. Indeed, even
though GCs have a relatively small angular size on the sky, we still often
need hundreds of beams to cover their half-mass radii, which are the regions
where most pulsars are likely to be located. This is especially true in the
case of $\omega$-Cen, which is the most massive and has the largest projected
size among all known GCs.

In this paper, we report the result of the searches for the new pulsars in two
observations of $\omega$-Cen made with MeerKAT. Section 2 describes the
observational parameters and analysis of the recorded data. Section 3 presents
the new discoveries, while Section 4 discusses the uncovered pulsar population
and its implications in our understanding of the cluster.

##  2 Observations and data analysis

###  2.1 Observations

We carried out two multibeam observations of $\omega$-Cen with MeerKAT on 2021
March 21 and 26, as part of the TRAPUM GC survey. The cluster was observed for
4 hours in each observation using the L-band receivers, which cover the
856–1712 MHz frequency range. The coherent beams were synthesised using the
Filterbanking BeamFormer User-Supplied Equipment (FBFUSE, Barr, 2018). A
tiling of 704 coherent beams was generated by
mosaic333https://github.com/wchenastro/Mosaic to cover the cluster, as shown
in Figure 1. This was centred at the nominal centre of $\omega$-Cen, at
equatorial coordinates (J2000) 13:26:47.24, $-$47:28:46.5 and Galactic
coordinates 309.10, 14.97 (Harris, 2010), covering a circular region of 7.53
arcmin in radius, i.e. roughly twice the size of the half-light radius of the
cluster. During both observations, 60 antennas were used. To create a
hexagonal tiling, the corresponding synthesised beam shape was approximated by
an ellipse whose major and minor axis were 20.46 arcseconds and 9.32
arcseconds at the middle of the observation on March 21. The method and
performance of such approximation is discussed in Chen et al. (2021). The
March 26 observation started at almost the same hour angle, so the beam shape
is very similar.

For each beam, the observing band was split into 2048 channels and recorded
every 153.12 $\mu$s as filterbank search-mode files by Accelerated Pulsar
Search User Supplied Equipment (APSUSE) computing cluster. Once the
observation was over, the frequency channels were summed in groups of 16 after
being incoherently de-dispersed at the dispersion measure (DM) of 97 pc cm-3,
about the average of the previously known 5 pulsars in $\omega$-Cen. The
resulting sub-banded data preserved 128 channels, significantly reducing the
total data volume. With this resolution, the intra-channel smearing time is
1.57 ms, at the lowest channel and the upper bound of the search DM.

###  2.2 Sensitivity

We calculate a minimum detectable flux density of S${}_{\textrm{min}}$ = 10
$\mu$Jy for our search on a 4-hour observation, following the modified
radiometer equation given by Dewey et al. (1985):

$\displaystyle
S_{\textrm{min}}=\frac{\textrm{S/N}\,\beta\,T_{\textrm{sys}}}{\mathcal{E}_{\textrm{FFT}}\,G\sqrt{n_{\textrm{pol}}\,B\,\Delta\,t_{\textrm{obs}}}}\sqrt{\frac{\zeta}{1-\zeta}}$
(1)

where S/N is the minimal signal-to-noise of 10 for a valid detection; $\beta$
is the correction factor of 1.01 to compensate the loss of sensitivity during
the digitization process; $T_{\textrm{sys}}$ is the system temperature of
about 26
K444https://skaafrica.atlassian.net/rest/servicedesk/knowledgebase/latest/articles/view/277315585,
including the contribution from the receiver temperature of 18 K, the sky
temperature of about 3.5 K at L-band and the combination of the spillover
noise and the atmosphere of about 4.5 K; $\mathcal{E}_{\mathrm{FFT}}$ is the
FFT search efficiency which is 0.7 according to Morello et al. (2020); G is
the combined gain of the array, which is 2.65 K Jy-1 when 60 antennas are
used; $n_{\textrm{pol}}$ is the number of polarizations, which is 2 while
using the total power beamformer in TRAPUM observations; $B$ is the effective
receiver bandwidth of about 640 MHz, after removing the channels affected by
RFI; $\Delta t_{\textrm{obs}}$ is the integration time which was 4 hours for
both observations; $\zeta$ is the pulse’s apparent duty cycle, which we
consider to be 8% with the broadening effects from interstellar medium and the
instrument.

###  2.3 Data reduction and search pipeline

The data were transferred to Germany from South Africa in hard drives.
Searches for periodic signals were carried out using the Hercules computing
cluster555https://docs.mpcdf.mpg.de/doc/computing/clusters/systems/Radioastronomy.html.
For this work, we restricted our search to the beams located within the half-
light radius of $\omega$-Cen. The search of an individual beam started with an
RFI mitigation procedure using filtool666https://github.com/ypmen/PulsarX. It
applies two major operations on the filterbank data according to the
statistical analysis. Assuming that the noise follows a Gaussian distribution,
it calculates the kurtosis of samples in certain time unit across the full
band, to obtain the deviation of each channel then replace the outlying ones
with the mean. Similarly, it calculates the skewness of the samples to
identify RFIs because they are often non-Gaussian. After this, it normalizes
the data across channels so that each channel has an average of 0 and a
variance of 1.

The resulting filterbank data were fed to a pulsar search pipeline supervised
by pulsar_miner777https://github.com/alex88ridolfi/PULSAR_MINER. The pipeline
is based on various programs and utilities from
presto888https://github.com/scottransom/presto (Ransom, 2011), a pulsar search
and analysis toolkit. It first uses rfifind to examine the data and flag the
narrow and wide band interference. The operation outputs a report and a mask
for later use. Free electrons in the interstellar medium engaging with passing
electromagnetic waves leads to arrival delays across channels, which are
quantified as Dispersion Measure (DM). Thus, signals with 0 DM are terrestrial
which were identified by prepdata and zapped. After that, the data were split
into segments of different lengths, to be sensitive to binary pulsars with
different orbital periods (Ransom et al., 2003). The chosen lengths of the
segments were 10, 20, 30, 60, 120 minutes, in addition to the full length (4
h) of each observation. Each of these segments were de-dispersed into time
series using a number of trial DMs. The step between these DMs were set to 0.1
pc cm-3 which were determined using DDplan.py, a program that generates
suitable de-dispersion schemes considering the DM range, central frequency,
sampling time, number of channels and other parameters. The final range of DMs
searched was 85–115 pc cm-3, which is slightly wider than the range of 90–110
used in Dai et al. (2020). Each time series was then searched for periodic
signals in the Fourier domain using accelsearch (Ransom et al., 2002). In
order to account for the drifting of the frequency of a signal due to orbital
motion, accelsearch searches for signals that drift linearly over multiple
Fourier bins. The maximum number of bins drifted is set by the zmax parameter,
which was chosen to be 200 in our search. The candidates from the search went
through a sifting process, where only the ones with S/Ns higher than 4 and
within the interested DM and period ranges were delivered to the next stage.
In the final step, the sifted candidates were folded using prepfold with an
extra optimization on the period and DM. The plots of these candidates were
inspected by eye.

##  3 Discoveries

Our search resulted in the discovery of 13 new pulsars in the beams that were
analyzed in this work (shown with blue edges in Figure 1). Seven of them are
in binary systems and we are able to place orbital constraints based on two
observations. The results show that their orbital periods fall into two
groups: below 4 hours and around 1 day. Notably, all except one have light
companions. Other than the new discoveries, we also re-detected all the
previously known pulsars (Dai et al., 2020) in the cluster. The properties of
all the new pulsars and which segment they were discovered, are listed in
Table 1, while their integrated profiles are shown in Figure 2. The
distribution of their spin periods are presented in Figure 4. In the remainder
of the section, we discuss their characteristics in more detail.

Figure 1: Tiling and detections. Shown is the beam tiling pattern at the start
of the observation generated by mosaic on sky, centred at the optical centre
(denote with a cross) of $\omega$-Cen. The radii of the core and the half
light (Harris, 2010) are denoted using dashed line and dash-dotted lines. The
beams with blue edges were searched for pulsars in this work, roughly cover
the region within the half light radius. Lower right is a zoom-out view of the
tiling. The position indicated by green with black edges are known sources
published in (Dai et al., 2020). The positions indicated by black are the new
discoveries in this paper. Bold face indicates that they have timing positions
(A and B) or their positions have been constrained using multibeam detections.
Others were placed at the centers of the beams where they have the brightest
detections. Table 1: List of discoveries from this work and their properties.
The DMs were give by prefold, the orbital parameters were derived using tempo,
the positions were localized using SeeKAT and the numbers in parentheses of
the positions represent 2-$\sigma$ uncertainty of the last digit. The search
lengths are the length of segments which the pulsars were discovered (see
section 2.3).

| Pulsar | Type | $P$ | DM | $P_{b}$ | $x_{p}$ | $M^{\textrm{min}}_{c}$ | $\alpha$ | $\delta$ | Search length
---|---|---|---|---|---|---|---|---|---|---
|  |  | (ms) | (pc cm-3) | (d) | (lt-s) | M⊙ | J2000 | J2000 | (h)
| F | Isolated | 2.27 | 98.29 | - | - | - | $13^{\textrm{h}}26^{\textrm{m}}53^{\textrm{s}}(1)$ | $-47\degr 28\arcmin 28\arcsec(6)^{{\ddagger}}$ | 4
| G | Binary | 3.30 | 99.69 | 0.1087597(1) | 0.032203(7) | 0.018 | $13^{\textrm{h}}26^{\textrm{m}}37^{\textrm{s}}.1(2)$ | $-47\degr 29\arcmin 41\arcsec(1)$ | 0.5
| H | Binary | 2.52 | 98.09 | 0.1356948(5) | 0.021844(9) | 0.011 | $13^{\textrm{h}}26^{\textrm{m}}44^{\textrm{s}}(1)$ | $-47\degr 28\arcmin 55\arcsec(4)$ | 0.5
| I | Binary | 18.95 | 102.2 | 1.113(1)♭ | 0.165(1) | 0.020 | $13^{\textrm{h}}26^{\textrm{m}}29^{\textrm{s}}.0(1)$ | $-47\degr 30\arcmin 24\arcsec(1)$ | 4
| J | Isolated | 1.84 | 97.28 | - | - | - |  $13^{\textrm{h}}26^{\textrm{m}}51^{\textrm{s}}.7$(1) | $-47\degr 27\arcmin 09\arcsec(1)$ | 4
| K | Binary | 4.72 | 94.73 | 0.09387146(2) | 0.067945(5) | 0.043 | $13^{\textrm{h}}26^{\textrm{m}}38^{\textrm{s}}.1(1)$ | $-47\degr 27\arcmin 39\arcsec(2)$ | 0.5
| L | Binary | 3.54 | 101.5 | 0.1589282(4) | 0.061809(9) | 0.027 | $13^{\textrm{h}}27^{\textrm{m}}02^{\textrm{s}}.8(1)$ | $-47\degr 26\arcmin 49\arcsec(2)$ | 0.5
| M | Isolated | 4.60 | 101.4 | - | - | - | $13^{\textrm{h}}26^{\textrm{m}}59^{\textrm{s}}(1)$ | $-47\degr 30\arcmin 09\arcsec(6)^{{\ddagger}}$ | 4
| N | Binary | 6.88 | 101.2 | 1.0816(3)♭ | 0.1250(3) | 0.015 | $13^{\textrm{h}}26^{\textrm{m}}49^{\textrm{s}}.8(1)$ | $-47\degr 31\arcmin 25\arcsec(1)$ | 4
| O | Isolated | 6.16 | 94.27 | - | - | - | $13^{\textrm{h}}26^{\textrm{m}}48^{\textrm{s}}(1)$ | $-47\degr 27\arcmin 19\arcsec(6)^{{\ddagger}}$ | 4
| P | Isolated | 2.79 | 102.1 | - | - | - | $13^{\textrm{h}}26^{\textrm{m}}45^{\textrm{s}}(1)$ | $-47\degr 29\arcmin 42\arcsec(6)^{{\ddagger}}$ | 4
| Q | Binary | 4.13 | 95.92 | 1.18(8)♭ | 1.1(1) | 0.138 | $13^{\textrm{h}}26^{\textrm{m}}35^{\textrm{s}}(1)$ | $-47\degr 27\arcmin 54\arcsec(6)^{{\ddagger}}$ | 2
| R | Isolated | 10.29 | 102.1 | - | - | - | $13^{\textrm{h}}27^{\textrm{m}}10^{\textrm{s}}(1)$ | $-47\degr 29\arcmin 02\arcsec(6)^{{\ddagger}}$ | 4

Note:

♭The orbits of these pulsars are longer than the observation, the values shown
here are derived from two observations, hence they are subject to change when
more data are available.

‡ The position of these pulsars was denoted by the centers of the beams where
they are detected with highest S/N. Because of the lack of or faint detection
in their neighbouring beams, further localization is not practical with
SeeKAT.

|  |  |  |  |  |
---|---|---|---|---|---|---
|  |  |  |  |  |
Figure 2: Profiles of the new pulsars folded with 4 hours of data taken on 26
March 2021. The Y-axis is intensity, frequency and time from top to bottom
panels, and X-axis is the phase window from 0 to 1.

###  3.1 Isolated pulsars

From our two observations, PSRs J1326$-$4728F, J, M, O, P and R appear to be
isolated, because they do not show changes to their barycentric period to
within the limits of the uncertainty. Most of them are within the core radius
of the cluster (located at 2.37 arcmin from the nominal center of
$\omega$-Cen), while pulsar M is near the edge of the core and pulsar R is
between the core and half light radius (at 5.00 arcmin). Among them,
J1326$-$4728J has a period of 1.84 ms, making it the fastest spinning pulsar
so far discovered in this cluster. It also has a large duty-cycle of 43.7%.
For J1326$-$4728O and P, there were harmonic detections with similar S/N, it
is ambiguous as to the number of peaks in the profile and thus the spin period
is also ambiguous, further observations with polarimetric capability should
resolve these ambiguities. The isolated pulsars discovered in these two
observations are relatively weak compared to the re-detections of the known
pulsars discovered in Dai et al. (2020), suggesting that the sensitivity could
be one of the reasons why the new pulsars were not detected in the previous
observations with Parkes.

###  3.2 Binary pulsars

More than half of the discoveries in this work reside in binaries. To obtain
their orbital parameters, the data were first split into segments, from which
time-dependent barycentric period and period derivative estimates could be
made. An example of such practice for J1326$-$4728G is shown in Figure 3. With
these segmental sets of values, the orbital parameters were fit assuming a
circular orbit, using the fit_circular_orbit.py Python script from presto.
These solutions were then used as the initial guess for the orbital parameter
when constructing an initial ephemeris for each binary pulsar. Subsequently,
these ephemerides were iteratively improved by pulsar timing using the
tempo999https://sourceforge.net/projects/tempo software package. According to
the current ephemerides, pulsar J1326$-$4728I, N and Q have orbits longer than
the observation, wherefore, the solution of these orbits are not unique and
subject to change where more data are available.

Figure 3: Changes of periods of J1326$-$4728G during the first four-hour
observation. The data points are 30- minute segments and the curve is an
orbital fit of the period changes given by fit_circular_orbit.py. Figure 4:
Distribution of spin period of pulsars in $\omega$-Cen.

J1326$-$4728G and H have periods of 3.30 and 2.52 ms, respectively. The former
was found near the edge of the core and was detected brightly in both epochs
and in neighbouring beams, the latter was found near the centre of the
cluster. Preliminary orbital solutions suggest that they have orbits of
approx. 2.61 hours and 2.36 hours, respectively. According to the mass
function and assuming the mass of the pulsar is 1.4 M⊙ (also for the following
paragraphs), the mass range of the companion of J1326$-$4728G is [0.018,
0.042] M⊙ when the inclination angle between the orbital plane and the light
of sight is $0\degr$ and $64\degr$. The upper bound was chosen such that the
range covers 90% of the cases (Lorimer & Kramer, 2004), assuming the orbital
planes are isotropically distributed. Similarly, the mass range of the
companion of J1326$-$4728H is [0.011, 0.024] M⊙. The orbital periods and
minimum companion masses fall within the typical range of expected values for
“black widow” systems (Roberts, 2013). Though, no eclipses, which are often
seen in such systems, are observed in these two observations. But considering
the tight orbits, low mass companions and rapid spin periods, it is very
likely they are black widows.

J1326$-$4728K and L have periods of 4.71 and 3.53 ms, respectively. The former
was discovered close to the edge of the core and was brightly detected in both
epochs and in neighbouring beams, the latter was discovered between the radii
of the core and the half light. An interpulse can be clearly observed in the
profile of J1326$-$4728K. Fits to the detected spin period and spin period
derivatives suggest orbits of about 2.25 and 3.81 hours, with the ranges of
companion masses of [0.043, 0.101] M⊙ and [0.027, 0.063] M⊙, respectively.
Again, these are values typical of black widow systems. For pulsars K and L,
however, eclipses are clearly seen in both observations.

J1326$-$4728I has a spin period of 18.95 ms with a relatively large duty cycle
of 29.7%. It is found between the radius of the core and of the half light of
the cluster. The profile can be phase-aligned well, assuming a constant line-
of-sight acceleration in each observation. However, the sign of the
acceleration changed between two observations. This suggests that the pulsar
is moving in a binary system with an orbit that is significant longer than the
length of our observations. The orbital solution suggests that it has an orbit
of 26.71 hours and its companion has a mass range of [0.020, 0.046] M⊙.
Additionally, the period of the pulsar indicates it might be a mildly recycled
pulsar. But without an accurate period derivative, other binary evolution
scenario cannot be ruled out. Further observations could shed light on the
nature of this pulsar.

J1326$-$4728N and Q have periods of 6.88 and 4.13 ms. They were both
discovered at the edge of the core. Their orbital parameters suggest orbits of
about 25.96 and 28.41 hours, with mass ranges of their companions of [0.015,
0.035] M⊙ and [0.138, 0.342] M⊙, respectively. It is difficult to constrain
the orbital parameters for J1326$-$4728N, which could be partially due to the
low and uneven orbital phase coverage during two observations. Pulsar Q was
brightly detected in one observation, but could not be re-detected in the
expected beam in the other observation. There are, however, very faint
detections in two of the neighbouring beams of that nearest beam. The
companion mass of Q indicates that it could be a helium white dwarf. Follow-up
observations are crucial for improving the orbital solutions of these two
pulsars.

###  3.3 Localization and correlation with X-ray emission

Figure 5: Localization and nearby X-ray sources. Top: localization result of
J1326$-$4728G, red cross is the best position given by SeeKAT, black cross is
the position of the nearest X-ray source 24f from Henleywillis et al. (2018).
Blue shades are the likelihood map of the position. The solid line and dashed
line are 1-$\sigma$ and 2-$\sigma$ confident levels of the result. Bottom:
localization result of J1326$-$4728K and the position of the nearest X-ray
source 21d.

About half of the new pulsars were detected in multiple neighbouring beams,
allowing their positions to be constrained to a precision better than the size
of the coherent beam. The localization was carried out using the
SeeKAT101010https://github.com/BezuidenhoutMC/SeeKAT multibeam localization
software (Bezuidenhout et al., submitted). It performs maximum-likelihood
estimation to obtain the best position. The likelihood is calculated by
testing the theoretical gain of the position given by a constant point spread
function (PSF), against the S/N’s of the neighbouring detections. The PSFs
were generated using mosaic with the observational parameters.

Our two observations lasted four hours each, and the PSF changed with time.
Significant changes of the PSF would lead to deterioration of localization
quality, therefore we tried to mitigate this effects by using the S/N’s and
PSFs in different segments of the observations. For faint pulsars, it is not
practical to obtain detections in short segments, in those cases, the PSFs
corresponding to the middle of the observation were used. The result of the
localization were shown in Table 1. Apart from the new pulsars, we also
perform localization with the known pulsars, because three of them have no
timing position at the moment. The results for the known pulsars are listed in
Table 2. Here, we present the localization of J1326$-$4728G and K plotted with
X-ray sources (24f and 21d) within the 2-$\sigma$ confident levels shown in
Figure 5. It is worth noting that multibeam localization relies strongly on
the accuracy of the PSF and the quality of the data, such as errors on the
S/N, severity of the RFI, coherence of the signal etc., all of which were not
considered in this estimation.

Previous observations made with the Parkes telescope (Dai et al., 2020)
associated one X-ray source from Henleywillis et al. (2018) to PSR
J1326$-$4728B. However, there are many X-ray sources in this cluster that are
still unassociated. Hence, we compared them with the new radio pulsar
discoveries and noticed some of the new pulsars have nearby X-ray sources
within a few arcseconds, such as the binaries pulsars J1326-4728G, H, K and L.
However, it is difficult to deduce a firm association until we obtain more
accurate radio timing position from follow-up observations. Since we now know
that this cluster hosts a considerable amount of pulsars, it is likely that
some of the X-ray emission from this cluster comes from pulsars.

Table 2: Positions of the known pulsars localized with SeeKAT. The numbers in
parentheses represent 2-$\sigma$ uncertainty of the last digit.

Pulsar | $\alpha$ | $\delta$
---|---|---
| J2000 | J2000
A | $13^{\textrm{h}}26^{\textrm{m}}39^{\textrm{s}}.7(1)$ | $-47\degr 30\arcmin 11\arcsec(2)^{{\ddagger}}$
B | $13^{\textrm{h}}26^{\textrm{m}}49^{\textrm{s}}.7(1)$ | $-47\degr 29\arcmin 26\arcsec(2)^{{\ddagger}}$
C | $13^{\textrm{h}}26^{\textrm{m}}55^{\textrm{s}}.5(1)$ | $-47\degr 30\arcmin 13\arcsec(2)$
D | $13^{\textrm{h}}26^{\textrm{m}}32^{\textrm{s}}.5(1)$ | $-47\degr 28\arcmin 39\arcsec(4)$
E | $13^{\textrm{h}}26^{\textrm{m}}42^{\textrm{s}}.6(1)$ | $-47\degr 27\arcmin 22\arcsec(1)$

Note:

‡ J1326$-$4728A and B have timing positions from Dai et al. (2020), which are
$13^{\textrm{h}}26^{\textrm{m}}39^{\textrm{s}}.670,-47\degr 30\arcmin
11\arcsec.64$ and $13^{\textrm{h}}26^{\textrm{m}}49^{\textrm{s}}.563,-47\degr
29\arcmin 24\arcsec.62$.

##  4 Discussion

###  4.1 The pulsar population in $\omega$-Cen

The number of pulsars known in $\omega$-Cen, which was 5 before this work, is
now 18. Taking the characteristics of the previously known pulsars, plus those
of Table 1 at face value, we see that the pulsar population of $\omega$-Cen
appears to be dominated by isolated pulsars. Indeed, the previously known
isolated pulsars A, C, D, E plus the newly discovered isolated pulsars, F, J,
M, O, P and R represent 10 out of a total of 18 pulsars. Furthermore, we can
also see that, with the exception of Q, all binary pulsars have very low-mass
companions, typical of what one finds in “black widow” systems. It has been
noted that black-widow pulsars are more easily formed in GCs (King et al.,
2003), but their very high faction in $\omega$-Cen is still surprising. Of
these, three systems (I, N and Q) seem to have orbital periods of just over 1
day, the other four have orbital periods smaller than the length of one
observation, 4 hours.

However, before we advance with a detailed analysis of the characteristics of
this population, it is important to keep in mind several caveats.

First, about the percentage of isolated pulsars, the low stellar density of
$\omega$-Cen means that wide binary systems with orbital periods of tens or
even hundreds of days might be stable in this cluster (as, for example, PSR
B1310+18A, a 255-day, low-eccentricity binary in the globular cluster M53, see
Kulkarni et al. 1991). This suggests that some of the pulsars that are
apparently isolated might, with additional observations, be found to be part
of wide binary systems, where the apparent changes of spin period caused by
the orbital motion are less pronounced.

Second, as mentioned above, we have analyzed in our search several stretches
of data going up to 4 hours. Finding binaries in such long integrations is
more difficult than finding isolated pulsars, because of the loss of
sensitivity caused by the orbital motion; this is especially true for cases
where the total integration time is of the same order as the orbital period.
This is the reason why we find, from the last column in Table 1, that while
all isolated pulsars and longer period binaries were discovered in 4-h
segments (with the exception of Q, which was found in a 2-h segment), no
short-period binaries were found in such long segments: they were instead
found, invariably, in 30-minute segments. This means that, in this survey, we
are $\sqrt{4/0.5}=2.83$ times more sensitive to very faint isolated MSPs (or
MSPs with very wide obits) than to short binaries. This means that the
isolated pulsars and binary pulsars with low accelerations are over-
represented in our sample relative to the short-period binaries.

Third, for normal MSP - WD binaries with orbital periods of a few days, no
orbits can be firmly determined with only two observations; to determine their
orbits additional observations will be necessary. We note in this regard that
the 1-day orbits for I, N and Q are preliminary, thus these binary MSPs could
in principle have more massive companions and longer orbital periods.

Summarizing, and taking these caveats into account, it does appear that a) the
number isolated MSPs (10 out of 18) is large, but represents an over-
estimation of the pulsar population of the cluster, because they are easier to
find and also because some of them could potentially be wide binaries; b) the
confirmed short-orbit systems (those with short orbital periods: B, G, H, K
and L) are more numerous than the other confirmed longer-period binaries (I, N
and Q); their numbers are likely to be under-estimated, since we are less
sensitive to these short-orbit binaries; c) The longer-period binaries might
not be fully characterized yet, they could have larger orbital periods and
more massive companions; two of them (I and N) could still have very light
companions.

From this, we can conclude firmly that the fraction of black widow systems in
$\omega$-Cen is unusual. The only comparable GC is M28, where they represent
half of the total binary population (Douglas et al., 2022). However, M28 has
vastly different properties, in particular a much denser core. In that GC,
almost all long-period binaries have been either disrupted, or show
significant eccentricities or even possible signs of having undergone
secondary exchange encounters. Only the BWs survive because their very short
orbital periods make them much more difficult to perturb: they are smaller
targets and require a more energetic encounter for the orbit to change
significantly.

###  4.2 Why this pulsar population is surprising

The total stellar encounter rate ($\Gamma$) and the stellar encounter rate per
binary ($\gamma$) are functions of the core radius ($r_{c}$) and central
density ($\rho_{c}$) of the GCs, with $\Gamma\propto\rho_{c}^{1.5}r_{c}^{2}$
and $\gamma\propto\rho_{c}^{0.5}r_{c}^{-1}$ (Verbunt & Hut, 1987; Verbunt &
Freire, 2014). The first parameter, $\Gamma$, gives an approximate prediction
of the size of the pulsar population, the second gives an approximate
prediction of how disturbed the binaries are: in GCs with low $\gamma$, the
population resembles the MSP population in the Galactic disk. As $\gamma$
increases, more frequent encounters have a higher chance of perturbing the
orbits of the binaries and on disrupting systems, producing more eccentric
binaries and more isolated pulsars overall. Finally, for core-collapsed
clusters with 5 or more pulsars known, such as Terzan 1, NGC 6517, NGC 6522,
NGC 6624, NGC 6752 and M15, the population is completely dominated by isolated
pulsars, while a large percentage of binaries results from secondary exchange
encounters (for recent discussions, see Ridolfi et al. 2021, 2022).

$\omega$-Cen has a low stellar density compared to most other GCs
($10^{3.15}\rm L_{\odot}pc^{-3}$), which results in rather low values for
$\Gamma$ and $\gamma$ of respectively 4.3 and 0.11 (these are normalized to
the values of the GC M4, as in Verbunt & Freire 2014, as in that work we have
used the values of $r_{c}$ and $\rho_{c}$ from Harris 2010). As a comparison,
47 Tuc has $\Gamma=29.2$ and $\gamma=6.6$.

These $\Gamma$ values mean that LMXBs and MSP binaries should form in
$\omega$-Cen at a rate $\sim 7$ times smaller than that of 47 Tuc. However,
the 18 pulsars in $\omega$-Cen are two thirds of the 27 pulsars detected in 47
Tuc with the same telescope (MeerKAT) and receivers (L-band, see Ridolfi et
al. 2021). This is partly a result of the fact that the latter searches were
done only with a single beam using 44 antennas within the 1-km core, not all
64 antennas and many hundreds of beams as the search described here; however
this difference is in part compensated by the larger distance of $\omega$-Cen
(5.2 kpc) compared to 47 Tuc (4.5 kpc). In any case, $\omega$-Cen is
surprisingly effective in producing MSPs.

The low $\gamma$ means that, once formed, there is not much chance of a
significant disruption of these systems: indeed, the interval between
successive interactions with other stars should be $\sim 60$ times larger in
$\omega$-Cen than in 47 Tuc. In the latter cluster, we see a MSP population
that resembles the pulsar population of the Galactic disk (except for the
absence of long-period binaries, see Ridolfi et al. 2016; Freire et al. 2017).
Therefore, the pulsar population of $\omega$-Cen should also be similar to the
MSP population in the Galactic disk, which is dominated by binaries, in an
approximate rate of 4 to 1. As discussed above, the fraction of isolated
pulsars in $\omega$-Cen appears to be large, but this might be mostly due to
selection effects, so it is still possible that the pulsar population of
$\omega$-Cen is dominated by binaries.

More difficult to explain is the small population of MSP-WD systems and long-
period binaries. As discussed above, the predominance of very tight systems is
likely real, and it is not something that can be expected from the dynamical
characteristics of this GC.

It is therefore clear that the properties of the pulsar populations of
$\omega$-Cen cannot be fully explained by two simple parameters like $\Gamma$
and $\gamma$: despite selection effects, we can already conclude that this
cluster has a larger pulsar population than expected from its $\Gamma$, too
many black widow systems and possibly too many isolated pulsars. That $\Gamma$
and $\gamma$ do not provide a full description of pulsar populations in GCs
was already highlighted by Verbunt & Freire (2014) when they compared the
pulsar populations of NGC 6440/1 with that of Terzan 5, which has a similar
$\Gamma$ and $\gamma$: the spin period distribution of the pulsars in both
populations is completely inconsistent. Therefore, the characteristics of the
pulsar populations in GCs must depend on additional factors, like the past
evolutionary history of the GCs and potentially their metallicity.

In the specific case of $\omega$-Cen, the past history of the system might be
of paramount importance. For instance, $\omega$-Cen could have been the
nucleus of a dwarf galaxy (Hilker & Richtler, 2000; Bekki & Freeman, 2003),
alternatively, it might have formed from the merger of several different GCs
(Calamida et al., 2020). Such explanations are motivated by the fact that this
GC has multiple stellar populations, with different metallicities and ages.
Given the typical ages of MSPs (many Gyr), the dramatic events in the history
of these GCs should be of paramount importance for an explanation of the
characteristics of its pulsar population today.

We note in this regard that there are other GCs in the Galaxy that are known
to have multiple stellar populations, and are likely associated with dwarf
galaxy systems or are the results of GC mergers. Two of the most prominent are
Terzan 5, where three distinct stellar populations have been found (Ferraro et
al., 2009; Origlia et al., 2013), and NGC 1851 (Carretta et al., 2011). All
have very abundant pulsar populations with about 50% and 40% of isolated
pulsars respectively (e.g., Ransom et al. 2005; Ridolfi et al. 2022). However,
their cores are so dense that the orbital characteristics of these pulsar
populations have likely been significantly altered by exchange encounters.
This is not the case for $\omega$-Cen, where the low $\Gamma$ and $\gamma$
mean that the orbital characteristics of the pulsar population have been
preserved for a long time; they should therefore reflect the earlier
evolutionary history of the cluster.

A detailed evaluation of these possibilities is beyond the scope of this work,
but it will be a profitable exercise, especially after the pulsar population
in $\omega$-Cen is better characterized.

##  5 Summary and future prospects

In this paper, we presented the discovery of 13 new pulsars in $\omega$-Cen,
which more than tripled the population of known pulsars in this cluster. They
are found within the core and also between the core and half light radius of
the cluster. Among them, six are isolated pulsars and the other seven are
binaries. More than half of the binaries have orbits less than 4 hours, which
is the length of the observations; three other binaries have orbital periods
of about 1 day, but confirming this will require additional observations. All
but one of the binaries have very light companions and two of them have
apparent eclipses.

Follow-up observations are crucial to improve the orbital parameters of the
wide binaries (I, N and Q) and help estimate their companion masses, which is
a first step to an accurate characterization of those systems. Additional
observations will also be important for deriving phase-connected timing
solutions. Thanks to the many beams systhesised in each observation, we were
able to constrain the positions of several of the new pulsars and compare them
with the position of X-ray unassociated sources; for some binaries, there is
an X-ray source nearby, within a few arcseconds. With timing solutions, these
positions will become orders of magnitude more precise, this will either
confirm or rule out some of our preliminary associations. Multi-wavelength
observations should be carried out to check if these sources are still
emitting X-rays, in order to establish whether the emission comes from the
pulsars themselves or from ongoing accretion. This kind of observations can
also be used to investigate other X-ray sources that have no radio signals
detected because it is possible that there is a LMXB system there and the
pulsar is still accreting. The timing solutions, with precise estimates of
acceleration and proper motions, will also be important for characterizing the
gravitational field of $\omega$-Cen (Prager et al., 2017; Freire et al., 2017;
Abbate et al., 2018).

The large pulsar population, the large number of isolated pulsars and the
fraction of black widow systems are surprising, considering the small
encounter rate and the low encounter rate per binary of this GC. These
parameters, although useful for an approximate characterization of the pulsar
population of GCs, clearly do not tell the full story; it is very likely, for
instance, that the past dynamical history of the GCs and the stellar evolution
in binaries play important roles. This implies that the accurate
characterization of the pulsar populations in GCs in general, and $\omega$-Cen
in particular, provides valuable material for the study of stellar and cluster
evolution. A particularly interesting possibility is that the pulsar
population in $\omega$-Cen came from different smaller clusters that might
have merged to form it (Calamida et al., 2020).

There are still more than half of the total beams outside the half light
radius (5′) that have not been searched. Searching them will be important for
finding out how centrally condensed the pulsar population of $\omega$-Cen is
compared to other clusters. The Parkes survey by Dai et al. (2020), which at
L-band has a beam radius of 7.5′, has only found pulsars within, or very near
the core, as discovered by our recent MeerKAT localisations and their
subsequent pulsar timing. Our discoveries are also mostly within the core,
with only five pulsars between 1 and 2 core radii from the centre of the
cluster. The number of X-ray sources in $\omega$-Cen decreases significantly
beyond 2 core radii, but still presents a detectable excess compared to the
background beyond 3 core radii (Henleywillis et al., 2018). This suggests that
additional pulsars might be detectable outside the half-light radius, but
likely in significantly smaller numbers.

The ‘dynamical relaxation time” in the core of $\omega$-Cen is 4 Gyr, while
the median relaxation time for the cluster as a whole is 12 Gyr. For 47 Tuc,
these numbers are 0.07 and 3.5 Gyr respectively, for Terzan 5, they are 0.037
and 0.34 Gyr respectively (Harris, 2010). What this means is that, in 47 Tuc
and Terzan 5, enough time has elapsed for mass segregation to occur, all
pulsars in these two clusters (with the exception of 47 Tuc X, Ridolfi et al.
2016) have moved to within 2′from their centres, and are likely in dynamical
equilibrium with the remaining stars of the cluster (i.e., they are a
“relaxed" population, see e.g., Heinke et al. 2005). In $\omega$-Cen, this
process takes much longer. This means that the current pulsar distribution,
especially outside the core, likely reflects the “original" dynamics of the
pulsars within the cluster (either where they formed, or where they were
placed by previous interactions of the cluster). A detailed dynamical study of
this distribution could thus provide additional clues on the origin of this
unusual pulsar population.

We also note that future TRAPUM observations with UHF-band (550-1100 MHz) and
S-band (1750-3500 MHz) receivers will very likely further increase the
population of known pulsars in $\omega$-Cen in all regions by probing
different spectral windows.

## Acknowledgements

TRAPUM observations used the FBFUSE and APSUSE computing clusters for data
acquisition, storage and analysis. These clusters were funded and installed by
the Max-Planck-Institut für Radioastronomie and the Max-Planck-Gesellschaft.
WC, AR and FA acknowledge continuing valuable support from the Max-Planck
Society. LV acknowledges financial support from the Dean’s Doctoral Scholar
Award from the University of Manchester. APo, AR and MBu gratefully
acknowledge financial support by the research grant “iPeska” (P.I. Andrea
Possenti) funded under the INAF national call Prin-SKA/CTA approved with the
Presidential Decree 70/2016. APo, AR, MBu also acknowledge support from the
Ministero degli Affari Esteri e della Cooperazione Internazionale - Direzione
Generale per la Promozione del Sistema Paese - Progetto di Grande Rilevanza
ZA18GR02. 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. The National Radio Astronomy
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. RPB acknowledges support ERC Starter Grant ‘Spiders’ under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement number 715051).

## Data availability

The data underlying this article will be shared on reasonable request to the
TRAPUM collaboration.

## References

  * Abbate et al. (2018) Abbate F., Possenti A., Ridolfi A., Freire P. C. C., Camilo F., Manchester R. N., D’Amico N., 2018, MNRAS, 481, 627
  * Abbate et al. (2022) Abbate F., et al., 2022, MNRAS, 513, 2292
  * Abdo et al. (2009) Abdo A. A., et al., 2009, Science, 325, 845
  * Abdo et al. (2010) Abdo A. A., et al., 2010, ApJS, 187, 460
  * Barr (2018) Barr E. D., 2018, in Weltevrede P., Perera B. B. P., Preston L. L., Sanidas S., eds, IAU symposium and colloquium proceedings series Vol. 337, Pulsar Astrophysics the Next Fifty Years. pp 175–178, doi:10.1017/S1743921317009036
  * Bedin et al. (2004) Bedin L. R., Piotto G., Anderson J., Cassisi S., King I. R., Momany Y., Carraro G., 2004, ApJ, 605, L125
  * Bekki & Freeman (2003) Bekki K., Freeman K. C., 2003, MNRAS, 346, L11
  * Benacquista & Downing (2013) Benacquista M. J., Downing J. M. B., 2013, Living Reviews in Relativity, 16, 4
  * Bhattacharya & van den Heuvel (1991) Bhattacharya D., van den Heuvel E. P. J., 1991, Phys. Rep., 203, 1
  * Calamida et al. (2020) Calamida A., et al., 2020, ApJ, 891, 167
  * Camilo (2018) Camilo F., 2018, Nature Astronomy, 2, 594
  * Camilo et al. (2015) Camilo F., et al., 2015, ApJ, 810, 85
  * Carretta et al. (2011) Carretta E., Lucatello S., Gratton R. G., Bragaglia A., D’Orazi V., 2011, A&A, 533, A69
  * Chen et al. (2021) Chen W., Barr E., Karuppusamy R., Kramer M., Stappers B., 2021, Journal of Astronomical Instrumentation, 10, 2150013
  * Clark (1975) Clark G. W., 1975, ApJ, 199, L143
  * Dai et al. (2020) Dai S., Johnston S., Kerr M., Camilo F., Cameron A., Toomey L., Kumamoto H., 2020, ApJ, 888, L18
  * DeCesar et al. (2015) DeCesar M. E., Ransom S. M., Kaplan D. L., Ray P. S., Geller A. M., 2015, ApJ, 807, L23
  * Dewey et al. (1985) Dewey R. J., Taylor J. H., Weisberg J. M., Stokes G. H., 1985, ApJ, 294, L25
  * Douglas et al. (2022) Douglas A., et al., 2022, ApJ, 927, 126
  * Edwards et al. (2001) Edwards R. T., van Straten W., Bailes M., 2001, ApJ, 560, 365
  * Ferraro et al. (2009) Ferraro F. R., et al., 2009, Nature, 462, 483
  * Forbes & Bridges (2010) Forbes D. A., Bridges T., 2010, MNRAS, 404, 1203
  * Freire et al. (2004) Freire P. C., Gupta Y., Ransom S. M., Ishwara-Chandra C. H., 2004, ApJ, 606, L53
  * Freire et al. (2017) Freire P. C. C., et al., 2017, MNRAS, 471, 857
  * Harris (2010) Harris W. E., 2010, arXiv e-prints, p. arXiv:1012.3224
  * Heinke et al. (2005) Heinke C. O., Grindlay J. E., Edmonds P. D., Cohn H. N., Lugger P. M., Camilo F., Bogdanov S., Freire P. C., 2005, ApJ, 625, 796
  * Henleywillis et al. (2018) Henleywillis S., Cool A. M., Haggard D., Heinke C., Callanan P., Zhao Y., 2018, MNRAS, 479, 2834
  * Hilker & Richtler (2000) Hilker M., Richtler T., 2000, A&A, 362, 895
  * Hobbs et al. (2020) Hobbs G., et al., 2020, Publ. Astron. Soc. Australia, 37, e012
  * Ibata et al. (1994) Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature, 370, 194
  * Ibata et al. (2019) Ibata R. A., Bellazzini M., Malhan K., Martin N., Bianchini P., 2019, Nature Astronomy, 3, 667
  * Jonas & MeerKAT Team (2016) Jonas J., MeerKAT Team 2016, in MeerKAT Science: On the Pathway to the SKA. p. 1
  * King et al. (2003) King A. R., Davies M. B., Beer M. E., 2003, MNRAS, 345, 678
  * Kulkarni et al. (1991) Kulkarni S. R., Anderson S. B., Prince T. A., Wolszczan A., 1991, Nature, 349, 47
  * Lee et al. (2009) Lee J.-W., Kang Y.-W., Lee J., Lee Y.-W., 2009, Nature, 462, 480
  * Liu et al. (2014) Liu K., Eatough R. P., Wex N., Kramer M., 2014, MNRAS, 445, 3115
  * Lorimer & Kramer (2004) Lorimer D. R., Kramer M., 2004, Handbook of Pulsar Astronomy. Vol. 4
  * Lynch et al. (2012) Lynch R. S., Freire P. C. C., Ransom S. M., Jacoby B. A., 2012, ApJ, 745, 109
  * Morello et al. (2020) Morello V., Barr E. D., Stappers B. W., Keane E. F., Lyne A. G., 2020, MNRAS, 497, 4654
  * Origlia et al. (2013) Origlia L., Massari D., Rich R. M., Mucciarelli A., Ferraro F. R., Dalessandro E., Lanzoni B., 2013, ApJ, 779, L5
  * Possenti et al. (2005) Possenti A., D’Amico N., Corongiu A., Manchester D., Sarkissian J., Camilo F., Lyne A., 2005, in Rasio F. A., Stairs I. H., eds, Astronomical Society of the Pacific Conference Series Vol. 328, Binary Radio Pulsars. p. 189
  * Prager et al. (2017) Prager B. J., Ransom S. M., Freire P. C. C., Hessels J. W. T., Stairs I. H., Arras P., Cadelano M., 2017, ApJ, 845, 148
  * Ransom (2011) Ransom S., 2011, PRESTO: PulsaR Exploration and Search TOolkit, Astrophysics Source Code Library, record ascl:1107.017 (ascl:1107.017)
  * Ransom et al. (2002) Ransom S. M., Eikenberry S. S., Middleditch J., 2002, AJ, 124, 1788
  * Ransom et al. (2003) Ransom S. M., Cordes J. M., Eikenberry S. S., 2003, ApJ, 589, 911
  * Ransom et al. (2005) Ransom S. M., Hessels J. W. T., Stairs I. H., Freire P. C. C., Camilo F., Kaspi V. M., Kaplan D. L., 2005, Science, 307, 892
  * Ridolfi et al. (2016) Ridolfi A., et al., 2016, MNRAS, 462, 2918
  * Ridolfi et al. (2021) Ridolfi A., et al., 2021, MNRAS,
  * Ridolfi et al. (2022) Ridolfi A., et al., 2022, A&A, 664, A27
  * Roberts (2013) Roberts M. S. E., 2013, in van Leeuwen J., ed., IAU Symposium Vol. 291, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years. pp 127–132 (arXiv:1210.6903), doi:10.1017/S174392131202337X
  * Stappers & Kramer (2016) Stappers B. W., Kramer M., 2016, in MeerKAT Science: On the Pathway to the SKA. p. 9
  * Venter et al. (2009) Venter C., De Jager O. C., Clapson A. C., 2009, ApJ, 696, L52
  * Verbunt & Freire (2014) Verbunt F., Freire P. C. C., 2014, A&A, 561, A11
  * Verbunt & Hut (1987) Verbunt F., Hut P., 1987, in Helfand D. J., Huang J. H., eds, Vol. 125, The Origin and Evolution of Neutron Stars. p. 187
  * Vleeschower et al. (2022) Vleeschower L., et al., 2022, MNRAS, 513, 1386


# Optical soliton formation controlled by angle-twisting in photonic moiré
lattices

Qidong Fu,1 Peng Wang,1 Changming Huang,2 Yaroslav V. Kartashov,3,4 Lluis
Torner,3,5 Vladimir V. Konotop,6,7 Fangwei Ye1∗
_1 School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China_
_2 Department of Electronic Information and Physics, Changzhi University,
Shanxi 046011, China_
_3 ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science
and Technology, 08860 Castelldefels (Barcelona), Spain_
_4 Institute of Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow,
108840, Russia_
_5 Universitat Politecnica de Catalunya, 08034 Barcelona, Spain_
_6 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa,
Campo Grande, Ed. C8, Lisboa 1749-016, Portugal_
_7 Centro de Física Teórica e Computacional, Universidade de Lisboa, Campo
Grande, Ed. C8, Lisboa 1749-016, Portugal_
_∗ Corresponding author<EMAIL_ADDRESS>
Exploration of the impact of synthetic material landscapes featuring tunable
geometrical properties on physical processes is a research direction of the
highest current interest because of the outstanding phenomena that are
continuously uncovered. Twistronics and the properties of wave excitations in
moiré lattices are salient examples. Moiré patterns bridge the gap between
aperiodic structures and perfect crystals, thus opening the door to the
exploration of effects accompanying the transition from commensurate to
incommensurate phases. Moiré patterns have revealed profound effects in
graphene-based systems GrapheneMoire ; Woods14 ; BandFlatGraphene ;
UnconvenSupercond ; Science , they are used to manipulate ultracold atoms
atomic1 ; atomic2  and to create gauge potentials Gauge , and are observed in
colloidal clusters colloids . Recently, it was shown that photonic moiré
lattices enable the observation of the two-dimensional localization-to-
delocalization transition of light in purely linear systems LDTscirep ;
LDTnature . Here we employ moiré lattices optically-induced in photorefractive
nonlinear media Efremidis02 ; Fleischer03 ; Freedman06  to elucidate the
formation of optical solitons under different geometrical conditions
controlled by the twisting angle between the constitutive sublattices. We
observe the formation of solitons in lattices that smoothly transit from fully
periodic geometries to aperiodic ones, with threshold properties that are a
pristine direct manifestation of flat-band physics LDTnature .

By and large, the linear transport and localization properties of excitations
in a material are intimately determined by its inner symmetry and geometrical
properties, including a periodic or aperiodic nature, as it occurs in
electronic electron01 , atomic matter01 ; matter02 , optical optical01 ;
optical02 , or two-dimensional material graphene01  systems. This feature is
directly related to the nature of the eigenstates of the system, which can be
extended or localized. When an underlying material exhibits nonlinear
response, formation of self-sustained excitation, alias solitons, becomes
possible solitons01 ; solitons02 ; solitons03 . Nevertheless, the properties
of solitons are still strongly impacted by the linear spectrum of the system.
Optical media offer a unique laboratory for the investigation of solitons in
different environments. Thus, two-dimensional self-trapping and soliton
formation have been investigated in fully periodic optical lattices
Fleischer03 ; Yang03 ; Efremidis03 ; Neshev03 , as well as in quasicrystals,
which are characterized by broken translational invariance Freedman06 ;
Ablowitz06 ; Law10 ; Ablowitz12 ; Denz10 . However, formation of solitons at
the transition from aperiodic to periodic systems has been never explored
experimentally because of lacking of a suitable setting. Optical moiré
lattices LDTnature  offer a powerful platform enabling such study. Created
with incommensurate geometries, moiré patterns may allow localization of light
even in the linear limit due to the existence of a large number of extremely
flat bands in their spectra, a property that therefore strongly impacts the
diffraction of beams in such media. Thus, since soliton states can emerge due
to a balance between diffraction and self-phase-modulation induced by
nonlinearity, moiré lattices allow the investigation of the formation of
solitons controlled by the geometry of the induced optical potential. Here we
provide the experimental evidence of such possibility, by reporting
qualitative differences in soliton excitation dynamics in commensurate and
incommensurate moiré lattices.

A photonic moiré lattice is created in a photorefractive crystal by a shallow
modulation of the refractive index in the $(x,y)$ plane induced by two
mutually rotated, or twisted, periodic square sublattices generated by light
with ordinary polarization. The latter is chosen so that the corresponding
beams practically do not experience self-action in the crystal, thus
propagating undistorted. In contrast, light with extraordinary polarization
experiences a strong nonlinear response. Propagation along the $z$-direction
of a signal beam in such polarization is governed by the nonlinear Schrödinger
equation for the dimensionless amplitude $\Psi({\bm{r}}_{\perp},z)$
Efremidis02 ; Fleischer03 :

$\centering i\frac{\partial\Psi}{\partial
z}=-\frac{1}{2}\nabla_{\perp}^{2}\Psi+\frac{V_{0}}{1+I(\bm{r}_{\perp})+|\Psi|^{2}}\Psi.\@add@centering$
(1)

Here $\nabla_{\perp}=(\partial/\partial x,\partial/\partial y)$; ${\bf
r_{\perp}}=(x,y)$ is the radius-vector in the transverse plane; $z$ is the
longitudinal coordinate scaled to the characteristic length $2\pi
n_{e}\lambda$, where $\lambda$ is the wavelength (in our experiments
$\lambda=632.8$ nm); $n_{e}$ is the unperturbed refractive index of the
crystal experienced by the extraordinary-polarized light; and $V_{0}>0$ is the
dimensionless applied dc field. Here we set $V_{0}=5$, which corresponds to
$5.7\times 10^{4}$ V/m dc electric field applied to the crystal;
$I({\bm{r}}_{\perp})=\left|p_{1}V({\bm{r}}_{\perp})+p_{2}V(S{\bm{r}}_{\perp})\right|^{2}$
is the moiré pattern composed by two ordinary-polarized periodic sublattices
$V({\bm{r}}_{\perp})$ and $V(S{\bm{r}}_{\perp})$ interfering in the $(x,y)$
plane [here $S=S(\theta)$ is the matrix of rotation in the $(x,y)$ plane by
the angle $\theta$]; $p_{1}$ and $p_{2}$ are the amplitudes of the first and
second square sublattices, respectively. Each square sublattice
$V({\bm{r}}_{\perp})$ is formed by the interference of four plane waves
LDTnature . In the following we set the amplitude of the first sublattice to
$p_{1}=0.5$, which corresponds to an average intensity $I_{\text{av}}\approx
1.9~{}\text{mW}/\text{cm}^{2}$, and tune the amplitude $p_{2}$ of the second
sublattice. For such parameters, the actual refractive index modulation depth
in the moiré pattern illustrated in Fig. 1(a-c) is of the order of $\delta
n\sim 10^{-4}$.

Figure 1:  Moiré patterns and properties of their linear eigenmodes. Example
of periodic (a,c) and aperiodic (b) moiré lattices $I({\bm{r}}_{\perp})$
produced by two superimposed square sublattices with $p_{1}=0.5$ and
$p_{2}=0.3$ rotated by the angles $\theta=\textrm{arctan}(5/12)$ (a),
$\theta\approx 0.167\pi$ (b), $\theta=\textrm{arctan}(3/4)$ (c). Blue arrows
in (a),(c) indicate the primitive lattice vectors. (d) Form-factor (inverse
width) of the most localized linear eigenmode of the lattice versus rotation
angle $\theta$ and depth $p_{2}$ of the second sublattice at $p_{1}=0.5$.
Green arrows indicate angles corresponding to the (3,4,5) and (5,12,13)
Pythagorean triples. The top lattice bands for two Pythagorean angles
$\theta=\arctan(3/4)$ and $\theta=\arctan(5/12)$ superimposed in one plot for
$p_{2}=0.1$ (e) and $p_{2}=0.3$ (f). Bloch momenta $k_{x,y}$ are normalized to
the width of the Brillouin zone equal to $K\approx 1.265$ for
$\theta=\arctan(3/4)$ and to $K\approx 0.785$ for $\theta=\arctan(5/12)$.

Two-dimensional Pythagorean moiré lattices composed of two square Bravais
sublattices $p_{1}V({\bm{r}}_{\perp})$ and $p_{2}V(S{\bm{r}}_{\perp})$ (the
point group $D_{4}$) rotated with respect to each other around common lattice
site, are periodic (commensurate) structures only when the rotation angle
$\theta$ satisfies $\cos\theta=a/c$, $\sin\theta=b/c$, where the positive
integers $(a,b,c)$ having no common divisors except $1$, constitute a
primitive Pythagorean triple, i.e., $a^{2}+b^{2}=c^{2}$ (Fig. 1a,c). We call
such angles Pythagorean. For other rotation angles the pattern is aperiodic
(incommensurate, or almost periodic in mathematical terms), as illustrated in
Fig. 1b. The linear spectrum of the lattices LDTscirep ; LDTnature  can be
obtained by omitting the nonlinear term $|\Psi|^{2}$ in Eq. (1) and searching
for the corresponding linear eigenmodes in the form
$\Psi({\bm{r}}_{\perp},z)=\psi({\bm{r}}_{\perp})e^{ibz}$, where $b$ is the
linear propagation constant and $\psi({\bm{r}}_{\perp})$ is the transverse
field distribution. To characterize the mode localization we use the integral
form-factor $\chi=(\iint|\Psi|^{4}d^{2}{\bm{r}}_{\perp})^{1/2}/U$, with $U$
being the mode power $U=\iint|\Psi|^{2}d^{2}{\bm{r}}_{\perp}$. The larger
$\chi$ the stronger the localization.

The dependence of the form-factor of the mode with the largest $b$ (this is
the most localized mode) on $\theta$ and $p_{2}$ is shown in Fig. 1d. For
Pythagorean angles the mode is delocalized for any depth $p_{2}$ of the second
sublattice because in this case the moiré pattern is periodic, but for non-
Pythagorean angles the mode becomes localized if $p_{2}$ exceeds some critical
value, $p_{2}^{\rm cr}\approx 0.18$ corresponding to the linear localization-
delocalization threshold. The physical origin of this phenomenon is the
suppressed diffraction due to flatness of the allowed bands of the effective
Pythagorean lattice approximating real incommensurate moiré pattern at
$p_{2}>p_{2}^{\rm cr}$ LDTnature . This observation comes from a general rule:
the higher the order of the primitive Pythagorean triple (determined by the
integer $c$) the larger the area of the respective primitive cell of the
lattice (see the blue arrows in Fig. 1a,c) and the smaller the width (in $b$)
of the allowed bands, thus indicating a reduced diffraction strength. This
behavior is illustrated in Fig. 1e,f that compares top bands of the Floquet-
Bloch spectra calculated for two different Pythagorean angles and two
different sublattice amplitudes $p_{2}$. As discussed below, the angular-
dependent band flattening has direct implications for soliton formation in the
lattices, because it exposes that the diffraction strength experienced by
narrow linear inputs in the lattices notably decreases with increase of the
order of the primitive Pythagorean triple. This is particularly clear at
$p_{2}>p_{2}^{\rm cr}$, where patterns akin to discrete diffraction are
observed in the linear limit and is less pronounced at $p_{2}<p_{2}^{\rm cr}$,
when such patterns co-exist with a rapidly expanding broader background.

Figure 2:  Families of two-dimensional solitons in moiré lattices. Soliton
power $U$ and peak amplitude $\psi_{\text{max}}$ versus propagation constant
$b$ (left column a, b, c) and representative soliton profiles (right column)
are shown for incommensurate moiré lattice with the rotation angle
$\theta\approx 0.139\pi$ below ($p_{2}=0.1$, the upper panels) and above
($p_{2}=0.3$, the middle panels) the critical value $p_{2}^{\rm cr}$, as well
as for a lattice with the Pythagorean rotation angle $\theta=\arctan(3/4)$
with $p_{2}=0.3$ (bottom panels). The profiles shown in the right column
correspond to red dots on the $\psi_{\text{max}}(b)$ curves.

Turning now to the nonlinear regime, we look for soliton solutions of Eq. (1)
in the form $\Psi({\bm{r}}_{\perp},z)=\psi({\bm{r}}_{\perp})e^{ibz}$, where
$b$ is the nonlinear propagation constant, which for $V_{0}>0$ (focusing
nonlinearity) exceeds the propagation constants of the linear eigenmodes.
Solitons form families characterized by the dependencies of the peak amplitude
$\psi_{\text{max}}=\text{max}{|\psi|}$ and power $U$ on $b$, as shown in Fig.
2. Since incommensurate moiré lattices have either delocalized (at
$p_{2}<p_{2}^{\rm cr}$) or localized (at $p_{2}>p_{2}^{\rm cr}$) linear modes,
solitons in such lattices show a completely different behaviour in the low-
amplitude limit depending on the $p_{2}$ value. When linear localized modes do
not exist ($p_{2}<p_{2}^{\rm cr}$), then consistent with the behavior of
homogeneous media Townes , solitons in incommensurate moiré lattices can exist
only if they carry a power $U$ that exceeds a certain threshold value
$U_{\text{th}}$ (Fig. 2a) below which they quickly diffract. However, when
$p_{2}>p_{2}^{\rm cr}$ and linear localized states exist, the soliton family
bifurcates from the respective linear localized mode, remaining well-localized
at any power $U$. Under such conditions, solitons do not feature a power
existence threshold (Fig. 2b). In contrast, in commensurate moiré lattices
corresponding to the Pythagorean angles, solitons always exhibit a nonzero
existence power threshold for any depth of the second sublattice $p_{2}$ (Fig.
2c), and strongly expand at low amplitudes. In all the cases shown in Fig. 2,
the peak soliton amplitude vanishes in the linear limit (see red curves in
Fig. 2a-c), while far from it the solitons become strongly localized. All
shown soliton families are either completely stable, when power threshold is
zero, or unstable at low amplitudes and get stabilized at high amplitudes,
when the power threshold does not vanish.

Figure 3:  Thresholds for soliton formation in moiré lattices. Comparison of
theoretically calculated (a,c) and experimentally measured (b,d) dependencies
of threshold for soliton formation on rotation angle $\theta$ for the depth of
the second sublattice above localization-delocalization threshold $p_{2}=0.3$
(a,b) and below localization-delocalization threshold $p_{2}=0.1$ (c,d).

Here we report the experimental observation of all such scenarios using
lattices with readily tunable twisting angle. More specifically, the
transition between commensurate and incommensurate moiré lattices or between
the regimes, where delocalization and localization takes place in the linear
limit, can be explored by adjusting the phase mask used for the creation of
the sublattices or by tuning their relative amplitudes. The comparison between
the theoretical predictions and the experimental observations of the threshold
power $U_{\text{th}}$ as a function of the rotation angle $\theta$ between the
sublattices is shown in Fig. 3. For $p_{2}$ above the critical value (Fig.
3a), $U_{\text{th}}$ vanishes for non-Pythagorean angles while it exhibits
narrow peaks around the Pythagorean angles. As visible in Fig. 1e,f, the band
curvature decreases with the increase of the order of the Pythagorean triple.
Thus, the nonlinearity required to balance the curvature-induced diffraction
decreases too, in agreement with the observations depicted in Fig. 3a,b. The
first, second, and third highest peaks occur near the Pythagorean angles
corresponding to the primitive triples $(3,4,5)$, $(5,12,13)$, and
$(8,15,17)$, respectively, in agreement with the properties of the bands at
such angles. Also, the heights of the peaks are observed to diminish as areas
of the primitive cells of the respective moiré patterns increases. Notice that
the area of a Pythagorean cell, that is defined by $c$, non-monotonically
depends on the twist angle. The smallest area of the primitive (square) cell
of a Pythagorean moiré pattern corresponds to the primitive triple $(3,4,5)$,
that explains the locations of the absolute maxima in all panels of Fig. 3.
Because the bands become flatter when the areas of the primitive cells
increase, other $U_{\textrm{th}}$ maxima are almost undetectable. This is
consistent with the approximation of incommensurate Pythagorean moiré patterns
by commensurate ones established in LDTnature . In all experimental results
reported, as well as in all numerical simulations showing dynamical soliton
excitation, we observed stability of two-dimensional beams.

Remarkably, the above relation between the order of the Pythagorean triple
associated with the lattice and the threshold for soliton formation, even
though less pronounced, was still observed in the regime $p_{2}<p_{2}^{\rm
cr}$, where localization in the linear system is impossible even in
incommensurate lattices (Fig. 3c,d). Under these conditions the soliton
formation threshold $U_{\textrm{th}}$ does not vanish for any rotation angles,
although it remains sensitive to the diffraction properties of the lattice.
Our observations show that it achieves the maximal value around
$\theta=\arctan(3/4)$ angle that corresponds to the commensurate Pythagorean
lattice associated with the $(3,4,5)$ triple and having the smallest possible
primitive cell. In all cases we observed an excellent agreement between the
theoretical predictions (Fig. 3a,c) and the experimental (Fig. 3b,d) results.

The threshold power $U_{\text{th}}$ for soliton formation was measured by
evaluating the soliton content $C=U_{\text{out}}/U_{\text{in}}$, where
$U_{\text{in}}$ is the power of the input beam measured at $z=0$ and
$U_{\text{out}}$ is the power of the beam that remains at the output after $2$
cm of propagation, within a pinhole of radius of $27~{}\mu$m approximately
equal to one period of the sublattice forming the moiré pattern. In the
numerical modeling of the dynamical soliton excitation within the frames of
the model (1), we used the same criterion calculating the output power for
Gaussian inputs $\Psi=A\,\textrm{exp}(-{\bm{r}}_{\perp}^{2}/r_{0}^{2})$ within
a circle of radius equal to the sublattice period, albeit at large propagation
distances ($z=1000$) to avoid the presence of transient effects.

In Fig. 4 we show the dependence of the soliton content $C$ on the input power
$U$, together with the low-power and high-power output beams, for a moiré
lattice with $p_{2}>p_{2}^{\rm cr}$ (i.e., above the linear localization-
delocalization threshold), for different rotation angles. In lattices
corresponding to the Pythagorean angles $\arctan(3/4)$ and $\arctan(5/12)$
(black and red lines in Fig. 4a,b) one observes a sharp jump from diffraction
to a soliton content close to unity upon increase of the input power. This
means that practically all input power goes into a soliton. The power
corresponding to the abrupt increase in $C$ is the dynamical threshold for
soliton formation $U_{\text{th}}$ \- it is this quantity that is plotted in
Fig. 3b and d. For non-Pythagorean angles, i.e. for incommensurate moiré
lattices, soliton content $C$ remains very high irrespective of the input
power (green lines in Fig. 4a,b) because such lattice supports linear modes
that are effectively excited by the input beam. Figures 4c-e illustrate that
for Pythagorean angles the output beam becomes localized only at sufficiently
high power, while for non-Pythagorean angles the output is always localized,
consistent with theoretical predictions and with expectations based on the
flatness of the bands at such twisting angles.

Figure 4:  Soliton formation above the linear localization-delocalization
threshold. Experimentally measured (a) and numerical (b) soliton content
versus input power for Pythagorean rotation angles $\theta=\arctan(3/4)$
(curve 1), $\arctan(5/12)$ (curve 2), and for non-Pythagorean one $0.139\pi$
(curve 3). The second sublattice depth is $p_{2}$=0.3. Lattice profiles and
corresponding low- and high-power output distributions for experimentally
measured signal beam for (c) $\theta=\arctan(5/12)$, (e)
$\theta=\arctan(3/4)$, and (d) $\theta\approx 0.139\pi$ angles. Power levels
are indicated on the plots.

Figure 5 depicts the observed behavior of the soliton content for lattices
with $p_{2}<p_{2}^{\rm cr}$ that thus cannot support linear localized states.
Now we observed step-like behavior of the soliton content $C$ for all
Pythagorean and non-Pythagorean angles. Namely, solitons are excited only
above a threshold power (see outputs in Fig. 5a-e). Notice the higher values
of the threshold power relative to the previous cases, consistent with Fig. 3c
(in all the cases the used input powers are well below the saturation limit).
Above threshold, the values of $C$ are similar for different rotation angles,
indicating that the excited states are strongly localized on the scale of
about one sublattice period.

Figure 5:  Soliton formation below linear localization-delocalization
threshold. Experimentally measured (a) and theoretically calculated (b)
soliton content versus input power for rotation angles $\theta$=arctan(5/12)
(curve 1), 0.139$\pi$ (curve 2), 0.167$\pi$ (curve 3), arctan(3/4) (curve 4),
0.222$\pi$ (curve 5), and 0.250$\pi$ (curve 6) for second sublattice depth
$p_{2}$=0.1. Lattice profiles and corresponding low- and high-power output
distributions for signal beam for Pythagorean (c) $\theta$=arctan(5/12), (e)
$\theta$=arctan(3/4), and non-Pythagorean (d) $\theta\approx 0.139\pi$ angles.

Summarizing, we have observed the excitation of two-dimensional solitons in
Pythagorean moiré patterns and showed that their excitation dynamics is
dictated by the fundamental transition from commensurate to incommensurate
lattice geometries. Incommensurability and the relation between the depths of
the two sublattices forming the moiré pattern has been shown to be crucial for
soliton excitation. In all cases, the behavior of the soliton formation
threshold is a direct manifestation of the band structure of the moiré
lattices resulting from the different twisting angles of the sublattices and,
in particular, of the general band-flattening associated to the geometry of
the moiré lattices. We anticipate that similar phenomena may occur in moiré
patterns composed of sublattices of other crystallographic symmetries, such as
twisted honeycomb lattices, and in other physical systems where flat-bands
induced by geometry arise.

## References

  * (1) Decker, R., Wang, Y., Brar, V. W., Regan, W., Tsai, H.-Z., Wu, Q., Gannett, W., Zettl, A. & Crommie, M. F. Local Electronic Properties of Graphene on a BN Substrate via Scanning Tunneling Microscopy, Nano Lett. 11, 2291-2295 (2011).
  * (2) Woods, C. R., Britnell, L., Eckmann, A., Ma, R. S., Lu, J. C., Guo, H. M., Lin, X., Yu, G. L., Cao, Y., Gorbachev, R. V., Kretinin, A. V., Park, J., Ponomarenko, L. A., Katsnelson, M. I., Gornostyrev, Yu. N., Watanabe, K., Taniguchi, T., Casiraghi, C., Gao, H-J., Geim, A. K. & Novoselov, K. S. Commensurate–incommensurate transition in graphene on hexagonal boron nitride, Nat. Phys. 10, 451 (2014).
  * (3) Cao, Y., Fatemi, V., Demir, A., Fang, S., Tomarken, S. L., Luo, J. Y., Sanchez-Yamagishi, J. D., Watanabe, K., Taniguchi, T., Kaxiras, E., Ashoori, R. C. & Jarillo-Herrero, P. Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018).
  * (4) Cao, Y., Fatemi, V., Fang, S., Watanabe, K., Taniguchi, T., Kaxiras, E. & Jarillo-Herrero, P. Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018).
  * (5) Ahn, S. J., Moon, P., Kim, T.-H., Kim, H.-W., Shin, H.-C., Kim, E. H., Cha, H. W., Kahng, S.-J., Kim, P., Koshino, M., Son, Y.-W., Yang, C.-W. & Ahn, J. R. Dirac electrons in a dodecagonal graphene quasicrystal, Science 786, 782-786 (2018).
  * (6) González-Tudela, A. and Cirac, J. I., Cold atoms in twisted-bilayer optical potentials, Phys. Rev. A 100 053604 (2019)
  * (7) Salamon, T., Celi, A., Chhajlany, R. W., Frrot, I., Lewenstein, M., Tarruell, L., and Rakshit, D., Simulating twistronics without a twist, (2019), arXiv:1912.12736 [cond-mat.quant-gas].
  * (8) San-Jose, P., González, J. & Guinea, F., Non-Abelian Gauge Potentials in Graphene Bilayers, Phys. Rev. Lett. 108, 216802 (2012).
  * (9) Cao, X., Panizon, E., Vanossi, A., Manini, N., & Bechinger, C., Orientational and directional locking of colloidal clusters driven across periodic surfaces, Nat. Phys. 15, 774 (2019).
  * (10) Huang, C., Ye, F., Chen X., Kartashov, Y. V., Konotop, V. V., & Torner, L., Localization-delocalization wavepacket transition in Pythagorean aperiodic potentials, Sci. Rep. 6, 32546 (2016).
  * (11) Wang, P., Zheng, Y., Chen, X., Huang, C., Kartashov, Y. V., Torner, L., Konotop, V. V., & Ye, F., Localization and delocalization of light in photonic moiré lattices, Nature 577, 42 (2020).
  * (12) Efremidis, N. K., Sears, S., Christodoulides, D. N., Fleischer, J. W. & Segev, M. Discrete solitons in photorefractive optically induced photonic lattices. Phys. Rev. E 66, 046602 (2002).
  * (13) Fleischer J. W., Segev M., Efremidis N. K. & Christodoulides D. N., Observation of two-dimensional discrete solitons in optically induced nonlinear photonic lattices, Nature, 422, 147 (2003).
  * (14) Freedman, B., Bartal, G., Segev, M., Lifshitz, R., Christodoulides, D. N. & Fleischer, J. W. Wave and defect dynamics in nonlinear photonic quasicrystals, Nature 440, 1166 (2006).
  * (15) Brandes, T. & Kettemann, S., The Anderson transition and its ramifications: Localization, quantum interference, and interactions (Springer, 2003).
  * (16) Morsch, O. & Oberthaler, M., Dynamics of Bose-Einstein condensates in optical lattices. Rev. Mod. Phys. 78, 179 (2006).
  * (17) Billy, J., Sanchez-Palencia, L., Bouyer, P. & Aspect, A., Direct observation of Anderson localization of matter waves in a controlled disorder, Nature 453, 891?894 (2008).
  * (18) Wiersma, D. S., Disordered photonics, Nat. Photon. 7, 188-196 (2013).
  * (19) Segev, M., Silberberg, Y., & Christodoulides, D. N., Anderson localization of light, Nat. Photon. 7, 197-204 (2013).
  * (20) Das Sarma, S., Adam, S., Hwang, E. H., & Rossi, E. Electronic transport in two-dimensional graphene, Rev. Mod. Phys. 83, 407-470 (2011).
  * (21) Lederer, F., Stegeman, G. I., Christodoulides, D. N., Assanto, G., Segev, M., & Silberberg, Y., Discrete solitons in optics, Phys. Rep. 463, 1-126 (2008).
  * (22) Kartashov, Y. V., Astrakharchik, G., Malomed, B., & Torner, L., Frontiers in multidimensional self-trapping of nonlinear fields and matter, Nat. Rev. Phys. 1, 185 (2019).
  * (23) Chen, Z., Segev, M., & Christodoulides, D. N., Optical spatial solitons: historical overview and recent advances, Rep. Prog. Phys. 75, 086401 (2012).
  * (24) Yang J., & Musslimani, Z. H., Fundamental and vortex solitons in a two-dimensional optical lattice, Opt. Lett. 28, 2094 (2003).
  * (25) Efremidis, N. K., Hudock, J., Christodoulides, D. N., Fleischer, J. M., Cohen, O., & Segev M., Two-dimensional optical lattice solitons, Phys. Rev. Lett. 91, 213906 (2003).
  * (26) Neshev, D., Ostrovskaya, E., Kivshar, Y., & Krolikowski, W., Spatial solitons in optically induced gratings, Opt. Lett. 28, 710 (2003).
  * (27) Ablowitz, M.J., Ilan, B., Schonbrun, E. & Piestun R., Solitons in two-dimensional lattices possessing defects, dislocations, and quasicrystal structures, Phys. Rev. E 74, 035601 (R) (2006).
  * (28) Law, K. J. H., Saxena, A., Kevrekidis, P. G. & Bishop, A. R. Stable structures with high topological charge in nonlinear photonic quasicrystals. Phys. Rev. A 82, 035802 (2010).
  * (29) Ablowitz, M. J., Antar, N., Bakirta?, ?. & Ilan, B. Vortex and dipole solitons in complex two-dimensional nonlinear lattices. Phys. Rev. A 86, 033804 (2012).
  * (30) Xavier, J., Boguslawski, M., Rose, P., Joseph, J., & Denz, C., Reconfigurable optically induced quasicrystallographic three-dimensional complex nonlinear photonic lattice structures, Adv. Materials 22, 356 (2010).
  * (31) Chiao, R. Y., Garmire, E., & Townes, C. H., Self-trapping of optical beams, Phys. Rev. Lett. 13, 479 (1964).

## METHODS

### Experimental setup.

The experiment is carried out in a biased SBN: 61 photorefractive crystal
whose dimensions are 5$\times$5$\times$20 mm3. The experimental setup is
similar to the one used in LDTnature . A continuous-wave laser, with
wavelength $\lambda$=532 nm and ordinary polarization, is used to write the
desirable moiré lattices into the sample by superposing two square lattice
patterns (sublattices) that have tunable amplitudes $p_{1,2}$ and twist angles
$\theta$. To probe the induced moiré lattice, an extraordinarily polarized
beam from He-Ne laser ($\lambda$=632.8 nm) is launched into the sample, and
its propagation through the sample is monitored. The power of the probe beam
can be varied in the range from several nano-Watts (nW) to micro-Watts
($\mu$W) by using a variable attenuator to control nonlinear self-action due
to photorefractive nonlinearity. The beam power in pre-set pinholes with
desired radius $R$, connected to the soliton content, is acquired using laser
beam profiler.

### Data Availability

The data that support the findings of this study are available from the
corresponding author F. Y. upon reasonable request.

### Acknowledgments

Q. F., P. W and F. Y acknowledge the support of the National Natural Science
Foundation of China under Grant Nos. No. 61475101 and No. 11690033. Y.V.K. and
L.T. acknowledge support from the Severo Ochoa Excellence Programme, Fundacio
Privada Cellex, Fundacio Privada Mir-Puig, and CERCA/Generalitat de Catalunya.
V. V. K. acknowledges financial support from the Portuguese Foundation for
Science and Technology (FCT) under Contract no. UIDB/00618/2020.

### Author contributions

Q. F and P. W contribute equally to this work. All authors contribute
significantly to the work.

### Competing interests

The authors declare no competing interests.