#  Mapping atomic orbitals with the transmission electron microscope: Images
of defective graphene predicted from first-principles theory

Lorenzo Pardini  Physics Department and IRIS Adlershof, Humboldt-Universität
zu Berlin, 12489 Berlin, Germany  Stefan Löffler  Institute of Solid State
Physics, Vienna University of Technology, Austria  University Service Centre
for Transmission Electron Microscopy, Vienna University of Technology, Austria
Department of Materials Science and Engineering, McMaster University,
Hamilton, Ontario, Canada  Giulio Biddau  Physics Department and IRIS
Adlershof, Humboldt-Universität zu Berlin, 12489 Berlin, Germany  Ralf Hambach
Electron Microscopy Group of Materials Science, University of Ulm, 89081 Ulm,
Germany  Fraunhofer Institute for Applied Optics and precision engineering
IOF, Jena, Germany  Ute Kaiser  Electron Microscopy Group of Materials
Science, University of Ulm, 89081 Ulm, Germany  Claudia Draxl  Physics
Department and IRIS Adlershof, Humboldt-Universität zu Berlin, 12489 Berlin,
Germany  European Theoretical Spectroscopy Facility (ETSF)  Peter
Schattschneider  Institute of Solid State Physics, Vienna University of
Technology, Austria  University Service Centre for Transmission Electron
Microscopy, Vienna University of Technology, Austria

###### Abstract

Transmission electron microscopy has been a promising candidate for mapping
atomic orbitals for a long time. Here, we explore its capabilities by a first
principles approach. For the example of defected graphene, exhibiting either
an isolated vacancy or a substitutional nitrogen atom, we show that three
different kinds of images are to be expected, depending on the orbital
character. To judge the feasibility of visualizing orbitals in a real
microscope, the effect of the optics’ aberrations is simulated. We demonstrate
that, by making use of energy-filtering, it should indeed be possible to map
atomic orbitals in a state-of-the-art transmission electron microscope.

###### pacs:

07.05.Tp

The possibility to see atomic orbitals has always attracted great scientific
interest. At the same time, however, the real meaning of ”measuring orbitals”
has been a subject that scientists have long and much dwelt upon (see, e.g.,
Schwarz (2006) and references therein). In the past, significant efforts have
been devoted to the development of experimental approaches and theoretical
models that allow for orbital reconstruction from experimental data Schwarz
(2006). Based on the generation of higher harmonics by femtoseconds laser
pulses, a tomographic reconstruction of the highest occupied molecular orbital
(HOMO) for simple diatomic molecules in the gas phase was proposed Itatani
_et al._ (2004). Direct imaging of the HOMO and the lowest-occupied molecular
orbital (LUMO) of pentacene on a metallic substrate was theoretically
predicted and experimentally verified with scanning-tunnelling microscopy
(STM) Repp  _et al._ (2005). More recently, real-space reconstruction of
molecular orbitals from angle-resolved photoemission data has been
demonstrated Puschnig  _et al._ (2009). This method has been subsequently
further developed to retrieve both the spatial distribution Puschnig  _et al._
(2011) and the phase of electron wavefunctions of pentacene and
perylene-3,4,9,10-tetracarboxylic dianhydride (PTCDA) adsorbed on silver
Lüftner  _et al._ (2013).

The reconstruction of charge densities and chemical bonds using transmission
electron microscopy (TEM) has been considered Zuo  _et al._ (1999); Meyer  _et
al._ (2011); Haruta  _et al._ (2013); Neish  _et al._ (2013); Oxley  _et al._
(2014), but only recently the possibility of probing selected transitions to
specific unoccupied orbitals by using energy-filtered TEM (EFTEM) was
demonstrated theoretically. A first example for the capability of this
approach was provided with the oxygen K-edge of rutile TiO2 Löffler  _et al._
(2013). However the interpretation of experimental TEM images for systems like
rutile would be complicated because of the multiple elastic scattering of
electrons that occurs in thick samples.

In this Letter, we suggest defective graphene Geim and Novoselov (2007);
Lehtinen  _et al._ (2004); Sammalkorpi  _et al._ (2004); Carlsson and
Scheffler (2006); Yazyev and Helm (2007) as the prototypical two dimensional
(2D) material to demonstrate the possibility of mapping atomic orbitals using
EFTEM . We break the ideal $sp^{2}$ hybridization by introducing two different
kinds of defects, namely a single isolated vacancy and a substitutional
nitrogen atom. This lifts the degeneracy of the $p$-states, inducing strong
modifications to the electronic properties compared to the pristine lattice
Sammalkorpi  _et al._ (2004); Novoselov  _et al._ (2005); Son  _et al._
(2006); Heersche  _et al._ (2007); Lee  _et al._ (2005); Hou  _et al._ (2013);
Hansson  _et al._ (2000); Ewels  _et al._ (2002). By selecting certain
scattering angles, dipole-allowed transitions dominate the electron energy
loss spectroscopy (EELS) signal Jorissen (2007). A single-particle description
can be safely adopted, since many-body effects do not play a major role in the
excitation process. Overall, TEM images of these systems can be interpreted in
terms of bare $s-p$ transitions.

In an EFTEM experiment, an incoming beam of high-energy electrons (of the
order of 100 keV) is directed to the target where it scatters at the atoms
either elastically or inelastically. The outgoing electron beam is detected
and analyzed. State of the art image simulations generally only include
elastic scattering of the electrons using the multi-slice approach Kirkland
(2010). In the case of EFTEM for a thin sample, the influence of elastic
scattering becomes negligible, and inelastic scattering gives the dominant
contribution to the formation of the images.

The key quantity to describe the inelastic scattering of electrons, which is
probed by EELS, is the mixed dynamic form factor (MDFF). It can be interpreted
as a weighted sum of transition matrix elements between initial and final
states $\phi_{i}$ and $\phi_{j}$ of the target electron Schattschneider  _et
al._ (2000); Nelhiebel  _et al._ (1999); Kohl and Rose (1985):

$S\left(\bf{q},\bf{q^{\prime}};\it{E}\right)=\sum_{i,j}\left\langle\phi_{i}\left|e^{i\bf{q}\cdot\bf{r}}\right|\phi_{j}\right\rangle\left\langle\phi_{j}\left|e^{-i\bf{q^{\prime}}\cdot\bf{r}}\right|\phi_{i}\right\rangle\delta\left(E_{j}-E_{i}-E\right)$
(1)

with energies $E_{i}$ and $E_{j}$. $E$ is the energy-loss of the fast electron
of the incident beam, $\bf{q}$ and $\bf{q}^{\prime}$ are the wavevectors of
the perturbing and induced density fluctuations, respectively. If many-body
effects can be neglected, this picture can be simplified for dipole-allowed
transitions. In this case, using the spherical harmonics as basis for the
target states and referring to transitions originating from a single state (as
in $s-p$ excitations), the MDFF is Löffler  _et al._ (2013)

$S\left(\bf{q},\bf{q^{\prime}};\it{E}\right)\propto\sum_{\mu
LM,\mu^{\prime}L^{\prime}M^{\prime}}{\left\langle
j_{\lambda=1}\left(q\right)\right\rangle_{LE}\left\langle
j_{\lambda=1}\left(q^{\prime}\right)\right\rangle_{L^{\prime}E}}Y^{\mu}_{\lambda=1}\left(\bf{q}\right)^{*}Y^{\mu^{\prime}}_{\lambda=1}\left(\bf{q^{\prime}}\right)\Xi_{\lambda=1\mu
LM,\lambda=1\mu^{\prime}L^{\prime}M^{\prime}}\left(E\right)$ (2)

where $Y^{\mu}_{\lambda=1}\left(\bf{q}\right)$ are spherical harmonics,
$\left\langle j_{\lambda=1}\left(q\right)\right\rangle_{LE}$ is an integral of
the spherical Bessel function $j_{\lambda=1}\left(q\right)$ weighted over the
initial and final radial wavefunctions. $L$ and $M$ indicate the azimuthal and
magnetic quantum number of the final state of the target electron, and
$\lambda$ and $\mu$ are the angular momenta transferred during the transition.
$\Xi_{\lambda\mu LM,\lambda^{\prime}\mu^{\prime}LM}\left(E\right)$ is a
quantity that describes crystal-field effects and is proportional to the
cross-density of states (XDOS)

$\sum_{n\bf{k}}D^{n\bf{k}}_{LM}\left(D^{n\bf{k}}_{L^{\prime}M^{\prime}}\right)^{*}\delta\left(E_{n}\left(\bf{k}\right)-E\right)$
(3)

where $D^{n\bf{k}}_{LM}$ is the angular part of the final wave function, $n$
is the band index, and $\bf{k}$ is a k-point in the first Brillouin zone.
Compared to the density of states (DOS), the XDOS includes also non-diagonal
terms connecting states with different angular momenta. As $\Xi_{\lambda=1\mu
LM,\lambda=1\mu^{\prime}LM}\left(E\right)$ is a hermitian matrix Löffler  _et
al._ (2013), the MDFF can be diagonalized. Therefore, assuming that the
target’s final states of the $s-p$ excitation are not degenerate, the
transition matrix elements reflect the azimuthal shape of the final single-
particle states and can thus be separated by using energy-filtering.

Ground-state calculations are performed using density-functional theory and
the full-potential augmented planewave plus local-orbital method, as
implemented in exciting Gulans  _et al._ (2014). Introducing a vacancy or a
substitutional atom, a 5$\times$5 supercell is set up, hosting 49 and 50
atoms, respectively (Fig. 1 and 2 in the Supplemental Material sup ). The
space group and thus the number of inequivalent carbon atoms (13) is the same
in both cases. We adopt a lattice parameter of $a$=4.648 bohr, corresponding
to a bond length of 2.683 bohr, while the cell size perpendicular to the
graphene plane is set to c=37.794 bohr in order to prevent interactions
between the periodically repeated layers. Exchange-correlation effects are
treated by the PBE functional Perdew  _et al._ (1996). The Brillouin zone is
sampled with an 8$\times$8$\times$1 k-point grid. The structures are relaxed
down to a residual force lower than 0.0005 Ha/bohr acting on each atom.
Interatomic distances between atoms of the relaxed structures, up to the
seventh nearest neighbor, are given in Table 1. Upon relaxation, the atoms
surrounding the vacancy move slightly away from it, thus shortening the bond
lengths with the next nearest neighbors, $d_{1-2}$, compared to the
unperturbed system. The effect of the vacancy extends up to the fourth
neighbours, whereas it is almost negligible for more distant atoms (more
information about the relaxed structures can be found in the Supplemental
Material sup ). In the case of nitrogen doping, the substitutional atom does
not strongly influence the atomic configuration of the system. This happens
because the nitrogen-carbon bond length is just slightly shortened with
respect to the carbon-carbon bond length in pristine graphene. For all the
systems, we have investigated dipole-allowed transitions at the K-edge of
carbon, assuming an incoming electron beam perpendicular to the graphene
plane.

