# Disc Tearing in a Be Star: Predicted 3D Observations

M.W. Suffak,1 C.E. Jones,1 A.C. Carciofi,2
1Department of Physics and Astronomy, Western University, London, ON N6A 3K7,
Canada
2Instituto de Astronomia, Geofísica e Ciéncias Atmosféricas, Universidade de
São Paulo, Brazil
E-mail<EMAIL_ADDRESS>
(Accepted 2023 November 23. Received 2023 November 20; in original form 2023
September 20)

###### Abstract

We build on our previous work involving smoothed particle hydrodynamic
simulations of Be stars, by using the model that exhibited disc tearing as
input into the three-dimensional Monte Carlo radiative transfer code hdust to
predict observables from a variety of viewing angles throughout the disc
tearing process. We run one simulation at the start of each orbital period
from 20 to 72 orbital periods, which covers two complete disc tearing events.
The resulting trends in observables are found to be dependent on the relative
position of the observer and the tearing disc. The H$\rm\alpha$ equivalent
width, $V$ magnitude, and polarization can all increase or decrease in any
combination depending on the viewpoint of the observer. The H$\rm\alpha$ line
profile also displays changes in strength and peak separation throughout the
tearing process. We show how the outer disc of the torn system can have a
large effect on the H$\rm\alpha$ line profile, and also contributes to a
wavelength-dependent polarization position angle, resulting in a similar
sawtooth shape to the polarization percentage. Finally, we compare our
predictions to Pleione (28 Tau) where evidence has suggested that a disc
tearing event has occurred in the past. We find that our tearing disc model
can broadly match the trends seen in Pleione’s observables, as well as produce
the two-component H$\rm\alpha$ lines observed in Pleione. This is the
strongest evidence, thus far, of Pleione’s disc having indeed experienced a
tearing event.

###### keywords:

binaries: general – circumstellar matter – stars: emission-line, Be – stars:
individual, Pleione

††pubyear: 2023††pagerange: Disc Tearing in a Be Star: Predicted 3D
Observations–References

##  1 Introduction

Classical Be stars are rapidly rotating B-type stars that, at some point, have
had Balmer emission lines in their spectrum due to a gaseous disc that has
emanated from the star itself (Collins, 1987; Porter & Rivinius, 2003). The
exact process of this disc feeding is uncertain, however non-radial pulsations
(Baade et al., 2016) coupled with rapid rotation is the mechanism thought to
be likely responsible (for the most recent review, see Rivinius et al., 2013).
The evolution of classical Be star discs is currently best-described by the
viscous decretion disc (VDD) model of Lee et al. (1991) where material is
launched from the star via some unknown mechanism, and then is transported
outwards via viscous effects.

In nearly all modelling efforts involving the VDD model, the Be star disc has
been assumed to be in the equatorial plane of the central star. However, a
number of recent works (Cyr et al., 2017; Brown et al., 2019; Suffak et al.,
2022) have used hydrodynamic models to investigate Be star discs that are
tilted out of the equatorial plane due to the influence of a binary companion
who’s orbit is misaligned to the equatorial plane of the Be star. The most
recent work of Suffak et al. (2023) examined the temperature structure and
observables of Be star discs that are tilted away from the star’s equatorial
plane, and was the first to study tilted Be star discs using radiative
transfer models. Such a scenario may be caused by supernova kicks, leading the
resulting neutron star to have a misaligned orbit to the Be star spin axis
(Martin et al., 2009). The misaligned binary star then provides a torque on
the Be star disc, causing the disc to tilt out of the equatorial plane and
undergo nodal precession about the binary angular momentum vector (Martin et
al., 2011). Additionally, if the Be star is actively losing mass to the disc,
then the disc feels a competing viscous torque trying to keep the disc aligned
with the Be star equator. In special cases, the two competing torques can
cause the disc to undergo tearing (Nixon et al., 2013; Doğan et al., 2015;
Suffak et al., 2022).

Suffak et al. (2022) computed six simulations where the binary companion orbit
was misaligned to the Be star equatorial plane, one of which exhibited the
phenomenon of disc tearing. Over the course of 100 binary orbital periods, the
disc underwent three separate disc tearing events, spaced roughly 30 orbital
periods apart. After the disc tears, the inner and outer disc nodally precess
independently, while the inner disc is still being injected with material from
the Be star. This results in the inner disc growing radially and tilting back
towards the Be star equatorial plane, while the outer disc undergoes large
inclination changes during precession, and gradually dissipates. This process
allows the disc tearing phenomenon to repeat multiple times.

The Be star Pleione (28 Tau, HD23862) has been a focus of recent publications
(Marr et al., 2022; Martin & Lepp, 2022) due to evidence suggesting the disc
of Pleione undergoes episodes of disc tearing, where the disc splits into
separate components: an inner disc and an outer precessing disc. In their
detailed modelling of Pleione, Marr et al. (2022) postulated that the disc
tearing scenario of Suffak et al. (2022) is the best explanation for drops in
the H$\alpha$ equivalent width (EW) and $V$ magnitude of Pleione. Martin &
Lepp (2022) then followed this analytically, finding that a broken disc is the
only way to match the precession period and disc size of Pleione. A two-disc
scenario was also proposed earlier for Pleione by Tanaka et al. (2007). For
other studies of Pleione, see the introduction of Marr et al. (2022).

In this work, we use the Smoothed Particle Hydrodynamics (sph) simulation from
Suffak et al. (2022) that showed disc tearing, as input into the 3-dimensional
non-local thermodynamic equilibrium (3D nLTE) Monte Carlo radiative transfer
code hdust (Carciofi & Bjorkman, 2006). This will produce fully 3D predictions
of observables of the disc before, during, and after a disc tearing episode,
and include a full velocity prescription from the sph code, which is a first
for Be star research. Not only do we compare this to the data from Pleione,
but we show the general trends in observables that disc tearing may produce
when viewed from a complete variety of orientations.

Section 2 of this paper briefly summarizes the details of the disc tearing
model of Suffak et al. (2022), Section 3 discusses the general trends seen in
the disc tearing observables from a variety of vantage points, Section 4 shows
the comparison between the disc tearing observables and Pleione’s variability,
and Section 5 provides our discussion and conclusions.

##  2 Simulation Details

The details about the sph code used by Suffak et al. (2022), as well as more
specifics on the simulation setup, can be found in section 2.1 of their paper.
The simulation initially starts with no disc, and at each subsequent timestep,
equal-mass particles are injected into the Be star’s equatorial plane with
Keplerian velocities to grow the disc. Most of the injected mass falls back
onto the Be star immediately, and a small fraction of the material stays to
form the disc. The orbit of the binary companion is misaligned by $None$ to
this equatorial plane. Both the binary’s orbital plane and the spin axis of
the Be star are fixed throughout the simulation. For convenience, parameters
used in the simulation from Suffak et al. (2022) are listed in Table 1. We
note that the orbital and stellar parameters are not the same as Pleione, as
the goal of this paper is to produce observables of a broken disc system in
general for broad applicability, but later we show similarities to the Pleione
system.

Table 1: Simulation parameters for the sph model from Suffak et al. (2022).

[c]  Parameter Value Be Star Mass 8 $\rm M_{\odot}$ Be Star Radius 5 $\rm
R_{\odot}$ Be Star Effective Temperature 20000 K Binary Star Mass 8 $\rm
M_{\odot}$ Binary Star Radius 5 $\rm R_{\odot}$ Disc Temperature 12000 K a
Viscosity Parameter $\alpha$ 0.1 Mass Injection Rate $10^{-8}\,\rm
M_{\odot}yr^{-1}$ Particle Injection Radius $1.04\,\rm R_{p}$ Binary Orbit
Misalignment Angle $None$ Binary Orbital Period 30 days Binary Orbital Radius
20.5 $\rm R_{p}$

  * a

The disc temperature is set to 60% of the primary star’s effective
temperature. This value was found by Carciofi & Bjorkman (2006) to be the
average temperature in the isothermal regions of the disc.

To adopt input from the sph simulation for hdust, we use an interface that
takes the sph particles at a desired timestep, and creates a computational
grid in a format that hdust accepts. We created the grid in spherical
coordinates with 50 radial cells, 50 azimuthal cells, and 50 polar cells. We
explored a larger number of cells for better resolution, but the differences
were negligible. In a standard grid, cells are spaced logarithmically in the
radial and polar direction, since a typical disc is concentrated nearest the
equator and near the central star. However, because of the complexity of the
disc configurations in this simulation, we have designed the grid such that
the cells are more closely spaced when the disc density is highest in each
vertical column.

For the stellar parameters that hdust requires, we input the source (primary
star) as having a rotation critical fraction, $W$, of 0.7, which is average
for Be stars (Rivinius et al., 2013), a mass of 8 $\rm M_{\odot}$, polar
radius of 4 $\rm R_{\odot}$ (consistent with $W\,=\,0.7$ to give an equatorial
radius of 5 $\rm R_{\odot}$), and a luminosity of 2300 $\rm L_{\odot}$
(consistent with the primary’s temperature in Table 1).

In Suffak et al. (2022), this simulation was run for 100 orbital periods
before the disc was allowed to dissipate. The simulation starts with a single
disc that tilts progressively out of the equatorial plane as it becomes more
massive. When tearing occurs, the inner disc tilts back to the equatorial
plane, while the outer disc starts precessing. What follows is inner disc
growth into a new full disc, while the outer, detached disc, partially
dissipates and perhaps recombines with the inner disc. During these 100
periods, the disc tearing process occurred three times, with the first episode
starting at 25 orbital periods, and then repeating approximately every 30
orbital periods.

##  3 Analyzing the Tearing Disc

###  3.1 Disc Geometry

