# NICER Detection of Thermal X-ray Pulsations from the Massive Millisecond
Pulsars PSR J0740$+$6620 and PSR J1614$-$2230

Michael T. Wolff Space Science Division, U.S. Naval Research Laboratory,
Washington, DC 20375-5352, USA  Sebastien Guillot IRAP, CNRS, 9 avenue du
Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France  Université de
Toulouse, CNES, UPS-OMP, F-31028 Toulouse, France  Slavko Bogdanov Columbia
Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York,
NY 10027, USA  Paul S. Ray Space Science Division, U.S. Naval Research
Laboratory, Washington, DC 20375-5352, USA  Matthew Kerr Space Science
Division, U.S. Naval Research Laboratory, Washington, DC 20375-5352, USA
Zaven Arzoumanian  Astrophysics Science Division, NASA’s Goddard Space Flight
Center, Greenbelt, MD 20771, USA  Keith C. Gendreau  Astrophysics Science
Division, NASA’s Goddard Space Flight Center, Greenbelt, MD 20771, USA  M.
Coleman Miller Department of Astronomy and Joint Space-Science Institute,
University of Maryland, College Park, MD 20742-2421, USA  Alexander J.
Dittmann Department of Astronomy and Joint Space-Science Institute, University
of Maryland, College Park, MD 20742-2421, USA  Wynn C. G. Ho Department of
Physics and Astronomy, Haverford College, 370 Lancaster Avenue, Haverford, PA
19041, USA  Lucas Guillemot  Laboratoire de Physique et Chimie de
l’Environnement et de l’Espace, LPC2E, CNRS-Université d’Orléans, F-45071
Orléans, France  Station de Radioastronomie de Nançay, Observatoire de Paris,
CNRS/INSU, F-18330 Nançay, France  Ismael Cognard  Laboratoire de Physique et
Chimie de l’Environnement et de l’Espace, LPC2E, CNRS-Université d’Orléans,
F-45071 Orléans, France  Station de Radioastronomie de Nançay, Observatoire de
Paris, CNRS/INSU, F-18330 Nançay, France  Gilles Theureau  Laboratoire de
Physique et Chimie de l’Environnement et de l’Espace, LPC2E, CNRS-Université
d’Orléans, F-45071 Orléans, France  Station de Radioastronomie de Nançay,
Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France  LUTH, Observatoire
de Paris, PSL Research University, CNRS, Université Paris Diderot, Sorbonne
Paris Cité, F-92195 Meudon, France  Kent S. Wood  TSC, Alexandria, VA, USA

###### Abstract

We report the detection of X-ray pulsations from the rotation-powered
millisecond-period pulsars PSR J0740+6620 and PSR J1614$-$2230, two of the
most massive neutron stars known, using observations with the Neutron Star
Interior Composition Explorer (NICER). We also analyze X-ray Multi-Mirror
Mission (XMM-Newton) data for both pulsars to obtain their time-averaged
fluxes and study their respective X-ray fields. PSR J0740+6620 exhibits a
broad double-peaked profile with a separation of $\sim 0.4$ in phase. PSR
J1614$-$2230, on the other hand, has a broad single-peak profile. The broad
modulations with soft X-ray spectra of both pulsars are indicative of thermal
radiation from one or more small regions of the stellar surface. We show the
NICER detections of X-ray pulsations for both pulsars and also discuss the
phase relationship to their radio pulsations. In the case of PSR J0740+6620,
this paper documents the data reduction performed to obtain the pulsation
detection and prepare for pulse profile modeling analysis.

pulsars: general — pulsars: individual (PSR J1614$-$2230, PSR J0740+6620) —
stars: neutron — X-rays: stars