Table 1:  Bond lengths between atoms up to the seventh nearest neighbor, $d_{0-1}$, $d_{1-2}$, $d_{2-3}$, $d_{2-4}$, $d_{4-6}$, and $d_{6-7}$, for graphene doped with nitrogen (top row) and with a vacancy (bottom row). $\Delta d$ are the relative deviations from those of pristine graphene. $0$ indicates the defect site. ( System |  | $d_{0-1}$ | $d_{1-2}$ | $d_{2-3}$ | $d_{2-4}$ | $d_{4-6}$ | $d_{4-7}$
---|---|---|---|---|---|---|---
N-doped |  $d$ [bohr] | 2.675 | 2.675 | 2.683 | 2.689 | 2.683 | 2.692
| $\Delta d$ | -0.3% | -0.3% | - | +0.2% | - | +0.3%
Vacancy |  $d$ [bohr] | 2.689 | 2.665 | 2.678 | 2.712 | 2.676 | 2.687
| $\Delta d$ | +0.2% | -0.7% | -0.2% | +1.1% | -0.3% | +0.1%

In Fig. 1, the projected density of states (PDOS) of pristine graphene (upper
panel) and of the first nearest-neighbor atom for nitrogen-doped graphene
(middle panel) and graphene with a single vacancy (bottom), respectively, is
plotted for empty states up to 12 eV above the Fermi energy. Here, $x$, $y$,
and $z$ represent the local Cartesian coordinates at the individual atomic
sites as determined by the point-group symmetry. In particular, $z$ is the
axis perpendicular to the graphene, i.e., ($x$, $y$) plane. All the other
atoms of the defective systems exhibit a PDOS with very similar character as
in pristine graphene, besides the second and third nearest neighbors which are
slightly affected by the defect sup .

Figure 1:  Local projected density of states (PDOS) of carbon in pristine
graphene (upper panel), and of the first nearest-neighbor atom for nitrogen-
doped graphene (middle panel) and graphene with a single vacancy (bottom).
$p_{x}$, $p_{y}$, and $p_{z}$ states are indicated by the red, green, and
black lines, respectively. Note that, in the lower panel, the peak close to
the Fermi energy exceeds the boundary of the box, with a maximum at about 2
states/eV. In the case of pristine graphene, the red line indicates the sum of
$p_{x}$ and $p_{y}$. The three colored areas labeled as A, B, and C indicate
the energy ranges for which the TEM images have been calculated.

In pristine graphene, antibonding $\pi^{*}$ and $\sigma^{*}$ states are
clearly recognizable at about 2 eV and 9 eV, respectively. As already reported
in literature Amara  _et al._ (2007); Hou  _et al._ (2012); El-Barbary  _et
al._ (2003); Lambin  _et al._ (2012), the introduction of a vacancy or a
substitutional nitrogen has a significant influence on the electronic
structure. A consequence of the doping atom is lifting the degeneracy of
$p_{x}$ and $p_{y}$ that is significant for the first nearest neighbors
(middle panel in Fig. 1). This effect is particularly evidenced by the
appearance of bands at about 5 eV, which exhibit $p_{y}$ character. Here,
three different regions can be easily identified: a) From 0 eV to 4 eV, the
bands have only $p_{z}$ character; energy ranges that present such DOS
character will be referred to as $T_{p_{z}}$. b) For energies higher than 6
eV, there are contributions from $p_{x}$, $p_{y}$, and $p_{z}$. The only
difference to ideal graphene is the lifted degeneracy of $p_{x}$ and $p_{y}$.
This defines a new kind of region, named $T_{p_{x,y,z}}$. c) Between 4 eV and
6 eV, the $p_{x}$ character of the first nearest neighbor is much less
pronounced than that of $p_{y}$, while all the other atoms have only $p_{z}$
character; this region will be referred to as $T_{p_{y,z}}$. In the case of
graphene with a vacancy, the same kinds of regions can be identified, but
corresponding to different energy ranges. Here, the $T_{p_{z}}$ type is found
between 0.5 eV and 7 eV; $T_{p_{x,y,z}}$ encompasses energies above 7 eV;
$T_{p_{y,z}}$ is a small energy window, just few tenths of eV close to the
Fermi energy. We find similar kinds of DOS characters also for damaged
nitrogen-doped graphene, i.e., graphene with a substitutional nitrogen and a
vacancy located near it; such defects have beed recently reported in TEM
measurements of nitrogen-doped graphene Susi  _et al._ (2012). Details of this
calculation and the corresponding simulated TEM images can be found in the
Supplemental Material sup .

To investigate the impact of the local electronic structure (PDOS) on the
EFTEM images, we first consider the ideal case of a perfect microscope with an
acceleration voltage of 300 keV. In this case, the recorded images correspond
to the intensity of the exit wavefunction in the multislice simulation (Cowley
and Moodie, 1957). The finite resolution of the spectrometer is taken into
account by simulating images every 0.05 eV in 2 eV-broad energy ranges
(regions A, B, and C in Fig. 2) and then summing them up to get the final
images. Each image is shown in contrast-optimized grayscale.

Figure 2:  Simulated real-space intensity of the electron’s exit wave function
after propagation of an incident planewave through a graphene layer in
presence of a nitrogen substitutional atom (upper panels) and a vacancy
(bottom panels). The colored lines indicate the energy regions highlighted in
Fig. 1.

First, we analyze graphene with nitrogen doping (Fig. 2, upper panels). Here,
in the region close to the Fermi level (A in Fig. 1) there are only
contributions from $p_{z}$ orbitals. The image is then formed by disk-like
features where their arrangement clearly visualizes the missing atoms (upper
left panel in Fig. 2). At an energy loss between 4 and 6 eV above the carbon
K-edge (region B in Fig. 1), there is a $T_{p_{y,z}}$-like region. We expect
to see contributions from $p_{y}$ of the atom closest to the nitrogen, but no
(or very little) signal coming from the other atoms. This happens because
$p_{y}$ lies on a plane perpendicular to the incoming electron beam and its
magnitude is more intense than the one of $p_{z}$; thus its contribution to
the final signal overcomes the one from $p_{z}$ states. Consequently, only the
$p_{y}$ orbitals of the three atoms surrounding the nitrogen are visible,
which are pointing towards the defect, as imposed by the local $D_{3h}$
symmetry (upper middle panel in Fig. 2). At an excitation energy between 8 and
10 eV above the K-edge (region C in Fig. 1), instead, there are contributions
from all the $p$ states. Since, however, $p_{x}$ and $p_{y}$ lie in a plane
perpendicular to the beam axis, their contribution to the final signal
dominates over the one from $p_{z}$ states. As a consequence, the image is
composed of ring-like features, stemming solely from $p_{x}$ and $p_{y}$
states, arranged in hexagons (upper right panel in Fig. 2). Due to symmetry
breaking, the intensity is not uniform, neither along a ring (since $p_{x}$
and $p_{y}$ states are non-degenerate) nor among different rings (due to non-
equivalent atomic sites).

The corresponding images for the system with a vacancy (bottom panels in Fig.
2) appear nearly identical to the ones above, but at different energy ranges.
This can be understood by comparing the PDOS of the two systems. Between 0 and
2 eV, for instance, we have a $T_{p_{z}}$-like region in the case of nitrogen-
doped graphene, and both $T_{p_{y,z}}$ and $T_{p_{z}}$ in the case of graphene
with a vacancy. Due to the similarity of the two systems, we will, in the
following, focus on doped graphene and show the corresponding analysis for
graphene with a single vacancy in the Supplemental Material sup .

In order to predict the outcome of real experiments, we now visualize the
effect of the optics’ aberrations and of a more realistic acceleration voltage
on these images. We have simulated an electron-beam acceleration voltage of 80
keV and the operating parameters of two different kinds of microscopes, the
FEI Tecnai G2 F20 and FEI Titan G2 60-300. The former has a spherical
aberration $C_{S}$ = 1.2 mm, corresponding to a an extended Scherzer defocus
of 849 Å, while the latter is a last-generation aberration-corrected
microscope, i.e., exhibiting no spherical and chromatic aberrations. In view
of that, chromatic aberrations are not included in the calculations. The
images corresponding to the energy regions A, B, and C are shown in Fig. 3.

Figure 3:  Simulated TEM image of nitrogen-doped graphene. The coloured lines
indicate the energy regions highlighted in Fig. 1. 80 keV incident beam energy
and lenses as in a Titan (upper panels) and a Tecnai microscope (bottom
panels) were assumed.

Because of the lower resolution of the Tecnai microscope, all the features are
blurred (lower panels) compared to the ideal images. Therefore, neither the
atomic positions, nor the orbital shapes can be retrieved from them. On the
other hand, images simulated by taking into account the aberration-corrected
optics of the Titan microscope are very sharp and let us identify all the
features already observed for the idealized situation previously described.
This can be easily seen, comparing the upper panels of Figs. 3 and 2. In
particular, at an energy loss of 5 eV (region B), the $p_{y}$ orbitals are
visible, as in the ideal images. This clearly demonstrates the potential
ability of aberration-corrected microscopes to visualize atomic orbitals with
EFTEM, especially in a system like graphene. This conclusion also holds when
considering noise caused by the finite electron dose (see Supplemental
Material sup  for corresponding images).

In summary, we have predicted the possibility of performing orbital mapping in
low-dimension systems using EFTEM and we have demonstrated it with the
prototypical example of defective graphene. In particular, we have shown that,
as far as the optics is concerned, reasonable image resolution may already
nowadays be experimentally achievable with last generation aberration-
corrected microscopes like a FEI Titan G2 60-300 and even more with improved
instruments of the next generation. However, additional work is necessary to
reduce artifacts such as noise, drift, instabilities and damage. The inelastic
cross section for the carbon K-edge ionisation is about a factor of 10 smaller
than the elastic scattering cross section on a carbon atom Reimer (1993). The
intensity collected within an energy window of 2 eV as in Fig. 2 is $<$5% of
the total K-edge intensity. So, order-of-magnitude-wise, in order to obtain
the same SNR as in elastic imaging, we need at least 200 times more incident
dose which means a dwell time of the order of several minutes for last
generation TEMs in EFTEM mode. There is no fundamental law that would forbid
such an experiment with today’s equipment; however, it is hampered by drift
(which must be well below the interatomic distance during the exposure time),
instabilities, and radiation damage. A new route to circumvent radiation
damage based on an EFTEM low-dose technique was proposed recently Meyer  _et
al._ (2014). This may solve the problem in future.

We have identified three different kinds of images that are expected to be
acquired in an EFTEM experiment, depending on the character of the DOS: When
only $p_{z}$ states are present in the electronic structure, the corresponding
images are composed of disk-like features. When the DOS is characterized by
contributions from all $p$ states, ring-like features are seen that, however,
only originate from a convolution of $p_{x}$ and $p_{y}$ states, while the
$p_{z}$ character is not visible. When the $p_{y}$ character strongly exceeds
the one of $p_{x}$, only a single orbital is recorded. We expect this work to
trigger new experiments on defective graphene and similar systems.

##  I Acknowledgments

Financial support by the Austrian Science Fund (projects I543-N20 and
J3732-N27) and the German Research Foundation within the DACH framework is
acknowledged.

##  II Bibliography styles