Figure 1: Schematic of the tearing disc system, showing the inner disc (dark
grey) separated from the outer disc (light grey), both centered around the
primary star (blue). In this frame, the $z$-axis is the spin axis of the Be
star. The azimuthal angle ($\phi$, used for defining the obseving location) is
also shown defined with $\phi=$$$ along the positive $x$-axis. Note that the
relative sizes and inclination of the inner and outer discs to each other and
to the central star would vary over time.

Before discussing observables of the tearing disc produced by hdust, we first
describe our coordinate orientation. A schematic of a torn disc is shown in
Figure 1, where we show an example of the inner and outer components of the
disc when tearing occurs, although it is important to remember that this is
just a schematic, and the inclinations of the two discs change over time. We
define our coordinate system with two angles: the azimuthal angle $\phi$,
which is in the range [$None$,$None$) as shown in Figure 1, and the polar
angle $\theta$, defined from $None$ at the pole of the star, to $None$ at the
equator of the star.

Overall, we compute observables at 152 different observing positions using ten
$\theta$ values, every ten degrees from $None$ (pole-on with the top of the
star) to $None$ (pole-on with the bottom of the star), and eight $\phi$
angles, every 45 degrees between $None$ and $None$. The observables are
computed at the start of every period between 20 and 72 orbital periods, which
covers two full disc tearing episodes, consisting of the disc tearing and then
regrowing. During a disc tearing episode, the disc exists as two independent
pieces: an inner disc whose inner boundary is still anchored at the Be star
equator due to ongoing mass injection, and a detached outer ring, each
inclined at different angles and nodally precessing independently. This means
the projected area and projected velocity the observer sees from both disc
pieces are very different at a given time, and can also change significantly
over time. Visually, we determine the radius at which the disc tears is
approximately 3 $R_{*}$. While mass injection is ongoing, the inner disc will
feel competing torques from the binary as well as the addition of material,
and the outer disc will only feel the binary torque and precess more quickly
than the inner disc. Through analyzing the angular momentum vector of the
inner and outer disc, we can calculate the longitude of the ascending node of
each disc section throughout the simulation. The two discs do not undergo a
full precession period as two separate pieces, so to find the precession
periods of the inner and outer disc, we can look at the slope of the ascending
node versus time when the two discs are completely separated (between 30 and
35 orbital periods), and compare what the period would be if the discs kept
this rate of precession without change. Through this, we find the outer disc
precesses approximately five times faster than the inner disc when the two
discs are independently precessing.

It is not practical to go through the resulting observables for all 80
observing positions that we simulated, so we will focus on a few
representative situations that best describe and explain the broad trends seen
for all observing angles. In our disc tearing models, we chose to examine the
H$\alpha$ line, $V$-band magnitude, $V$-band polarization percentage, and the
corresponding polarization position angle (PA).

###  3.2 Photometry, Polarization, and Equivalent Width

Throughout the simulation, as the disc cycles between tearing and regrowing,
the observables also oscillate. The H$\alpha$ equivalent width (EW) can vary
as much as 1 to 2 nm, the change in $V$ magnitude can be on order of
hundredths or tenths of a magnitude, percent polarization amplitude can be as
much as 1%, and the polarization PA changes by tens of degrees. The exact
amount of change in each observable is dependent on the orientation of the
observer to the star and the tearing disc, as does the way in which the
observables change with respect to each other.

Figure 2 shows examples of this variability where we have plotted the
H$\rm\alpha$ EW, photometry, polarization, and inclination of the inner
($<3\,R_{*}$) and outer ($>3\,R_{*}$) disc relative to the observer looking
pole-on with the star $(\theta\,=\,$$$). Immediately before the disc tears it
is tilted a small amount from the equatorial plane, accounting for the small,
but non-zero, polarization signal. When the disc tears, the projected area of
the disc that an observer would see from this pole-on angle decreases; as the
optically thick H$\rm\alpha$ line emission strength is directly proportional
to this area, the EW drops as disc tearing begins. After the disc tears, as
the outer disc undergoes precession and is viewed more egde-on, the inner disc
slowly realigns with the equatorial plane, which means pole-on observing
angles will see the inner disc more face-on. As the $V$-band excess and
polarization of a Be star disc originate from the inner portion of the disc
(Carciofi, 2011), seeing the inner disc more face-on accounts for the
brightening in $V$ magnitude, and the decrease in percent polarization as the
polarization vectors would largely cancel due to symmetry. The variation in
the PA depends on the direction of the tilt and orientation of the
polarization axes in the simulation. While the inner disc is expected to cause
much of the polarized flux, the outer disc also contributes to this flux,
otherwise we would not see such large PA changes. In this case, the PA varies
by about $None$ as the outer disc undergoes precession.

Figure 2: Top to bottom, plots of the computed H$\rm\alpha$ EW, $V$ magnitude,
percent polarization, polarization PA, and average disc inclination relative
to the observer of the inner ($<3\,R_{*}$) and outer ($>3\,R_{*}$) disc, for a
pole-on observing angle of $\theta\,=\,$$$ and $\phi\,=\,$$$. The $x$-axis is
in units of binary orbital periods. The vertical lines indicate times where
the disc tearing begins (red, dot-dashed line) and where the disc is whole
again (green, dotted line).

As the observer moves to more equator-on views, tearing in the disc can lead
to the opposite trends that are seen for pole-on angles. In many cases, the
$V$ magnitude now decreases along with the H$\alpha$ equivalent width as the
disc tears. The polarization percentage can both increase or decrease
depending on how the inner disc orientation changes with respect to the
observer. In some other cases, the EW or $V$ magnitude can go through multiple
increases and decreases between the disc tearing and regrowing. An example of
the trends when tearing is viewed equator-on is shown in Figure 3.

Figure 3: Same format as Figure 2, but for an observing angle of
$\theta\,=\,$$$ and $\phi\,=\,$$$.

####  3.2.1 Spectropolarimetry

Figure 4: Polarization PA (top), and net polarization percentage (bottom),
versus wavelength, for the tearing disc at 31 $\rm P_{orb}$ (shortly after
disc tearing). Viewed from $\theta\,=\,$$$, $\phi\,=\,$$$.

Figure 5: Left to right, the disc image coloured by pixel polarization angle,
disc polarized flux, and polarization PA versus wavelength, with the blue
highlighted region corresponding to the wavelength the two other panels are
sampled from. The top plot is just before the Balmer jump, and the bottom plot
is just after the Balmer jump. The viewing angle is $\theta\,=\,$$$,
$\phi\,=\,$$$.

The polarization values presented in Figures 2 and 3 are averaged over the $V$
band, so the wavelength-dependent nature of Be star polarization is lost when
computing these values. It is well known that the amount of polarization from
a Be star disc takes on a wavelength dependent sawtooth-shaped pattern, due to
the imprint of varying hydrogen opacity on the polarized light (see Haubois et
al., 2014, for example), with jumps in the polarization occurring at the
various hydrogen series boundaries. In examining spectropolarimetry from the
tearing disc model, we find this behaviour in the polarization percentage,
however, for some observing angles, we also find a sawtooth pattern in the
polarization PA right after disc tearing, when the detached disc maintains its
maximum mass (Figure 4). This behaviour is unexpected and would never be seen
in a normal flat Be star disc since all of the disc material would be in the
same plane, and the large majority of the polarization comes from the very
innermost disc. However, it is not a surprising feature when considering both
inner and outer discs contribute to the polarization.

In Figure 5, the left panel shows the disc image produced by hdust, with each
pixel coloured according to its PA, calculated by

${\rm PA}=\frac{1}{2}\arctan\frac{U}{Q}.$ (1)

The middle panel shows this same disc image, but with the pixels coloured by
polarized flux, and the rightmost panel highlights the wavelength range used
in the left and middle images. The middle panel is essentially a map of the
disc density. In comparing the left and middle panel, we can see that the
outer disc has a much different polarization PA than the inner disc due to
their opposing orientations. The middle panel also allows us to see that right
before the Balmer jump (top row in Figure 5), when the PA is at a negative
maximum, there is much less scattered light from the inner disc than right
after the Balmer jump (bottom row in Figure 5), while the amount of scattered
light from the outer disc remains the same. This means the outer disc is
contributing a larger proportion of the total scattered light before the
Balmer jump, enough to significantly change the PA from its value when the
inner disc dominates. Incidentally, Figure 5 demonstrates the statement made
above that the inner disc contributes more to the polarized flux than the
outer disc, but also shows that the contribution of the latter is not
negligible.

###  3.3 The H$\alpha$ Line

We also examined the H$\alpha$ line profile, with respect to how the disc
tearing process can change its shape. Overall the H$\alpha$ line shape keeps a
largely “normal" double-peaked profile, however the line strength and peak
separation can change as the disc tilts and precesses before/during/after the
tearing process, depending on the orientation relative to the observer. Figure
6 shows examples of this for a pole-on and equator-on viewing angle, by
plotting the dynamical H$\alpha$ spectrum across an entire episode of disc
tearing and recombination. This Figure shows how not only the line strength,
but also the peak separation and strength/width of the central reversal, can
vary throughout the disc tearing process as different portions of the Doppler
shifted disc material can be viewed from the observer’s vantage point.

Figure 6: (Top) Dynamical H$\rm\alpha$ spectra of the tearing disc model from
55 to 72 orbital periods.(Bottom) H$\rm\alpha$ line profiles used to make the
top plot. The colourbar indicates the normalized flux at each velocity along
the x-axis. A velocity of zero indicates H$\rm\alpha$ line center. (Left) The
observing angle is $\theta\,=\,$$$ and $\phi\,=\,$$$. (Right) The observing
angle is $\theta\,=\,$$$ and $\phi\,=\,$$$.