††facilities: NICER, XMM-Newton††software: astropy (Astropy Collaboration et
al., 2013, ascl:1304.002), PINT (ascl:1902.007,
https://github.com/nanograv/pint), HEAsoft (ascl:1408.004,
https://heasarc.nasa.gov/lheasoft/)

##  1 Introduction

Rotation powered millisecond pulsars (MSPs) are usually old neutron stars
(NSs) with spin periods in the $\sim$1–30 ms range. Their fast rotation is
thought to result from the transfer of angular momentum of accreted matter
from a close binary companion during a so-called “recycling” phase (Alpar et
al., 1982; Radhakrishnan & Srinivasan, 1982). MSPs can produce emission from
radio to gamma-rays through particle acceleration and radiation processes in
the pulsar magnetosphere. However, in the X-ray band substantial emission can
originate from hot spots on the NS surface, and such emission can be an
important probe of NS properties. This paper will focus on the X-ray emission
from two pulsars found to be massive NSs.

Radio pulse timing observations of MSPs have resulted in precisely measured NS
masses near or above 2 M⊙ for several objects to date. For PSR J1614$-$2230,
Demorest et al. (2010) found $M_{\rm NS}=1.97\pm 0.04\,M_{\odot}$ and
Arzoumanian et al. (2018a) found $M_{\rm NS}=1.908\pm 0.016\,M_{\odot}$; for
PSR J0348+0432, Antoniadis et al. (2013) found $M_{\rm NS}=2.01\pm
0.04\,M_{\odot}$; for PSR J0740+6620, Cromartie et al. (2020) found $M_{\rm
NS}=2.14_{-0.09}^{+0.10}\,M_{\odot}$ and Fonseca et al. (2021) found $M_{\rm
NS}=2.08\pm 0.07\,\mbox{$M_{\odot}$}$. A potential fourth pulsar with high
mass is PSR J1811$-$2405 for which Ng et al. (2020) find $M_{\rm
NS}=2.0_{-0.5}^{+0.8}\,M_{\odot}$. Finally, radial velocity measurements in
the optical band have identified potential high NS mass candidates as well.
Romani et al. (2021) measured $M_{\rm NS}=2.13\pm 0.04\,\mbox{$M_{\odot}$}$
for the MSP PSR J1810$+$1744, and Kandel & Romani (2020) report $M_{\rm
NS}=2.28_{-0.09}^{+0.10}\,M_{\odot}$ for the redback binary system PSR
J2215$+$5135 . Indeed, Romani et al. (2021) explore the possibility that a
substantial fraction of the NS population may exist with high masses near 2.0
$M_{\odot}$ and above. Depending on the evolutionary pathway of a given MSP,
an appreciable amount of mass can be accreted from the companion star: for
example, Tauris et al. (2011) find that PSR J1614$-$2230 could have accreted
between $\sim$0.05$\,M_{\odot}$ and $\sim$0.45$\,M_{\odot}$.

Such massive NSs are of great interest for understanding the NS mass-radius
relation and the equation of state (EOS) of cold, dense matter. High-mass NSs
enable us to constrain the very uncertain EOS at the highest densities and
connect with the EOS at lower densities obtained from less massive NSs.
Furthermore, each theoretical model of NS structure predicts a particular
maximum mass (Oppenheimer & Volkoff 1939; see also, e.g., Lattimer 2012),
above which black hole formation occurs. For example, GW190814 indicates the
merger of a black hole and a second compact object with mass 2.6$\,M_{\odot}$
(Abbott et al., 2020). If the second object is also a black hole, then the
maximum mass of a (slowly-rotating) NS is between 2$\,M_{\odot}$ and
2.6$\,M_{\odot}$, and there is either a small gap or no gap between the
highest mass NS and lowest mass black hole. On the other hand, if knowledge
gained from 2$\,M_{\odot}$ NSs indicates an EOS that can support much higher
mass NSs, then objects such as the secondary of GW190814 could be a NS (see
Abbott et al. 2018).

Even more stringent constraints on the dense matter EOS can be obtained if the
mass measurement is combined with a radius estimate from modeling of the
pulsed surface thermal X-ray radiation (see, e.g., Pavlov & Zavlin, 1997;
Miller & Lamb, 1998; Bogdanov et al., 2007; Psaltis & Özel, 2014). Obtaining
NS radius measurements using this technique is one of the primary scientific
goals of the Neutron star Interior Composition Explorer (NICER) mission. The
first such constraints were presented in Miller et al. (2019) and Riley et al.
(2019) based on a deep NICER exposure of the nearby isolated MSP PSR
J0030$+$0451 (Bogdanov et al., 2019).

NICER observations of massive MSPs have the potential to provide unique
constraints on the properties of cold, catalyzed matter beyond nuclear density
even when very X-ray faint, as discussed in depth by Miller (2016) for the
case of PSR J1614$-$2230\. Because more compact stars tend to have lower
fractional modulation amplitudes, all else being equal, measurement of the
X-ray pulsation amplitude will provide an upper limit on the compactness
$GM/Rc^{2}$ and thus a lower limit on the radius given the precisely-
determined mass of this star. A lower limit to the radius of a high-mass NS,
if that lower limit exceeds 10–11 km, would provide important input to nuclear
theories. An upper limit on the stellar radius cannot be obtained in the same
way. For example, moving the spot center toward a rotational pole can reduce
the amplitude by an arbitrary factor regardless of how large the star is. A
robust minimum radius can be inferred, using assumptions that otherwise
maximize the overall amplitude (including the assumptions that the spot is
point like and centered on the rotational equator).

PSR J0740+6620 is an MSP with $P_{\rm{spin}}=2.89$ ms discovered in the course
of the Green Bank Northern Celestial Cap (GBNCC) pulsar survey (Lynch et al.,
2018). It is in a nearly circular 4.8 day binary orbit, possibly with an
ultra-cool white dwarf companion (Beronya et al., 2019). Recently, Fonseca et
al. (2021) obtained an updated mass measurement of $M_{\rm NS}=2.08\pm 0.07$
$M_{\odot}$ (at 68.3% confidence) for PSR J0740+6620, from the observed radio
Shapiro delay, making it the NS with the largest precisely established mass
(first reported by Cromartie et al. 2020). Fonseca et al. (2021) also reported
an updated distance measurement from timing parallax and the Shklovskii effect
of $d=1.14^{+0.17}_{-0.15}$ kpc (68.3% confidence).

PSR J1614$-$2230 was discovered in a Parkes radio search targeting
unidentified EGRET sources (Hessels et al., 2005; Crawford et al., 2006). It
is an MSP with $P_{\rm{spin}}=3.15$ ms bound to a massive white dwarf
companion in a $P_{\rm{b}}=8.7$ d orbit. The pulsar is particularly important
for the NICER mission because it is one of the most massive NSs known, with
the most recent measurement of $M=1.908\pm 0.016\mbox{$\,M_{\odot}$}$
(Arzoumanian et al., 2018b). This precise measurement resulted from the
detection of a strong signature of Shapiro delay in the binary (Demorest et
al., 2010). PSR J1614$-$2230 has a parallax of 1.5$\pm$0.1 mas, so it lies at
a parallax distance of $d=670\mbox{${}^{+50}_{-40}$}$ pc (Arzoumanian et al.,
2018b). X-ray pulsations from this pulsar were first claimed at the $\sim
4\sigma$ confidence level based on XMM-Newton observations (Pancrazi et al.,
2012).

Initial exploratory X-ray observations of both pulsars with NICER, XMM-Newton,
and Swift showed that they are very faint sources. We have undertaken a
systematic survey with NICER of the sample of nearby rotation-powered MSPs to
detect and characterize their pulsed X-ray radiation, with a focus on finding
thermal X-ray pulsations that are desirable for constraints on the NS mass-
radius relation and the dense matter EOS. The NICER observations of eight
other MSPs were presented in Ray et al. (2019) and Guillot et al. (2019).

In the present paper, we describe NICER observations of the two nearby binary
MSP systems hosting massive neutron stars, PSR J0740+6620 and PSR
J1614$-$2230, and report the detection of thermal X-ray pulsations from both.
The paper is organized as follows. In Section 2 we describe the X-ray
observations and the data reduction procedures. Because the photon event data
for PSR J0740$+$6620 have been made available for the NS radius inference
analyses of Miller et al. (2021) and Riley et al. (2021), the data extraction
procedure is more stringent for that pulsar than for PSR J1614$-$2230\. In
Section 3 we present our findings from the X-ray observations together with
the known radio properties of these two pulsars. In Section 4, we discuss
issues that are raised either by the observations themselves or the follow-on
inference analyses. Finally, in Section 5 we summarize our principal
conclusions.

##  2 Observations

The Neutron Star Interior Composition Explorer (NICER; see Gendreau &
Arzoumanian, 2017) is a NASA astrophysics explorer mission that has been
operating on the International Space Station (ISS) since 2017 June. The NICER
X-ray Timing Instrument (XTI) is an array of 52 active silicon drift detectors
each paired with concentrator optics that are nominally sensitive to photons
in the 0.2$-$12 keV range and have absolute timing precision of $\sim$100
nanoseconds rms for tagging of photon detection times (Prigozhin et al.,
2016). The event data for both sources were processed using HEASOFT v6.28 and
the NICER-specific package NICERDAS v6, together with the NICER calibration
files v20200202. However, somewhat different criteria for extracting NICER
event data were used for the two pulsars. Furthermore, the XMM-Newton data
were obtained under a number of different circumstances, as described below.

Table 1: NICER observation summary for PSR J0740+6620 ObsID | Dates | Exposure | Exposure
---|---|---|---
Range | Range | Raw (ks) | Filtered (ks)
1031020101 – 1031020128 | 2018-09-21 – 2018-10-29 | 322.1 | 302.8
2031020101 – 2031020257 | 2019-04-05 – 2020-02-29 | 1142.2 | 983.5
3031020201 – 3031020244 | 2020-03-01 – 2020-04-17 | 403.2 | 360.1
Table 2: NICER observation summary for PSR J1614$-$2230 ObsID | Dates | Exposure | Exposure
---|---|---|---
Range | Range | Raw (ks) | Filtered (ks)
0060310101 – 0060310109 | 2017-07-04 – 2017-07-16 | 17.2 | 15.2
1060310101 – 1060310257 | 2017-07-18 – 2019-03-01 | 410.3 | 332.3
2060310201 – 2060310264 | 2019-03-06 – 2019-09-12 | 143.7 | 126.8
3060310201 – 3060310298 | 2020-04-02 – 2021-03-01 | 281.3 | 260.0

###  2.1 NICER

PSR J0740+6620 was observed by NICER for 1867.5 ks between 2018 September 21
and 2020 April 17, while PSR J1614$-$2230 was observed by NICER for 852.5 ks
over the period spanning 2017 July 4 to 2021 March 1. Both pulsars are
challenging to observe because of their faintness, although PSR J1614$-$2230
is somewhat brighter than PSR J0740+6620. Because both are faint, everything
possible has been done to minimize backgrounds, including making observations
primarily during ISS orbital night. Special care is warranted for PSR
J0740+6620 because data sets are being used to obtain a NS radius constraint
from the light curve modeling of the pulsed component (which is derived from
NICER) relative to the total phase-averaged flux (which is derived from XMM-
Newton imaging). Furthermore, there are nearby sources of comparable
brightness that the NICER concentrating detector optics do not exclude and
their contributions must be estimated. A summary of NICER observations of our
target pulsars can be found in Tables 1 and 2.

For the pulsar PSR J0740+6620, all events included in our and subsequent
analyses were obtained when the pulsar was at least 80∘ from the Sun. This
minimized the effects of scattered solar X-rays and detector noise resulting
from photons whose energies are just below the NICER energy range. We extract
the event data in an ObsID by ObsID manner. We first filter out event data
outside the 0.25$-$3.0$\,{\rm keV}$ energy range. We also filter the
observations keeping only those times that have a Cosmic Ray Cutoff Rigidity
greater than 2.0 GeV/$c$ and a space weather Kp index less than 5.0. We
further require that the observations by NICER include the full complement of
52 detectors active during the observations. We exclude the detector with
DET_ID 34 since this detector tends to be noisy. This ensures that the event
data passed to the light curve analysts have an internally consistent set of
51 out of 52 available detectors and thus utilize a constant effective area
configuration.

We also demand that NICER is always pointed to better than 15 arcseconds of
the source in order to include the data. The NICER vignetting function is
known to be relatively flat-topped within 2 arcminutes of the instrument
boresight so the pointing accuracy filter removes rare, large pointing errors
as a possible significant source of brightness variation. In August 2020
maintenance on the individual detectors was begun and this involved cycling
individual detectors off for a period of annealing. Thus, for our initial
investigation of PSR J0740+6620, utilizing event data obtained after the
detector maintenance was started would have created a NICER XTI effective area
inconsistency across this date boundary and we choose to avoid this
complication. The pulsar phase for each event is calculated with the
photonphase tool of the PINT software package (Jing et al., 2020) using the
timing solutions for each pulsar described in §2.3. The end result here is a
set of good time intervals (GTIs) with event data and low background noise.

For PSR J1614$-$2230, because we are not yet at the stage of performing a full
light curve analysis to obtain an estimated NS radius, we slightly relaxed our
event acceptance criteria. Namely, we employ a Kp index constraint of 4.0 and
a Sun angle constraint of 60∘. Furthermore, some of the NICER detectors
suffered from a brief single event upset (SEU) on 8 July 2019, that caused
them to develop anomalous time stamps until those detectors could be reset. We
exclude the three days of NICER event data between the SEU occurrence and the
reset. Because PSR J1614$-$2230 is somewhat brighter than PSR J0740+6620 this
change in event acceptance criteria was found to be the best trade-off between
integration time accumulation and background noise.

The NICER non-imaging field of view (FOV) with a $\sim 3\hbox{${}^{\prime}$}$
radius does not resolve X-ray sources near the targets and thus such sources
can contribute to the effective background counts for the target source.
Therefore, the true X-ray phase-averaged count rate for faint pulsars may be
difficult to accurately determine. In order to partially mitigate this issue,
we make use of XMM-Newton observations to better determine the X-ray fluxes of
the pulsars and investigate the presence of other X-ray sources in the NICER
FOV.

###  2.2 XMM-Newton

Table 3: XMM-Newton observation summary for PSR J0740+6620

Instrument | ObsID | Date | Exposure | Exposure
---|---|---|---|---
|  |  | Raw (ks) | Filtered (ks)
XMM-MOS1 | 0851181601 | 2019-10-26 | 6.55 | 5.74
| 0851181401 | 2019-10-28 | 8.69 | 4.94
| 0851181501 | 2019-11-01 | 11.45 | 7.28
XMM-MOS2 | 0851181601 | 2019-10-26 | 9.33 | 7.49
| 0851181401 | 2019-10-28 | 8.68 | 3.33
| 0851181501 | 2019-11-01 | 11.44 | 7.86
XMM-pn | 0851181601 | 2019-10-26 | 6.42 | 4.11
| 0851181401 | 2019-10-28 | 6.00 | 0.91
| 0851181501 | 2019-11-01 | 8.42 | 1.78

Note. — All XMM-Newton exposures were obtained in ‘Full-frame’ mode. The “Raw”
exposure times report the value of the “LIVETIME” FITS keyword. The “Filtered”
exposure times report the sum of the GTIs after the filtering described in
Section 2.2.

Table 4: XMM-Newton observation summary for PSR J1614$-$2230

Instrument | ObsID | Date | Exposure | Exposure
---|---|---|---|---
|  |  | Raw (ks) | Filtered (ks)
XMM-MOS1 | 0404790101 | 2007-02-08 | 45.6 | 37.3
XMM-MOS2 | 0404790101 | 2007-02-08 | 45.7 | 37.0
XMM-pn | 0404790101 | 2007-02-08 | 38.7 | 21.9

Note. — All XMM-Newton exposures were obtained in ‘Full-frame’ mode. The “Raw”
exposure times report the value of the “LIVETIME” FITS keyword. The “Filtered”
exposure times report the sum of the GTIs after the filtering described in
Section 2.2.

For PSR J0740+6620, we acquired three separate XMM-Newton data sets on 2019
October 26 (ObsID 0851181601), 2019 October 28 (ObsID 0851181401), and 2019
November 1 (ObsID 0851181501) through the Director’s Discretionary Time
program. The three EPIC instruments were operated in ‘Full Frame’ mode with
the ‘Thin’ optical blocking filters in place. Due to the long instrument event
readout times in ‘Full Frame’ mode (73.4 ms for EPIC-pn and 2.6 s for EPIC-
MOS), the data cannot be used for a pulse timing analysis. As these
observations occurred close to the telescope perigee passage, a large fraction
of the exposures were affected by intense particle flaring. These are filtered
out using the standard procedure, with special care to exclude times of high
background flaring. For the EPIC-MOS cameras, we extracted events with
PATTERN<=12, and for the EPIC-pn camera we extracted events with PATTERN<=4.
For both instruments we excluded time ranges with significantly higher than
average count rates in both the full 0.2–12 keV band and the high energy bands
(10–12 keV for EPIC-MOS and 10–15 keV for EPIC-pn). This procedure resulted in
the exposure times listed in Table 3.

For PSR J1614$-$2230, there are four separate observations of the pulsar in
the XMM-Newton HEASARC archive. We utilize archival XMM-Newton observations
(ObsID: 0404790101) originally reported by Pancrazi et al. (2012).

###  2.3 Radio pulsar timing

To search for pulsations, we require pulsar timing models that provide
sufficient precision over the full NICER observation span to assign pulse
phases to the NICER data without needing any additional trials over folding
periods. These models are provided by ground-based pulsar timing programs at
radio frequencies.

PSR J0740+6620 is regularly observed as part of the NANOGrav pulsar timing
program (Alam et al., 2021) with both the Green Bank Telescope (GBT) and
Canadian Hydrogen Intensity Mapping Experiment (CHIME). Fonseca et al. (2021)
have performed an analysis of these data, extending the work done by Cromartie
et al. (2020). Since this covers the entire span of the NICER data, we use
that timing model to assign pulse phases for the present work.

For PSR J1614$-$2230, we constructed a timing model using a merged dataset
from the NANOGrav 12.5yr release (Alam et al., 2021) and the Nançay radio
telescope (NRT). The NRT data were obtained from MJD 57424 (2016 Feb 6)
through 58921 (2020 Mar 13). The NRT backends and data processing procedure
are described in Guillemot et al. (2016). We analyzed this combined data set
using Tempo2 and fit a pulsar timing model containing astrometric, rotational,
and binary parameters for the pulsar. The model included a time-variable
dispersion measure (DM) using a second-order polynomial plus a solar wind
model. The resulting pulse timing solution is presented in Table 5.

Table 5: Ephemeris of PSR J1614$-$2230\. Digits in parentheses represent the $1\sigma$ uncertainty on the last quoted digit of a parameter value.  Parameter | Value
---|---
Right Ascension (J2000) | 16:14:36.507995(7)
Declination (J2000) |  $-22$:30:31.3329(5)
Proper motion in R.A. (mas yr-1) | 3.87(3)
Proper motion in Decl. (mas yr-1) | $-32.3(2)$
Parallax, $\pi$ (mas) | 1.58(4)
Epoch of position (MJD) | 56823
Spin frequency, $\nu$ (Hz) | 317.3789418919128(3)
Spin frequency derivative, $\dot{\nu}$ (s-2) | $-9.69490(4)\times 10^{-16}$
Epoch of period (MJD) | 56823
Dispersion measure, DM (pc cm-3) | 34.48759(3)
Dispersion measure derivative, DM1 | $-0.000302(6)$
Dispersion measure 2nd derivativeaaThe default definition of DM2 (and higher terms) in Tempo2 changed in 2020 June. We use the corrected Taylor series representation here., DM2 | $-1.75(4)\times 10^{-4}$
Epoch of DM (MJD) | 56823
Solar wind density at 1 AU (NE_SW; cm-3) | 7.33(16)
Binary model | ELL1
Binary orbital period, $P_{\rm{b}}$ (days) | 8.6866194173(2)
First derivative of orbital period, $\dot{P_{b}}$ | 1.28(5)$\times 10^{-12}$
Projected semi-major axis of orbit, $x$ (lt-s) | 11.29119744(4)
Epoch of ascending node passage, Tasc (MJD) | 52331.17011020(5)
First Laplace parameter, $\sin\,\omega$ | $9.7(6)\times 10^{-8}$
Second Laplace parameter $\cos\,\omega$ | $-1.345(4)\times 10^{-6}$
Companion mass, $M_{c}$ ($M_{\odot}$) | 0.4943(19)
Sine of orbital inclination, $\sin i$ | 0.999896(4)
Terrestrial time standard (CLK) | TT(BIPM2017)
Units of barycentric time | TDB
Solar System ephemeris | DE421

##  3 Pulsation Detection

Figure 1: NICER pulse profiles of (a) PSR J0740+6620 in the energy range
0.30–1.25 keV and of (b) PSR J1614$-$2230 in the energy range 0.59–2.20 keV,
folded with their respective radio spin ephemerides. Two rotational cycles are
shown for clarity. The blue dashed lines correspond to the radio profiles
obtained with the NRT in the 1.4 GHz band. Figure 2:  (Left) XMM-Newton
combined MOS1 and MOS2 0.5$-$8.0 keV image of PSR J0740$+$6620 and nearby
sources. The green cross shows the radio position of the pulsar and the inner
and outer white circles are 2 and 4 arcminutes away from the pulsar where the
NICER vignetting function drops quickly from 0.90 to 0.11, respectively. The
AGN source SDSS J074115.14+662234.9 is 3.57 arcminutes to the north-east of
PSR J0740$+$6620\. At this position SDSS J074115.14+662234.9 falls at the 21%
position on the NICER vignetting function. (Right) XMM-Newton combined MOS1
and MOS2 0.5$-$8.0 keV image of PSR J1614$-$2230 and nearby sources. The green
cross shows the radio position of the pulsar and the inner and outer white
circles are 2 and 4 arcminutes away from the pulsar, respectively. The blue
crosses are nearby sources from the Chandra source catalog (Evans et al.,
2010).

After standard processing and filtering of the NICER data as described above,
we isolate the events from the 0.25 to 3.0 keV range and search for
pulsations. Because for both pulsars we have established radio ephemerides we
ignore trials factors and simply assign pulse phases as described above. For
PSR J0740+6620 a simple Z${}_{2}^{2}$ test yields a test statistic of
Z${}_{2}^{2}=251.19$, corresponding to a detection significance of
$15.35\sigma$ for 574405 events in the 0.31 to 1.22 keV range. We take
advantage of the fact that the filtered NICER count rate (about 0.6–1.0 s-1)
is dominated by the background, with the pulsar contributing only a few
percent to the total count rate. We can use the total count rate in each GTI
interval as a proxy for the effective background count rate and employ the GTI
sorting method described in Guillot et al. (2019). The GTI sorting method
breaks the original GTI intervals into shorter segments and adds each segment
into the evaluated time series, making sure that the overall pulsation
significance increases with each addition. The lowest count rate segments are
included first, and we add in segments that have increasing count rates (thus
higher background) so long as the detection significance continues to
increase. Once the background noise becomes sufficiently large that adding
additional GTIs from the original distribution decreases rather than increases
the detection significance, we stop the procedure. When we apply this method
to PSR J0740+6620 we increase the overall pulsation detection significance to
$\sigma=15.46$ for 521004 events in the 0.31 to 1.18 keV range. The GTI
selection is particularly effective for faint pulsars with count rates $\sim
0.01$ c/s such as PSR J0740+6620. The event lists and GTIs are then saved in
both text and FITS file format so that they can be utilized in subsequent
processing, in particular the light curve modeling of Miller et al. (2021) and
Riley et al. (2021) .

For PSR J1614$-$2230 we confirm the pulsations originally found by Pancrazi et
al. (2012) using XMM-Newton timing-mode data. Searching the originally
extracted 734.3 ks of event data for this pulsar in the 0.25 to 3.00 keV range
detects the pulsations with significance Z${}_{2}^{2}=72.89$ giving
$\sigma=7.81$ for 335614 events in the 0.59 to 2.20 keV range. When we
optimize the good time intervals as above, our detection significance rises to
$\sigma=8.31$ but the total time is cut somewhat to 644.9 ks, the number of
included events is likewise reduced to 267122, but the optimal energy range is
the same at 0.59 to 2.20 keV.

Figure 1 shows the folded NICER pulse profiles of PSRs J0740$+$6620 and
J1614$-$2230 in the energy bands where the maximum pulse detection
significance is found as determined by the Z${}_{2}^{2}$-test (de Jager et
al., 1989; de Jager & Büsching, 2010), together with their radio pulse
profiles aligned in rotational phase according to the absolute pulse arrival
times provided by the radio timing model of each pulsar.

##  4 Discussion

The NICER observations presented here reveal PSR J0740$+$6620 and PSR
J1614$-$2230 as pulsed millisecond-variable X-ray sources with phase-broadened
thermal pulsations. Both NSs have masses near 2 $M_{\odot}$ and thus are
important X-ray sources to elucidate the cold matter EOS near the maximum NS
mass. PSR J0740$+$6620 shows a double pulse light curve in the NICER energy
range while PSR J1614$-$2230 shows one broad pulse. The relative phasing of
the radio and X-rays shows that for PSR 0740+6620 the radio pulses are aligned
with the thermal X-ray pulses and that the X-ray pulses are separated by
roughly 0.4 in phase. In the case of PSR J1614-2230 the dominant radio pulse
occurs near the end of the broad single thermal pulse although there is a very
small second pulse in the radio light curve. Both pulsars are established as
emitting from radio to $\gamma$-ray energies, indicating that non-thermal
emission from the magnetosphere is present. We found no evidence, however, of
non-thermal X-ray emission above 2.5 keV, almost certainly because such
emission is too faint to be detected in these pulsars.

The goals of these NICER and XMM-Newton observations are to establish the
existence of X-ray pulsations from both pulsars and to estimate the time-
averaged flux of each pulsar. The XMM-Newton measurements help constrain the
total flux from the NS, enabling a better estimate of the pulsed emission
fraction in the NICER data, and thus the estimate of the NS compactness. A
correct pulsed NICER flux from the pulsar, if paired with an erroneously high
pulse-averaged flux from XMM-Newton, would yield a decreased pulsed fraction
and thus an incorrect lower NS radius. Conversely, a spuriously low XMM-Newton
flux would yield a larger NS radius. Because of the high NS masses in these
systems, even a solid lower limit in the radius of each pulsar would be
valuable for constraining the NS EOS because the cores of such massive NSs
contain higher-density matter. An upper limit on the compactness, along with
the known NS masses, can yield a lower limit on the pulsar radius.

These NICER data sets for PSR J0740$+$6620 and PSR J1614-2230 are very quiet
in that the background count rates are strongly minimized. Also, for PSR
J0740$+$6620 the effective area as a function of time should be very constant
because we always utilized the events from 51 out of 52 active detectors.

The XMM-Newton observations of the field of PSR J0740$+$6620 allow us to
estimate the expected count rates of the likely AGN in the field that is 3.57
arcminutes north and east from the pulsar. In fact, investigation of this
X-ray point source shows that there are really two extragalactic sources
within 3.157 arcseconds of one another [SDSS J074114.62+662235.2 and SDSS
J074115.14+662234.9, from the Sloan Digital Sky Survey Data Release 9 (Ahn et
al., 2012)] so the XMM-Newton observations we obtained can not easily
distinguish if the X-rays are coming from one or both of these extragalactic
sources. Using the SAS utility eregionanalyse we estimate a For the position
of SDSS J074115.14+662234.9 eregionanalyse gives a count rate of
$\sim$3.1$\times 10^{-2}$ for this AGN and we assume a non-thermal spectrum
with power law spectral index $\Gamma\sim 2.71$. However, since the AGN (at
3.57 arcminutes away) is at the 21% level of the NICER vignetting curve
(Takahashi, 2013) as a function of angular radius from the look direction of
the XTI, this implies the AGN will only add $\sim 0.7$ when NICER is pointed
directly at PSR J0740$+$6620\. Weak sources in the field will also contribute
some counts because the NICER concentrator optics will cause some photons from
these sources to fall sufficiently close to the center of the detectors to be
included in the background.

##  5 Conclusions

The event data from NICER and XMM-Newton we described here for PSR J0740+6620
were used by Miller et al. (2021); Raaijmakers et al. (2021); Riley et al.
(2021) to investigate the radius and EOS of the pulsar in this system. The
event list we derive for this pulsar is available for download on Zenodo.org
website111https://zenodo.org/. For PSR J0740+6620 the NICER observations have
shown X-ray pulsations with two broad thermal X-ray peaks separated in phase
by about $\sim$0.4. In the case of PSR J1614$-$2230, we have confirmed the
X-ray pulsations originally reported by Pancrazi et al. (2012) and find one
broad thermal X-ray peak in the phase diagram. We show the fields of these
pulsars in XMM-Newton images and find that both pulsars are faint and nearby
X-ray sources will contribute significantly to their respective NICER
concentrator count rates.

This work was supported in part by NASA through the NICER mission and the
Astrophysics Explorers Program. NICER science team members at NRL are
supported by NASA under Interagency Agreement NNG200808A. S.B. was funded in
part by NASA grants NNX17AC28G and 80NSSC20K0275. M.C.M. thanks the Radboud
Excellence Initiative for supporting his stay at Radboud University. W.C.G.H.
acknowledges support through grant 80NSSC20K0278 from NASA. This research has
made use of data and/or software provided by the High Energy Astrophysics
Science Archive Research Center (HEASARC), which is a service of the
Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics
Division of the Smithsonian Astrophysical Observatory. We acknowledge use of
NASA’s Astrophysics Data System (ADS) bibliographic services and the ArXiv.
The authors thank the XMM-Newton observatory, specifically the help of Felix
Fuerst, for the Director’s Discretionary Time observations of PSR J0740+6620.
, DS9 ( Joye & Mandel (2003), ascl:0003.002), emcee (ascl:1303.002,
https://github.com/dfm/emcee), and the XMM-Newton Scientific Analysis System
(SAS, ascl:1404.004, https://www.cosmos.esa.int/web/xmm-newton/sas).

## References

  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101
  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJ, 896, L44
  * Ahn et al. (2012) Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2012, ApJS, 203, 21
  * Alam et al. (2021) Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021, ApJS, 252, 4
  * Alpar et al. (1982) Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Nature, 300, 728
  * Antoniadis et al. (2013) Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448
  * Arzoumanian et al. (2018a) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018a, ApJS, 235, 37
  * Arzoumanian et al. (2018b) —. 2018b, The Astrophysical Journal Supplement Series, 235, 37
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Beronya et al. (2019) Beronya, D. M., Karpova, A. V., Kirichenko, A. Y., et al. 2019, MNRAS, 485, 3715
  * Bogdanov et al. (2019) Bogdanov, S., Guillot, S., Ray, P., et al. 2019, ApJ, 887, L25
  * Bogdanov et al. (2007) Bogdanov, S., Rybicki, G. B., & Grindlay, J. E. 2007, ApJ, 670, 668
  * Cognard et al. (2017) Cognard, I., Freire, P. C. C., Guillemot, L., et al. 2017, ApJ, 844, 128
  * Crawford et al. (2006) Crawford, F., Roberts, M. S. E., Hessels, J. W. T., et al. 2006, ApJ, 652, 1499
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nature Astronomy, 4, 72
  * de Jager & Büsching (2010) de Jager, O. C., & Büsching, I. 2010, A&A, 517, L9
  * de Jager et al. (1989) de Jager, O. C., Raubenheimer, B. C., & Swanepoel, J. W. H. 1989, A&A, 221, 180
  * Demorest et al. (2010) Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081
  * Evans et al. (2010) Evans, I. N., Primini, F. A., Glotfelty, K. J., et al. 2010, ApJS, 189, 37
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, arXiv e-prints, arXiv:2104.00880
  * Gendreau & Arzoumanian (2017) Gendreau, K., & Arzoumanian, Z. 2017, Nature Astronomy, 1, 895
  * Guillemot et al. (2016) Guillemot, L., Smith, D. A., Laffon, H., et al. 2016, A&A, 587, A109
  * Guillot et al. (2019) Guillot, S., Kerr, M., Ray, P. S., et al. 2019, ApJ, 887, L27
  * Guillot et al. (2016) Guillot, S., Kaspi, V. M., Archibald, R. F., et al. 2016, MNRAS, 463, 2612
  * He et al. (2013) He, C., Ng, C. Y., & Kaspi, V. M. 2013, ApJ, 768, 64
  * Hessels et al. (2005) Hessels, J., Ransom, S., Roberts, M., et al. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 328, Binary Radio Pulsars, ed. F. A. Rasio & I. H. Stairs, 395
  * Jing et al. (2020) Jing, L., et al. 2020, ApJ, submitted
  * Joye & Mandel (2003) Joye, W. A., & Mandel, E. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook, 489
  * Kandel & Romani (2020) Kandel, D., & Romani, R. W. 2020, ApJ, 892, 101
  * Lattimer (2012) Lattimer, J. M. 2012, Annual Review of Nuclear and Particle Science, 62, 485
  * Linares et al. (2018) Linares, M., Shahbaz, T., & Casares, J. 2018, ApJ, 859, 54
  * Lynch et al. (2018) Lynch, R. S., Swiggum, J. K., Kondratiev, V. I., et al. 2018, ApJ, 859, 93
  * Miller (2016) Miller, M. C. 2016, ApJ, 822, 27
  * Miller & Lamb (1998) Miller, M. C., & Lamb, F. K. 1998, ApJ, 499, L37
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24
  * Miller et al. (2021) —. 2021, ApJ, submitted
  * Ng et al. (2020) Ng, C., Guillemot, L., Freire, P. C. C., et al. 2020, MNRAS, 493, 1261
  * Oppenheimer & Volkoff (1939) Oppenheimer, J. R., & Volkoff, G. M. 1939, Physical Review, 55, 374
  * Pancrazi et al. (2012) Pancrazi, B., Webb, N. A., Becker, W., et al. 2012, A&A, 544, A108
  * Pavlov & Zavlin (1997) Pavlov, G. G., & Zavlin, V. E. 1997, ApJ, 490, L91
  * Prigozhin et al. (2016) Prigozhin, G., Gendreau, K., Doty, J. P., et al. 2016, in Proc. SPIE, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, 99051I
  * Psaltis & Özel (2014) Psaltis, D., & Özel, F. 2014, ApJ, 792, 87
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Riley, T. E., et al. 2021, ApJ, submitted
  * Radhakrishnan & Srinivasan (1982) Radhakrishnan, V., & Srinivasan, G. 1982, Current Science, 51, 1096
  * Ray et al. (2019) Ray, P. S., Guillot, S., Ransom, S. M., et al. 2019, ApJ, 878, L22
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, submitted
  * Romani et al. (2021) Romani, R. W., Kandel, D., Filippenko, A. V., Brink, T. G., & Zheng, W. 2021, ApJ, 908, L46
  * Takahashi (2013) Takahashi, T. 2013, Mem. Soc. Astron. Italiana, 84, 776
  * Tauris et al. (2011) Tauris, T. M., Langer, N., & Kramer, M. 2011, MNRAS, 416, 2130


# Arbitrage Equilibrium, Invariance, and the Emergence of Spontaneous Order
in the Dynamics of Birds Flocking

Abhishek Sivaram  Complex Resilient Intelligent Systems Laboratory, Department
of Chemical Engineering, Columbia University, New York, NY 10027, U.S.A.
Venkat Venkatasubramanian<EMAIL_ADDRESS>Complex Resilient Intelligent
Systems Laboratory, Department of Chemical Engineering, Columbia University,
New York, NY 10027, U.S.A.

###### Abstract

The physics of active biological matter, such as bacterial colonies and bird
flocks, exhibiting interesting self-organizing dynamical behavior has gained
considerable importance in recent years. Recent theoretical advances use
techniques from hydrodynamics, kinetic theory, and non-equilibrium statistical
physics. However, for biological agents, these don’t seem to recognize
explicitly their critical feature, namely, the role of survival-driven purpose
and the attendant pursuit of maximum utility. Here, we propose a novel game-
theoretic framework and show a surprising result that the bird-like agents
self-organize dynamically into flocks to approach a stable arbitrage
equilibrium of equal effective utilities. While it has been well-known for
three centuries that there are constants of motion for passive matter, it
comes as a surprise to discover that the dynamics of active matter populations
could also have an invariant. What we demonstrate is for ideal systems,
similar to the ideal gas or Ising model in thermodynamics. The next steps
would involve examining and learning how real swarms behave compared to their
ideal versions. Our theory is not limited to just birds flocking but can be
adapted for the self-organizing dynamics of other active matter systems.

Active matter, Arbitrage equilibrium, Invariance, Self-organization,
Statistical teleodynamics, Flocking, Emergence

##  I Introduction

Active matter describes systems composed of large numbers of self-actualizing
agents, which consume and dissipate energy resulting in interesting dynamical
behavior [1, 2, 3, 4, 5, 6, 7, 8, 9]. Biological examples of such systems
include self-organizing bio-polymers, bacteria, schools of fish, flocks of
birds, and so on. This study is focused on birds flocking.

Flocking has been studied extensively from dynamical systems and statistical
mechanics perspectives [10, 11, 12, 13, 14]. Such analyses have contributed
substantially to our evolving understanding of interesting emergent properties
such as phase segregation, flock stability, etc. However, these approaches
don’t seem to model explicitly the critical feature of active biological
agents, namely, the role of purpose and its attendant pursuit of maximum
utility. Being biological agents, birds are innately purposeful, driven by the
goal to survive and thrive in challenging environments as Darwin explained. We
believe any comprehensive theory of active biological matter has to overtly
account for this defining characteristic of the agents.

We address this need by using a novel game-theoretic framework, which we call
statistical teleodynamics [15, 16, 17, 18, 19, 20]. The name comes from the
Greek word telos, which means goal. Just as the dynamical behavior of gas
molecules is driven by thermal agitation (hence, thermodynamics), the dynamics
of purposeful agents is driven by the pursuit of their goals and, hence,
teleodynamics. Statistical teleodynamics may be considered as the natural
generalization of statistical thermodynamics for purpose-driven agents in
active matter. It is a synthesis of the central concepts and techniques of
potential games theory with those of statistical mechanics towards a unified
theory of emergent equilibrium phenomena in active and passive matter [20].

In this paper, we study a model of birds flocking that is inspired by the
well-known Reynolds’s Boids model, analytically and computationally. In the
next two sections, we briefly review the dynamical systems and statistical
mechanics perspectives, respectively. This is followed by our statistical
teleodynamics formulation, which starts with a quick introduction to potential
games. We conclude with a discussion of the main results and their
implications.

##  II Dynamical Models of Flocking

Flocking has been studied using dynamical models, which describe the time
evolution of the position $\mathbf{r}_{i}$ of the $i^{\text{th}}$ agent and
its velocity $\mathbf{v}_{i}$ using pre-specified rules for change. The state
of the flock at any given time is specified by specifying $\mathbf{r}_{i}$ and
$\mathbf{v}_{i}$ for all agents. When the agents move at a constant speed
$v_{0}$, the state of the system is then determined by the set of agents’
positions and velocity directions or orientations
$\\{\mathbf{r}_{i},\mathbf{s}_{i}\\}_{i=1}^{N}$, where $N$ is the number of
agents. Both the Reynolds’s model [10] and the Vicsek model [11] describe the
time evolution of an agent’s velocity, but using different force models.

An agent $i$ is said to be affected an agent $j$, if $j$ is in the
neighborhood of $i$, $\mathcal{N}^{i}$. The neighborhood $\mathcal{N}^{i}$ of
$i$ is defined by a matrix whose elements are $n_{ij}$, where

$\displaystyle n_{ij}=\begin{cases}1&\text{if $j$ is a neighbor of $i$}\\\
0&\text{otherwise}\\\ \end{cases}$ (1)

The span of the neighborhood is specified in terms of the absolute distance
between $i$ and $j$, and a size parameter $r_{0}$ [11, 21], such that

$\displaystyle
n_{ij}=\begin{cases}1&|\mathbf{r}_{i}^{(t)}-\mathbf{r}_{j}^{(t)}|<r_{0}\\\
0&\text{otherwise}\\\ \end{cases}$ (2)

As seen from Eq. 2, we consider agent $i$ to be its own neighbor. One can also
define a neighborhood in terms of a fixed topology of nearest neighbors [12,
13, 14], but we don’t use this specification in this study. It follows that
the number of neighbors of an agent $i$ is given by $n_{i}=\sum_{j}n_{ij}$.

In the Reynolds’s model, the agents (called boids) obey the following three
rules as they fly around:

  1. 1.

Rule of cohesion: A boid steers to move towards the average position of local
flockmates

  2. 2.

Rule of separation: A boid steers to avoid collision and crowding of local
flockmates

  3. 3.

Rule of alignment: A boid steers towards the average heading of local
flockmates

In general, we can write the net effect of these forces on the $i^{\text{th}}$
boid by the equation (Eq. 3),

$\displaystyle\mathbf{v}_{i}(t+\Delta
t)=\mathbf{v}_{i}(t)+\left(\frac{a}{n_{i}}\sum_{j}n_{ij}(\mathbf{r}_{j}-\mathbf{r}_{i})+b\sum_{j}n_{ij}(\mathbf{r}_{i}-\mathbf{r}_{j})\frac{c}{n_{i}}\sum_{j}n_{ij}(\mathbf{v}_{j}-\mathbf{v}_{i})+\bm{\eta}(t)\right)\Delta
t$ (3)

where $a,b,$ and $c$ are parameters corresponding to the rules of cohesion,
separation, and alignment, respectively. Parameter $\bm{\eta}$ is the
uncorrelated noise in the agent’s velocity. The time-scale, $\Delta t$, in Eq.
3 can be subsumed in $a,b,c$. The Vicsek model, similarly, updates velocity
purely as a function of the alignment of the agent with its neighbors, though
modifications have been proposed to include pairwise forces [13, 21] (see
also, Supplementary Information, Section 2).

It’s important to note that in both the Reynolds model and the Viscek model,
there is no concept of the “final” state, or an equilibrium state, of the
system as time tends to infinity. In this regard, they are like molecular
dynamics simulations of molecules where there is no concept of a final
equilibrium state in the equations. They can only determine the immediate next
move of the molecules, at any given time, not their final configurations. The
final outcome is determined, a posteriori, after the simulations are run for a
long time. We highlight this important point here as our theory differs
conceptually in this regard.

While dynamical models of flocking have been studied extensively in
literature, work has also been reported using statistical mechanics to analyze
the flocking behavior. These methods typically use a maximum entropy
formulation via an Ising-model inspired Hamiltonian of the boids’ interaction
[13, 14].

##  III Statistical Teleodynamics of Flocking – A Game-theoretic Formulation

The Reynolds and Vicsek models specify bottom-up agent-level dynamical
behavior, but they don’t provide an analytical framework to predict the
dynamics of the entire flock. This is determined only computationally via
agent-based simulations. That is, there is no analytical framework to derive
the behavior of the whole from the behaviors of the parts.

On the other hand, the statistical mechanics formulation [13] is a top-down
approach that starts with the specification of the Hamiltonian of the flock
and then imposes the maximum entropy distribution on it. It is not clear why
maximum entropy, which is obviously relevant for passive matter systems, would
be applicable for survival-driven birds. The typical statistical mechanics
approach uses the superficial similarity between spins in magnetic systems
(e.g., the Ising model) and the orientation of the birds to apply maximum
entropy methods. The deeper question of why this is conceptually relevant for
the birds is not addressed. Most importantly, all these approaches don’t seem
to recognize explicitly that active agents such as birds act instinctively to
improve their survival prospects.

We address these challenges using our statistical teleodynamics framework [15,
16, 17, 18, 19, 20]. In this theory, the fundamental quantity is an agent’s
effective utility which is a measure of the net benefits pursued by the agent.
Every agent behaves strategically to increase its effective utility by
switching states and exploiting arbitrage opportunities.

In our theory of flocking, we propose that birds are arbitrageurs that always
maneuver to increase their effective utilities, which determine their survival
prospects, dynamically in flight. The effective utility of a bird depends on
its position, speed, and alignment with the rest of the members in its
neighborhood.

Thus, we interpret the three rules of engagement for Reynolds’s boids not as
externally imposed forces on the agents but as innate, self-actualizing,
properties of the agents acquired over millions of years of Darwinian
evolution. These instinctive characteristics enable the agents to incessantly
search for better effective utilities in order to improve their survival
chances.

Hence, we believe that the proper formulation of flocking ought to start with
a model of effective utility that a bird uses to make such decisions
dynamically in flight. Seen from this perspective, we suggest that birds do
not fly randomly (as statistical mechanics-based formulations implicitly
assume), but maneuver strategically to improve their utilities. We exploit
this critical insight to model the dynamical behavior of birds in flight by
using the concepts and techniques from potential games.

In potential games, there exists a single scalar-valued global function,
called a potential ($\phi(\mathbf{x})$) that has the necessary information
about the payoffs or the utilities of the agents. The gradient of the
potential is the utility, $h_{i}$, of the $i$th agent [22, 23, 24, 25, 15].
Therefore, we have

${h}_{i}(\mathbf{x})\equiv{\partial\phi(\mathbf{x})}/{\partial x_{i}}$ (4)

where $x_{i}=N_{i}/N$ and $\mathbf{x}$ is the population vector. A potential
game reaches equilibrium, called Nash equilibrium (NE), when the potential
$\phi(\mathbf{x})$ is maximized. Furthermore, this Nash equilibrium is unique
if $\phi(\mathbf{x})$ is strictly concave [23]. At Nash equilibrium, all
agents enjoy the same effective utility, i.e., $h_{i}=h^{*}$. In fact, the
equality of effective utilities in all states is the fundamental criterion of
game-theoretic equilibrium for active matter. It is an arbitrage equilibrium
[19] where the agents don’t have any incentive to switch states anymore as all
states provide the same effective utility $h^{*}$. Thus, the maximization of
$\phi$ and $h_{i}=h^{*}$ are exactly equivalent criteria, and both specify the
same outcome, namely, an arbitrage equilibrium.

There is a deep and beautiful connection between potential game theory and
statistical mechanics as discussed by Venkatasubramanian [17, 20]. Since an
elaborate discussion about this would take us afar from the objectives of this
paper, we refer interested readers to [17, 20] .

###  III.1 Garud’s Utility: Position Dependence

Our goal here is to develop a simple model of the effective utility ($h_{i}$)
of our boid-like agent, called garud (after the legendary king of birds,
Garuda, in Indian mythology).

We want the model to be an appropriate coarse-grained description of the
system that can make useful predictions not restricted by system-specific
nuances. We have tried, deliberately, to keep the model as simple as possible
without losing key insights and relevance to empirical phenomena. One can add
more complexity as and when desired later on. What we are aiming for is the
equivalent of the ideal gas model or the Ising model for birds flocking.

We develop our teleodynamical model using Reynolds’s model as the start, but
our approach is not restricted to this example alone; it is applicable to
other models as well.

We consider a garud’s position in the frame of reference of the center of mass
of its neighborhood. We then apply the rule of cohesion and rule of separation
to formulate the model for utility.

The rule of cohesion requires the garuds to come together and hence an
$i^{\text{th}}$ garud’s utility increases as it has more neighbors, $n_{i}$.
However, the increased utility comes at the cost of congestion, the disutility
of congestion (corresponding to the rule of separation). The trade-off between
the two terms, the benefit-cost trade-off, results in an inverted-U profile,
which, following Venkatasubramanian [17], can be parameterized as,

$h_{r}^{(i)}=\alpha n_{i}-\beta n_{i}^{2}$ (5)

where $h_{r}^{(i)}$ is the position component of the utility for the
$i^{\text{th}}$ $\alpha,\beta>0$. Note that the positional dependence is
accounted for in the computation of $n_{i}$. Given a configuration of
$\\{\mathbf{r}_{i}\\}$, the neighborhood of the $i^{\text{th}}$ garud is
defined by the parameter $r_{0}$, where if $j^{\text{th}}$ garud is within
this radius then it is considered a neighbor.

This in turn identifies a direction of increased utility, given by,

$\frac{\partial h_{r}^{(i)}}{\partial\mathbf{r}_{i}}=\alpha\frac{\partial
n_{i}}{\partial\mathbf{r}_{i}}-2\beta n_{i}\frac{\partial
n_{i}}{\partial\mathbf{r}_{i}}$ (6)

$\dfrac{\partial n_{i}}{\partial\mathbf{r}_{i}}$ is dependent on the garuds in
the perimeter of the neighborhood of the reference garud $i$.

###  III.2 Garud’s Utility: Velocity Dependence

The utility of a garud is also dependent on the velocity of its neighbors in
that the garud attempts to match the orientation with its neighboring garuds.
This utility component ($h_{v}^{(i)}$) can be written as,

$h_{v}^{(i)}={\gamma}\sum_{j}n_{ij}\frac{\mathbf{v}_{i}}{|\mathbf{v}_{i}|}\cdot\frac{\mathbf{v}_{j}}{|\mathbf{v}_{j}|}$
(7)

The utility of the $i^{\text{th}}$ garud, then, depends on the orientation of
the other garuds in its neighborhood, i.e.,
$\mathbf{s}_{i}\cdot\mathbf{s}_{j}$ where $j$ is a neighbor of garud $i$. This
gives the alignment utility for the $i^{\text{th}}$ garud as,

$h_{v}^{(i)}=\gamma\sum_{j}n_{ij}\mathbf{s}_{i}\cdot\mathbf{s}_{j}$ (8)

where $n_{ij}$ shows if the $j^{\text{th}}$ garud is a neighbor of garud $i$,
i.e.,

$\displaystyle
n_{ij}=\begin{cases}1&|\mathbf{r}_{i}{(t)}-\mathbf{r}_{j}{(t)}|<r_{0}\\\
0&\text{otherwise}\\\ \end{cases}$ (9)

If each garud is perfectly aligned with its neighbors, this utility component
is maximal, whereas if they are oriented in the opposite direction this is
minimal. Therefore, the garuds prefer to be aligned. This gives the
$i^{\text{th}}$ garud an arbitrage opportunity to adjust its velocity vector
towards this direction to increase its utility. This opportunity for
increasing its utility generates a self-propelled force on the $i^{\text{th}}$
garud. If the $i^{\text{th}}$ garud is not aligned with its neighbors, this
direction of increased utility is given by,

$\frac{\partial
h_{v}^{(i)}}{\partial\mathbf{s}_{i}}=\gamma\sum_{j}n_{ij}\mathbf{s}_{j}$ (10)

###  III.3 Garud’s Effective Utility

There is one other utility component remaining to be considered. This is not
stated explicitly in the three rules of the boids. However, it is implied
because it is assumed that the boids have to be moving constantly.

So, as a garud incessantly moves and jockeys for better positions and
orientations, its ability to do so is hampered by the competition from other
garuds in the neighborhood that are also trying to do the same. As
Venkatasubramanian explains [17], this disutility of competition can be
modeled as $-\delta\ln n_{i}$.

This term, when integrated to obtain the potential $\phi(\mathbf{x})$, leads
to entropy in statistical mechanics. Thus, maximizing the potential
$\phi(\mathbf{x})$ is equivalent to maximizing entropy under certain
conditions. For more details, the reader is referred to Venkatasubramanian
[17, 20].

Now, by combining all these components, we arrive at the effective utility for
the $i^{\text{th}}$ garud given by

$\displaystyle h_{i}$ $\displaystyle=\alpha n_{i}-\beta
n_{i}^{2}+{\gamma}\sum_{j}n_{ij}\mathbf{s}_{i}\cdot\mathbf{s}_{j}-\delta\ln
n_{i}$ (11) $\displaystyle=\alpha n_{i}-\beta n_{i}^{2}+\gamma
n_{i}l_{i}-\delta\ln n_{i}$ (12)

where $l_{i}$ =
$\frac{1}{n_{i}}\sum_{j}n_{ij}\mathbf{s}_{i}\cdot\mathbf{s}_{j}$ is the
average alignment of agent $i$. Without any loss of generality, $\delta$ can
be assumed to be 1, and will be done so for the rest of this paper.

When $\alpha,\beta,\gamma=0$, the garuds don’t have any preferences and hence
fly randomly. This is what is captured by the remaining $-\ln n_{i}$ term,
which we call entropic restlessness.

Statistical teleodynamics, via potential game theory, proves that the self-
organizing dynamics of the garuds would eventually result in an arbitrage
equilibrium where the effective utilities of all the garuds are the same,
i.e., $h_{i}=h^{*}$.

In the next section, we discuss our simulation results that confirm this
prediction.

##  IV Results

Figure 1: Effective utility and its derivative as a function of the number of
neighbors $n_{i}$, for different values of alignment $l_{i}$
($\alpha,\beta,\gamma,\delta=0.5,0.005,0.25,1$). There are two locations where
the derivative of the effective utility is zero for different alignments.

For the simulation details, the reader is referred to the methodology section
VII below. If the garuds are flying randomly, without any rules of behavior,
then this base case corresponds to $\alpha,\beta,\gamma,=0;\delta=1$. This
result is discussed in Section S4.

For the other cases, the effective utility function in Eq. 12 is plotted in
Fig. 1 in terms of the number of neighbors of garud $i$ ($n_{i}$), for
different alignments ($l_{i}$) for a given set of $\alpha,\beta,\gamma$. We
see that there are two values of $n_{i}$ ($n_{-}$ and $n_{+}$) where the
gradient of utility, for a given value of alignment, is zero. These values are
determined by $\left(\dfrac{(\alpha+\gamma l_{i})\mp\sqrt{(\alpha+\gamma
l_{i})^{2}-8\beta\delta}}{4\beta}\right)$ (see, Supplementary Information,
Section 1).

In Fig. 1, at the lower value ($n_{-}$), any deviation in $n_{i}$ increases
the utility of the garud, and hence leads to an unstable point. However, for
the higher $n_{i}$ point ($n_{+}$), we see that any deviation reduces the
garud’s utility. Therefore, this leads to a stable point as any deviation
would bring a garud back to the higher utility state. Therefore, this is the
point a garud will try to reach to maximize its utility. For example, for the
red curve, this would correspond to the point where $n_{i}=73.6$.

However, despite this point’s stability, a garud will not be able to stay
there indefinitely as the other garuds in its neighborhood are constantly
changing their positions and orientations in their flights. Therefore, the
$i^{\text{th}}$ garud would be fluctuating around this point.

In Fig. 2 and 3, we show the simulation results of both the Reynolds’s Boids
model and our utility-driven Garud model dynamics. We show the results of the
Reynolds’s model for illustrative purposes only as our objective is not to
mimic the Reynolds’s model exactly. We just want to show that the utility-
driven model’s collective behavior is very similar to that of the Reynolds
model.

The results are shown for different parameter values of ($a,b,c$) and
($\alpha,\beta,\gamma$). From the simulations, we obtain a set of position and
velocity values $\\{\mathbf{r}_{i},\mathbf{v}_{i}\\}$ of each agent $i$ at
every time step. Once this is obtained, we extract the features $n_{i}$ and
$l_{i}$ for all the agents. This in turn, is used to compute the average
number of neighbors $\bar{n}$ and average alignment $\bar{l}$ for the entire
population for all time points (Fig. 2).

Figure 2 shows the snapshots at different time points of the evolution in
3D-space and the $n_{i}$-$l_{i}$ phase-space.

Figure 2: Trajectory of the agents and the corresponding phase portrait in
$n_{i}$-$l_{i}$ space for the average number of neighbors of each individual
agent during the course of the simulation, corresponding to (a) Reynolds’s
boids for $a=0.5,b=0.01,c=0.5$ and (b) Utility-driven garuds for
$\alpha=0.5,\beta=0.01,\gamma=0.25,\delta=1$ Figure 3: Trajectory of the
average of number of neighbors of each agent and average alignment of the
agents in the $n_{i}$, $l_{i}$ phase-space, and corresponding estimated
averages for the (a) Reynolds’s model (b) Utility-driven model ($\delta=1$)

While the exact dynamics, the exact configuration of the population, and the
time-scale of evolution of these two models cannot be the same (and that is
not the aim, either), we observe that the qualitative patterns of collective
behavior are very similar in both cases. In particular, we see that all agents
gravitate towards a certain region in the $n_{i}-l_{i}$ phase space for both
models. They both start at the lower far right point (where the average
alignment is near zero as the agents are all randomly oriented initially, and
closely packed) and evolve towards the upper center-right region in black. We
notice a qualitative match of both the trajectories.

Furthermore, we can also see a quantitative similarity between the two models
for specific parameters (Fig. 3). Fig. 3a shows the phase-space for the
Reynolds’s boids model, and Fig. 3b shows the same for the utility-driven
model.

In both models, we notice similar features of evolution towards the arbitrage
equilibrium states, starting from the lower right point at time $t=0$ to
ending in the colored regions, where the average number of neighbors and the
average alignment fall in similar corresponding regions. The plots in the
right show the average alignment and average number of neighbors of the agents
in the last 100 time steps.

While both models exhibit similar collective behaviors, it is not apparent,
however, from the three rules of the Reynolds’s model, that its dynamics would
result in an equilibrium state in the $n_{i}-l_{i}$ phase space.

This is different for the utility-driven model, however. Since its potential
game formulation predicts an arbitrage equilibrium outcome, it is clear right
from the beginning where in the phase space the system is going to end up in.
We can make a quantitative prediction about the average $n_{i}^{*}$ and
$h^{*}$ values at equilibrium.

This ability to predict the final outcome of the collective behavior of the
population, given the individual agent-level properties captured in the
utility function $h_{i}$, is an important defining feature and strength of the
statistical teleodynamics framework. An additional characteristic is the
ability to prove the stability of the final outcome.

We also ran the simulations for different time-step sizes of $\Delta
t=0.01,0.1,0.5$ in Eq. 3 to understand the dynamics of the evolution better.
Note in Figure 4, at the start ($t=0$), the utilities of all the garuds are
spread out, with many having negative utility values, and the average utility
($\bar{h}$) low.

But as the dynamics evolves, every garud tries to increase its utility by
maneuvering to a better neighborhood and better orientation, the distribution
becomes narrower, the average utility keeps increasing, and reaches a near-
maximum value ($\bar{h}=22.17\pm 2.90$ in Fig. 4a) and fluctuates around it.
Note that this is around the maximum theoretical value of about 23.8 (given by
Eq. 12), where the histogram peaks. This suggests that nearly all the garuds
have similar effective utilities asymptotically, approaching the maximum.
This, of course, is the arbitrage equilibrium outcome predicted by the theory
(see also, Supplementary information, Fig. S3).

The garuds do not converge exactly on $h^{*}$ but fluctuate around it because
of the stochastic dynamics. This is also seen in Table 1 where nearly the top
10 $\%$ of the garuds at a particular time-step are very close to the maximum
utility value. In fact, the top $50\%$ of the garuds have an average utility
of greater than 23.

This arbitrage equilibrium state is unique only if the potential
$\phi(\mathbf{x})$ is strictly concave [23]. For garuds, this is not the case
in general as the concavity would depend on $\alpha$, $\beta$, and $\gamma$
having some particular values. So, for the typical case where
$\phi(\mathbf{x})$ is not concave, there could be multiple equilibrium
configurations of the garuds. Thus, instead of an equilibrium point in the
$n_{i}-l_{i}$ phase space, one has an equilibrium region, in general. In other
words, invoking a terminology from chaos and nonlinear dynamics, there would
be a basin of attraction in the phase space where the garuds finally settle in
and fly around. This is what we see in Fig. 2 and 3 in the colored regions.

Figure 4: Histogram of the utility ($h_{i}$) of the garuds for $\alpha,\beta,\gamma,\delta=0.5,0.005,0.25,1$ corresponding to (a) $\Delta t=0.01$, (b) $\Delta t=0.1$ and (c) $\Delta t=0.5$. Dashed line shows the average utility at a particular time. Table 1: Utility of different percentiles of the garuds at the $1000^{\text{th}}$ time step corresponding to Figure 4 Population |  Time step size, $\Delta t$ | Average utility
---|---|---
Top $1\%$ | 0.01 | 23.79
| 0.1 | 23.80
| 0.5 | 23.76
Top $1-10\%$ | 0.01 | 23.77
| 0.1 | 23.78
| 0.5 | 23.68
Top $10-50\%$ | 0.01 | 23.60
| 0.1 | 23.51
| 0.5 | 23.35
Top $50-75\%$ | 0.01 | 22.83
| 0.1 | 22.64
| 0.5 | 22.58
Bottom $50\%$ | 0.01 | 18.59
| 0.1 | 15.59
| 0.5 | 17.91

###  IV.1 Stability of the Arbitrage Equilibrium

We can ascertain the stability of this equilibrium by performing a Lyapunov
stability analysis [17]. A Lyapunov function $V$ is a continuously
differentiable function that takes positive values everywhere except at the
equilibrium point (i.e., $V$ is positive definite), and decreases (or is
nonincreasing) along every trajectory traversed by the dynamical system
($\dot{V}$ is negative definite or negative semidefinite). A dynamical system
is locally stable at equilibrium if $\dot{V}$ is negative semidefinite and is
asymptotically stable if $\dot{V}$ is negative definite.

Following Venkatasubramanian [17], we identify our Lyapunov function
$V(n_{i})$

$\displaystyle V(n_{i})=\phi^{*}(n_{i}^{*})-\phi(n_{i})$ (13)

where $\phi^{*}$ is the potential at the Nash equilibrium (recall that
$\phi^{*}$ is at its maximum at NE) and $\phi(n_{i})$ is the potential at any
other state. Note that $V(n_{i})$ has the desirable properties we seek: (i)
$V(n^{*})$ = 0 at NE and $V(n_{i})$ $>$ 0 elsewhere, i.e., $V(n_{i})$ is
positive definite; (ii) since $\phi(n_{i})$ increases as it approaches the
maximum, $V(n_{i})$ decreases with time, and hence it is easy to see that
$\dot{V}$ is negative definite. Therefore, the arbitrage equilibrium is not
only stable but also asymptotically stable.

Our simulation results confirm this theoretical prediction (see Figure 5).
After the garuds population reached equilibrium, we disturbed the equilibrium
state by randomly changing the positions and/or velocities of the garuds. The
simulation is then continued from the new disturbed far-from-equilibrium
state. We conducted experiments with three kinds of disturbances:

  * •

Disturbance 1:Velocity disturbance, where each garud’s velocity is changed to
a random orientation and magnitude.

  * •

Disturbance 2: Position disturbance, where each garud’s position is randomly
changed.

  * •

Disturbance 3: Position and velocity disturbance, where both position and
velocity vectors are changed.

Figure 5: The trajectory of the garuds as a function of time and the
corresponding phase-portraits for
$\alpha,\beta,\gamma,\delta=0.5,0.005,0.25,1$ for stability analysis case
studies (a) Disturbance 1, (b) Disturbance 2, and (c) Disturbance 3. The
disturbances occur after equilibrium at time step 101.

As seen in Figure 5, after the 100${}^{\text{th}}$ time step when the
population had reached equilibrium, we introduced these disturbances. The
101${}^{\text{st}}$ time step shows in red color the new far-from equilibrium
states, where the average utility has dropped considerably.

In all cases, the population recovers quickly, typically in another 100 time
steps or so to reach the original equilibrium region (shown in green). The
figure shows the disturbance (red) and recovery (green) in both the 3-D space
and the phase space.

In Fig. 5a, as the velocities are randomized at the 101${}^{\text{st}}$ time
step, the alignment goes down to 0, but recovers to the original equilibrium
quickly. In Fig. 5b, as the new configuration corresponds to a similar value
of average number of neighbors as before, the disturbance is not that much.
Note that the drop in average utility is small. In Fig. 5c, we see that this
disturbance is huge, pushing the configuration close to the original random
state. But still, the population is able to recover to the arbitrage
equilibrium quickly.

This shows that the arbitrage equilibrium region is not only stable, but
asymptotically stable. That is, the garuds flocking configuration is resilient
and self-healing. Given the speed of the recovery, it could possibly be
exponentially stable but we have not proved this analytically here.

The asymptotic stability of this arbitrage equilibrium is similar to that of
the income-game dynamics as discussed by Venkatasubramanian [15, 17] using a
similar Lyapunov stability analysis.

##  V Conclusion

For three centuries we have known that there are constants of motion, such as
energy and momentum, for passive matter. Nevertheless, it comes as a surprise
to discover that the dynamics of active matter populations could also have an
invariant, namely, the effective utility. However, the role of invariance here
is different from its role in dynamics. The constants of motion such as energy
and momentum are conserved, but effective utility is not.

The role of this invariance is more like that of set point tracking and
disturbance rejection in feedback control systems. These are called the
regulation and servo problems, respectively, in control theory [26, 27]. The
system, i.e., the garuds population, adjusts itself dynamically and
continually, in a feedback control-like manner, to maintain its overall
effective utility.

It is important to emphasize, however, that this control action is
decentralized as opposed to the typical centralized control system in many
engineering applications. The agents individually self-organize, adapt, and
dynamically course-correct to offset the negative impact on their effective
utilities by other agents or other external sources of disturbance. The
population as a whole stochastically evolves towards the stable basin of
attraction in the phase space in a self-organized and distributed-control
fashion.

This is essentially Adam Smith’s invisible hand mechanism of economics. As
Smith observed[28], ”It is not from the benevolence of the butcher, the
brewer, or the baker, that we expect our dinner, but from their regard to
their own interest. We address ourselves, not to their humanity but to their
self-love, and never talk to them of our own necessities but of their
advantages.” Thus, every garud is pursuing its own self-interest, to increase
its own $h_{i}$, and a stable collective order spontaneously emerges via such
self-organization.

Invariants are quite rare in physics, rarer still in biology and economics.
That’s why it is exciting to see them as their presence usually signals
something deep, something fundamental. In physics, their existence has
revealed deeply fundamental symmetries of the cosmos, as Emmy Noether showed.
Therefore, it is important to understand the implications of our discovery in
sufficient depth and breadth.

We do realize, of course, that our bird-like agents are not real birds. Our
model and simulations are not real biological systems. Nevertheless, our
results suggest intriguing possibilities for real biological entities that
need to be explored carefully.

Therefore, the interesting and surprising result, seen both analytically and
in simulations, that the emergent arbitrage equilibrium is asymptotically
stable is an important one with potentially far-reaching consequences,
particularly in biological, ecological, and economic contexts. For example,
this could be an important mechanism of pattern formation and pattern
stability in biological systems. Populations of cells could self-organize,
under different spatial and temporal conditions and constraints, driven by
their incessant and instinctive hunt for better utilities, to settle into
various stable basins of attraction - i.e., into different types of stable
emergent order - to form stable organized structures. Their asymptotic
stability property bestows upon them the resilient self-healing feature found
so commonly in many biological systems. This process could be a core mechanism
behind the design, control, and optimization of stable biological systems via
self-organization.

This mechanism is applicable to different length and time scales, from
molecular to macroscopic to planetary scales. These results raise several
interesting questions about populations of biological active matter competing
with one another. For instance, consider all the different kinds of microbial
populations in the human body, or for that matter, in any living organism. Not
only are all the individual microbes competing with one other strategically
for resources to increase their effective utilities, to improve their survival
and growth fitness, different microbial populations are also competing with
one another at a higher-scale.

So, is there a hierarchy of arbitrage equilibria? That is, the microbes in a
given population at some arbitrage equilibrium among themselves (say, level-1
equilibrium), and such populations themselves are in equilibrium with one
another at a higher level (say, level-2), and so on. Is there a planet-scale
equilibrium? That is, are all the living species on our planet, along with the
environment, at some arbitrage equilibrium? Or are we evolving towards one,
just like the garuds population did in this study. What happens when this
equilibrium is upset by either internal disturbances (such as climate change)
or external shocks (such as asteroid impact)?

As one can see, our theory is not limited to just birds flocking. It is also
applicable to the self-organizing dynamics and evolution of a wide variety of
systems in physics, biology, sociology, and economics. As Venkatasubramanian
et al.[20] showed the emergence of the exponential energy (i.e., Boltzmann)
distribution for gas molecules can be modeled by the effective utility

$h_{i}=-\beta\ln E_{i}-\ln N_{i}.$ (14)

Similarly, they showed,[20] as examples of biological systems, bacterial
chemotaxis can be modeled by

$h_{i}=-\alpha c_{i}-\ln N_{i}$ (15)

and the emergence of ant craters by

$h_{i}=b-\frac{\omega r_{i}^{a}}{a}-\ln N_{i}.$ (16)

The same study showed how the Schelling game-like segregation dynamics in
sociology can be modeled by

$\displaystyle h_{i}=\eta N_{i}+\ln(H-N_{i})-\xi N_{i}^{2}-\ln N_{i}$ (17)

and the income game in economics by

$h_{i}=\alpha\ln S_{i}-\beta\left(\ln S_{i}\right)^{2}-\ln N_{i}.$ (18)

What we have is for ideal systems, similar to the ideal gas or Ising model in
thermodynamics. Just as real gases and liquids don’t behave exactly like their
ideal versions in statistical thermodynamics, we don’t expect real biological
systems (or economic or ecological systems) to behave like their ideal
counterparts in statistical teleodynamics. Nevertheless, the ideal versions
serve as useful starting and reference points as we develop more comprehensive
models of active matter systems. The next steps would involve examining and
learning how real-world biological systems behave compared to their ideal
versions. This would, of course, necessitate several modifications to the
ideal models.

We note, from equations 14-18, a certain pattern in the structure of the
effective utility functions in different domains. Thus, we see that the same
mathematical and conceptual framework is able to predict and explain the
emergence of spontaneous order via self-organization to reach arbitrage
equilibrium in dynamical systems in physics, biology, sociology, and
economics.

This kind of universality is particularly striking, prompting us to conclude
with a quote a from the inimitable Richard Feynman that seems apropos here:
”Nature uses only the longest threads to weave her patterns, so that each
small piece of her fabric reveals the organization of the entire tapestry.” It
appears that the emergence of spontaneous order via self-organizing stable
arbitrage equilibria is such a thread.

##  VI Acknowledgements

This work was supported in part by the Center for the Management of Systemic
Risk (CMSR), Columbia University. This manuscript was written when the
corresponding author (VV) was the Otto Monsted Distinguished Visiting
Professor at the Danish Technical University (DTU) as well as a resident of
Nyhavn 18, Copenhagen, as a guest of the Danish Central Bank in the summer of
2022. It is with great pleasure that VV acknowledges the generous support of
these institutions and the warm hospitality of his colleagues in the Chemical
Engineering Department at DTU.

##  VII Authorship contribution statement

Venkat Venkatasubramanian: Conceptualization, Formal analysis, Funding
acquisition, Methodology, Investigation, Supervision, Writing – original
draft, Writing –review & editing.

Abhishek Sivaram: Software, Formal analysis, Writing –review & editing, Data
curation, Validation.

##  VIII Methodology

We created a simulation of 1000 garuds in a periodic-box of dimensions
$20\times 20\times 20$, where each garud’s neighborhood is a sphere with
radius $r_{0}=3$. Each garud starts at a random location and orientation
inside a $10\times 10\times 10$ block. Speed of each garud is limited between
$0.5$ and $1$. The update algorithm works similar to the Reynolds’s garuds
update, except the force is driven by the numerical estimates of direction of
increased utility (additively based on position and velocity). An additional
noise is also added to the velocity update strategy similar to the Reynolds’s
model to capture the erroneous strategies of velocity update for each garud.
This is given by a noise parameter (0.01, unless specified) times the
magnitude of the velocity. The noise indicates that a garud does not make
perfect choices in updating its velocity.

## References

  * Toner  _et al._ [2005] J. Toner, Y. Tu, and S. Ramaswamy, Annals of Physics 318, 170 (2005).
  * Narayan  _et al._ [2007] V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 105 (2007).
  * Ramaswamy [2010] S. Ramaswamy, Annu. Rev. Condens. Matter Phys. 1, 323 (2010).
  * Marchetti  _et al._ [2013] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Reviews of Modern Physics 85, 1143 (2013).
  * Cates and Tailleur [2013] M. E. Cates and J. Tailleur, EPL (Europhysics Letters) 101, 20010 (2013).
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Gonnella  _et al._ [2015] G. Gonnella, D. Marenduzzo, A. Suma, and A. Tiribocchi, Comptes Rendus Physique 16, 316 (2015).
  * Berkowitz [2020] R. Berkowitz, Physics 13, s146 (2020).
  * O’Byrne and Tailleur [2020] J. O’Byrne and J. Tailleur, Physical Review Letters 125, 208003 (2020).
  * Reynolds [1987] C. W. Reynolds, in  _Proceedings of the 14th annual conference on Computer graphics and interactive techniques_ (1987) pp. 25–34.
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Physical review letters 75, 1226 (1995).
  * Ballerini  _et al._ [2008] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, _et al._ , Proceedings of the national academy of sciences 105, 1232 (2008).
  * Bialek  _et al._ [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Proceedings of the National Academy of Sciences 109, 4786 (2012).
  * Bialek  _et al._ [2014] W. Bialek, A. Cavagna, I. Giardina, T. Mora, O. Pohl, E. Silvestri, M. Viale, and A. M. Walczak, Proceedings of the National Academy of Sciences 111, 7212 (2014).
  * Venkatasubramanian  _et al._ [2015] V. Venkatasubramanian, Y. Luo, and J. Sethuraman, Physica A: Statistical Mechanics and its Applications 435, 120 (2015).
  * Venkatasubramanian [2017a] V. Venkatasubramanian, Langmuir 33, 11703 (2017a).
  * Venkatasubramanian [2017b] V. Venkatasubramanian,  _How Much Inequality Is Fair?: Mathematical Principles of a Moral, Optimal, and Stable Capitalist Society_ (Columbia University Press, 2017).
  * Venkatasubramanian and Luo [2019] V. Venkatasubramanian and Y. Luo, in  _Network Theory and Agent-Based Modeling in Economics and Finance_ (Springer, 2019) pp. 159–174.
  * Kanbur and Venkatasubramanian [2020] R. Kanbur and V. Venkatasubramanian, The European Physical Journal Special Topics 229, 1661 (2020).
  * Venkatasubramanian  _et al._ [2022] V. Venkatasubramanian, A. Sivaram, and L. Das, Computers & Chemical Engineering , 107887 (2022).
  * Grégoire and Chaté [2004] G. Grégoire and H. Chaté, Physical review letters 92, 025702 (2004).
  * Rosenthal [1973] R. W. Rosenthal, International Journal of Game Theory 2, 65 (1973).
  * Sandholm [2010] W. H. Sandholm,  _Population games and evolutionary dynamics_ (MIT press, 2010).
  * Easley  _et al._ [2010] D. Easley, J. Kleinberg, _et al._ ,  _Networks, crowds, and markets_ , Vol. 8 (Cambridge university press Cambridge, 2010).
  * Monderer and Shapley [1996] D. Monderer and L. S. Shapley, Games and economic behavior 14, 124 (1996).
  * Seborg  _et al._ [2016] D. E. Seborg, T. F. Edgar, D. A. Mellichamp, and F. J. Doyle III,  _Process dynamics and control_ (John Wiley & Sons, 2016).
  * Åström and Murray [2021] K. J. Åström and R. M. Murray,  _Feedback systems: an introduction for scientists and engineers_ (Princeton university press, 2021).
  * Smith [1776] A. Smith,  _An Inquiry into the Nature and Causes of the Wealth of Nations_ , Vol. 1 (W. Strahan and T. Cadell, London, 1776).