## References

  * Schwarz (2006) W. H. E. Schwarz, Angewandte Chemie International Edition 45, 1508 (2006).
  * Itatani  _et al._ (2004) L. Itatani, J. Levesque, D. Zeidler, H. Niikura, H. Pepin, J. C. Kieffer, P. B. Corkum, and D. M. Villeneuve, Nature 432, 867 (2004).
  * Repp  _et al._ (2005) J. Repp, G. Meyer, S. M. Stojković, A. Gourdon, and C. Joachim, Phys. Rev. Lett. 94, 026803 (2005).
  * Puschnig  _et al._ (2009) P. Puschnig, S. Berkebile, A. J. Fleming, G. Koller, K. Emtsev, T. Seyller, J. D. Riley, C. Ambrosch-Draxl, F. P. Netzer, and M. G. Ramsey, Science 326, 702 (2009).
  * Puschnig  _et al._ (2011) P. Puschnig, E.-M. Reinisch, T. Ules, G. Koller, S. Soubatch, M. Ostler, L. Romaner, F. S. Tautz, C. Ambrosch-Draxl, and M. G. Ramsey, Phys. Rev. B 84, 235427 (2011).
  * Lüftner  _et al._ (2013) D. Lüftner, T. Ules, E. M. Reinisch, G. Koller, S. Soubatch, F. S. Tautz, M. G. Ramsey, and P. Puschnig, Proceedings of the National Academy of Sciences (2013).
  * Zuo  _et al._ (1999) J. M. Zuo, M. Kim, M. O’Keeffe, and J. C. H. Spence, Nature 401, 49–52 (1999), 10.1038/43403.
  * Meyer  _et al._ (2011) J. C. Meyer, S. Kurasch, H. J. Park, V. Skakalova, D. Künzel, A. Groß, A. Chuvilin, G. Algara-Siller, S. Roth, T. Iwasaki, U. Starke, J. H. Smet, and U. Kaiser, Nature Materials 10, 209 (2011).
  * Haruta  _et al._ (2013) M. Haruta, T. Nagai, N. R. Lugg, M. J. Neish, M. Nagao, K. Kurashima, L. J. Allen, T. Mizoguchi, and K. Kimoto, Journal of Applied Physics 114, 083712 (2013).
  * Neish  _et al._ (2013) M. J. Neish, N. R. Lugg, S. D. Findlay, M. Haruta, K. Kimoto, and L. J. Allen, Phys. Rev. B 88, 115120 (2013).
  * Oxley  _et al._ (2014) M. P. Oxley, M. D. Kapetanakis, M. P. Prange, M. Varela, S. J. Pennycook, and S. T. Pantelides, Microscopy and Microanalysis 20, 784 (2014).
  * Löffler  _et al._ (2013) S. Löffler, V. Motsch, and P. Schattschneider, Ultramicroscopy 131, 39 (2013).
  * Geim and Novoselov (2007) A. K. Geim and K. S. Novoselov, Nat Mater 6, 183 (2007).
  * Lehtinen  _et al._ (2004) P. O. Lehtinen, A. S. Foster, Y. Ma, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. Lett. 93, 187202 (2004).
  * Sammalkorpi  _et al._ (2004) M. Sammalkorpi, A. Krasheninnikov, A. Kuronen, K. Nordlund, and K. Kaski, Phys. Rev. B 70, 245416 (2004).
  * Carlsson and Scheffler (2006) J. M. Carlsson and M. Scheffler, Phys. Rev. Lett. 96, 046806 (2006).
  * Yazyev and Helm (2007) O. V. Yazyev and L. Helm, Phys. Rev. B 75, 125408 (2007).
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * Son  _et al._ (2006) Y.-W. Son, M. L. Cohen, and S. G. Louie, Nature 444, 347 (2006).
  * Heersche  _et al._ (2007) H. B. Heersche, P. Jarillo-Herrero, J. B. Oostinga, L. M. K. Vandersypen, and A. F. Morpurgo, Nature 446, 56 (2007).
  * Lee  _et al._ (2005) G.-D. Lee, C. Z. Wang, E. Yoon, N.-M. Hwang, D.-Y. Kim, and K. M. Ho, Phys. Rev. Lett. 95, 205501 (2005).
  * Hou  _et al._ (2013) Z. Hou, X. Wang, T. Ikeda, K. Terakura, M. Oshima, and M.-a. Kakimoto, Phys. Rev. B 87, 165401 (2013).
  * Hansson  _et al._ (2000) A. Hansson, M. Paulsson, and S. Stafström, Phys. Rev. B 62, 7639 (2000).
  * Ewels  _et al._ (2002) C. Ewels, M. Heggie, and P. Briddon, Chemical Physics Letters 351, 178 (2002).
  * Jorissen (2007) K. Jorissen,  _The Ab Initio Calculation of Relativistic Electron Energy Loss Spectra_ , Ph.D. thesis, Universiteit Antwerpen (2007).
  * Kirkland (2010) E. J. Kirkland,  _Advanced computing in electron microscopy_ (Springer US, 2010).
  * Schattschneider  _et al._ (2000) P. Schattschneider, M. Nelhiebel, H. Souchay, and B. Jouffrey, Micron 31, 333 (2000).
  * Nelhiebel  _et al._ (1999) M. Nelhiebel, N. Luchier, P. Schorsch, P. Schattschneider, and B. Jouffrey, Philosophical Magazine Part B 79, 941 (1999), http://dx.doi.org/10.1080/13642819908214851 .
  * Kohl and Rose (1985) H. Kohl and H. Rose,  _Theory of image formation by inelastically scattered electrons in the electron microscopy_ , Advences in electronics and electron physics (Academic press, 1985).
  * Gulans  _et al._ (2014) A. Gulans, S. Kontur, C. Meisenbichler, D. Nabok, P. Pavone, S. Rigamonti, S. Sagmeister, U. Werner, and C. Draxl, J. Phys.: Condens. Matter 26, 363202 (2014).
  * (31) See Supplemental Material at [http link] for more information about the DOS of the second nearest neighbors, the DOS of the atom most distant from the defect, the simulated EFTEM images of graphene with a single vacancy with lenses as in Tecnai and Titan microscopes, details of the calculation of damaged nitrogen-doped graphene and the corresponding images, and the simulated images with finite electron dose.
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Amara  _et al._ (2007) H. Amara, S. Latil, V. Meunier, P. Lambin, and J.-C. Charlier, Phys. Rev. B 76, 115423 (2007).
  * Hou  _et al._ (2012) Z. Hou, X. Wang, T. Ikeda, K. Terakura, M. Oshima, M.-a. Kakimoto, and S. Miyata, Phys. Rev. B 85, 165439 (2012).
  * El-Barbary  _et al._ (2003) A. A. El-Barbary, R. H. Telling, C. P. Ewels, M. I. Heggie, and P. R. Briddon, Phys. Rev. B 68, 144107 (2003).
  * Lambin  _et al._ (2012) P. Lambin, H. Amara, F. Ducastelle, and L. Henrard, Phys. Rev. B 86, 045448 (2012).
  * Susi  _et al._ (2012) T. Susi, J. Kotakoski, R. Arenal, S. Kurasch, H. Jiang, V. Skakalova, O. Stephan, A. V. Krasheninnikov, E. I. Kauppinen, U. Kaiser, and J. C. Meyer, ACS Nano 6, 8837 (2012).
  * Cowley and Moodie (1957) J. M. Cowley and A. F. Moodie, Acta Crystallographica 10, 609 (1957).
  * Reimer (1993) L. Reimer,  _Transmission Electron Microscopy_ (Springer US, 1993).
  * Meyer  _et al._ (2014) J. Meyer, J. Kotakoski, and C. Mangler, Ultramicroscopy 145, 13 (2014).