For particular observing orientations, at certain times, the H$\alpha$ line
presents a “two-component" profile, where the line is clearly composed of a
broad and narrow component, that wouldn’t be seen in an ordinary disc. Figure
7 shows an example of one of these two-component lines on the left, with the
corresponding disc image produced by hdust in the H$\alpha$ velocity bin that
is highlighted in blue on the line profile in the left panel of the Figure.
From this, we can see that the inner double peaks in the line originate from
the inner disc, while fast moving wings of the two-component line are due to
the outer ring that has separated from the inner disc during the tearing
process. This is consistent with the geometry of the disc from this observing
angle, as the inner disc is close to face-on and thus has lower projected
velocities than the outer ring which is edge-on with the observer. So we can
say that the line profile is in fact a combination of two profiles, one
produced by the inner disc, and one from the outer disc. The visible
separation between the broad and narrow components in the line profile are
strongest just after the disc tears, and tend to vanish as the disc regrows
and the outer disc partially dissipates and perhaps recombines with the rest
of the disc.

Figure 7: A two-component H$\rm\alpha$ line (left) along with the
corresponding hdust image (right), showing the emission in the velocity bin
highlighted in blue on the line profile. The top plot is for the red peak of
the line, while the bottom plot is for the red wing of the line as indicated
on the left panels.

###  3.4 Role of the Outer Disc

The conclusion in the previous paragraph naturally raises the question: how
much does the outer disc contribute to the H$\rm\alpha$ emission line, as well
as the rest of the resulting spectrum? To test this, we compute an hdust
simulation at 65 orbital periods - when the disc is completely torn into two
parts - and remove the outer disc from the computational grid, meaning this
simulation will produce observables from only the star and inner disc.

Figure 8: H$\rm\alpha$ line computed with hdust, viewed directly pole-on with
the star and at 65 orbital periods, of the entire system (blue) and with the
outer disc removed (orange). Figure 9: Visible and near-infrared spectrum
computed with hdust, viewed directly pole-on with the star and at 65 orbital
periods, of the entire system (blue) and with the outer disc removed (light
blue).

The results of this test are shown in Figures 8 and 9, where we show, for a
pole-on view, a comparison of the H$\rm\alpha$ line, as well as the visible
and near-infrared spectrum, between the simulation with the entire torn disc
system, and simulation containing just the inner disc and central star.

When analyzing the results of the H$\rm\alpha$ line, one must remember that
the shape of the line is largely due to the relative velocities between the
disc and the observer, which will vary with changing inclination of the disc
to the observer. In Figure 8, we can see that removing the outer disc has a
large effect on the emission line, removing the wings of the line entirely,
while the inner double-peaks of the line remain approximately unchanged. This
reaffirms our previous findings that, from this observing angle, the outer
disc was responsible for the higher velocity wings of the line, while the
inner disc produces the inner double-peaked line structure. It is important to
understand that for other observing angles, the reverse may be true and the
outer disc may be responsible for the lower velocity portions of the line, and
the inner disc responsible for the wings, depending on the vantage point.

In Figure 9, we see that removing the outer disc has very little effect on the
visible and near-infrared spectrum. The same is seen in the polarization
signature as well, with the exception that the PA no longer has the sawtooth
shaped wavelength dependence. This is expected as the visible and near-
infrared excess of Be star discs, as well as the polarization, is produced
within the first few stellar radii of the disc, nearest the star and where the
disc is densest (Rivinius et al., 2013; Haubois et al., 2012; Haubois et al.,
2014).

##  4 Comparison of the Disc Tearing Model to Pleione

The data on Pleione collected by Marr et al. (2022) spans a timeframe of
decades and includes $V$-band photometry, H$\rm\alpha$ EW, percent
polarization, as well as polarization PA. The authors also published
H$\rm\alpha$ line profiles from 2005 to 2019. To facilitate comparison between
our simulation results and the data of Pleione, we obtained the H$\rm\alpha$
spectra from the online journal version of Marr et al. (2022), and secured the
data presented in their Figure 14 (private communication).

Pleione shows a sharp drop in H$\rm\alpha$ EW and $V$ magnitude, while the
polarization percentage is beginning to increase, around 2007 (MJD $\approx$
54000). This is the behaviour that Marr et al. (2022) was unable to match with
conventional disc models, and then hypothesized that disc tearing was
occurring.

As stated in Section 3, for some observing angles close to equator-on with the
star, we see these exact trends of EW and $V$ magnitude decreasing while
polarization increases. Figure 10 shows a side-by-side comparison of Pleione’s
data from Marr et al. (2022) with the observables of our disc tearing model
computed at $\theta\,=\,$$$ and $\phi\,=\,$$$. While the drops in H$\rm\alpha$
EW and $V$ magnitude from the disc tearing model are not as steep as those in
Pleione, the overall trends are the same, and in the case of EW, the trend is
remarkably similar before this sharp drop as well.

Other observing angles in our simulation are able to match some features in
Pleione as well. We show a second example in Figure 11, viewed from the same
polar angle as Figure 10, but from $\phi\,=\,$$$. While the decrease in $V$
magnitude and H$\rm\alpha$ EW do not occur simultanesouly here, the system
shows a much steeper drop in photometry and polarization. The trends in EW and
PA do not differ significantly from Figure 10 and still match those seen in
Pleione. The viewing angle of $\theta\,=\,$$$ corresponds to the same stellar
inclination that Marr et al. (2022) identified for Pleione. It is worth noting
that these trends arealso observed for polar angles greater than $None$,
spanning multiple azimuthal angles.

Figure 10: H$\rm\alpha$ EW, $V$ magnitude, percent polarization, and
polarization PA data of Pleione from Marr et al. (2022) (left, top to bottom),
and predictions from our disc tearing model from 20 to 72 orbital periods
viewed from $\theta\,=\,$$$ and $\phi\,=\,$$$. Figure 11: Same as Figure 10,
but the model is viewed from $\theta\,=\,$$$ and $\phi\,=\,$$$.

From 2005 to 2019, the H$\rm\alpha$ line in Pleione transitioned from strongly
double-peaked in 2005 to a shell-line post-2007. In the time between these two
profiles, the line displays the same two-component behaviour that we
discovered in our simulations, as discussed in Section 3.3. Figure 12 shows
the H$\rm\alpha$ line in our simulation seen from $\theta\,=\,$$$,
$\phi\,=\,$$$, at a time when the disc is beginning to tear, overlayed with
Pleione’s H$\rm\alpha$ line from 2006 from Marr et al. (2022). We see that the
shape of the lines are qualitatively the same, with slight differences in the
peak heights and width of the line. In our simulation, the lines are expected
to be wider than those of Pleione due to a more massive star than Pleione
adopted in this work. This means the material around the star is rotating much
faster in our simulation than it would around Pleione.

Figure 12: Comparison between one of Pleione’s H$\rm\alpha$ lines from Marr et
al. (2022) that shows the two-component line profile, and one such line from
our disc tearing model, taken at 60 $\rm P_{orb}$, viewed from $\theta\,=\,$$$
and $\phi\,=\,$$$.

##  5 Discussion and Conclusions

In this work, we used the disc tearing sph simulation from Suffak et al.
(2022) as input into the radiative transfer code hdust to examine the trends
in observables that a tearing disc would produce. This is the first time a
true 3D sph and radiative transfer model has been used to examine observations
while a disc tears. We produced simulated observables for the disc tearing
model at the start of every orbital period, from 20 to 72 orbital periods and
calculated the $V$-band magnitude, percent polarization, polarization PA, and
the H$\rm\alpha$ EW, in addition to analyzing the H$\rm\alpha$ line shape.

Much like the findings of Suffak et al. (2023), who examined the temperature
and observables of flat (non-warped) tilted discs, the variation in the
observables of a tearing disc are entirely dependent on the orientation of the
disc with respect to the observer. When the disc tears, any of the observables
we examined could increase or decrease with some correlated or anti-correlated
with other observables (see Figures 2 and 3). However, we do find that usually
the start of the disc tearing process corresponds to a maximum or minimum in
all observables we examined. Furthermore, for a majority of the observing
angles, the H$\rm\alpha$ EW seems to peak near the start of disc tearing, and
then decreases. This is due to most observing angles seeing a smaller
projected area of the disc when the disc tears. With respect to the
H$\rm\alpha$ line, we find throughout the disc tearing process the line
strength and peak separation can change significantly. This is expected, as
tearing in the disc changes the projected area and relative velocities of the
disc no matter what angle the disc is being viewed from.

We also find that the H$\rm\alpha$ line can present a “two-component" line
profile, where there is a clear separation between the broad line wings and
the narrow inner double-peaked structure. We find that the outer disc can
significantly contributes to the wings of the line in these cases (see Figure
8), which shows how mapping each line component to a particular disc portion
is difficult without a good grasp of the ensuing disc geometry. The outer disc
is can also affect the polarization PA near the hydrogen series limits,
causing the PA to have a sawtooth shape wavelength dependence, however it
minimally contributes to the continuum flux and polarization magnitude. This
is consistent with past publications of Be stars showing that the inner most
area of a Be star disc is responsible for continuum and polarized emission,
since this is usually where there is a large amount of UV radiation and
electrons for scattering due to this region having the largest density.

The outer disc also does not contribute equally to all emission lines produced
by the disc. It is well known that emission lines of Be stars have different
emitting regions (see figure 3 of Wisniewski et al., 2007, for example), with
the H$\rm\alpha$ line having one of the largest emitting volumes (Rivinius et
al., 2013). Because many emission lines form close to the star (especially IR
lines), disc tearing is not expected to effect all lines equally, and many
will not show the same “two component" profile seen in H$\rm\alpha$ . This
could be useful information, allowing the radius at which the disc is tearing
to be predicted. The geometry of disc warping and tearing in combination with
different emitting volumes may also cause some emission lines to be in a shell
phase while others may be in a normal Be-phase. This is certainly something
that deserves a thorough investigation, but is beyond the scope of this study.

We then examined the similarities between our predictions and observations of
the Be star Pleione (28 Tau). Marr et al. (2022) proposed that the long-term
behaviour in Pleione could be due to disc tearing, but at the time did not
have the 3D tools developed to investigate this hypothesis. In comparing our
model to Pleione, we find that the observables produced by our disc tearing
model viewed from $\theta\,=\,$$$, $\phi\,=\,$$$, can roughly match the trends
that have been observed in Pleione. From Figure 10, we can see that the
H$\rm\alpha$ EW and $V$ magnitude drop at approximately the same time in our
models, as they do in Pleione, while the percent polarization reaches a
minimum and begins to increase, and the PA has peaked, beginning to decrease.
The drops that we see from this perspective in the EW and photometry are not
seemingly as sharp as those observed in Pleione, however the observing angle
$\theta\,=\,$$$ and $\phi\,=\,$$$, seen in Figure 11, can produce these
steeper drops in $V$ magnitude and polarization, despite the photometry and EW
not decreasing at the same time. In examining the H$\rm\alpha$ line, we also
find a similar two-component emission line that Pleione has (Figure 12), at
$\theta\,=\,$$$ and $\phi\,=\,$$$. This inclination is significant as Marr et
al. (2022) determined the stellar inclination of Pleione to be $None$ through
modelling Pleione when it did not have a disc.

We caution that the sph simulation parameters and stellar parameters are not
the same as the Pleione system; however, the general qualitative agreement
does add additional support to the hypothesis that Pleione undergoes dis
tearing events. The H$\rm\alpha$ line and other observables will change with
different model parameters such as the mass ratio, mass-loss rate, the
viscosity parameter $\alpha$, and the binary orbit inclination, all of which
also change the tearing radius in the disk (see equation 8 of Doğan et al.,
2015). So there are likely many combinations of parameters that could produce
similar observable trends and line shape as Pleione at a consistent
inclination. We also simulated systems with the same parameters but a lower
mass ratio than the one presented here, and find disc tearing does occur for a
mass ratio of a half, with even lower mass ratios still to be investigated.

The results of this paper present the first study at the observables of a
tearing Be star disc, using a 3D sph and radiative transfer model in
combination. The trends presented here will be invaluable in providing a
baseline for evidence of disc tearing in binary systems with discs. Our
example of this using the Be star Pleione provides the most conclusive
evidence to date that disc tearing is occurring in Pleione, despite different
simulation parameters used in this work. Our study opens the door to
understand a broad range of astrophysical objects with discs for the first
time, allowing complex disk systems to be investigated with fully 3D modelling
techniques, and is applicable to other investigations where circumstellar
discs tilt, become unstable, and tear.

## Acknowledgements

The authors thank the anonymous referee for their comments and suggestions
which helped improve the paper. C.E.J. acknowledges support through the
National Science and Engineering Research Council of Canada. A. C. C.
acknowledges support from CNPq (grant 311446/2019-1) and FAPESP (grants
2018/04055-8 and 2019/13354-1). This work was made possible through the use of
the Shared Hierarchical Academic Research Computing Network (SHARCNET).

## Data Availability

No new data was generated in this work.

## References

  * Baade et al. (2016) Baade D., et al., 2016, A&A, 588, A56
  * Brown et al. (2019) Brown R. O., Coe M. J., Ho W. C. G., Okazaki A. T., 2019, MNRAS, 488, 387
  * Carciofi (2011) Carciofi A. C., 2011, in Neiner C., Wade G., Meynet G., Peters G., eds, IAU Symposium Vol. 272, Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits. pp 325–336 (arXiv:1009.3969), doi:10.1017/S1743921311010738
  * Carciofi & Bjorkman (2006) Carciofi A. C., Bjorkman J. E., 2006, ApJ, 639, 1081
  * Collins (1987) Collins George W. I., 1987, in Slettebak A., Snow T. P., eds, IAU Colloq. 92: Physics of Be Stars. p. 3
  * Cyr et al. (2017) Cyr I. H., Jones C. E., Panoglou D., Carciofi A. C., Okazaki A. T., 2017, MNRAS, 471, 596
  * Doğan et al. (2015) Doğan S., Nixon C., King A., Price D. J., 2015, MNRAS, 449, 1251
  * Haubois et al. (2012) Haubois X., Carciofi A. C., Rivinius T., Okazaki A. T., Bjorkman J. E., 2012, ApJ, 756, 156
  * Haubois et al. (2014) Haubois X., Mota B. C., Carciofi A. C., Draper Z. H., Wisniewski J. P., Bednarski D., Rivinius T., 2014, ApJ, 785, 12
  * Lee et al. (1991) Lee U., Osaki Y., Saio H., 1991, MNRAS, 250, 432
  * Marr et al. (2022) Marr K. C., Jones C. E., Tycner C., Carciofi A. C., Silva A. C. F., 2022, ApJ, 928, 145
  * Martin & Lepp (2022) Martin R. G., Lepp S., 2022, MNRAS, 516, L86
  * Martin et al. (2009) Martin R. G., Tout C. A., Pringle J. E., 2009, MNRAS, 397, 1563
  * Martin et al. (2011) Martin R. G., Pringle J. E., Tout C. A., Lubow S. H., 2011, MNRAS, 416, 2827
  * Nixon et al. (2013) Nixon C., King A., Price D., 2013, MNRAS, 434, 1946
  * Porter & Rivinius (2003) Porter J. M., Rivinius T., 2003, PASP, 115, 1153
  * Rivinius et al. (2013) Rivinius T., Carciofi A. C., Martayan C., 2013, A&ARv, 21, 69
  * Suffak et al. (2022) Suffak M., Jones C. E., Carciofi A. C., 2022, MNRAS, 509, 931
  * Suffak et al. (2023) Suffak M., Jones C. E., Carciofi A. C., de Amorim T. H., 2023, Non-LTE Monte Carlo Radiative Transfer. III. The thermal properties of Tilted and Warped Be Star Discs, Submitted
  * Tanaka et al. (2007) Tanaka K., Sadakane K., Narusawa S.-Y., Naito H., Kambe E., Katahira J.-I., Hirata R., 2007, PASJ, 59, L35
  * Wisniewski et al. (2007) Wisniewski J. P., Kowalski A. F., Bjorkman K. S., Bjorkman J. E., Carciofi A. C., 2007, ApJ, 656, L21


# Responsibility-associated Multi-agent Collision Avoidance
with Social Preferences

Yiwei Lyu1, Wenhao Luo2 and John M. Dolan3 ∗This work was supported in part by
the CMU Argo AI Center for Autonomous Vehicle Research and the Faculty
Research Grant award at UNC Charlotte.1The author is with the Department of
Electrical and Computer Engineering, Carnegie Mellon University. Email:
<EMAIL_ADDRESS>author is with the Department of Computer Science,
University of North Carolina at Charlotte. Email<EMAIL_ADDRESS>author is with the Robotics Institute, Carnegie Mellon University. Email:
<EMAIL_ADDRESS>
###### Abstract

This paper introduces a novel social preference-aware decentralized safe
control framework to address the responsibility allocation problem in multi-
agent collision avoidance. Considering that agents do not necessarily
cooperate in symmetric ways, this paper focuses on semi-cooperative behavior
among heterogeneous agents with varying cooperation levels. Drawing upon the
idea of Social Value Orientation (SVO) for quantifying the individual
selfishness, we propose a novel concept of Responsibility-associated Social
Value Orientation (R-SVO) to express the intended relative social implications
between pairwise agents. This is used to redefine each agent’s social
preferences or personalities in terms of corresponding responsibility shares
in contributing to the coordination scenario, such as semi-cooperative
collision avoidance where all agents interact in an asymmetric way. By
incorporating such relative social implications through proposed Local
Pairwise Responsibility Weights, we develop a Responsibility-associated
Control Barrier Function-based safe control framework for individual agents,
and multi-agent collision avoidance is achieved with formally provable safety
guarantees. Simulations are provided to demonstrate the effectiveness and
efficiency of the proposed framework in several multi-agent navigation tasks,
such as a position-swapping game, a self-driving car highway ramp merging
scenario, and a circular position swapping game.

##  I Introduction

Multi-agent systems have been extensively studied in recent years given its
potential in modelling and handling complex tasks through cooperative
behaviors among agents, e.g. search and rescue [1], environmental exploration
[2], and coordinated autonomous driving [3, 4]. One critical challenge for the
success of multi-agent interaction often lies in the safety assurance when
producing the inter-agent behaviors. While this has been studied in terms of
safe control and collision avoidance for autonomous systems such as multi-
robot systems [5, 6, 7, 8, 4], it is often assumed that the interacting agents
in a collision avoidance scenario such as that shown in Fig. 1, are either
fully cooperative in a symmetric manner or non-cooperative as passive moving
obstacles [9, 8, 10, 11, 4]. Considering the possible varying coordination
levels among agents in a real-world scenario, e.g. multi-vehicle traffic with
egoistic and altruistic drivers showing various social preferences in
participating in the collision avoidance, it is desired to (a) model such
heterogeneity in the agents’ collision avoidance behaviors and (b) design the
decentralized safe control framework accordingly for autonomous agents to
accommodate the possible heterogeneity and ensure safety assurance during
interaction.

Figure 1:  A multi-agent semi-cooperative task where all agents share
different social preferences or personalities. How can they achieve efficient
and safe interaction?