* Andersson and Kokkotas [2001] N. Andersson and K. D. Kokkotas, The r-mode instability in rotating neutron stars, International Journal of Modern Physics D 10, 381 (2001), https://doi.org/10.1142/S0218271801001062 .
  * Özel and Freire [2016] F. Özel and P. Freire, Masses, radii, and the equation of state of neutron stars, Annual Review of Astronomy and Astrophysics 54, 401 (2016), https://doi.org/10.1146/annurev-astro-081915-023322 .
  * Chatziioannou [2020] K. Chatziioannou, Neutron-star tidal deformability and equation-of-state constraints, Gen. Relativ. Gravit. 52, 109 (2020).
  * Moustakidis [2015] C. C. Moustakidis, Effects of the nuclear equation of state on the $r$-mode instability and evolution of neutron stars, Phys. Rev. C 91, 035804 (2015).
  * Akmal  _et al._ [1998] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Equation of state of nucleon matter and neutron star structure, Phys. Rev. C 58, 1804 (1998).
  * Margaritis  _et al._ [2020] C. Margaritis, P. S. Koliogiannis, and C. C. Moustakidis, Speed of sound constraints on maximally rotating neutron stars, Phys. Rev. D 101, 043023 (2020).
  * Rhoades and Ruffini [1974] C. E. Rhoades and R. Ruffini, Maximum mass of a neutron star, Phys. Rev. Lett. 32, 324 (1974).
  * Kalogera and Baym [1996] V. Kalogera and G. Baym, The maximum mass of a neutron star, Astrophys. J. 470, L61 (1996).
  * Koranda  _et al._ [1997] S. Koranda, N. Stergioulas, and J. L. Friedman, Upper Limits Set by Causality on the Rotation and Mass of Uniformly Rotating Relativistic Stars, Astrophys. J. 488, 799 (1997).
  * Chamel  _et al._ [2013] N. Chamel, P. Haensel, J. Zdunik, and A. Fantina, On the Maximum Mass of Neutron Stars, Int. J. Mod. Phys. E 22, 1330018 (2013).
  * Podkowka  _et al._ [2018] D. M. Podkowka, R. F. P. Mendes, and E. Poisson, Trace of the energy-momentum tensor and macroscopic properties of neutron stars, Phys. Rev. D 98, 064057 (2018).
  * Xia  _et al._ [2021] C. Xia, Z. Zhu, X. Zhou, and A. Li, Sound velocity in dense stellar matter with strangeness and compact stars, Chin. Phys. C 45, 055104 (2021).
  * Tews  _et al._ [2018a] I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Constraining the speed of sound inside neutron stars with chiral effective field theory interactions and observations, Astrophys. J. 860, 149 (2018a).
  * Feynman  _et al._ [1949] R. Feynman, N. Metropolis, and E. Teller, Equations of State of Elements Based on the Generalized Fermi-Thomas Theory, Phys. Rev. 75, 1561 (1949).
  * Baym  _et al._ [1971] G. Baym, C. Pethik, and P. Sutherland, The Ground State of Matter at High Densities: Equation of State and Stellar Models, Astrophys. J. 170, 299 (1971).
  * Moustakidis [2017] C. C. Moustakidis, The stability of relativistic stars and the role of the adiabatic index, Gen. Rel. Grav. 49, 68 (2017).
  * Lattimer and Prakash [2001] J. M. Lattimer and M. Prakash, Neutron star structure and the equation of state, The Astrophys. J. 550, 426 (2001).
  * Postnikov  _et al._ [2010] S. Postnikov, M. Prakash, and J. M. Lattimer, Tidal love numbers of neutron and self-bound quark stars, Phys. Rev. D 82, 024016 (2010).
  * Lattimer and Prakash [2005] J. M. Lattimer and M. Prakash, Ultimate energy density of observable cold baryonic matter, Phys. Rev. Lett. 94, 111101 (2005).
  * Lattimer and Prakash [2000] J. M. Lattimer and M. Prakash, Nuclear matter and its role in supernovae, neutron stars and compact object binary mergers, Physics Reports 333-334, 121 (2000).
  * Baiotti [2019] L. Baiotti, Gravitational waves from neutron star mergers and their relation to the nuclear equation of state, Prog. Part. Nucl. Phys. 109, 103714 (2019).
  * Flanagan and Hinderer [2008] E. E. Flanagan and T. Hinderer, Constraining neutron-star tidal love numbers with gravitational-wave detectors, Phys. Rev. D 77, 021502 (2008).
  * Hinderer [2008] T. Hinderer, Tidal love numbers of neutron stars, Astrophys. J. 677, 1216 (2008).
  * Damour and Nagar [2009] T. Damour and A. Nagar, Relativistic tidal properties of neutron stars, Phys. Rev. D 80, 084035 (2009).
  * Hinderer  _et al._ [2010] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral, Phys. Rev. D 81, 123016 (2010).
  * Fattoyev  _et al._ [2013] F. J. Fattoyev, J. Carvajal, W. G. Newton, and B.-A. Li, Constraining the high-density behavior of the nuclear symmetry energy with the tidal polarizability of neutron stars, Phys. Rev. C 87, 015806 (2013).
  * Lackey and Wade [2015] B. D. Lackey and L. Wade, Reconstructing the neutron-star equation of state with gravitational-wave detectors from a realistic population of inspiralling binary neutron stars, Phys. Rev. D 91, 043002 (2015).
  * Takátsy and Kovács [2020] J. Takátsy and P. Kovács, Comment on “tidal love numbers of neutron and self-bound quark stars”, Phys. Rev. D 102, 028501 (2020).
  * Thorne [1998] K. S. Thorne, Tidal stabilization of rigidly rotating, fully relativistic neutron stars, Phys. Rev. D 58, 124031 (1998).
  * Abbott  _et al._ [2019] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Properties of the binary neutron star merger gw170817, Phys. Rev. X 9, 011001 (2019).
  * Friedman and Stergioulas [2013] J. L. Friedman and N. Stergioulas,  _Rotating Relativistic Stars_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2013).
  * Takami  _et al._ [2011] K. Takami, L. Rezzolla, and S. Yoshida, A quasi-radial stability criterion for rotating relativistic stars, Mon. Not. R. Astron. Soc. Lett. 416, L1 (2011), https://academic.oup.com/mnrasl/article-pdf/416/1/L1/6140636/416-1-L1.pdf .
  * Weih  _et al._ [2017] L. R. Weih, E. R. Most, and L. Rezzolla, On the stability and maximum mass of differentially rotating relativistic stars, Mon. Not. R. Astron. Soc. Lett. 473, L126 (2017), https://academic.oup.com/mnrasl/article-pdf/473/1/L126/24841378/slx178.pdf .
  * Stergioulas [1996] N. Stergioulas, http://www.gravity.phys.uwm.edu/rns/ (1996).
  * Stergioulas and Friedman [1995] N. Stergioulas and J. L. Friedman, Comparing models of rapidly rotating relativistic stars constructed by two numerical methods, Astrophys. J. 444, 306 (1995).
  * Komatsu  _et al._ [1989] H. Komatsu, Y. Eriguchi, and I. Hachisu, Rapidly rotating general relativistic stars. I - Numerical method and its application to uniformly rotating polytropes, Mon. Not. R. Astron. Soc. 237, 355 (1989).
  * Cook  _et al._ [1994a] G. Cook, S. Shapiro, and S. Teukolsky, Rapidly rotating polytropes in general relativity, Astrophys. J. 422, 227 (1994a).
  * Goodstein [1985] D. L. Goodstein,  _States of Matter_ (Dove, New York, 1985).
  * Fetter and Walecka [2003] A. L. Fetter and J. D. Walecka,  _Quantum Theory of Many-Particle Systems_ (Dover,Mineola,New York, 2003).
  * Prakash [1994] M. Prakash,  _The Equation of State and Neutron Star_ (unpublished, 1994).
  * Burgio  _et al._ [2007] G. F. Burgio, M. Baldo, O. E. Nicotra, and H.-J. Schulze, A microscopic equation of state for protoneutron stars, Astrophys. J. Suppl. 308, 387 (2007).
  * Xu  _et al._ [2007] J. Xu, L.-W. Chen, B.-A. Li, and H.-R. Ma, Effects of isospin and momentum dependent interactions on liquid–gas phase transition in hot asymmetric nuclear matter, Phys. Lett. B 650, 348 (2007).
  * Moustakidis [2008] C. C. Moustakidis, Temperature and momentum dependence of single-particle properties in hot asymmetric nuclear matter, Phys. Rev. C 78, 054323 (2008).
  * Moustakidis and Panos [2009] C. C. Moustakidis and C. P. Panos, Equation of state for $\beta$-stable hot nuclear matter, Phys. Rev. C 79, 045806 (2009).
  * Tan  _et al._ [2016] N. H. Tan, D. T. Loan, D. T. Khoa, and J. Margueron, Mean-field study of hot $\beta$-stable protoneutron star matter: Impact of the symmetry energy and nucleon effective mass, Phys. Rev. C 93, 035806 (2016).
  * Takatsuka  _et al._ [1994] T. Takatsuka, S. Nishizaki, and J. Hiura, Particular Properties of Dense Supernova Matter, Prog. Theor. Phys. 92, 779 (1994), https://academic.oup.com/ptp/article-pdf/92/4/779/5358408/92-4-779.pdf .
  * Takatsuka [1996] T. Takatsuka, Hot Neutron Stars at Birth and Related Problems*), Prog. Theor. Phys. 95, 901 (1996), https://academic.oup.com/ptp/article-pdf/95/5/901/5280763/95-5-901.pdf .
  * Friedman  _et al._ [1988] J. L. Friedman, J. R. Ipser, and R. D. Sorkin, Turning point method for axisymmetric stability of rotating relativistic stars, Astrophys. J. 325, 722 (1988).
  * Sorkin [1982] R. D. Sorkin, A stability criterion for many parameter equilibrium families, Astrophys. J. 257, 847 (1982).
  * Goussard  _et al._ [1997] J. O. Goussard, P. Haensel, and J. L. Zdunik, Rapid uniform rotation of protoneutron stars, Astron. Astrophys. 321, 822 (1997).
  * Marques  _et al._ [2017] M. Marques, M. Oertel, M. Hempel, and J. Novak, New temperature dependent hyperonic equation of state: Application to rotating neutron star models and $i\text{$-$}q$ relations, Phys. Rev. C 96, 045806 (2017).
  * LORENE [1998] LORENE, Lorene: Langage objet pour la relativité numérique (1998).
  * Miller  _et al._ [2019] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Harding, W. C. G. Ho, J. M. Lattimer, R. M. Ludlam, S. Mahmoodifar, S. M. Morsink, P. S. Ray, T. E. Strohmayer, K. S. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin, and Y. Soong, PSR j0030+0451 mass and radius from NICER data and implications for the properties of neutron star matter, The Astrophysical Journal 887, L24 (2019).
  * Greif  _et al._ [2019] S. K. Greif, G. Raaijmakers, K. Hebeler, A. Schwenk, and A. L. Watts, Equation of state sensitivities when inferring neutron star and dense matter properties, Monthly Notices of the Royal Astronomical Society 485, 5363 (2019), https://academic.oup.com/mnras/article-pdf/485/4/5363/28262089/stz654.pdf .
  * Radice  _et al._ [2018] D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, GW170817: Joint constraint on the neutron star equation of state from multimessenger observations, The Astrophysical Journal 852, L29 (2018).
  * Tews  _et al._ [2018b] I. Tews, J. Margueron, and S. Reddy, Critical examination of constraints on the equation of state of dense matter obtained from gw170817, Phys. Rev. C 98, 045804 (2018b).
  * Radice and Dai [2019] D. Radice and L. Dai, Multimessenger parameter estimation of gw170817, The European Physical Journal A 55, 50 (2019), https://doi.org/10.1140/epja/i2019-12716-4 .
  * Kiuchi  _et al._ [2019] K. Kiuchi, K. Kyutoku, M. Shibata, and K. Taniguchi, Revisiting the lower bound on tidal deformability derived by AT 2017gfo, The Astrophysical Journal 876, L31 (2019).
  * Coughlin  _et al._ [2018] M. W. Coughlin, T. Dietrich, Z. Doctor, D. Kasen, S. Coughlin, A. Jerkstrand, G. Leloudas, O. McBrien, B. D. Metzger, R. O’Shaughnessy, and S. J. Smartt, Constraints on the neutron star equation of state from AT2017gfo using radiative transfer simulations, Monthly Notices of the Royal Astronomical Society 480, 3871 (2018), https://academic.oup.com/mnras/article-pdf/480/3/3871/33043575/sty2174.pdf .
  * Coughlin  _et al._ [2019] M. W. Coughlin, T. Dietrich, B. Margalit, and B. D. Metzger, Multimessenger Bayesian parameter inference of a binary neutron star merger, Monthly Notices of the Royal Astronomical Society: Letters 489, L91 (2019), https://academic.oup.com/mnrasl/article-pdf/489/1/L91/30032497/slz133.pdf .
  * Most  _et al._ [2018] E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, New constraints on radii and tidal deformabilities of neutron stars from gw170817, Phys. Rev. Lett. 120, 261103 (2018).
  * Foley  _et al._ [2020] R. J. Foley, D. A. Coulter, C. D. Kilpatrick, A. L. Piro, E. Ramirez-Ruiz, and J. Schwab, Updated parameter estimates for GW190425 using astrophysical arguments and implications for the electromagnetic counterpart, Monthly Notices of the Royal Astronomical Society 494, 190 (2020), https://academic.oup.com/mnras/article-pdf/494/1/190/32982848/staa725.pdf .
  * Rezzolla  _et al._ [2018] L. Rezzolla, E. R. Most, and L. R. Weih, Using gravitational-wave observations and quasi-universal relations to constrain the maximum mass of neutron stars, The Astrophysical Journal 852, L25 (2018).
  * Raithel  _et al._ [2018] C. A. Raithel, F. Özel, and D. Psaltis, Tidal deformability from GW170817 as a direct probe of the neutron star radius, Astrophys. J. 857, L23 (2018).
  * Zhao and Lattimer [2018] T. Zhao and J. M. Lattimer, Tidal deformabilities and neutron star mergers, Phys. Rev. D 98, 063020 (2018).
  * De  _et al._ [2018] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, Tidal deformabilities and radii of neutron stars from the observation of gw170817, Phys. Rev. Lett. 121, 091102 (2018).
  * Annala  _et al._ [2018] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Gravitational-wave constraints on the neutron-star-matter equation of state, Phys. Rev. Lett. 120, 172703 (2018).
  * Fasano  _et al._ [2019] M. Fasano, T. Abdelsalhin, A. Maselli, and V. Ferrari, Constraining the neutron star equation of state using multiband independent measurements of radii and tidal deformabilities, Phys. Rev. Lett. 123, 141101 (2019).
  * Fattoyev  _et al._ [2018] F. J. Fattoyev, J. Piekarewicz, and C. J. Horowitz, Neutron skins and neutron stars in the multimessenger era, Phys. Rev. Lett. 120, 172702 (2018).
  * Landry  _et al._ [2020] P. Landry, R. Essick, and K. Chatziioannou, Nonparametric constraints on neutron star matter with existing and upcoming gravitational wave and pulsar observations, Phys. Rev. D 101, 123007 (2020).
  * Most  _et al._ [2020] E. R. Most, L. J. Papenfort, L. R. Weih, and L. Rezzolla, A lower bound on the maximum mass if the secondary in GW190814 was once a rapidly spinning neutron star, Mon. Not. R. Astron. Soc. Lett. 499, L82 (2020).
  * Cook  _et al._ [1994b] G. Cook, S. Shapiro, and S. Teukolsky, Rapidly rotating neutron stars in general relativity: Realistic equations of state, Astrophys. J. 424, 823 (1994b).
  * Salgado  _et al._ [1994] M. Salgado, S. Bonazzola, E. Gourgoulhon, and P. Haensel, High precision rotating neutron star models 1: Analysis of neutron star properties, Astron. Astrophys. 291, 155 (1994).
  * Tsang  _et al._ [2019] C. Tsang, M. Tsang, P. Danielewicz, F. Fattoyev, and W. Lynch, Insights on skyrme parameters from gw170817, Phys. Lett. B 796, 1 (2019).
  * Weih  _et al._ [2020] J.-B. Weih, J.-J. Lu, G. F. Burgio, Z.-H. Li, and H.-J. Schulze, Are nuclear matter properties correlated to neutron star observables?, Europ. Phys. J. A 56, 63 (2020).
  * Thorne [1974] K. S. Thorne, Disk-accretion onto a black hole. ii. evolution of the hole, Astrophy. J. 191, 507 (1974).
  * Morsink  _et al._ [1999] S. M. Morsink, N. Stergioulas, and S. R. Blattnig, Quasi-normal modes of rotating relativistic stars: Neutral modes for realistic equations of state, Astrophys. J. 510, 854 (1999).


# Data-Driven Contract Design for Multi-Agent Systems with Collusion Detection

Nayara Aguiar, Parv Venkitasubramaniam and Vijay Gupta Nayara Aguiar and Vijay
Gupta are with the Department of Electrical Engineering, University of Notre
Dame, IN, USA<EMAIL_ADDRESS>Vekitasubramaniam is with the
Department of Electrical Engineering, Lehigh University, PA, USA
<EMAIL_ADDRESS>
###### Abstract