Collision avoidance for distributed multi-agent systems often relies on
individual robot controller to constrain their own motion based on the assumed
behavior models of the surrounding robots and environments. For example,
reactive collision avoidance methods such as reciprocal velocity obstacles
[12, 13], safety barrier certificates [14, 11], and Buffered Voronoi Cells
[15] have been extensively used to model the collision-free multi-robot
behaviors by characterizing the admissible control space or local collision-
free position space for each robot. The pre-computed safety guarantee is
achieved with the assumption of fully cooperative neighboring robots under the
same control design or passive obstacles moving with piece-wise constant
velocity known to each robot. Although the heterogeneity of individual agent’s
behavior can be described by varying capability-related parameters such as
safety radius, actuation limits, etc. As discussed in [16], it is still
challenging to model heterogeneous behavior patterns such as varying social
preferences/personalities that generalize different levels of conservativeness
and coordination among interacting robots in collision avoidance scenario.

To incorporate the social preferences or personalities information into a
multi-agent control framework, some researchers draw ideas from sociology and
psychology, for instance, the concept of Social Value Orientation (SVO) [17,
18] as quantified selfishness to characterize egoistic and altruistic agents.
These are promising potential solutions, as concepts like SVO provide a
quantitative way to measure the agent choice over weighting its own rewards
against the others in social dilemmas. There has been growing interest in
leveraging these concepts to design heterogeneous multi-agent interaction
mechanisms, e.g. [19, 20, 21].

Motivated by these considerations, we focus on safe interaction in
heterogeneous multi-agent systems. In order to model the influence of
individual social preference or personality on the interactive behaviors, e.g.
collision avoidance between pairwise agents, we propose Responsibility-
associated Social Value Orientation (R-SVO) as a novel variant of traditional
SVO. This enables a unified representation to describe the relative social
implications among agents during interaction. For example, when operating in
close proximity, how likely are agents with various absolute individual
selfishness measures to yield to each other for collision avoidance? With
R-SVO defining such relative responsibility allocation among agents, we then
propose a Control Barrier Function (CBF)-based computational framework for
decentralized safe control that encodes R-SVO in terms of shared
responsibility into the collision-free behavior design. Such a framework is
able to accommodate the various quantified social preferences of individual
agents by allocating responsibility shares accordingly when enforcing
reciprocal safe constraints. The framework is proved to be valid in ensuring
overall safety for the multi-agent system.

Our main contributions are: 1) Drawing upon the idea of Social Value
Orientation (SVO), we propose a novel concept of Responsibility-associated
Social Value Orientation (R-SVO) to redefine the relative difference in social
preference or personality among agents while preserving the intended social
implications and other properties inherited from SVO; 2) we use Local Pairwise
Responsibility Weights to encode agents’ social preference in responsibility
allocation in multi-agent collision avoidance; and 3) we present a
decentralized safe control framework for individual agents to ensure provably
correct collision-free interaction in multi-agent systems. Simulations on
multi-agent position-swapping games are provided to demonstrate the advantages
of the proposed framework in terms of saving overall task completion time and
deadlock resolution efficiency.

##  II Related Work

Figure 2:  Concept of Social Value Orientation in angular representation [19].
The angular value ranges from zero to $\frac{\pi}{2}$, representing the
changes in weight of how an agent value its own reward against the rewards of
others.

Safe control in terms of collision avoidance is critical for multi-agent
systems. To minimize the deviation from robots’ primary task execution due to
safety considerations, reactive collision avoidance methods such as reciprocal
velocity obstacles [5, 6], safety barrier certificates [11, 16, 8, 22] and
buffered Voronoi cells [7, 23] have been presented to minimally revise the
robot’s task-related controller subject to collision avoidance constraints.
Among these methods, Control Barrier Function (CBF) [10] based approaches such
as safety barrier certificates [11, 8] become increasingly popular due to the
advantages of rendering larger admissible control space and addressing
nonlinear agent dynamics with formally provable safety guarantees. In multi-
robot collision avoidance, the safety barrier certificates define an
admissible safety region around individual robot’s controller from which the
resulting multi-robot behaviors are collision-free. To consider collision
avoidance with heterogeneous agents, [16] proposed a distributed safe control
framework using CBF to address robots with varying dynamics in terms of safety
radius, acceleration limits, etc, yet it is difficult to design individual
robot controller reflecting the varying degrees of coordination in multi-robot
interaction, e.g. how to determine the proper control space of an altruistic
agent when interacting with an egoistic agent for safety guarantee.

To explore interaction among heterogeneous agents with different social
preferences or personalities, some works borrow ideas from social psychology
like Social Value Orientation (SVO) [17, 18, 24] to quantify the effect of
agent personalities on behavior modeling. SVO is usually represented in
angular notation $\phi$ with different values corresponding to various
personalities, as shown in Fig. 2. The angular value $\phi$ ranges from zero
to $\frac{\pi}{2}$, corresponding to different personalities ranging from
egoistic ($\phi\in[0,\frac{\pi}{4}]$) to prosocial ($\phi=\frac{\pi}{4}$) and
then to altruistic ($\phi\in(\frac{\pi}{4},\frac{\pi}{2}]$). [20] presents a
weighted buffered Voronoi Cell method which incorporates SVO as a linear
weight into the definition of the Voronoi tesselation. This enables egoistic
agents to have a larger action space than less egoistic agents to avoid
collisions. However, instead of analyzing how to decide the appropriate SVO
value for different agents, SVO is simply assigned to agents directly for use
without any justification. A more principled way to construct the appropriate
SVO for heterogeneous agents in multi-agent systems is needed. [19] estimates
the SVO values of human drivers online as the coefficients in a human decision
model to predict their future trajectories. However, the SVO values estimated
from human drivers are only used as a way to express human intentions to the
autonomous vehicle. In this paper, we introduce a novel concept called
Responsibility-associated SVO (R-SVO) to measure the relative degree of
personality differences among agents, and demonstrate how it can be
incorporated into a CBF-based safe controller for semi-cooperative multi-agent
interaction. To the best of our knowledge, this is the first paper
incorporating agent social preferences in CBF-based safe control frameworks
and leveraging the measurement of relative personality difference among agents
to contribute towards the responsibility allocation in multi-agent collision
avoidance.

The remainder of the paper is organized as follows: Section III reviews the
background on Control Barrier Functions. Section IV presents the proposed
Responsibility-associated Control Barrier Function-based decentralized control
framework. Simulations and discussion of the proposed algorithm in three
examples are provided in Section V.

##  III Background on Control Barrier Functions

Control Barrier Functions (CBF) [10] is a method used to define an admissible
control space for safety assurance of dynamical systems. One of its important
properties is its forward-invariance guarantee of a desired safety set.
Consider the following nonlinear system in control affine form:

$\dot{x}=f(x)+g(x)u$ (1)

where $x\in\mathcal{X}\subset\mathbb{R}^{n}$ and
$u\in\mathcal{U}\subset\mathbb{R}^{m}$ are the system state and control input
with $f$ and $g$ assumed to be locally Lipschitz continuous. A desired safety
set $x\in\mathcal{H}$ can be denoted by a safety function $h(x)$:

$\mathcal{H}=\\{x\in\mathbb{R}^{n}:h(x)\geq 0\\}$ (2)

Thus the control barrier function for the system to remain in the safety set
can be defined as follows [10]:

######  Definition 1.

(Control Barrier Function) Given a dynamical system (1) and the set
$\mathcal{H}$ defined in (2) with a continuously differentiable function
$h:\mathbb{R}^{n}\rightarrow\mathbb{R}$, then $h$ is a control barrier
function (CBF) if there exists a class $\mathcal{K}$ function for all
$x\in\mathcal{X}$ such that

$\sup_{u\in\mathcal{U}}\
\\{L_{f}h(x)+L_{g}h(x)u\\}\geq-\kappa\big{(}h(x)\big{)}$ (3)

where $\dot{h}(x,u)=L_{f}h(x)+L_{g}h(x)u$ with $L_{f}h,L_{g}h$ as the Lie
derivatives of $h$ along the vector fields $f$ and $g$.

A commonly selected class $\mathcal{K}$ function is $\kappa(h(x))=\gamma h(x)$
[10, 25, 26], where $\gamma\in\mathbb{R}^{\geq 0}$ is a CBF design parameter
controlling system behaviors near the boundary of $h(x)=0$. Hence, the
admissible control space in (3) can be redefined as

$\mathcal{B}(x)=\\{u\in\mathcal{U}:\dot{h}(x,u)+\gamma h(x)\geq 0\;\\}$ (4)

It is proved in [10] that any controller $u\in\mathcal{B}(x)$ will render the
safe state set $\mathcal{H}$ forward-invariant, i.e., if the system (1) starts
inside the set $\mathcal{H}$ with $x(t=0)\in\mathcal{H}$, then it implies
$x(t)\in\mathcal{H}$ for all $t>0$ under controller $u\in\mathcal{B}(x)$.

##  IV Method

In this section, we start by introducing the definition and meaning of
Responsibility-associated Social Value Orientation (R-SVO), and its connection
with tradition SVO in terms of properties and social implications. This is the
first step to show how to construct the approriate R-SVO for agents
considering the relative difference in personalities between themselves and
their pairwise agents. Then we present the proposed Responsibility-associated
CBF-based safe control framework, where Local Pairwise Responsibility Weight
is defined to connect from R-SVO to responsibility allocation weights used in
the control framework. Then detailed proofs are provided showing that the
proposed framework can always formally guarantee interaction safety among all
agents. The pipeline of the proposed framework is described in Fig. 3.

Figure 3:  Proposed framework. We propose the concept of Responsibility-
associated Social Value Orientation (R-SVO) to measure the relative difference
in agent personalities. Local Pairwise Responsibility Weight uses the
information from (R-SVO) and is applied as a weight in the CBF-based
decentralized safe controller to solve the responsibility allocation problem
in multi-agent collision avoidance.

###  IV-A Social Value Orientation-based Responsibility Allocation