In applications such as participatory sensing and crowd sensing, self-
interested agents exert costly effort towards achieving an objective for the
system operator. We study such a setup where a principal incentivizes multiple
agents of different types who can collude with each other to derive rent. The
principal cannot observe the efforts exerted directly, but only the outcome of
the task, which is a noisy function of the effort. The type of each agent
influences the effort cost and task output. For a duopoly in which agents are
coupled in their payments, we show that if the principal and the agents
interact finitely many times, the agents can derive rent by colluding even if
the principal knows the types of the agents. However, if the principal and the
agents interact infinitely often, the principal can disincentivize agent
collusion through a suitable data-driven contract.

###### Index Terms:

Agent collusion, game theory, dynamic contract.

##  I Introduction

Many scenarios in smart infrastructure systems require a system operator to
incentivize self-interested agents to exert costly effort to make decisions
that align with the goal of the operator. For instance, in participatory
sensing, a system operator requires many autonomous sensors to take
measurements to allow estimation of a global quantity. The operator cannot
observe directly the effort of each sensing agent (possibly for privacy
reasons) and the agent might not benefit directly from the goal of the
operator and thus needs to be compensated based on noisy outputs. Incentive
design in this setting has been analyzed, particularly in the context of
participatory sensing and crowd sensing, for classical estimation as well as
learning based tasks (see [1, 2] and the references therein).

One possible abstraction of the problem is as a Stackelberg game (e.g. as
adopted in the context of jamming [3, 4]). A principal seeks to incentivize
self-interested rational agents to exert costly effort to perform a task. The
principal acting as the leader issues a contract that specifies the payment
structure. The agents as the followers respond by deciding whether to accept
the contract or not. If the contract is accepted, the agents further decide
how much effort to exert. If the agents and the principal interact multiple
times, then the entire process is repeated. In this paper, we focus on a
duopoly, where each agent $i\in\mathcal{N}:=\\{1,2\\}$ can be defined as a
type $\theta_{i}$ which influences his utility and the quality of his task
output. Further, we assume that the principal can observe only the task
output, but not the effort exerted by each agent, leading to a problem of
moral hazard.

Single-agent, single-principal problems with such (or similar) abstractions
have been widely explored in the literature. Most such works that consider
multiple agents being present assume that either the agents are decoupled or
that collusion is prohibited, thus reducing the analysis to pairwise
interactions between each individual agent and the principal (possibly with an
additional optimization problem at the principal on which agents to select for
the contract). We focus on the effect of collusion, specifically seeking to
answer the question if it is possible for the agents to derive rent through
collusion. We show that there is a difference between the principal and the
agents interacting finitely or infinitely many times. Rent is possible in the
former situation, but the principal can propose a dynamic contract with a
collusion detection mechanism disincentivizing collusion in the latter.
Interestingly, this result holds even if the principal knows the types of the
agents.

Three lines of literature are particularly relevant here. Dynamic contract
design for multi-agent systems with information assymetry has been studied in
game theory (see, e.g., [5, 6]. In signal processing applications, a repeated
game is analyzed in [7], where the principal performs random verification and
assigns a reputation score to each sensor. A dynamic estimation problem is
studied in [8], leading to an equilibrium in which strategic sensors use a
memory-less policy. In [9], a data-driven contract is analyzed for the
situation when the quality of the output produced by participating agents
depends on their true types. As opposed to these works, we consider multiple
agents that are coupled both due to their compensation and possibility of
collusion.

Another line of relevant literature is that of payment schemes which couple
agents. Competition among data aggregators is studied in [10], where
aggregators are coupled in their costs and data sources are coupled in their
compensation. Payment coupling is also considered in [11] in a demand response
problem. However, the issue of collusion has received much less attention in
the literature.

The third line of literature is one which studies collusion. In [12], the
formation of a coalition is analyzed for a class of static estimation problems
with strategic sensors. In a static setting, the principal does not have
information about whether agents are communicating to devise a joint strategy,
which may lead to agents being able to derive rent by colluding. In contrast,
our work considers the problem when the agents and the principal interact
repeatedly and have the ability to change their strategies. We show that, in
the limit of infinitely many interactions, the principal can utilize the
previous task outputs to perform a hypothesis test on whether collusion is
underway and update the contract offered to disincentivize it.

The remainder of this paper is organized as follows. Section II introduces the
problem formulation proposed. The static problem is solved in Section III,
where we derive conditions for collusion to benefit the agents. Section IV
presents the dynamic contract which enables the principal to verify if
collusion is present. Lastly, Section V concludes this work.

##  II Problem Formulation

Consider a Stackelberg game in which a principal, as a leader, establishes a
contract to incentivize two agents to exert costly effort to perform a task.
The task output of each agent $i\in\mathcal{N}:=\\{1,2\\}$ is represented by
the random variable $X_{i}$, which can be observed by the principal. We let
$X_{i}$ be defined by a function
$f_{i}(a_{i},\theta_{i},\eta_{i}):\mathcal{A}\times\Theta\times\mathbb{R}\mapsto\mathbb{R}$,
which depends on the agent’s effort $a_{i}$ and type $\theta_{i}$, and on the
random noise $\eta_{i}$. The noise variables $\\{\eta_{i}\\}$ are assumed to
be independent and identically distributed across the agents, and bounded in
the interval $\eta_{i}\in[\underline{\eta}_{i},\overline{\eta}_{i}]$.

The timeline for a specified time horizon $K$ is as follows. Initialize $k=0$.

  1. 1.

$t=4k+1$: The principal issues a contract defined by the payment
$w_{i}(\mathbf{X})$, where $\mathbf{X}:=\\{X_{1},X_{2}\\}$.

  2. 2.

$t=4k+2$: Each agent decides whether or not to (i) accept this contract, and
(ii) collude with the other agent.

  3. 3.

$t=4k+3$: Non-colluding agents exert an effort $a_{i}$ which maximizes their
individual utility, and colluding agents exert the effort $a_{i}$ that
increases both of their utilities as compared to the non-colluding solution.

  4. 4.

$t=4k+4$: For each participating agent, the effort $a_{i}$ produces an output
$X_{i}$, which is compensated by the principal according to
$w_{i}(\mathbf{X})$.

  5. 5.

$k=k+1$ and repeat steps 1-5 while $k<K$.

In this problem, the utility of the principal is given by

$U=\mathbb{E}\left[S(\mathbf{X})-\sum_{i=1}^{2}w_{i}(\mathbf{X})\right],$ (1)

where $S(\mathbf{X})=X_{1}+X_{2}$ is the benefit derived due to the outputs
produced, and the expectation is taken over the output noise. The utility of
each agent is defined as

$u_{i}=\mathbb{E}\left[w_{i}(\mathbf{X})-h(a_{i},\theta_{i})\right],$ (2)

where, for concreteness, the effort cost is defined as
$h(a_{i},\theta_{i})=\frac{a_{i}^{2}}{2\theta_{i}}$, so that it is increasing
in the effort, and decreasing in the type. Finally, we assume that
$X_{i}=a_{i}+Q(\theta_{i})\eta_{i}$, where $Q(\theta_{i})$ is a increasing
function of the type.

##  III One-Shot Contract

Having a limited budget, the principal uses a Cournot-like payment function
which discourages large aggregated efforts:

$w_{i}(\mathbf{X})=X_{i}\left(\lambda-\sum_{i=1}^{2}X_{i}\right),$ (3)

where $\lambda$ is a design variable. Given the utility function (1) and the
payments (3), the problem of the principal is given by

$\displaystyle\underset{\lambda}{\max}~{}U$ (4)
$\displaystyle\text{s.t.}~{}\mathbb{E}[w_{i}(\mathbf{X})-h$
$\displaystyle(a_{i}^{*},\theta_{i})]\geq 0~{}\forall i$ $\displaystyle
a_{i}^{*}=\text{argmax}~{}$
$\displaystyle\mathbb{E}\left[w_{i}(\mathbf{X})-h(a_{i},\theta_{i})\right]~{}\forall
i,$

where the first constraint is the individual rationality (IR) constraint which
reflects that the agents will only accept the contract if their expected
utility is non-negative, and the second constraint is the incentive
compatibility constraint which enforces that the efforts considered in the
problem are such that they maximize the agents’ utilities.

After solving his problem, the principal issues the contract to the agents.
The payment (3) couples the outputs from both agents, and therefore
competition is induced at the agents. Each agent solves the problem
$\underset{a_{i}}{\max}~{}u_{i}(X_{i},X_{-i}),$ where we make explicit the
dependence of agent $i$’s utility on his own output $X_{i}$ and the output of
the other agent $X_{-i}$.

We begin by analyzing this duopoly in a scenario where the agents play their
Cournot strategies, and then we establish conditions for which the agents can
derive rent by colluding.

### Cournot Duopoly

We use backwards induction to find the equilibrium of this Cournot
competition.

######  Proposition 1.

Given payments $w_{i}(\mathbf{X})$ as expressed in (3), and in the absence of
collusion, the principal maximizes his own utility by choosing

$\lambda^{*}=\max\left\\{\frac{(2\theta_{1}+1)(2\theta_{2}+1)-\theta_{1}\theta_{2}}{2(\theta_{1}+1)(\theta_{2}+1)},\underline{\lambda}\right\\}$
(5)

where $\underline{$\lambda$}$ is the minimum value that guarantees individual
rationality for all agents, and is given by

$\underline{\lambda}=\max\left\\{\frac{\sqrt{2}Q(\theta_{i})\sigma[(2\theta_{1}+1)(2\theta_{2}+1)-\theta_{1}\theta_{2}]}{(\theta_{-i}+1)\sqrt{\theta_{i}(2\theta_{i}+1)}}\right\\}.$
(6)

In response, the optimal efforts of agents $i\in\\{1,2\\}$ is

$a_{i}^{*}=\frac{\lambda^{*}\theta_{i}(\theta_{-i}+1)}{(2\theta_{i}+1)(2\theta_{-i}+1)-\theta_{i}\theta_{-i}},$
(7)

where $\theta_{-i}$ is the other agent in the duopoly.

###### Proof.

For a given $\lambda$, the utility of agent $i=1$ is

$\displaystyle u_{1}$ $\displaystyle=a_{1}(\lambda-
a_{1}-a_{2})-Q(\theta_{1})^{2}\sigma^{2}-a_{1}^{2}/2\theta_{1},$ (8)

where we used the fact that the outputs $X_{i}$ are independent,
$\mathbb{E}[X_{i}]=a_{i}~{}\forall i$, and
$\mathbb{E}[X_{i}^{2}]=a_{i}^{2}+Q(\theta_{1})^{2}\sigma^{2}~{}\forall i$.
This function is concave in the decision variable $a_{1}$, and thus we use the
first order condition to derive the equilibrium effort

$\frac{\partial u_{1}}{\partial
a_{1}}=\lambda-2a_{1}-a_{2}-\frac{a_{1}}{\theta_{1}}=0\implies
a_{1}^{*}=\frac{\theta_{1}(\lambda-a_{2}^{*})}{2\theta_{1}+1}.$ (9)

Proceeding similarly for agent $i=2$, we can find an expression for the
equilibrium effort $a_{2}^{*}$ as a function of $a_{1}^{*}$, then substitute
it in (9) to find the expression (7). For individual rationality to hold, the
agents must have a non-negative utility. For the optimal efforts derived, this
condition reduces to

$u_{i}=\frac{\lambda^{2}\theta_{i}(\theta_{-i}+1)^{2}(2\theta_{i}+1)}{2\left[(2\theta_{i}+1)(2\theta_{-i}+1)-\theta_{i}\theta_{-i}\right]^{2}}-Q(\theta_{i})^{2}\sigma^{2}\geq
0~{}\forall i,$ (10)

which establishes a lower bound on $\lambda$. With these optimal efforts, the
principal’s utility can be written as

$\displaystyle U=\lambda$ $\displaystyle
g(\theta_{1},\theta_{2})-\lambda^{2}(g(\theta_{1},\theta_{2})-g(\theta_{1},\theta_{2})^{2})$
$\displaystyle+[Q(\theta_{1})^{2}+Q(\theta_{2})^{2}]\sigma^{2},$

where $g(\theta_{1},\theta_{2})=(a_{1}^{*}+a_{2}^{*})/\lambda$ and is given by

$g(\theta_{1},\theta_{2})=\frac{\theta_{1}(\theta_{2}+1)+\theta_{2}(\theta_{1}+1)}{(2\theta_{1}+1)(2\theta_{2}+1)-\theta_{1}\theta_{2}}.$
(11)

The principal’s utility can be shown to be concave in $\lambda$. Thus, we can
use the first order condition to find the optimal decision as in the first
expression in (5). This value of $\lambda^{*}$ satisfies the individual
rationality condition (10) for all agents if

$Q(\theta_{i})^{2}\sigma^{2}\leq\frac{\theta_{i}(2\theta_{i}+1)}{8(\theta_{i}+1)^{2}}~{}\forall
i.$ (12)

If the inequality above does not hold, then $\lambda^{*}$ must be equal to
maximum of the lower bounds established through (10). ∎

The optimal effort (7) can be shown to be increasing in agent $i$’s own type
for a fixed $\theta_{-i}$, and thus agents with higher type values will exert
more effort. Conversely, if agent $i$’s type is fixed, then his effort is
decreasing with $\theta_{-i}$. This is due to the trade-off established
through the payment function. We can also observe that agents with high output
variance will require higher payments to participate in the contract, as their
output may be very low even if a high effort is exerted.

### Collusion Strategy

Collusion becomes individually rational for both agents if the utility they
derive by exerting the optimal effort (7) can be surpassed by picking off-
equilibrium efforts that will yield a surplus. In essence, the agents will
collude if there exists a Pareto improvement that makes them both better off.
Formally, this condition can be expressed as

$\exists~{}(\overline{a}_{i},\overline{a}_{-i},P)~{}\text{s.t.}~{}u_{i}(\overline{a}_{i},\overline{a}_{-i})+P>u_{i}(a_{i}^{*},a_{-i}^{*})~{}\forall
i,$ (13)

where $P$ is a possible side payment among the agents, which is positive
(negative) for the agent receiving (making) it.

Consider a scenario in which the two agents collude to exert efforts which,
combined, yield a certain $\theta_{i}$-monopoly effort. Such effort is the
optimal effort exerted by an agent of type $\theta_{i}$ when he is the only
agent providing outputs to the principal, and is defined as
$a^{M}_{i}=\lambda\theta_{i}/(2\theta_{i}+1)$.

For simplicity, we assume the following collusion strategies:

  * •

if agents are of the same type, they equally split the amount of effort
exerted to reach monopoly level, and no side payment is performed, and

  * •

if agents are of different types, the agent with higher type exerts all the
monopoly-level effort and makes a side payment to the lower-type agent, who
exerts no effort.

The next result identifies conditions for the agents to have an incentive to
collude.

######  Theorem 1.

If the agents are of the same type, that is $\theta_{1}=\theta_{2}=\theta$,
then collusion to mimic the $\theta_{1}$-monopoly effort by choosing the
efforts $\overline{a}_{1}=\overline{a}_{2}=a_{1}^{M}/2$ is individual rational
for both agents if $\theta>(1+\sqrt{17})/8\approx 0.64$. Further, if the
agents are of different types
$\theta_{1}=\theta_{h}>\theta_{\ell}=\theta_{2}$, and satisfy one of
conditions

$\displaystyle\text{C1:}~{}\theta_{\ell}\geq 1+\sqrt{2}$ (14) or
$\displaystyle\text{C2:}\begin{cases}0\leq\theta_{\ell}<1+\sqrt{2}~{}\text{and}~{}\\\
\theta_{h}\geq\frac{\sqrt{(\theta_{\ell}+1)(17\theta_{\ell}+9)}+3\theta_{\ell}+1}{2(\theta_{\ell}+2)}\end{cases}$
(15)

then there exists a side payment $P$ from the $\theta_{h}$ agent to the
$\theta_{\ell}$ agent that makes the colluding efforts
$\overline{a}_{1}=a_{h}^{M},~{}\overline{a}_{2}=0$ to mimic a
$\theta_{h}$-monopoly individual rational.

###### Proof.

Agents of the same type $\theta$ collude by exerting the same effort and
without making a side payment. Then, their IR constraint (13) is the same and
can be written as

$\overline{a}_{i}(\lambda-\overline{a}_{i}-\overline{a}_{-i})-\overline{a}_{i}^{2}/2\theta>a_{i}^{*}(\lambda-
a_{i}^{*}-a_{-i}^{*})-{a_{i}^{*}}^{2}/2\theta.$ (16)

Substituting the proposed collusion efforts
$\overline{a}_{1}=\overline{a}_{2}=a_{1}^{M}/2$ and the Cournot strategies
(7), we can find the inequality

$\displaystyle\frac{\lambda^{2}\theta(4\theta+3)}{8(2\theta+1)^{2}}$
$\displaystyle>\frac{\lambda^{2}\theta(\theta+1)^{2}(2\theta+1)}{2\left[(2\theta+1)^{2}-\theta^{2}\right]^{2}}$
(17) $\displaystyle\implies 4\theta^{5}+11\theta^{4}+8\theta^{3}$
$\displaystyle>2\theta^{2}+4\theta+1,$ (18)

which holds for $\theta>(1+\sqrt{17})/8$.

If agents are of different types $\theta_{\ell}$ and
$\theta_{h}>\theta_{\ell}$, let the side payment from agent $\theta_{h}$ to
agent $\theta_{\ell}$ be a fraction of agent $\theta_{h}$’s expected payment,
that is $|P|=\alpha\mathbb{E}[w_{h}(\mathbf{X})]$, where $0\leq\alpha\leq 1$.
Define
$D=2\left[(2\theta_{h}+1)(2\theta_{\ell}+1)-\theta_{h}\theta_{\ell}\right]^{2},$
$R_{h}=\frac{Q(\theta_{h})^{2}\sigma^{2}}{\lambda^{2}\theta_{h}},$ and
$R_{\ell}=\frac{Q(\theta_{\ell})^{2}\sigma^{2}}{\lambda^{2}\theta_{h}}.$ For
the collusion efforts $\overline{a}_{1}=a_{h}^{M}$ and $\overline{a}_{2}=0$,
the constraint (13) for agent $\theta_{h}$ reduces to

$\frac{(\theta_{h}+1/2)-\alpha(\theta_{h}+1)}{(2\theta_{h}+1)^{2}}+\alpha
R_{h}>\frac{(\theta_{\ell}+1)^{2}(2\theta_{h}+1)}{D},$ (19)

from which we can find the following upper bound on $\alpha$

$\alpha<\frac{(2\theta_{h}+1)D/2-(\theta_{\ell}+1)^{2}(2\theta_{h}+1)^{3}}{D\left[(\theta_{h}+1)-R_{h}(2\theta_{h}+1)^{2}\right]}.$
(20)

For agent $\theta_{\ell}$, the individual rationality constraint (13) becomes

$\frac{\alpha(\theta_{h}+1)}{(2\theta_{h}+1)^{2}}-\alpha
R_{h}>\frac{\theta_{\ell}(\theta_{h}+1)^{2}(2\theta_{\ell}+1)}{D\theta_{h}}-R_{\ell},$
(21)

leading to the lower bound on the side payment fraction

$\alpha>\frac{[\theta_{\ell}(\theta_{h}+1)^{2}(2\theta_{\ell}+1)-D\theta_{h}R_{\ell}](2\theta_{h}+1)^{2}}{D\theta_{h}\left[(\theta_{h}+1)-R_{h}(2\theta_{h}+1)^{2}\right]}.$
(22)

For this collusion to be individual rational for both agents, the bounds
established for $\alpha$ must yield a non-empty interval. The conditions
presented in Theorem 1 are sufficient, and can be found for the case where
agent $\theta_{\ell}$’s output has zero variance, so that $R_{\ell}=0$ and the
lower bound (22) becomes higher. ∎

When agents are of the same type, the lower bound on $\theta$ presented in
Theorem 1 indicates that, when the type value of the agents is too low, they
cannot improve their utilities by colluding to mimic the $\theta_{1}$-monopoly
effort. Since the effort is increasing with type, this could be seen as a
limiting situation in which the agents cannot exert an effort that is high
enough to drive prices significantly. Similarly, when agents are of different
types, collusion becomes rational when the agent with the lowest type exceeds
a certain threshold; otherwise, the high-type agent must have a high enough
type value for collusion to occur. We note that, under a static contract, the
agents satisfying these conditions would still derive rent by colluding if
they interact finitely many times with the principal.

In the following result, we relax the assumption that agents of the same type
would mimic their own type’s monopoly to provide conditions for which they can
extract extra rent by choosing a combined effort that imitates a different
monopoly.

######  Corollary 1.

Let agents of the same type $\theta$ equally split the effort exerted, and not
make side payments when colluding. Then, for
$\theta<\hat{\theta}<\theta(4\theta+3)$, collusion to mimic a
$\hat{\theta}$-monopoly yields a higher utility to both agents than if they
mimic their own type’s monopoly.

###### Proof.

We check if it is individual rational for both agents to mimic a
$\hat{\theta}$-monopoly, for any $\hat{\theta}\neq\theta$. For this case, the
agents have the same individual rationality constraint, given by

$\frac{\lambda^{2}\hat{\theta}(4\theta\hat{\theta}+4\theta-\hat{\theta})}{8\theta(2\hat{\theta}+1)^{2}}>\frac{\lambda^{2}\theta(4\theta+3)}{8(2\theta+1)^{2}},$
(23)

which compares the agent’s utility in a collusion at a $\hat{\theta}$-monopoly
and at a $\theta$-monopoly. This inequality can be evaluated to show that it
only holds if $\theta<\hat{\theta}\leq\theta(4\theta+3)$. ∎

Corollary 1 shows that it is not individual rational for agents of the same
type $\theta$ to jointly exert a $\hat{\theta}$-monopoly effort for any
$\hat{\theta}<\theta$. In this case, the payments received would be too low
for extra rent to be derived. On the other hand, pretending to be a monopoly
of a type $\hat{\theta}>\theta$ can be beneficial to both agents, as long as
this higher type satisfies the upper bound presented. For any other
$\hat{\theta}>\theta(4\theta+3)$, the effort cost would dominate, nullifying
any possible increase in the payment received.

##  IV Dynamic Contract

If the principal and the agents interact continuously in an infinite horizon
setting, the principal is able to perform a noisy verification of the agents’
outputs, so as to learn about their behavior and adjust the contract if
necessary. For notational simplicity, in what follows, let

$\hat{\lambda}_{i}=\frac{\lambda^{*}(\theta_{-i}+1)}{(2\theta_{1}+1)(2\theta_{2}+1)-\theta_{1}\theta_{2}}.$
(24)

In the dynamic contract, the principal collects a sequence of $n$ outputs
$\mathcal{X}_{i}:=\\{X_{i}^{1},...,X_{i}^{n}\\}$ from each agent $i$, and
performs the following hypothesis test:

$\displaystyle\delta^{i}(\mathcal{X}_{i})$
$\displaystyle=\begin{cases}H_{0}^{i}&\text{if}~{}\big{|}\overline{X}_{i}-\hat{\lambda}_{i}\theta_{i}\big{|}<\epsilon\\\
H_{1}^{i}&\text{otherwise}\end{cases}$ (25)

where $\overline{X}_{i}=\frac{1}{n}\sum_{k=1}^{n}X_{i}^{k}$ is the sample
mean. The hypothesis $H_{0}^{i}$ indicates that agent $i$ is playing the
Cournot equilibrium, whereas $H_{1}^{i}$ flags that agent $i$ is deviating
from this strategy, which is indicative that collusion is underway. If both
agents are flagged, the principal updates the payment function so as to
penalize them for the remainder of the time, disincentivizing collusion. This
verification process is repeated every $n$ time instants, and the principal
utilizes the latest $n$ samples to perform the hypothesis test, as detailed in
Algorithm 1.

Algorithm 1 Dynamic Contract

1:$k\leftarrow 1$

2:for $k\leq n$ do

3: Principal issues contract (3) with $\lambda$ as in Eq. (5)

4: if Agent $i$ accepts contract and exerts effort $a_{i}$ then

5: $X_{i}^{k}\leftarrow a_{i}+Q(\theta_{i})\eta_{i}$

6: end if

7: $k\leftarrow k+1$

8:end for

9:repeat

10: if Any agent rejects any contract for $k\in[k-n+1,k]$ then

11: $\lambda\leftarrow 0$

12: else if $\delta^{i}(X_{i}^{k-n+1},...,X_{i}^{k})=H_{1}^{i}~{}\forall i$
then

13: $\lambda\leftarrow 0$

14: else

15: $\lambda$ unchanged

16: end if

17: $t\leftarrow 1$

18: for $t\leq n$ do

19: Principal issues contract (3) with most recent $\lambda$

20: if Agent $i$ accepts contract then

21: $X_{i}^{t}\leftarrow a_{i}+Q(\theta_{i})\eta_{i}$

22: end if

23: $k\leftarrow k+1$

24: $t\leftarrow t+1$

25: end for

26:until $k\leq K$

######  Theorem 2.

For the dynamic contract in Algorithm 1, if the time window during which data
is collected is $n\in o(K)$, then the principal can detect if the agents are
playing off-equilibrium strategies almost surely as $K\to\infty$.

###### Proof.

When agent $i$ plays Cournot, the probability of false alarm can be bounded as

$\text{Prob}\left\\{\left|\overline{X}_{i}-\hat{\lambda}_{i}\theta_{i}\right|\geq\epsilon\right\\}\leq
2\exp{\left(-\frac{2n\epsilon^{2}}{(\overline{\eta}_{i}-\underline{\eta}_{i})^{2}}\right)}$

using Hoeffding’s inequality. This probability falls exponentially to zero as
$n\to\infty$. If the agents collude, the probability of missed detection is
given by
$\text{Prob}\\{|\overline{X}_{i}-\hat{\lambda}_{i}\theta_{i}|<\epsilon\\}$. In
this case, the sample mean $\overline{X}_{i}\to a_{i}$ as $n\to\infty$ by the
law of large numbers. Since the efforts exerted by the colluding agents are
different from their Cournot equilibrium,
$\text{Prob}\\{|\overline{X}_{i}-\hat{\lambda}_{i}\theta_{i}|<\epsilon\\}\to
0$. Thus, as $K\to\infty$, the probability of incorrect detection falls to
zero as long as $n\in o(K)$. ∎

The last result indicates that the principal can eliminate the incentive for
collusion in the long run by dynamically updating the contract according to
the outputs received.

##  V Conclusion

We formulated and analyzed a contract design problem for a principal who seeks
to incentivize agents to exert costly effort to perform a task. We considered
agents to be coupled in their compensation through a Cournot-like payment, and
allowed them to collude. For a static scenario, we showed that agents can
extract rent by colluding. However, if the principal interacts infinitely
often with the agents, a dynamic contract can be implemented by the principal
to use historical data as a way to detect and disincnetivize such collusion.

## References

  * [1] H. Gao, C. H. Liu, W. Wang, J. Zhao, Z. Song, X. Su, J. Crowcroft, and K. K. Leung, “A survey of incentive mechanisms for participatory sensing,” _IEEE Communications Surveys Tutorials_ , vol. 17, no. 2, pp. 918–943, 2015.
  * [2] F. Restuccia, S. K. Das, and J. Payton, “Incentive mechanisms for participatory sensing: Survey and research challenges,” _ACM Trans. Sen. Netw._ , vol. 12, no. 2, Apr. 2016. [Online]. Available: https://doi-org.proxy.library.nd.edu/10.1145/2888398
  * [3] A. Garnaev, A. Petropulu, W. Trappe, and H. V. Poor, “A power control game with uncertainty on the type of the jammer,” in _2019 IEEE Global Conference on Signal and Information Processing (GlobalSIP)_ , 2019, pp. 1–5.
  * [4] A. Mukherjee and A. L. Swindlehurst, “Jamming games in the mimo wiretap channel with an active eavesdropper,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 1, pp. 82–91, 2013.
  * [5] B. Holmström, “Managerial incentive problems: A dynamic perspective,” _The review of Economic studies_ , vol. 66, no. 1, pp. 169–182, 1999.
  * [6] J. Prat and B. Jovanovic, “Dynamic contracts when the agent’s quality is unknown,” _Theoretical Economics_ , vol. 9, no. 3, pp. 865–914, 2014.
  * [7] D. G. Dobakhshari, P. Naghizadeh, M. Liu, and V. Gupta, “A reputation-based contract for repeated crowdsensing with costly verification,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 23, pp. 6092–6104, 2019.
  * [8] F. Farokhi, A. M. H. Teixeira, and C. Langbort, “Estimation with strategic sensors,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 2, pp. 724–739, 2017.
  * [9] P. Venkitasubramaniam and V. Gupta, “Data-driven contract design,” in _2019 American Control Conference (ACC)_ , 2019, pp. 2283–2288.
  * [10] T. Westenbroek, R. Dong, L. J. Ratliff, and S. S. Sastry, “Competitive statistical estimation with strategic data sources,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 4, pp. 1537–1551, 2020.
  * [11] D. G. Dobakhshari and V. Gupta, “A contract design approach for phantom demand response,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 5, pp. 1974–1988, 2019.
  * [12] F. Farokhi, I. Shames, and M. Cantoni, “Promoting truthful behavior in participatory-sensing mechanisms,” _IEEE Signal Processing Letters_ , vol. 22, no. 10, pp. 1538–1542, 2015.


# Maximum mass and stability of differentially rotating neutrons
stars††thanks: Presented at The 8th Conference of the Polish Society on
Relativity

Paweł Szewczyk  Dorota Gondek-Rosińska
Pablo Cerdá-Durán  Astronomical Observatory Warsaw University, 00-478 Warsaw,
Poland  Observatori Astronòmic, Universitat de València, E-46980, Paterna
(València), Spain

###### Abstract

We present our study of stability of differentially rotating, axisymmetric
neutron stars described by a polytropic equation of state with $\Gamma=2$. We
focus on quasi-toroidal solutions with a degree of differential rotation
$\widetilde{A}=1$. Our results show that for a wide range of parameters
hypermassive, quasi-toroidal neutron stars are dynamically stable against
quasi-radial perturbations, which may have implications for newly born neutron
stars and binary neutron stars mergers.

##  1 Introduction

Differential rotation seems to appear naturally in many dynamical scenarios
involving neutron stars (NS), including the collapse of stellar cores (see
e.g. [23]) and binary neutron star (BNS) mergers (see e.g. [13]). Its
stabilizing effect may allow for configurations with masses significantly
higher than the mass limit for rigidly rotating neutron stars. Its study is
relevant for the understanding of black hole formation in those astrophysical
scenarios with consequences in observations of core-collapse supernovae (CCSN)
and BNS mergers, especially with current gravitational wave ground-based
observatories (LIGO, Virgo and Kagra [15, 1, 12]) and future ones (the
Einstein Telescope and Cosmic Explorer).

###  1.1 Equilibrium models of differentially rotating NS

The solution space of differentially rotating neutron stars in equilibrium was
already extensively studied by different authors. It was shown by [3] that
differentially rotating NS with masses significantly larger than non-rotating
or rigidly rotating NS can exist and be stable against radial collapse and bar
formation. Those with masses larger than the limit for rigidly rotating
objects are called hypermassive NS.

However, studying the whole solution space has proven to be numerically
challenging. The existence of different types of solutions of differentially
rotating neutron stars was for the first time found by [2] using relativistic
highly accurate and stable multi-domain spectral numerical code FlatStar. Most
importantly, for a given degree of differential rotation, the solution is not
uniquely determined by the maximal density and angular momentum of the NS (or
any other suitable pair of parameters), as is the case for rigid rotation.
Instead, different types of solutions may coexist for the same parameters. The
maximum mass for different degrees of differential rotation and different
solution types was presented by [10] and [19] for polytropes, showing that the
most massive configurations are obtained for modest degree of differential
rotation. Similar results were obtained for strange quark stars [20] and NS
with several realistic equations of state [7].

While many studies in the past were using a rotation law of [14], which is
mainly consistent with CCSN remnants [23], the rotation law observed in
simulations of BNS merger remnants departs significantly from that one.
Rotation laws better suited for BNS mergers have been proposed by [22], and
its impact on the solution space of equilibrium models studied by [11].

###  1.2 Stability properties of hypermassive NS

For non-rotating NS the limit for both secular and dynamical stability occurs
at the point of maximal mass ($M_{TOV}$). This criterion can be, to some
point, extended to rigidly rotating NS. The so-called turning point criterion
was presented by [9] and proven to be a sufficient criterion for instability.
It states that the point of maximal gravitational mass $M$ on a sequence of
configurations of fixed angular momentum $J$ ($J$-constant turning points),
or, alternatively, the point of minimal gravitational mass on a sequence of
fixed rest mass $M_{0}$ ($M_{0}$-constant turning points) marks the onset of
instability. This criterion, however, does not give the exact threshold to
collapse. The neutral-stability point where F-mode frequency vanishes differs
from the turning-point line [21]. Numerical simulations confirm that the
neutral-stability line marks the threshold to prompt collapse.

For rigidly rotating NS, the $J$-constant turning points coincide with the
$M_{0}$-constant turning points, but it is no longer the case for differential
rotation of a given degree. While other authors usually refer to the former,
in this paper we use the latter as we find it to be a closer estimate of
stability threshold.

On secular timescales, differential rotation transforms into rigid rotation
due to effects of viscosity and magnetic breaking [16, 6, 17]. By definition,
hypermassive NS have masses that cannot be supported by rigid rotation only.
This eventually may lead to a delayed collapse and delayed emission of
gravitational waves. There is no clear criterion of dynamical stability for
hypermassive NS to tell if the collapse will be prompt or delayed.

Various authors have studied the stability properties of differentially
rotating NS by means of numerical simulations. An example of hypermassive NS
dynamically stable to both radial instabilities and bar formation was
presented by [3]. In [24] the authors explore the limit of stability to quasi-
radial oscillations for differentially rotating NS, excluding quasi-toroidal
configurations. The threshold to collapse proves to be close to the
($J$-constant) turning-point line, which is still a valid sufficient criterion
of dynamical instability. A caveat for the large masses supported by many
these works is that they may be subject to non-axisymmetric corrotational
instabilities (usually known as low-$T/|W|$ instabilities, see e.g. [18]) that
are able to transport efficiently angular momentum and erase differential
rotation. Although there is no clear criterion for the onset of this
instability, all studied cases in the literature of NS with quasi-toroidal
shape (e.g. [8]) have shown the dynamical growth of these instabilities.

##  2 Equilibrium models

We consider axisymmetric, stationary configuration of rotating fluid in
cylindrical coordinates $(t,\rho,z,\phi)$. The configurations we study are
highly flattened and cylindrical coordinates are more practical than spherical
ones. The line element associated with such configuration may be written in
form:

$ds^{2}=-e^{2\nu}dt^{2}+e^{2\mu}(d\rho^{2}+dz^{2})+W^{2}e^{-2\nu}(d\phi-\omega
dt)^{2}\text{,}$ (1)

with four metric potentials $\mu,W,\nu,\omega$, being functions of $\rho$ and
$z$ only, due to symmetry.

In general, the properties of matter are defined by the equation of state
(EOS). Here we use a polytropic EOS with $\Gamma=2$ (or $N=1$ in alternate
notation) which yields a relation between total mass-energy density $\epsilon$
and pressure $p$:

$\epsilon(p)=p+\sqrt{\frac{p}{K}}\text{,}$ (2)

where $K$ is the polytropic constant.

We often use a dimensionless value of relativistic enthalpy as a main
thermodynamical parameter, which in case of polytrope with $\Gamma=2$ can be
expressed as:

$H=\text{log}(1+2K\epsilon_{B})\text{,}$ (3)

where $\epsilon_{B}$ is the rest-mass density.

To describe the differential rotation profile, one need to specify the
rotation law. For the fluid four-velocity $u^{\alpha}$ and angular velocity
$\Omega$, it can be defined as function $u^{t}u_{\phi}=F(\Omega)$. Here we use
the J-const rotation law [9]:

$F(\Omega)=A^{2}(\Omega_{c}-\Omega)\text{,}$ (4)

with $\Omega_{c}$ being the angular velocity on the rotation axis and $A$ a
constant describing the steepness of the profile. This produces a rotation
profile consistent with remnants of CCSN. As a dimensionless measure of the
degree of differential rotation, we use the value of
$\widetilde{A}=\frac{r_{e}}{A}$ with $r_{e}$ being the star radius on
equatorial plane.

To construct the initial data we solve the relativistic field equations for
four metric potentials on the $\rho,z$ plane using the highly accurate code
FlatStar. It uses an efficient multi-domain spectral method to construct
equilibrium models of rotating compact objects. For more technical details see
[2] and appendix of [10].

In this paper, we focus on the degree of differential rotation
$\widetilde{A}=1$. According to classification of [2], all configurations
studied here are of type C. Our main interest lies in quasi-toroidal
configurations, which produce the largest masses and are not extensively
studied by other authors. Figure 1 shows the relativistic enthalpy profile of
one of such configurations. The maximal value, $H_{\max}$, is not in the
center of star, and hence the classification as quasi-toroidal. We select 5
sequences with constant rest mass of configurations close to the stability
limit estimated by the turning point criterion.

Figure 1: Example of a quasi-toroidal initial configuration ($H_{\max}=0.26$,
$M_{0}=0.33$). Color coded, the relativistic enthalpy $H$ in a meridional
cross section. The maximal value $H_{\max}$ is found far off the rotation axis
($y=0$).

##  3 Stability against quasi-radial instabilities

To test the stability of the selected configurations, we perform relativistic
axisymmetric hydrodynamical simulations using the CoCoNuT code [4], which uses
the conformally flat approximation (CFC). In the 3+1 split, the line element
reads:

$ds^{2}=-\alpha^{2}dt^{2}+\gamma_{ij}(dx^{i}+\beta^{i}dt)(dx^{j}+\beta^{j}dt),$
(5)

where $\gamma_{ij}=\Phi^{4}\delta_{ij}$ in the CFC approximation, with the
conformal factor $\Phi^{6}=e^{2\mu-u}$ (using variables from equation 1). The
accuracy of CFC was tested, for example, by [5], showing that the
astrophysical properties (such as mass) are reconstructed with only a small
discrepancy less than 5%.

To induce the collapse in unstable configurations, we introduce a radial
perturbation in form of an additional velocity component. The amplitude of
this perturbation is carefully chosen to make the mean value of the maximal
density match the initial value in stable solutions (for unstable solutions we
extrapolate the amplitude from the stable region). We test our results by
comparing them with the results of [24] and [21], finding them to be in
agreement.

For all our models we inspect the evolution of maximal density in time. For
stable configurations we see the value of density oscillating around the
initial value. Unstable configurations show an exponential growth of the
maximal density in the first few ms, marking the prompt collapse to a BH.
Figure 2 shows stable and unstable configurations on $H_{\max}-M$ plane.

Figure 2: Simulated configurations divided into stable (green marks) and
unstable (red marks) to quasi-radial perturbations. Blue dashed line shows the
line of ($M_{0}$-constant) turning-points, being the first estimate of
stability. The orange dashed line marks the boundary between spheroidal and
quasi-toroidal configurations. Limit of mass for this degree of differential
rotation, limit for rigid rotation and sequence of non-rotating NS are
presented for reference.

##  4 Summary

We have selected a sample of quasi-toroidal configurations of neutron stars
with polytropic equation of state. We used a polytrope with $\Gamma=2$ and
j-constant rotation law. By performing a numerical relativistic hydrodynamical
evolution, we tested stability of the selected equilibria against
axisymmetrical (quasi-radial) perturbations. We show that differential
rotation allows the existence of dynamically stable models with masses almost
twice as massive as $M_{\rm TOV}$. These stable configurations, if formed
during CCSN or BNS mergers, may undergo a significantly delayed collapse to a
black hole. Further study is needed to inspect the stability properties
against non-axisymmetrical perturbations on dynamical timescales.

##  5 Acknowledgments

This work was partially supported by the Polish National Science Centre grants
No. 2017/26/M/ST9/00978 and 2022/45/N/ST9/04115, by POMOST/2012-6/11 Program
of Foundation for Polish Science co-financed by the European Union within the
European Regional Development Fund, by the Spanish Agencia Estatal de
Investigación (Grants No. PGC2018-095984-B-I00 and PID2021-125485NB-C21)
funded by MCIN/AEI/10.13039/501100011033 and ERDF A way of making Europe, by
the Generalitat Valenciana (PROMETEO/2019/071), and by COST Actions CA16104
and CA16214.

## References

  * [1] F. Acernese and M. Agathos  “Advanced Virgo: a second-generation interferometric gravitational wave detector”  In _Classical and Quantum Gravity_ 32.2, 2015, pp. 024001  DOI: 10.1088/0264-9381/32/2/024001
  * [2] M. Ansorg, D. Gondek-Rosińska and L. Villain  “On the Solution Space of Differentially Rotating Neutron Stars in General Relativity”  In _MNRAS_ 396, 2008  DOI: 10.1111/j.1365-2966.2009.14904.x
  * [3] Thomas W. Baumgarte, Stuart L. Shapiro and Masaru Shibata  “On the Maximum Mass of Differentially Rotating Neutron Stars”  In _ApJ_ 528.1, 2000, pp. L29–L32  DOI: 10.1086/312425
  * [4] P. Cerdá-Durán, J.. Font, L. Antón and E. Müller  “A new general relativistic magnetohydrodynamics code for dynamical spacetimes”  In _A &A_ 492.3, 2008, pp. 937–953  DOI: 10.1051/0004-6361:200810086
  * [5] Gregory B. Cook, Stuart L. Shapiro and Saul A. Teukolsky  “Testing a simplified version of Einstein’s equations for numerical relativity”  In _Phys. Rev. D_ 53.10, 1996, pp. 5533–5540  DOI: 10.1103/PhysRevD.53.5533
  * [6] Matthew D. Duez, Yuk Tung Liu, Stuart L. Shapiro and Branson C. Stephens  “General relativistic hydrodynamics with viscosity: Contraction, catastrophic collapse, and disk formation in hypermassive neutron stars”  In _Phys. Rev. D_ 69.10, 2004, pp. 104030  DOI: 10.1103/PhysRevD.69.104030
  * [7] P.. Espino and V. Paschalidis  “Revisiting the maximum mass of differentially rotating neutron stars in general relativity with realistic equations of state”  In _Phys. Rev. D_ 99.8, 2019, pp. 083017  DOI: 10.1103/PhysRevD.99.083017
  * [8] Pedro L. Espino, Vasileios Paschalidis, Thomas W. Baumgarte and Stuart L. Shapiro  “Dynamical stability of quasitoroidal differentially rotating neutron stars”  In _Phys. Rev. D_ 100.4, 2019, pp. 043014  DOI: 10.1103/PhysRevD.100.043014
  * [9] J.. Friedman, J.. Ipser and R.. Sorkin  “Turning-point method for axisymmetric stability of rotating relativistic stars”  In _ApJ_ 325, 1988, pp. 722–724  DOI: 10.1086/166043
  * [10] D. Gondek-Rosińska et al.  “A New View on the Maximum Mass of Differentially Rotating Neutron Stars”  In _ApJ_ 837, 2017, pp. 58  DOI: 10.3847/1538-4357/aa56c1
  * [11] Panagiotis Iosif and Nikolaos Stergioulas  “Equilibrium sequences of differentially rotating stars with post-merger-like rotational profiles”  In _MNRAS_ 503.1, 2021, pp. 850–866  DOI: 10.1093/mnras/stab392
  * [12] Kagra Collaboration, T. Akutsu and M. Ando  “KAGRA: 2.5 generation interferometric gravitational wave detector”  In _Nature Astronomy_ 3, 2019, pp. 35–40  DOI: 10.1038/s41550-018-0658-y
  * [13] Wolfgang Kastaun and Filippo Galeazzi  “Properties of hypermassive neutron stars formed in mergers of spinning binaries”  In _Phys. Rev. D_ 91 American Physical Society, 2015, pp. 064027  DOI: 10.1103/PhysRevD.91.064027
  * [14] H. Komatsu, Y. Eriguchi and I. Hachisu  “Rapidly rotating general relativistic stars. II - Differentially rotating polytropes”  In _MNRAS_ 239, 1989, pp. 153–171  DOI: 10.1093/mnras/239.1.153
  * [15] LIGO Scientific Collaboration, J. Aasi and B.. Abbott  “Advanced LIGO”  In _Classical and Quantum Gravity_ 32.7, 2015, pp. 074001  DOI: 10.1088/0264-9381/32/7/074001
  * [16] Stuart L. Shapiro  “Differential Rotation in Neutron Stars: Magnetic Braking and Viscous Damping”  In _ApJ_ 544.1, 2000, pp. 397–408  DOI: 10.1086/317209
  * [17] M. Shibata, K. Taniguchi and K. Uryū  “Merger of binary neutron stars with realistic equations of state in full general relativity”  In _Phys. Rev. D_ 71.8, 2005, pp. 084021  DOI: 10.1103/PhysRevD.71.084021
  * [18] Masaru Shibata, Shigeyuki Karino and Yoshiharu Eriguchi  “Dynamical bar-mode instability of differentially rotating stars: effects of equations of state and velocity profiles”  In _MNRAS_ 343.2, 2003, pp. 619–626  DOI: 10.1046/j.1365-8711.2003.06699.x
  * [19] A.. Studzińska et al.  “Effect of the equation of state on the maximum mass of differentially rotating neutron stars”  In _MNRAS_ 463, 2016, pp. 2667–2679  DOI: 10.1093/mnras/stw2152
  * [20] M. Szkudlarek, D. Gondek-Rosińska, L. Villain and M. Ansorg  “Maximum Mass of Differentially Rotating Strange Quark Stars”  In _ApJ_ 879, 2019, pp. 44  DOI: 10.3847/1538-4357/ab1752
  * [21] K. Takami, L. Rezzolla and S. Yoshida  “A quasi-radial stability criterion for rotating relativistic stars”  In _MNRAS: Letters_ 416.1, 2011, pp. L1–L5  DOI: 10.1111/j.1745-3933.2011.01085.x
  * [22] K. Uryu et al.  “Modeling differential rotations of compact stars in equilibriums”  In _Phys. Rev. D_ 96 American Physical Society, 2017, pp. 103011  DOI: 10.1103/PhysRevD.96.103011
  * [23] L. Villain, J.. Pons, P. Cerdá-Durán and E. Gourgoulhon  “Evolutionary sequences of rotating protoneutron stars”  In _A &A_ 418, 2004, pp. 283–294  DOI: 10.1051/0004-6361:20035619
  * [24] L.. Weih, E.. Most and L. Rezzolla  “On the stability and maximum mass of differentially rotating relativistic stars”  In _MNRAS_ 473.1, 2018, pp. L126–L130  DOI: 10.1093/mnrasl/slx178