Consider a multi-agent system with a total number of agents $N\in\mathcal{N}$.
$\theta_{i}\in\mathbb{R}^{\geq 0}$ for $i=\\{1,...,N\\}$ is the global
personality score of the agent $i$, representing the degree of agent egoism.
The smaller $\theta_{i}$ is, more egoistically the agent $i$ behaves, and
similarly, the larger $\theta_{i}$ is, more altruistic the agent $i$ is.

Problem Statement: How to ensure the aforementioned multi-agent system is
collision-free, while taking agent personalities into consideration?

While addressing collision avoidance within a multi-agent system, pairwise
agent safety is considered. To this end, for any pair of agents $i$ and $j$,
it is important to evaluate the appropriate responsibility given their
relative degree of personalities. The motivation is that, although the global
personality scores $\theta_{i}$ and $\theta_{j}$ are available, suppose both
of them share the same large value, meaning both of them are altruistic
agents, it doesn’t add much information when we consider how to allocation
responsibilities between agents in collision avoidance.

In this work, inspired by SVO, we propose to use the Responsibility-associated
Social Value Orientation (R-SVO) to define the relative difference in
personalities between pairwise agents.

######  Definition 2.

For any pair of agents $i$ and $j$, given their global personality scores
$\theta_{i}$ and $\theta_{j}$, the  Responsibility-associated Social Value
Orientation (R-SVO) of each agent $\phi$ is defined as:

$\begin{split}\phi_{i}=\frac{\theta_{i}}{\theta_{i}+\theta_{j}}\cdot\frac{\pi}{2}\quad
s.t.\quad\theta_{i}+\theta_{j}>0\end{split}$ (5)

$\phi_{i}$ measures how the agent $i$ weights its reward against the reward of
agent $j$.

Here, in order to verify that R-SVO $\phi_{i}$ is a legitimate representation
of SVO, we prove 1) $\phi_{i}\in[0,\frac{\pi}{2}]$, and 2) the more egoistic
the agent $i$ is, the smaller $\phi_{i}$ is.

###### Proof.

1) To show R-SVO shares the same bounds of SVO,
$\phi_{i}\in[0,\frac{\pi}{2}]$: $\theta_{i},\theta_{j}\in\mathbb{R}^{\geq 0}$
indicates that $\phi_{i}\in\mathbb{R}^{\geq 0}$. The lower bound of
$\phi_{i}=0$ is achieved when $\theta_{i}=0$. The upper bound of $\phi_{i}=1$
is achieved when $\theta_{j}=0$.

For any $\theta_{i},\theta_{j}\in\mathbb{R}^{>0}$,

$\begin{split}\phi_{i}=\frac{\theta_{i}}{\theta_{i}+\theta_{j}}\cdot\frac{\pi}{2}=\frac{1}{1+\frac{\theta_{j}}{\theta_{i}}}\cdot\frac{\pi}{2}\\\
1+\frac{\theta_{j}}{\theta_{i}}\in[1,+\infty)\rightarrow\phi_{i}\in(0,\frac{\pi}{2}]\end{split}$
(6)

Therefore, we conclude the proof that $\phi_{i}\in[0,\frac{\pi}{2}]$.

2) The definition of R-SVO inherits the same social implication from SVO: The
smaller $\theta_{i}$ is, the larger $1+\frac{\theta_{j}}{\theta_{i}}$ is, the
smaller $\phi_{i}$ is, meaning the more egoistically agent $i$ behaves. If
$\theta_{i}=\theta_{j}$,

$\phi_{i}=\frac{\theta_{i}}{\theta_{i}+\theta_{j}}\cdot\frac{\pi}{2}=\frac{\pi}{4}$
(7)

[17, 18] also verify that $\phi=\frac{\pi}{4}$ represents the prosocial
personality, which is exactly the intended social implication to express. ∎

###  IV-B Responsibility-associated Control Barrier Function-based Safe
Control

In multi-agent systems, agents are usually modeled using a single integrator:
$\dot{x}=u$, where $x\in\mathbb{R}^{2}$ is the position of the agents, and
$u\in\mathbb{R}^{2}$ is the velocity as control input. For any pair of agents
$i$ and $j$ employing Control Barrier Function-based safe controllers, we
consider the particular choice of pairwise agent safety function
$h(x_{i},x_{j})$ and safety set $\mathcal{H}$ as follows.

$\begin{split}\mathcal{H}(x)&=\\{x\in\mathcal{X}:\;h(x_{i},x_{j})=||x_{i}-x_{j}||^{2}-R_{safe}^{2}\geq
0,\forall i\neq j\\}\end{split}$ (8)

where $x_{i},x_{j}\in\mathbb{R}^{2}$ are the positions of agents $i$ and $j$,
and $R_{safe}\in\mathbb{R}$ represents the predefined safety margin. The
admissible control space in Eq. 4 then becomes:

$\mathcal{B}(x)=\\{u_{i},u_{j}\in\mathcal{U}:\dot{h}(x_{i},x_{j},u_{i},u_{j})+\gamma
h(x_{i},x_{j})\geq 0\;\\}$ (9)

Given the agent dynamics equation, we have

$\begin{split}2(x_{i}-x_{j})&{}^{T}(u_{i}-u_{j})+\gamma h(x_{i},x_{j})\geq
0\\\ -2(x_{i}-x_{j})&{}^{T}u_{i}+2(x_{i}-x_{j})^{T}u_{j}\leq\gamma
h(x_{i},x_{j})\end{split}$ (10)

meaning that as long as $u_{i}$ and $u_{j}$ satisfy Eq. 10, then safety
between the agent pair $i$ and $j$ is guaranteed.

Incorporating the Responsibility-associated SVO introduced in the previous
section, we propose a decentralized safe control framework with responsibility
allocation:

######  Theorem 3.

In a multi-agent system, agent safety during an interaction is formally
guaranteed, if for any pair of agents $i$ and $j$, agent $i$ takes the Local
Pairwise Responsibility Weight $\omega_{i}=cos^{2}(\phi_{i})$, so that the
admissible control space in a centralized system in Eq. 10 is converted to:

$\begin{split}\mathcal{B}(x)=&\\{u_{i}\in\mathcal{U}_{i}:A_{i}u_{i}\leq\omega_{i}b_{i},\\\
A_{i}=-2(x_{i}-x_{j})^{T}&\in\mathbb{R}^{1\times 2},b_{i}=\gamma
h(x_{i},x_{j})\in\mathbb{R}\;\\}\end{split}$ (11)

Therefore, for any pair of agents in multi-agent interaction, the
Responsibility-associated CBF-based decentralized safe control framework is
formulated as a quadratic program:

$\begin{split}\min_{u_{i}\in\mathcal{U}_{i}}||u_{i}&-\bar{u}_{i}||^{2}\\\
s.t\quad u_{min}\leq&u_{i}\leq u_{max}\\\ A_{i}u_{i}\leq\omega_{i}b_{i},\quad
A_{i}=-2(x_{i}&-x_{j})^{T},\quad b_{i}=\gamma h(x_{i},x_{j})\end{split}$ (12)

where $\bar{u_{i}}\in\mathbb{R}^{2}$ is the nominal controller input, assumed
to be computed by a higher-level task-related planner, for example, a behavior
planner. $u_{min}$ and $u_{max}$ are the minimum and maximum allowed velocity.

###  IV-C Evaluation of Multi-agent System Safety

In this section, detailed proof is provided to show that the proposed
framework in Th. 11 is the necessary condition of formally provable safety
guarantees. In other words, we aim to show that for any pair of agents $i$ and
$j$, Th. 11 ensures the agents will not collide.

###### Proof.

Step 1: By applying Th. 11 to agent $j$, we get the safety constraint over
$u_{j}$ as:

$A_{j}u_{j}\leq\omega_{j}b_{j},\quad A_{j}=-2(x_{j}-x_{i})^{T},\quad
b_{j}=\gamma h(x_{j},x_{i})$ (13)

We know $A_{j}=-2(x_{j}-x_{i})^{T}=2(x_{i}-x_{j})^{T}$. The summation of the
left hand sides of the inequalities in Th. 11 and Eq. 13 is:

$\begin{split}A_{i}u_{i}+A_{j}u_{j}=-2(x_{i}-x_{j})^{T}u_{i}+2(x_{i}-x_{j})^{T}u_{j}\end{split}$
(14)

which is exactly the left hand side of the inequality in Eq. 10.

Step 2: By Eq. 8, we have $b_{j}=h(x_{j},x_{i})=h(x_{i},x_{j})=b_{i}$, and now
we sum the right hand sides of the two inequalities and get:

$\omega_{i}b_{i}+\omega_{j}b_{j}=(\omega_{i}+\omega_{j})b_{i}=(\cos^{2}(\phi_{i})+\cos^{2}(\phi_{j}))\gamma
h(x_{i},x_{j})$ (15)

As long as we can prove that $(\cos^{2}(\phi_{i})+\cos^{2}(\phi_{j}))\gamma
h(x_{i},x_{j})$ is an equal or a tighter bound compared to $\gamma
h(x_{i},x_{j})$, which is the right hand side of Eq. 10, safety between
pairwise agents during interaction can be formally guaranteed.

Step 3: By Def. 5, we have
$\phi_{i}=\frac{\theta_{i}}{\theta_{i}+\theta_{j}}\cdot\frac{\pi}{2}$ and
$\phi_{j}=\frac{\theta_{j}}{\theta_{i}+\theta_{j}}\cdot\frac{\pi}{2}$, so
$\phi_{j}=\frac{\pi}{2}-\phi_{i}$. Since $\phi_{i}\in[0,\frac{\pi}{2}]$,

$\begin{split}\cos^{2}(\phi_{j})&=\cos^{2}(\frac{\pi}{2}-\phi_{i})=\sin^{2}(\phi_{i})\\\
\Rightarrow\qquad\qquad(\cos^{2}&(\phi_{i})+\cos^{2}(\phi_{j}))\gamma
h(x_{i},x_{j})\\\ =(\cos^{2}(\phi_{i})+&\sin^{2}(\phi_{i}))\gamma
h(x_{i},x_{j})=\gamma h(x_{i},x_{j})\end{split}$ (16)

Thus the proof is concluded that the proposed Responsibility-associated
decentralized CBF in Th. 11 provides an equivalent safety guarantee to that of
the common centralized CBF in Eq. 10. ∎

The algorithm of the proposed framework is presented in Algorithm 1.

Algorithm 1 Responsibility-associated CBF-based
Decentralized Safe Control Framework

0: $\theta_{1,...,N},x^{0}_{1,...,N},\bar{u}^{t}_{1,...,N},\gamma,R_{safe}$

0: $u^{t}_{1,...,N}$

for $t=1:T$ do

for $i=1:N$ do

for $j=1:N$ except $i$ do

Compute Responsibility-associated Social Value Orientation $\phi_{i}$ (Eq. 5)

Compute Local Pairwise Responsibility Weight $\omega_{i}$ (Th. 11)

Calculate pairwise agent safety function $h(x_{i},x_{j})$ Eq. 8

Calculate $A_{i}^{t},b_{i}^{t}$ for safety constraint composition (Eq. 12)

Stack $A^{t}_{i},b^{t}_{i}$ for all surrounding agent $j$

end for

Minimum Deviation Control: $\min_{u^{t}_{i}}||u^{t}_{i}-\bar{u}^{t}_{i}||^{2}$

s.t $\quad u^{t}_{i}\in[u_{min},u_{max}],\quad
A^{t}_{i}u^{t}_{i}\leq\omega_{i}b^{t}_{i}$

end for

end for

The superscript $0$ represents the initial condition at $t=0$, and the
superscript $t$ represents variables at timestep $t$. At each timestep, we
consider all pairwise agents for safety constraint composition. For every
agent $i$ with all its pairwise agents $j$, R-SVO $\phi_{i}$ is calculated
based on the global personality score $\theta_{i}$ and $\theta_{j}$ to measure
the personality difference between agents. Then Local Pairwise Responsibility
Weight $\omega_{i}$ is computed based on R-SVO $\phi_{i}$. Using the Local
Pairwise Responsibility Weight, responsibility allocation is achieved between
the pairwise agents in collision avoidance, and the decentralized safety-
constraint over $u_{i}^{t}$ is composed with state (position) observations of
both agents $x_{i}^{t}$ and $x_{j}^{t}$. The proposed framework scales up well
with a larger number of agents and is highly generally applicable to other
real-time robotics applications.

##  V Simulation & Discussion

In this section, simulations using Matlab are provided to demonstrate the
validity and effectiveness of the proposed method. Three provided examples
include: a two-agent position swapping game, a self-driving car highway ramp
merging problem, and a multi-agent circular position swapping game.

###  V-A Pairwise-agent Interaction: Position Swap Game

A two-agent position swapping game is demonstrated to show how agent
personalities contribute to the pairwise agent responsibility allocation. Both
agents aim to navigate to their goal locations without colliding with each
other while employing a move-to-goal controller
$\bar{u}=-k\cdot(x-x_{target})$, where $k\in\mathbb{R}^{>0}$, and
$x_{target}\in\mathbb{R}^{2}$ is the goal position of each agent. A right hand
heuristic rule is designed to avoid deadlock situations as in [20]. The final
result is shown in Fig. 4. Two sets of trajectories of pairwise agents with
different predefined safety margins $R_{safe}=1,R_{safe}=2$ are shown in two
rows respectively, where the three columns represent three kinds of
responsibility-associated behavior corresponding to agents’ various
personalities. The dashed line represents the trajectory of agent $1$, which
starts from the lower left corner, and the solid line represents the
trajectory of agent $2$, which starts from the upper right corner.

In the first set of experiments (first column), agent $1$ is an egoistic agent
compared to agent $2$ ($\theta_{1}<\theta_{2}$), and the resulting Local
Pairwise Responsibility Weights are $\omega_{1}=0.8,\omega_{2}=0.2$, which
also means that $u_{2}$ is enforced by a tighter bound than $u_{1}$ when
considering collision avoidance. Therefore, the two agents first move straight
towards each others’ initial locations, with agent $1$ moving at a faster
velocity than agent $2$, since it enjoys a relatively flexible bound over
$u_{1}$. Once the two agents come to a point where deadlock is about to
happen, by applying the right hand heuristic rule, they take appropriate
actions to avoid collision.

In the second set of experiments (second column), agent $1$ shares the same
global personality as agent $2$ ($\theta_{1}=\theta_{2}$), leading to Local
Pairwise Responsibility Weights $\omega_{1}=0.5,\omega_{2}=0.5$, so that both
agents share the same portion of responsibility, leading to symmetric
behavior. In the third set of experiments (third column), agent $1$ is an
altruistic agent compared to agent $2$ ($\theta_{1}>\theta_{2}$), thus the
Local Pairwise Responsibility Weights are $\omega_{1}=0.2,\omega_{2}=0.8$.
Benefiting from a more flexible control constraint, Agent $2$ moves at a
higher velocity and proceeds a longer distance than agent $1$ before they
meet. Generally, the responsibility allocation between two agents is easier to
observe in examples with larger safety margin $R_{safe}$. The above examples
show how the proposed method generates various responsibility allocations in
producing safe interaction behaviors based on different agent personalities.

Figure 4:  Final positions of two agents in a position Swapping Game. The
dashed line represents the trajectory of agent $1$, which starts from the
lower left corner, and the solid line represents the trajectory of agent $2$,
which starts from the upper right corner. The colorful circles stands for
predefined safety margin for agents, with $R_{safe}=1$ for the first row and
$R_{safe}=2$ for the second row. The three columns corresponds to three
different Local Pairwise Responsibility Weight settings.

###  V-B Examples in Highway Ramp Merging Scenario

The proposed framework is generally applicable to different robotics
applications which may require coordination with agent social preferences or
personalities, such as autonomous driving, collaborative manufacturing robots,
and other general mobile robots. In Fig. 5, we show a ramp merging scenario
where two vehicles start from their initial location on each lane with the
same distance away from the fixed ramp merging point at $Y=70m$. The goal for
each vehicle is to pass the merging point as soon as possible without any
collision. Three cases are given as examples to demonstrate how the resulting
safe behaviors of the vehicles $V1,V2$ using our proposed decentralized method
align with vehicle’s individual personality/social preference (e.g. egoistic,
altruistic, etc) during interaction.

The system dynamics of a vehicle is described by double integrators as in [4]:

$\begin{split}\dot{X}&=\begin{bmatrix}\dot{x}\\\
\dot{v}\end{bmatrix}=\begin{bmatrix}0_{2\times 2}\;I_{2\times 2}\\\ 0_{2\times
2}\;0_{2\times 2}\end{bmatrix}\begin{bmatrix}x\\\
v\end{bmatrix}+\begin{bmatrix}0_{2\times 2}\;0_{2\times 2}\\\ I_{2\times
2}\;0_{2\times 2}\end{bmatrix}\begin{bmatrix}u\\\ 0\end{bmatrix}\\\
\end{split}$ (17)

where $x\in\mathbb{R}^{2},v\in\mathbb{R}^{2}$ are respectively the position
and linear velocity of each vehicle. $u\in\mathbb{R}^{2}$ represents the
acceleration control input that is determined by our proposed decentralized
responsibility-associated safe controller (Algorithm 1).

Figure 5:  Application of the proposed method in the highway ramp merging
scenario.

In case 1 when two vehicles $V1,V2$ have the same personality/social
preference with the corresponding global personality scores
$\theta_{1}=\theta_{2}$, our Algorithm 1 assigns the same Local Pairwise
Responsibility Weights value $\omega_{1}=\omega_{2}=0.5$ for computing the
decentralized vehicle controller. Hence it degenerates to the homogeneous safe
controller as in Fig. 4 (middle column) and as expected, deadlock situation
happens in this symmetric setting, no matter whether $V1$ and $V2$ are both
egoistic vehicles or altruistic vehicles.

Similar to the example in Fig. 4 (left column), in case 2 in Fig. 5 when $V2$
is altruistic compared to $V1$, $V1$ then has a less restrictive motion in
avoiding collisions due to the assigned smaller R-SVO by Algorithm 1, and thus
allowing it to take more aggressive actions to pass the merging point first
while ensuring safety, which matches the expectation based on the relative
personality between the two. In case 3, although $V1$ is an egoistic vehicle,
$V2$ is even more egoistic than $V1$. Therefore, $V2$ arrives the merging
point first, while $V1$ behaves more conservatively by enforcing a tighter
bound on its admissible action space.

To that end, from the three examples in Fig. 5 it is justified the safe
behavior generated by the proposed controller is well aligned with the agent
social preferences/personalities.

###  V-C Multi-agent Interaction: Circular Game

In this section, we apply the proposed framework in a multi-agent interaction
scenario, where three agent pairs, six agents in total, conduct a circular
game where all agents conduct position swapping with their paired agent while
ensuring no collisions with all other agents. The purpose of this example is
to demonstrate how the proposed method can efficiently solve the deadlock
situation in a crowded environment while remaining collision-free.

In the first example, a symmetric coordination is deployed, meaning no agent
personality information is taken into consideration, to serve as a comparison
baseline to the proposed method. The simulation result is shown in Fig. 6. It
is observed that at $t=140$, a deadlock situation is about to happen. All the
six agents forms a circle, and by applying the right hand heuristic rule, all
agents rotate slowly around the origin in the circular formation. The deadlock
resolution ends at $t=390$. The whole coordination task is completed at
$t=620$ and no collision happens.

Figure 6:  Symmetric example where no agent personality information is taken
into considerations. The dash lines retrospectively show the path of the
agents swapping their positions with the agents on their opposite direction,
e.g. agent 1 swaps its position with agent 4. The safety margin $R_{safe}$ is
set to be $1$m for all agents. The coordinate axes are in meters, and the
timestep is 0.01 seconds.

The result of the second example where the proposed asymmetric method is
deployed is shown in Fig. 7. All position swapping pairs consist of an
egoisitic agent and an altruistic agent. The deadlock situation occurs at
$t=142$. An interesting observation is that the way agents form the deadlock
resolution is quite different from the symmetric case in Fig. 6. Agents
$\\{1,3,5\\}$ form an outer triangle outside the inner triangle formed by
agents $\\{2,4,6\\}$. During the deadlock resolution period, the relatively
egoistic agents $\\{2,4,6\\}$ only travel half of the distance travelled by
agents $\\{1,3,5\\}$, who are relatively altruistic. With fewer agents
assigned to more than one rotation track, the rotation speed is greatly
improved as congestion is eliminated compared to the previous symmetric
example. The deadlock is resolved at $t=275$ and the whole circular position
swapping task is completed at $t=415$ with no collision happening among
agents.

A more general comparison of the proposed method and symmetric control is
provided in Fig. 8, where the average minimum distance among agents over time
is presented. Five sets of global personality scores are randomly initialized
in the range of 1 to 10 excluding repeated values, and tested by applying
symmetric control and the proposed method. The average minimum distance $D$ is
computed for all five sets of simulation. Overall, both symmetric and
asymmetric examples ensure collision-free interactions among agents. However,
the proposed method outperforms the symmetric example in terms of deadlock
resolution efficiency and the overall time it takes to complete the multi-
agent coordination task. The deadlock resolution efficiency is improved by
46.8% and the overall task completion time is improved by 33%.

Figure 7:  Asymmetric example where agent $\\{1,3,5\\}$ are designed to be
more altruistic with a larger $\theta$ and agent $\\{2,4,6\\}$ are designed to
be more egoistic with a smaller $\theta$. The dash lines retrospectively show
the path of the agents swapping their positions with the agents on their
opposite direction, e.g. agent 1 swaps its position with agent 4. Figure 8:
Comparison of the proposed method and symmetric control. The x-axis is
timestep and y-axis is the average minimum distance among all agents in the
multi-agent system.

###  V-D Who benefits the most?

An interesting question to ask is, who benefits the most from the multi-agent
coordination task: the egoistic agents or the altruistic agents?

We argue that both kinds of agents have their own advantages. Egoistic agents
gain a relatively flexible control space when considering avoiding potential
collision with others, and therefore their major tasks are less interrupted,
for example, there is less deviation from the nominal controller. Altruistic
agents have tighter safe control bounds in order to keep them safe, and as a
result, they can deviate from the nominal controller more frequently.

From another angle, the effect of different personalities on resulting
behavior can be viewed as a kind of balance between time and travel distance.
Egoistic agents stay close to the nominal control and the preferred path.
However, this property may sometimes lead to higher time cost. While the
altruistic agents are willing to actively accelerate deadlock resolutions, it
comes with a consequence of longer travelled distance.

##  VI CONCLUSIONS

In this paper, we present a novel social preference-aware decentralized safe
control framework. Responsibility allocation in collision avoidance is
conducted in a multi-agent system where heterogeneous agents are specified by
their individual personalities. Responsibility-associated Social Value
Orientation (R-SVO) is introduced to measure the relative difference in
personality among agents while preserving the intended social implications and
other inherited properties from traditional SVO. By leveraging the idea of
proposed Local Pairwise Responsibility Weights, a Responsibility-associated
Control Barrer Function-based safe control framework is presented along with
formal proofs on safety guarantees. The proposed method is generally
applicable to many real-time robotics applications where multi-agent
coordination is required in a safe and efficient manner. In future, we aim to
investigate how personalities could be assigned in the multi-agent system for
better task coordination.

## References

  * [1] J. Parker, E. Nunes, J. Godoy, and M. Gini, “Exploiting spatial locality and heterogeneity of agents for search and rescue teamwork,” _Journal of Field Robotics_ , vol. 33, no. 7, pp. 877–900, 2016.
  * [2] T. Rouček, M. Pecka, P. Čížek, T. Petříček, J. Bayer, V. Šalanskỳ, D. Heřt, M. Petrlík, T. Báča, V. Spurnỳ _et al._ , “Darpa subterranean challenge: Multi-robotic exploration of underground environments,” in _International Conference on Modelling and Simulation for Autonomous Systems_. Springer, 2019, pp. 274–290.
  * [3] C. Yu, X. Wang, X. Xu, M. Zhang, H. Ge, J. Ren, L. Sun, B. Chen, and G. Tan, “Distributed multiagent coordinated learning for autonomous driving in highways based on dynamic coordination graphs,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 21, no. 2, pp. 735–748, 2019.
  * [4] Y. Lyu, W. Luo, and J. M. Dolan, “Probabilistic safety-assured adaptive merging control for autonomous vehicles,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 10 764–10 770.
  * [5] J. Van den Berg, M. Lin, and D. Manocha, “Reciprocal velocity obstacles for real-time multi-agent navigation,” in _IEEE International Conference on Robotics and Automation_ , 2008, pp. 1928–1935.
  * [6] J. Alonso-Mora, A. Breitenmoser, P. Beardsley, and R. Siegwart, “Reciprocal collision avoidance for multiple car-like robots,” in _IEEE International Conference on Robotics and Automation_ , 2012.
  * [7] M. Wang, Z. Wang, J. Talbot, J. C. Gerdes, and M. Schwager, “Game theoretic planning for self-driving cars in competitive scenarios.” in _Robotics: Science and Systems_ , 2019.
  * [8] W. Luo, W. Sun, and A. Kapoor, “Multi-robot collision avoidance under uncertainty with probabilistic safety barrier certificates,” _Advances in Neural Information Processing Systems_ , vol. 33, 2020.
  * [9] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control Barrier Function based Quadratic Programs with Application to Adaptive Cruise Control,” in _53rd IEEE Conference on Decision and Control_ , 2014.
  * [10] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in _18th European Control Conference (ECC)_ , 2019.
  * [11] L. Wang, A. D. Ames, and M. Egerstedt, “Safety barrier certificates for collisions-free multirobot systems,” _IEEE Transactions on Robotics_ , vol. 33, no. 3, pp. 661–674, 2017.
  * [12] D. Bareiss and J. Van Den Berg, “Generalized reciprocal collision avoidance,” _The International Journal of Robotics Research_ , vol. 34, no. 12, pp. 1501–1514, 2015.
  * [13] A. Best, S. Narang, and D. Manocha, “Real-time reciprocal collision avoidance with elliptical agents,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2016, pp. 298–305.
  * [14] U. Borrmann, L. Wang, A. D. Ames, and M. Egerstedt, “Control barrier certificates for safe swarm behavior,” _IFAC-PapersOnLine_ , vol. 48, no. 27, pp. 68–73, 2015.
  * [15] D. Zhou, Z. Wang, S. Bandyopadhyay, and M. Schwager, “Fast, on-line collision avoidance for dynamic vehicles using buffered voronoi cells,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 2, pp. 1047–1054, 2017.
  * [16] L. Wang, A. D. Ames, and M. Egerstedt, “Safe certificate-based maneuvers for teams of quadrotors using differential flatness,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2017.
  * [17] W. B. G. Liebrand, “The effect of social motives, communication and group size on behaviour in an n-person multi-stage mixed-motive game,” _European Journal of Social Psychology_ , vol. 14, pp. 239–264, 1984.
  * [18] R. O. Murphy, K. A. Ackermann, and M. Handgraaf, “Measuring social value orientation,” _Judgment and Decision making_ , vol. 6, no. 8, pp. 771–781, 2011.
  * [19] W. Schwarting, A. Pierson, J. Alonso-Mora, S. Karaman, and D. Rus, “Social behavior for autonomous vehicles,” _Proceedings of the National Academy of Sciences_ , vol. 116, no. 50, pp. 24 972–24 978, 2019.
  * [20] A. Pierson, W. Schwarting, S. Karaman, and D. Rus, “Weighted buffered voronoi cells for distributed semi-cooperative behavior,” in _IEEE international conference on robotics and automation (ICRA)_ , 2020, pp. 5611–5617.
  * [21] B. Toghi, R. Valiente, D. Sadigh, R. Pedarsani, and Y. P. Fallah, “Social coordination and altruism in autonomous driving,” _arXiv preprint arXiv:2107.00200_ , 2021.
  * [22] L. Wang, A. Ames, and M. Egerstedt, “Safety barrier certificates for heterogeneous multi-robot systems,” in _American Control Conference (ACC)_ , 2016, pp. 5213–5218.
  * [23] G. Angeris, K. Shah, and M. Schwager, “Fast reciprocal collision avoidance under measurement uncertainty,” _International Symposium on Robotics Research (ISRR)_ , 2019.
  * [24] R. O. Murphy and K. A. Ackermann, “Social preferences, positive expectations, and trust based cooperation,” _Journal of Mathematical Psychology_ , vol. 67, pp. 45–50, 2015.
  * [25] J. Zeng, B. Zhang, and K. Sreenath, “Safety-critical model predictive control with discrete-time control barrier function,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 3882–3889.
  * [26] S. He, J. Zeng, B. Zhang, and K. Sreenath, “Rule-based safety-critical control design using control barrier functions with application to autonomous lane change,” 2021.