# Gravitational Redshift Detection from the Magnetic White Dwarf Harbored in
RX J1712.6$-$2414

Takayuki Hayashi  Center for Research and Exploration in Space Science and
Technology (CRESST II), Greenbelt, MD 20771, USA  Department of Physics,
University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD
21250, USA  NASA’s Goddard Space Flight Center, X-ray Astrophysics Division,
Greenbelt, MD 20771, USA  Hideyuki Mori  Institude of Space and Astronautical
Science (JAXA), Scientific Ballooning Research and Operation Group,
Sagamihara, Kanagawa, 252-5210, Japan  Koji Mukai  Center for Research and
Exploration in Space Science and Technology (CRESST II), Greenbelt, MD 20771,
USA  Department of Physics, University of Maryland, Baltimore County, 1000
Hilltop Circle, Baltimore, MD 21250, USA  NASA’s Goddard Space Flight Center,
X-ray Astrophysics Division, Greenbelt, MD 20771, USA  Yukikatsu Terada
Institude of Space and Astronautical Science (JAXA), Department of Space
Astronomy and Astrophysics, Sagamihara, Kanagawa, 252-5210, Japan  Graduate
School of Science and Engineering, Saitama University, 255 Shimo-Okubo,
Sakura-ku, Saitama City, Saitama 338-8570, Japan  Manabu Ishida  Institude of
Space and Astronautical Science (JAXA), Department of Space Astronomy and
Astrophysics, Sagamihara, Kanagawa, 252-5210, Japan  Department of Physics,
Tokyo Metropolitan University, 1-1 Minami-Osawa, Hachioji, Tokyo 192-0397,
Japan

###### Abstract

Gravitational redshift is a fundamental parameter that allows us to determine
the mass-to-radius ratio of compact stellar objects, such as black holes,
neutron stars, and white dwarfs (WDs). In the X-ray spectra of the close
binary system, RX J1712.6$-$2414, obtained from the Chandra High-Energy
Transmission Grating observation, we detected significant redshifts for
characteristic X-rays emitted from hydrogen-like magnesium, silicon ($\Delta
E/E_{\rm rest}\sim 7\times 10^{-4}$), and sulfur ($\Delta E/E_{\rm rest}\sim
15\times 10^{-4}$) ions, which are over the instrumental absolute energy
accuracy (${\Delta E/E_{\rm rest}\sim 3.3}\times 10^{-4}$). Considering some
possible factors, such as Doppler shifts associated with the plasma flow,
systemic velocity, and optical depth, we concluded that the major contributor
to the observed redshift is the gravitational redshift of the WD harbored in
the binary system, which is the first gravitational redshift detection from a
magnetic WD. Moreover, the gravitational redshift provides us with a new
method of the WD mass measurement by invoking the plasma-flow theory with
strong magnetic fields in close binaries. Regardless of large uncertainty, our
new method estimated the WD mass to be $M_{\rm WD}>0.9\,M_{\odot}$.

##  1 Introduction

Main-sequence (MS) stars having a mass of less than 8 $M_{\odot}$ will evolve
into a white dwarf (WD), a compact object with a radius of the order of
$10^{4}$ km. Besides, the WDs can be evolved into even by the MSs of more than
8 $M_{\odot}$ in binaries due to binary evolution. The WD is supported against
its gravity by electron degeneracy pressure. The WD in a close binary system
is fed by a mass accretion from its companion star. Hence, as the WD becomes
more massive, its radius shrinks to reinforce the degeneracy pressure.
However, the WD mass has an upper limit (Chandrasekhar mass $\sim$ $1.38$
$M_{\odot}$) where the degeneracy pressure no longer supports its gravity
(Chandrasekhar, 1931). A WD exceeding the mass limit is expected to bring
about an explosion called a type-Ia supernova or an accretion-induced collapse
into a neutron star. Therefore, the mass is a key parameter of the WDs. The
type-Ia supernovae are used to calculate the distances from our galaxy to
their host galaxies, demonstrating the accelerating expansion of the Universe
(Riess et al., 1998; Perlmutter et al., 1999).

A gravitational redshift enables us to directly measure the WD mass-radius
ratio. In weak gravitational regime, the gravitational redshift can be written
as

$v_{g}=cz=\frac{c\Delta E}{E_{\rm obs}}\simeq\frac{c\Delta E}{E_{\rm
rest}}=0.635\frac{M_{\odot}}{R_{\odot}}\,{\rm km\,s^{-1}},$ (1)

where $c$ is the speed of the light, $z$ is the redshift parameter, $E_{\rm
obs}$ and $E_{\rm rest}$ are observed and rest-frame energies, respectively,
and $\Delta E\equiv-(E_{\rm obs}-E_{\rm rest})$. Since Einstein proposed three
measurements to test the general relativity, one of which is to measure the
gravitational redshift from stars (Einstein, 1916), the gravitational redshift
has been employed to calculate the WD mass. For example, the gravitational
redshift of the first WD Sirius B is $v_{g}=cz=c\times(2.688\pm 0.026)\times
10^{-4}=80.65\pm 0.77$ km s-1 and thus the WD mass is $1.017\pm 0.025$ M⊙
(Joyce et al., 2018). This technique has been applied for the WDs in the
common proper motion binaries (Greenstein & Trimble, 1967; Koester, 1987;
Silvestri et al., 2001), in binaries in which the motion of both constituent
star were well determined (Sion et al., 1998; Long & Gilliland, 1999; Smith et
al., 2006; Steeghs et al., 2007; van Spaandonk et al., 2010; Parsons et al.,
2012, 2017; Joyce et al., 2018), or in an open cluster (Greenstein & Trimble,
1967; Wegner et al., 1989; Claver et al., 2001; Pasquini et al., 2019), where
the Doppler shift can be precisely measured to break the degeneracy with the
gravitational redshift. However, there is no detection report of the
gravitational redshift from highly magnetized ($B\gtrsim$ 0.1 MG), spun-up
($P_{\rm spin}\lesssim 10^{3}$ s) WDs, or from a WD with the X-ray.

Magnetic cataclysmic variables (mCVs) harbors a highly-magnetized WD, which is
often highly spun up via mass accretion. The cataclysmic variables (CVs) are
close binaries consisting of a WD and a companion of late-type star (Mukai,
2017). In the CVs, the accreting gas from the companion brings the angular
momentum to the WD, causing the WD spin-up and frequently forming the
accretion disk around the WD. Furthermore, as WD shrinks caused by mass gain,
its magnetic field should be strengthened (Das et al., 2013). Such a strong
magnetic field either prevents the accretion disk from reaching the WD surface
(intermediate polar: IP; Patterson 1994) or the prevents forming the accretion
disk itself (polar; Cropper 1990). The accreting gas is captured along the
magnetic field and accelerated by the WD potential to a velocity greater than
the sound velocity. Hence, a standing shock forms near the WD surface.
According to Rankine-Hugoniot relations, the post-shock gas is heated up to
$\sim 10^{8}$ K, and then the gas is highly ionized. While the post-shock gas
is descending to the WD, the plasma is cooled down by emitting X-rays. The
electrons and ions are gradually recombined so that the hydrogen- and helium-
like (H- and He-like) ions of various elements (e.g., neon (Ne), magnesium
(Mg), silicon (Si), sulfur (S), Argon (Ar), and iron (Fe)) are produced. The
recombined ions return to their ground states by emitting X-ray line
emissions, as is indicated by observed X-ray spectra. The ratio of the line
intensities is contingent upon the ion population in the plasma flow. The
features in the X-ray spectra, such as the cut-off energy of the continuum
spectra, line intensity ratio, and line energy shift caused by the Doppler
shift and/or the gravitational redshift, allow us to measure the temperature,
density, velocity, and gravitational potential in the post-shock plasma in
principle, all of which are linked to the WD potential, that is, the WD mass.

RX J1712.6$-$2414, also known as V2400 Oph, is an IP that was discovered from
the ROSAT All-Sky Survey (Buckley et al., 1995). Follow-up optical/near-
infrared observation detected a circular polarization with the spin period of
$927$ s, indicating WD’s relatively strong magnetic field of
$\bigl{(}9$–$27\bigr{)}\times 10^{6}$ G and that we always see only one of the
magnetic poles; the accretion flow onto this magnetic pole is nearly parallel
to our line of sight. This IP has been known as a diskless IP since it does
not show the $927$-s spin modulation in most cases but shows the synodic
$1003$-s modulation in the X-ray. On the other hand, the spin modulation was
observed in 2001 by XMM-Newton, and in 2005 and 2014 by Suzaku (Joshi et al.,
2019). However, the best-fit spectral model parameters of the 2001 XMM-Newton
observation are consistent with those of the 2000 XMM-Newton observation where
the spin modulation was not detected (Joshi et al., 2019), which means the
effect of whether the disk forms on the plasma flow structure is minor. We
observed RX J1712.6$-$2414 by High-Energy Transmission Grating (HETG) of the
Chandra observatory to investigate the velocity profile of the plasma in the
accretion flow. The HETG spectra potentially allow us to measure the Doppler
shift of $\sim 30$ km s-1 (Ishibashi et al., 2006).

##  2 Observation and data reduction

We carried out the X-ray observation of RX J1712.6$-$2414 with Chandra in May
2020. Table 1 shows the observation log of RX J1712.6$-$2414\. The observation
was divided into six intervals with the following observation IDs: 21274,
23038, 23039, 23244, 23267, and 23268. The total exposure time was $169.15$
ks. The HETG modules were inserted between the X-ray optics and the CCD chips.
The HETG consists of two independent gratings: the Medium Energy Grating (MEG)
covers the 0.4-5.0 keV with an energy resolution $\Delta E/E\sim$ 1/300 at 2
keV, and the High Energy Grating (HEG) covers the 0.8-10 keV with $\Delta
E/E\sim$ 200 at 6 keV. We chose the FAINT and Timed Exposure modes for the
instrumental setup. First, applying the latest calibration files of the
instruments (the 4.9.5 version), we reprocessed the observational data with
CIAO v4.12 (Fruscione et al., 2006). We did not apply any filters to the data.
After applying the barycentric correction, we then created the averaged X-ray
spectra, following the instructions of the Chandra
analysis111https://cxc.harvard.edu/ciao/.

Table 1: Chandra observations of RX J1712.6$-$2414 Observation ID | Exposure (ks) | Start date | Instrument
---|---|---|---
21274 | 24.7 | 2020-05-30 07:31:32 | HETGa
23038 | 22.9 | 2020-05-29 03:35:56 | HETGa
23039 | 37.6 | 2020-05-06 04:05:57 | HETGa
23244 | 37.6 | 2020-05-06 21:24:29 | HETGa
23267 | 24.7 | 2020-05-31 01:29:08 | HETGa
3268 | 21.7 | 2020-05-31 18:53:25 | HETGa
aHigh-Energy Transmission Grating

##  3 Analysis and result

Figure 1 shows X-ray HETG spectra around H-like Kα lines of Ne, Mg, Si, S, Ar,
and Fe. We focused on these emission lines to measure each energy shift
because they consist of only Kα1 and Kα2, whose intensity ratio is invariant
to the plasma temperature or density, and can be easily separated from lines
emitted by ions in the other states. To determine their energy centroids, we
extracted the spectra in the energy ranges of $0.98$–$1.07$, $1.44$–$1.51$,
$1.96$–$2.05$, $2.52$–$2.72$, $3.18$–$3.44$, and $6.85$–$7.10$ keV for Ne, Mg,
Si, S, Ar, and Fe, respectively. We fitted a power-law function to the
continuum and two Gaussians incorporating the redshift parameter (zgauss) to
the H-like Kα1 and Kα2 lines by using Xspec (version 12.11.0; Arnaud 1996).
The Gaussian centers were fixed at the rest-frame energies of Kα1 and Kα2
lines tabulated in Table 2 and their widths were fixed at 0. The redshift
parameter was common to the H-like Kα1 and Kα2 for each ions, and left free to
fit the model to the data. Their intensity ratio was fixed at the nominal of
$I_{{\rm K}_{\alpha 1}}/I_{{\rm K}_{\alpha 2}}=2$. The energy ranges we
selected were narrow; no photoelectric absorption was introduced.

While the energy shifts of the Kα lines from H-like Ne, Ar, and Fe ions are
marginal because of insufficient photon statistics, those from H-like Mg, Si,
and S ions are statistically constrained as shown in Table 2: $\Delta E/E_{\rm
rest}=6.9_{-0.2}^{+0.0}\times 10^{-4}$ for Mg, $7.4_{-0.7}^{+0.0}\times
10^{-4}$ for Si, and $15.4_{-4.6}^{+5.5}\times 10^{-4}$ for S, which
correspond to the line-of-sight velocities of $2.1_{-0.1}^{+0.0}\times
10^{2}$, $2.2_{-0.2}^{+0.0}\times 10^{2}$, and $4.6_{-1.4}^{+1.7}$ km
s${}^{-1}\times 10^{2}$, respectively. The errors represent the statistical
$90$% confidence level.

Figure 1:  X-ray spectra obtained from the Chandra High Energy Grating (HEG; black) and Medium Energy Grating (MEG; red). We show the spectra including emission lines from H-like ions of Ne, Mg, Si, S, Ar, and Fe. The thick solid lines are the best-fit models. Meanwhile, the thin solid lines are the components of the best-fit models, consisting of a power-law function and two Gaussians. The residuals from the best-fit model are shown at the bottom of each panel. We also show the position of the rest-frame energy of each emission line by a vertical dotted line and the energy shift value. Table 2: Summary of fitting to the H-like ion emission lines. | Ne | Mg | Si | S | Ar | Fe
---|---|---|---|---|---|---
Significance ($\sigma$) | 3.5 | 7.2 | 11.6 | 5.8 | 2.3 | 4.9
$E_{\rm rest}$ of Kα1 (keV)1 | 1.0220 | 1.4726 | 2.0061 | 2.6227 | 3.3230 | 6.9732
$E_{\rm rest}$ of Kα2 (keV)1 | 1.0215 | 1.4717 | 2.0043 | 2.6197 | 3.3182 | 6.9520
Energy centroid of Kα1 and Kα2 (keV)2 | 1.0218 | 1.4723 | 2.0055 | 2.6217 | 3.3214 | 6.9660
$z$ $\simeq\Delta E/E_{\rm rest}$ ($\times 10^{-4}$)2 | $1.0_{-\infty}^{+\infty}$ | $6.9_{-0.2}^{+0.0}$ | $7.4_{-0.7}^{+0.0}$ | $15.4_{-4.6}^{+5.5}$ | $10.5_{-12.4}^{+12.2}$ | $1.2_{-\infty}^{+\infty}$
$v\,(\times 10^{2}$ km s-1)2 | $0.3_{-\infty}^{+\infty}$ | $2.1_{-0.1}^{+0.0}$ | $2.2_{-0.2}^{+0.0}$ | $4.6_{-1.4}^{+1.7}$ | $3.1\pm 3.7$ | $0.4_{-\infty}^{+\infty}$
Energy centroid of Kα1 and Kα2 (keV)3 |  | 1.4722 | 2.0052 | 2.6212 |  |
$z$ $\simeq\Delta E/E_{\rm rest}$ ($\times 10^{-4}$)3 |  | 6.4${}_{-0.0}^{+0.4}$ | 3.4${}_{-0.4}^{+4.0}$ | 16.0${}_{-6.1}^{+0.0}$ |  |
$v\,(\times 10^{2}$ km s-1)3 |  | $1.9_{-0.0}^{+0.1}$ | $1.0_{-0.1}^{+1.1}$ | $4.8_{-1.8}^{+0.0}$ |  |
1 AtomDB: http://www.atomdb.org
2 Assuming the nominal intensity ratio $I_{{\rm K}_{\alpha 1}}/I_{{\rm
K}_{\alpha 2}}=2$
3 Assuming the optically-thick intensity ratio $I_{{\rm K}_{\alpha 1}}/I_{{\rm
K}_{\alpha 2}}=1$

The instrumental absolute energy accuracy is $\Delta E/E_{\rm rest}\simeq\pm
3.3\times 10^{-4}$ (i.e., $\pm 1\times 10^{2}$ km s-1)
222https://cxc.harvard.edu/proposer/POG/html/index.html. Even considering the
instrumental accuracy, the emission lines demonstrate the energy shift of
$\Delta E/E_{\rm rest}>3\times 10^{-4}$ and the line-of-light velocity of
$v>1\times 10^{2}$ km s-1.

##  4 Discussion

We found redshifts of $z\simeq\Delta E/E_{\rm rest}=7$ \- $15\times 10^{-4}$,
which correspond to the line-of-sight velocity of 200-450 km s-1 with the
H-like Kα lines of Mg, Si, and S from RX J1712.6$-$2414\. The detected
redshifts are statistically significant and surpass the instrumental accuracy
of $\Delta E/E_{\rm rest}$ $\simeq$ 3 $\times 10^{-4}$ (i.e., 100 km s-1 in
the line-of-light velocity).

We realized that the measured redshifts cannot be explained by the current
accretion models. The magnetically channeled accretion plasma flow is modeled
by energy and momentum equations combined with the equation of state of the
ideal gas, assuming that the plasma velocity is zero at the WD surface (Aizu,
1973) (see Appendix). The solutions determine the temperature, density, and
velocity profiles of the plasma along the flow. The maximum temperature of the
post-shock plasma in RX J1712.6$-$2414 was measured to be 23–26 keV (Yuasa et
al., 2010; Xu et al., 2016; Joshi et al., 2019). Therefore, the WD mass is
never less than 0.6 $M_{\odot}$ based on the jump conditions of the strong
shock (equation 5 and 8). Figure 2 shows the relations between the temperature
and $\Delta E/E_{\rm rest}$ (i.e., the line-of-sight velocity) along the flow
with the WD mass of 0.6 $M_{\odot}$. Note that we assume here the exact pole-
on geometry. Thus, the line-of-sight velocity equals to the actual one. A
simple analytic model with the isobaric approximation (Frank et al., 2002)
shows that the accreting plasma becomes faster with a lighter WD by comparing
at a certain temperature as is indicated by equation 14. The plasma velocity
measured with an emission line is the velocity of the local plasma whose
temperature is at the emissivity maximum of the corresponding line
(hereinafter called the line peak emissivity temperature). The local plasma
velocity can be measured with the centroid energy of the emission line which
has the emissivity peak there. The plasma velocity that the emission lines can
measure is that at the certain temperature (hereinafter called the line peak
emissivity temperature), where the corresponding species dominates the ion
population. We note that the velocity is independent of the cooling function
(Equation 14) and, therefore, whether the cyclotron cooling is significant
does not matter. In fact, the numerical calculation shows that the
temperature$-$velocity relations involving and not involving the cyclotron
cooling are almost identical for the high specific accretion rate (i.e.,
accretion rate per unit area) of $a$ = 1 g cm-2 cm-1 (Figure 2). Moreover, the
plasma is more quickly decelerated than the prediction of the isobaric model
because the pressure increases as the plasma descends. A small specific
accretion rate ($a$ = 0.01 g cm-2 cm-1 in Figure 2) enlarges the increase in
the pressure and makes the plasma velocity even slower (Hayashi & Ishida,
2014). In summary, the fastest flow is realized by the lightest WD mass ($0.6$
$M_{\odot}$ for RX J1712.6$-$2414) and a high enough specific accretion rate.
The observed line-of-sight velocities of
$\raisebox{-2.15277pt}{$\>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}\>$}1\times
10^{2}$ km s-1 are significantly faster than the theoretical fastest flow
($\simeq$ 30 km s-1 and 80 km s-1 at the line peak emissivity temperatures of
Mg and S, respectively).

Figure 2:  $\Delta E/E_{\rm rest}$ and the line-of-sight plasma velocity
measured with the emission lines of H-like Mg, Si, and S ions (squares). The
temperature of data points represents the line peak emissivity temperature of
the the corresponding ion (see the text). Each error bar shows the sum of the
statistical error for a $90$% confidence and the instrumental absolute energy
uncertainty. The lines represent theoretical temperature-velocity relations
(Hayashi & Ishida, 2014) of the plasma flow with the WD mass of $0.6$
$M_{\odot}$ assuming the exact pole-on geometry: thick solid lines are the
cases of $a$ = 1 g cm-2 cm-1 and $B$ = 0 G (black), $a$ = 1 g cm-2 cm-1 and
$B$ = 30M G (red), and $a$ = 0.01 g cm-2 cm-1 and $B$ = 0 G (blue). The black
dash-dotted line shows the isobaric accretion flow.

We investigated other possibilities that increase the $\Delta E/E_{\rm rest}$.
RX J1712.6$-$2414 is located at ($l$, $b$) = ($+359^{\circ}\\!\\!.87$,
$+8^{\circ}\\!\\!.74$) in the Galactic coordinate, and its distance and proper
motions ($\mu_{\rm ra}$, $\mu_{\rm dec}$) are 699.7${}_{-10.8}^{+9.7}$ pc and
(-1.765, 2.557) in mas yr-1 (Gaia Collaboration et al., 2022; Bailer-Jones et
al., 2021), respectively. Assuming the space motion along the galactic pole
$W=0$, the radial velocity is calculated at $-$2 km s-1. Indeed, the optical
spectra obtained from the South African Astronomical Observatory revealed that
the systemic velocity is less than $20$ km s-1(Buckley et al., 1995).
Furthermore, the line-of-sight velocity associated with the binary motion is
nullified in our phase-averaged spectra and does not affect the result.

Although significant optical depth shifts the energy centroids of the emission
lines (Del Zanna et al., 2002), it is not enough to explain the observed
redshifts. The Kα lines are ensembles of the Kα2 and Kα1 lines, whose
intensity ratio ($I_{{\rm K}_{\alpha 2}}/I_{{\rm K}_{\alpha 1}}$) affects the
energy centroid of the Kα lines. The optical depth at the Kα1 line is double
that at the Kα2 line. Therefore, the Kα1 line is more easily attenuated than
the Kα2 line; thus, the Kα line energy centroid shifts toward the red side.
This effect approximately halves $I_{{\rm K}_{\alpha 1}}/I_{{\rm K}_{\alpha
2}}$ and makes it unity at the optically thick limit (Kastner & Kastner, 1990;
Mathioudakis et al., 1999) (see Appendix). We fitted the power-law and
2-Gaussians model in the same manner in $\S$3 by assuming $I_{{\rm K}_{\alpha
1}}/I_{{\rm K}_{\alpha 2}}=1$. The computed $\Delta E/E_{\rm rest}$ and
velocity are also presented in Table 2, which still shows the velocities of
$\raisebox{-2.15277pt}{$\>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}\>$}1\times
10^{2}$ km s-1. Although all of them are consistent with the corresponding
result with $I_{{\rm K}_{\alpha 1}}/I_{{\rm K}_{\alpha 2}}=2$ within the
statistical error, the H-like Si Kα line gave us the best-fit $\Delta E/E_{\rm
rest}$, different from the corresponding by $4\times 10^{-4}$. However, this
line showed a fine structure in the fitting residual (see Figure 5), implying
that optically thick limit (i.e., $I_{{\rm K}_{\alpha 2}}/I_{{\rm K}_{\alpha
1}}=1$) is not a good assumption.

The H-like Kα lines from the pre-shock accreting gas do not contaminate those
lines from the post-shock plasma. The pre-shock accreting plasma close to the
shock is photoionized by the X-ray irradiation from the post-shock plasma and
emits H-like Kα lines (Luna et al., 2010). However, the velocity of the pre-
shock gas is 4.1$\times 10^{3}$ km s-1 and $\Delta E/E_{\rm rest}$ =
1.4$\times 10^{-2}$ with the pole-on geometry even if the WD is light as much
as possible i.e., 0.6 $M_{\odot}$. Such highly Doppler-shifted line
spectroscopically goes out of the line of the same ion from the post-shock
plasma.

Consequently, we conclude that the measured $\Delta E/E_{\rm rest}$ requires a
gravitational redshift caused by the WD. Considering the gravitational
redshift to our accretion-flow model (Hayashi & Ishida, 2014), as well as the
systemic velocity and the instrumental absolute energy accuracy, we estimated
the WD mass of RX J1712.6$-$2414 to be $>0.9$ $M_{\odot}$ (Figure 3). The WD
mass estimated by previous works (e.g.,
$0.62^{+0.06}_{-0.05}\,M_{\odot}$(Yuasa et al., 2010), $0.72\pm
0.05\,M_{\odot}$(Suleimanov et al., 2019), and
$0.67^{+0.06}_{-0.05}\,M_{\odot}$(Shaw et al., 2020)) is lighter than ours.
One plausible explanation for the discrepancy is the cyclotron cooling that
softens the X-ray spectrum and was not considered in the previous mass
estimations. The magnetic field of RX J1712.6$-$2414 is $(9-27)\times 10^{6}$
G and comparable to that of polars in which the cyclotron cooling is
significant (Wu et al., 1994). Another reasonable reason is the X-ray
reflection that is maximized at the pole-on geometry and causes complicated
systematic error (Hayashi et al., 2021). In previous studies, X-ray reflection
was not taken into account.

Lastly, we note that more precise spectral modeling reduces the contribution
of the plasma velocity to the redshift, which improves the accuracy of the
gravitational redshift estimation. We assume that the H-like Kα lines are
emitted at the corresponding line peak emissivity temperature. However, the
lines are indeed emitted from a temperature range specified for each ion.
Accreting plasma neighboring on the WD has higher density and thus emits more
intense X-rays. Meanwhile, the neighboring plasma is slower so that the
velocity average weighted over the X-ray intensity is slower than the velocity
at the line peak emissivity temperature we used in Figures 2 and 3. The more
precise spectral modeling may necessitate greater WD mass, but it is beyond
our main aim to report the gravitational redshift detection.

Figure 3:  Same as Figure 2, except for theoretical calculations involving the
binary system motion and the gravitational redshift. Thick solid lines
represent the calculations with the WD mass of $0.9\,M_{\odot}$ (red),
$1.3\,M_{\odot}$ (black), and $1.4\,M_{\odot}$ (blue). The dashes, dotted, and
thin lines show the components of the gravitational redshift, the plasma flow
velocity (Hayashi & Ishida, 2014), and the binary systemic velocity,
respectively. The pole-on geometry is assumed for the calculations.

##  5 Conclusion

We observed the diskless intermediate polar RX J1712.6-2414 with the High-
Energy Transmission Grating (HETG) of the Chandra Observatory to study the
velocity profile of the plasma in the accretion flow. We found significant
redshifts for the K$\alpha$ lines of hydrogen-like magnesium, silicon ($\Delta
E/E_{\rm rest}\sim 7\times 10^{-4}$), and sulfur ($\Delta E/E_{\rm rest}\sim
15\times 10^{-4}$) ions, which are above the instrumental absolute energy
accuracy (${\Delta E/E_{\rm rest}\sim 3.3}\times 10^{-4}$). We considered
several factors producing the redshift, such as the Doppler shift associated
with the plasma flow velocity and the systemic velocity, the optical depth,
and the gravitational redshift, and then concluded that the gravitational
redshift is the major contributor to the observed redshift. This is the first
gravitational redshift detection from a magnetic WD. The gravitational
redshift provides us with a new method of the WD mass measurement, which
estimates the WD mass to be $M_{\rm WD}>0.9\,M_{\odot}$.

### Plasma flow model

An accretion plasma flow channelled by a magnetic field is modeled as a
1-dimensional flow. Fundamental equations of a 1-dimensional flow are the mass
continuity equation:

$\displaystyle\frac{{\rm d}}{{\rm d}z}(\rho v)=0,$ (1)

the momentum equation:

$\displaystyle\rho v\frac{{\rm d}v}{{\rm d}z}+\frac{{\rm d}P}{{\rm d}z}=\rho
F,$ (2)

and the energy equation:

$\displaystyle\frac{\rm d}{{\rm d}z}\left[v\left(\frac{1}{2}\rho
v^{2}+\frac{\gamma P}{\gamma-1}\right)\right]=\rho vF-\varepsilon.$ (3)

Here, $\rho$ denotes mass density, $v$ presents flow velocity, $P$ denotes
pressure, $F$ is an external force, $\epsilon$ is a radiative cooling rate,
and $\gamma$ is an adiabatic index of $5/3$. The integral form of equation 1
is

$\displaystyle\rho v=a,$ (4)

where $a$ is called a specific accretion rate, that is, the accretion rate per
unit area. The simultaneous equations 2, 3 and 4 are resolved under initial
conditions derived from the strong-shock jump condition calculated by the
Rankine-Hugoniot relations with the free-fall velocity:

$\displaystyle v_{0}$ $\displaystyle=$ $\displaystyle 0.25\sqrt{2GM_{\rm
WD}/(R_{\rm WD}+h)},$ (5) $\displaystyle\rho_{0}$ $\displaystyle=$
$\displaystyle\frac{a}{v_{0}},$ (6) $\displaystyle P_{0}$ $\displaystyle=$
$\displaystyle 3av_{0},$ (7) $\displaystyle T_{0}$ $\displaystyle=$
$\displaystyle 3\frac{\mu m_{\rm H}}{k}v_{0}^{2},$ (8)

where $M_{\rm WD}$ is a WD mass, $R_{\rm WD}$ is a WD radius, $G$ is the
constant of gravitation, and $h$ is a shock height. The equation of state for
the ideal gas is used:

$\displaystyle P=\frac{\rho kT}{\mu m_{\rm H}}.$ (9)

A boundary condition of soft landing is also assumed at the WD surface:

$\displaystyle v_{\rm WD}=0.$ (10)

By assuming that the shock height is negligibly low compared with the WD
radius ($R_{\rm WD}\gg h$, i.e., the post-shock region is low enough), the
Bernoulli’s principle requires the total of the pressure and the ram pressure
is constant:

$\displaystyle P_{0}+\rho_{0}v_{0}^{2}=P_{\rm WD}={\rm constant},$ (11)

where $P_{\rm WD}$ is the pressure at the WD surface at which the velocity is
zero. From equations 4 and 7, we obtain

$\displaystyle P_{0}=3\rho_{0}v_{0}^{2}.$ (12)

In other words, the pressure is increased only by a factor of 4/3 through the
entire flow. Thus, an isobaric flow is a good approximation. With the constant
pressure, from equations 4 and 9, the following is derived:

$\displaystyle\frac{v}{v_{0}}=\frac{\rho_{0}}{\rho}=\frac{T}{T_{0}}.$ (13)

From equations 5, 8 and 13, the relation between the temperature and the
velocity is written as

$\displaystyle v=\frac{v_{0}}{T_{0}}T=\left(\frac{3\mu m_{\rm
H}}{4k}\sqrt{\frac{2GM_{\rm WD}}{R_{\rm WD}}}\right)^{-1}T.$ (14)

A result of this simple analytical model is shown in the dot-dash line
(labelled with “isobaric”) in Figures 2 and 3.

Assuming only the Bremsstrahlung for the radiation cooling process, equations
2, 3, 4, and 13 give us a simple thermal function of $z$:

$\displaystyle\frac{T}{T_{0}}=\left(\frac{z}{h}\right)^{2/5},$ (15)

where

$\displaystyle h=kT_{0}^{1/2}v_{0}/({\mu m_{\rm H}}\Lambda_{\rm
m,br}\rho_{0})$ (16)

and

$\displaystyle\Lambda_{\rm m,br}\sim 7\times 10^{20}\,{\rm
erg\,g^{-1}\,s^{-1}}.$ (17)

More realistic models have been numerically calculated. The cooling via line
emission is included into the cooling function according to plasma codes
(Yuasa et al., 2010; Hayashi & Ishida, 2014). Some models involve the
cyclotron cooling, which is usually important in the polars (this radiation
mainly appears in the infrared band), as well by assuming optically thick
radiation (Wada et al., 1980; Cropper et al., 1999; Belloni et al., 2021).
Note that the difference in the cooling function does not affect the relation
between the temperature and velocity, as shown in equation 14. Moreover, in a
real plasma flow with finite height, the gravitational force works. $F$ in
equations 2 and 3 is represented as follows to include this effect.

$\displaystyle F=\frac{GM_{\rm WD}}{(R_{\rm WD}+z)^{2}}.$ (18)

The calculation of the realistic model are shown in Figures 2 and 3.

### Optical depth

The observed intensity ratio of Kα1 and Kα2 is represented by

$\displaystyle\frac{I_{{\rm K}_{\alpha 1}}}{I_{{\rm K}_{\alpha
2}}}=\frac{I_{{\rm 0,K}_{\alpha 1}}}{I_{{\rm 0,K}_{\alpha 2}}}\frac{p_{{\rm
K}_{\alpha 1}}}{p_{{\rm K}_{\alpha 2}}}=2\frac{p_{{\rm K}_{\alpha 1}}}{p_{{\rm
K}_{\alpha 2}}}$ (19)

where $I_{0}$ is the nominal intensity, that is, the intensity at the moment
of emission, and $p$ is the probability that an X-ray photon escapes out from
the post-shock plasma. The optical depth ($\tau$) at the central energy of a
line is proportional to the oscillator strength ($f$) of the corresponding
transition (Jordan, 1978). Moreover, the oscillator strength of a Kα1 line is
double that of the Kα2 line as

$\displaystyle\frac{\tau_{{\rm K}_{\alpha 1}}}{\tau_{{\rm K}_{\alpha
2}}}=\frac{f_{{\rm K}_{\alpha 1}}}{f_{{\rm K}_{\alpha 2}}}=2.$ (20)

The ratio of the escape probability was calculated (Kastner & Kastner, 1990)
for a situation in which the emitters (i.e., excited H-like ions) and
absorbers (i.e., H-like ions at the ground state) are identically distributed
as in the post-shock plasma. Figure 4 depicts the escape probability ratio for
lines with an optical depth ratio of 2, i.e., Kα1 and Kα2. At the optically
thick limit, the escape probability ratio approximately converges on 0.5 and
$I_{{\rm K}_{\alpha 1}}/I_{{\rm K}_{\alpha 2}}\simeq 1$ at the optically thick
limit from equation 19.

Figure 4:  The ratio of the escape probability of K${}_{\alpha}1$ line to
K${}_{\alpha}2$ line as a function of optical depth of the K${}_{\alpha}1$
line (log$\tau_{\rm K_{\alpha 1}}$).

We fitted the power-law and 2-Gaussians model for Mg, Si, and S in the same
manner in the $\S$3 by using $I_{{\rm K}_{\alpha 1}}/I_{{\rm K}_{\alpha
2}}=1$. Figure 5 shows the best-fit models with the data and Table 2 shows the
best-fit energy shift and velocity.

Figure 5:  Same as the spectra of Mg, Si, and S in Figure 1, except for the
optically-thick limit assumption ($I_{{\rm K}_{\alpha 1}}/I_{{\rm K}_{\alpha
2}}$ = 1 for H-like ion Kα lines).

## Acknowledgement

The authors are grateful to all of the Chandra project members for developing
the instruments and their software, the spacecraft operations, and the
calibrations. We also thank Maruzen-Yushodo Co. Ltd. and Xtra Inc. for their
language editing service of our English. This research has made use of data
obtained from the Chandra Data Archive and the Chandra Source Catalog, and
software provided by the Chandra X-ray Center (CXC) in the application
packages CIAO and Sherpa. Support for this work was provided by the National
Aeronautics and Space Administration through Chandra Award Number GO9-20022A
issued by the Chandra X-ray Center, which is operated by the Smithsonian
Astrophysical Observatory for and on behalf of the National Aeronautics Space
Administration under contract NAS8-03060. This work was also supported by JSPS
Grant-in-Aid for Scientific Research(C) Grant Number JP21K03623.

## References

  * Aizu (1973) Aizu, K. 1973, Progress of Theoretical Physics, 50, 344, doi: 10.1143/PTP.50.344
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Belloni et al. (2021) Belloni, D., Rodrigues, C. V., Schreiber, M. R., et al. 2021, Astrophysical Journal Supplement Series, 256, 45, doi: 10.3847/1538-4365/ac141c
  * Buckley et al. (1995) Buckley, D. A. H., Sekiguchi, K., Motch, C., et al. 1995, Monthly Notices of the Royal Astronomical Society, 275, 1028, doi: 10.1093/mnras/275.4.1028
  * Chandrasekhar (1931) Chandrasekhar, S. 1931, Astrophysical Journal, 74, 81, doi: 10.1086/143324
  * Claver et al. (2001) Claver, C. F., Liebert, J., Bergeron, P., & Koester, D. 2001, ApJ, 563, 987, doi: 10.1086/323792
  * Cropper (1990) Cropper, M. 1990, Space Science Reviews, 54, 195, doi: 10.1007/BF00177799
  * Cropper et al. (1999) Cropper, M., Wu, K., Ramsay, G., & Kocabiyik, A. 1999, Monthly Notices of the Royal Astronomical Society, 306, 684, doi: 10.1046/j.1365-8711.1999.02570.x
  * Das et al. (2013) Das, U., Mukhopadhyay, B., & Rao, A. R. 2013, Astrophysical Journal, Letters, 767, L14, doi: 10.1088/2041-8205/767/1/L14
  * Del Zanna et al. (2002) Del Zanna, G., Landini, M., & Mason, H. E. 2002, Astronomy and Astrophysics, 385, 968, doi: 10.1051/0004-6361:20020164
  * Einstein (1916) Einstein, A. 1916, Annalen der Physik, 354, 769, doi: 10.1002/andp.19163540702
  * Frank et al. (2002) Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. D. R. Silva & R. E. Doxsey, 62701V, doi: 10.1117/12.671760
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A.G.A., Prusti, T., & et al. 2022, A&A, doi: 10.1051/0004-6361/202243940
  * Greenstein & Trimble (1967) Greenstein, J. L., & Trimble, V. L. 1967, ApJ, 149, 283, doi: 10.1086/149254
  * Hayashi & Ishida (2014) Hayashi, T., & Ishida, M. 2014, Monthly Notices of the Royal Astronomical Society, 438, 2267, doi: 10.1093/mnras/stt2342
  * Hayashi et al. (2021) Hayashi, T., Kitaguchi, T., & Ishida, M. 2021, Monthly Notices of the Royal Astronomical Society, 504, 3651, doi: 10.1093/mnras/stab809
  * Ishibashi et al. (2006) Ishibashi, K., Dewey, D., Huenemoerder, D. P., & Testa, P. 2006, Astrophysical Journal, Letters, 644, L117, doi: 10.1086/505702
  * Jordan (1978) Jordan, C. 1978, in Progress in Atomic Spectroscopy, part A, ed. W. Hanle & H. Kleinpoppen, 1453
  * Joshi et al. (2019) Joshi, A., Pandey, J. C., & Singh, H. P. 2019, The Astronomical Journal, 158, 11, doi: 10.3847/1538-3881/ab1ea6
  * Joyce et al. (2018) Joyce, S. R. G., Barstow, M. A., Holberg, J. B., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 2361, doi: 10.1093/mnras/sty2404
  * Kastner & Kastner (1990) Kastner, S. O., & Kastner, R. E. 1990, Journal of Quantitiative Spectroscopy and Radiative Transfer, 44, 275, doi: 10.1016/0022-4073(90)90033-3
  * Koester (1987) Koester, D. 1987, ApJ, 322, 852, doi: 10.1086/165779
  * Long & Gilliland (1999) Long, K. S., & Gilliland, R. L. 1999, Astrophysical Journal, 511, 916, doi: 10.1086/306721
  * Luna et al. (2010) Luna, G. J. M., Raymond, J. C., Brickhouse, N. S., et al. 2010, Astrophysical Journal, 711, 1333, doi: 10.1088/0004-637X/711/2/1333
  * Mathioudakis et al. (1999) Mathioudakis, M., McKenny, J., Keenan, F. P., Williams, D. R., & Phillips, K. J. H. 1999, Astronomy and Astrophysics, 351, L23
  * Mukai (2017) Mukai, K. 2017, Publications of the Astronomical Society of the Pacific, 129, 062001, doi: 10.1088/1538-3873/aa6736
  * Parsons et al. (2012) Parsons, S. G., Marsh, T. R., Gänsicke, B. T., et al. 2012, Monthly Notices of the Royal Astronomical Society, 420, 3281, doi: 10.1111/j.1365-2966.2011.20251.x
  * Parsons et al. (2017) Parsons, S. G., Gänsicke, B. T., Marsh, T. R., et al. 2017, Monthly Notices of the RAS, 470, 4473, doi: 10.1093/mnras/stx1522
  * Pasquini et al. (2019) Pasquini, L., Pala, A. F., Ludwig, H. G., et al. 2019, Astronomy and Astrophysics, 627, L8, doi: 10.1051/0004-6361/201935835
  * Patterson (1994) Patterson, J. 1994, Publications of the ASP, 106, 209, doi: 10.1086/133375
  * Perlmutter et al. (1999) Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, Astrophysical Journal, 517, 565, doi: 10.1086/307221
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, Astronomical Journal, 116, 1009, doi: 10.1086/300499
  * Shaw et al. (2020) Shaw, A. W., Heinke, C. O., Mukai, K., et al. 2020, Monthly Notices of the Royal Astronomical Society, 498, 3457, doi: 10.1093/mnras/staa2592
  * Silvestri et al. (2001) Silvestri, N. M., Oswalt, T. D., Wood, M. A., et al. 2001, AJ, 121, 503, doi: 10.1086/318005
  * Sion et al. (1998) Sion, E. M., Cheng, F. H., Szkody, P., et al. 1998, Astrophysical Journal, 496, 449, doi: 10.1086/305378
  * Smith et al. (2006) Smith, A. J., Haswell, C. A., & Hynes, R. I. 2006, Monthly Notices of the Royal Astronomical Society, 369, 1537, doi: 10.1111/j.1365-2966.2006.10409.x
  * Steeghs et al. (2007) Steeghs, D., Howell, S. B., Knigge, C., et al. 2007, Astrophysical Journal, 667, 442, doi: 10.1086/520702
  * Suleimanov et al. (2019) Suleimanov, V. F., Doroshenko, V., & Werner, K. 2019, Monthly Notices of the Royal Astronomical Society, 482, 3622, doi: 10.1093/mnras/sty2952
  * van Spaandonk et al. (2010) van Spaandonk, L., Steeghs, D., Marsh, T. R., & Parsons, S. G. 2010, Astrophysical Journal, Letters, 715, L109, doi: 10.1088/2041-8205/715/2/L109
  * Wada et al. (1980) Wada, T., Shimizu, A., Suzuki, M., Kato, M., & Hoshi, R. 1980, Progress of Theoretical Physics, 64, 1986, doi: 10.1143/PTP.64.1986
  * Wegner et al. (1989) Wegner, G., Reid, I. N., & McMahan, R. K. 1989, in IAU Colloq. 114: White Dwarfs, ed. G. Wegner, Vol. 328, 378, doi: 10.1007/3-540-51031-1_352
  * Wu et al. (1994) Wu, K., Chanmugam, G., & Shaviv, G. 1994, Astrophysical Journal, 426, 664, doi: 10.1086/174103
  * Xu et al. (2016) Xu, X.-j., Wang, Q. D., & Li, X.-D. 2016, Astrophysical Journal, 818, 136, doi: 10.3847/0004-637X/818/2/136
  * Yuasa et al. (2010) Yuasa, T., Nakazawa, K., Makishima, K., et al. 2010, Astronomy and Astrophysics, 520, A25, doi: 10.1051/0004-6361/201014542


# Boltzmann State-Dependent Rationality

Osher Lerner
University of California, Berkeley

###### Abstract

This paper expands on existing learned models of human behavior via a measured
step in structured irrationality. Specifically, by replacing the suboptimality
constant $\beta$ in a Boltzmann rationality model with a function over states
$\beta(s)$, we gain natural expressivity in a computationally tractable
manner. This paper discusses relevant mathematical theory, sets up several
experimental designs, presents limited preliminary results, and proposes
future investigations.

##  1 Introduction

With the eruption of deep learning we’ve created many strong ”intelligence”
engines. Yet despite impressive results in academia, there exists a great gap
in deploying systems in the world alongside humans. Currently, operational
frameworks for human-computer interaction are the missing link, and a key
component of these frameworks are models of human behavior.

In a collaboration setting, simple algorithms often overestimate human
rationality. The existence of a mathematically clear optimal sequence of
actions does not guarantee that a spatially and temporally imprecise operator
with a unique belief state will perform them. As collaborating humans, we
temper our expectations, communication, and planning according to task
difficulty, our partner’s comfort and acuity with the task, and the
uncertainty in their information and actions. In order to successfully
collaborate with us, robots must — at least implicitly — model these notions.

A standard method for modeling our suboptimality is ”Boltzmann Rationality”
[16], which models a noisy optimal trajectory based on exponentially weighting
the cost of available trajectories. However, this model still falls short in
accurately modeling human states. Research on models of systematic
suboptimality, such as the ”Boltzmann Policy Distribution” [8], if successful,
could open the flood-gates towards countless downstream competitive,
assistive, collaborative, and observational methods. Furthermore, we believe a
better understanding of human suboptimality could provide insight into how we
represent and plan tasks. In that manner, if we learn how humans are stupid,
along the way we may learn useful frameworks for how we are so smart.

##  2 Background

Under the moniker ”Industry 4.0”, the next industrial revolution is expected
to be the collaboration of humans and robots [3]. When we combine the
consistency, precision, and power of robots with the dexterity, intuition, and
general adaptive knowledge of humans, our industries are able to operate far
more effectively than using only one of the two. A hot area of research is
thus the efficient coordination of fleets of humans and robots. Much work has
been published on the management of dynamic new safety conditions [12],
fatigue [3], ergonomics, skill [10], supervision, cycle times, preference [9],
between humans, robots, and teams of humans and/or robots in assembly tasks.
However, current work is often rigidly fixed to hand-specified characteristics
and thus unable to adapt to general considerations like humans. The space of
strategies is vast, especially since collaboration can grow to great
industrious scales. [11] outlines the myriad space of considerations and
methods and concludes that in the 2020s we will soon reach levels beyond even
support, coordination, and cognition on to co-evolution. In this imagined
system, robots initiate proactive actions, agents are self organizing, and
cognition and empathy are shared amongst both humans and robots. This
structure requires accounting for numerous models of mutual prediction,
perception, uncertainty estimation, and well-being.

One key such model is human proficiency or task difficulty. This notion
appears to critically affect other aspects of collaboration tasks, such as the
importance of nonverbal communication [1] and the asymmetry of updates to
trust [13].

###  2.1 Models of Human Proficiency

In many supervised learning settings, machine learning models have reached
human ability, and now seek to go beyond. For this purpose, researchers have
modeled the suboptimality of the human labelers. In [14], they learn each
human annotator’s decision parameters inside a (simultaneously learned) latent
image representation space. This simultaneous optimization successfully
reveals human competence, difficulty of classification, and superhuman
annotations — all three of which are image-dependent.

In the context of RL, models of capability have been developed for robots,
since human operators often misunderstand joint limits of robotic arms. [7]
introduces a framework including a function $b_{h}$ mapping a state to a value
between 0 and 1 representing a human’s belief of how likely it is to be
reachable by the robot. This formalism is then flipped around to optimize
demonstrations to best update the human’s beliefs, and achieves better
performance than simply sampling from untraversed waypoints.

In [2], the authors use [5]’s annotator error maximum likelihood estimation
algorithm for imitation learning, surpassing rewards achieved by SOTA IL
algorithms in GridWorld, robotic manipulation, and chess endgame settings.
They embed demonstrators and states in a shared representation space and
calculate ”expertise level” as a similarity of the two, and theoretically show
they can recover the true optimal policy from suboptimal demonstrations.

###  2.2 Inverse RL

Given an environment with states $S$ (often partial observations), actions
$A$, and rewards $R$ for state transitions, the reinforcement learning problem
is to find a policy $\pi(a\mid s)$ optimizing expected returns of a trajectory
$\xi$ sampled from the environment.

$\displaystyle\pi^{*}=\arg\max_{\pi}\mathbb{E}_{\xi\sim\pi}[U(\xi)]$

IRL solves the inverse problem. Given an environment and rollout data of
trajectories of a certain policy, we want to calculate the cost function
$U(\xi)$ that the policy is maximizing. The space of cost functions is vast,
so this problem is under-determined given a finite number of samples. We wish
to find meaningful solutions, which depending on the context can mean
interpretable feature weights, well-conditioned for optimization, and/or
similarity to the true cost function (under some notion of similarity in cost
function space).

In practice, we observe data from policies other than the optimal policy
$\pi^{*}$. Most often for HRI tasks, the data we collect is of humans who are
optimizing an objective subject to context clues, internal beliefs, etc.
Relative to robot action spaces, human actions have imprecision in space and
time, and are roughly planned rather than exactly optimal. Algorithms solving
the IRL problem for HRI tasks must account for human suboptimality in its many
forms.

###  2.3 Boltzmann Rationality

The Boltzmann Rationality (BR) model over trajectories is formulated as the
solution to the near-optimal maximum entropy IRL problem.

That is,

$\displaystyle\max_{P}H(P)\,\,\,\,\,\,\,\,\,\,\text{ s.t.
}\,\mathbb{E}_{\xi_{D}\sim P}[U(\xi_{D})]\approx\min_{\xi}U(\xi)$

is solved by the following Boltzmann Rational distribution over trajectories:

$\displaystyle P(\xi)$ $\displaystyle=\frac{1}{Z}e^{-U(\xi)}$ (1)

where $Z$ is a constant normalization factor, and $U(\xi)$ is the cost
accumulated by the trajectory $\xi$. Typically, cost is calculated as a
function of each state, and summed over the trajectory. Although a discount
factor is sometimes used, for simplicity we will use $\gamma=1$. The
parametrization of the cost function is taken to be a vector of weights
weighing each feature of the state, computed by some map $\phi(s)$. These
features can be prescribed, learned offline, or updated online.

$\displaystyle P(\xi\mid\theta)$
$\displaystyle=\frac{1}{Z}e^{-\sum_{s\in\xi}\theta^{T}\phi(s)}$ (2)

In order to account for a range of possible ”distances” from optimality, an
”inverse temperature” parameter $\beta$ is introduced. Varying $\beta$ will
interpolate between the optimal policy $\beta=\infty$ and uniform policy
$\beta=0$ (though differently from action-space interpolation such as used in
[2]).

$\displaystyle P(\xi\mid\theta,\beta)$
$\displaystyle=\frac{1}{Z}e^{-\beta\sum_{s\in\xi}\theta^{T}\phi(s)}$ (3)

Note that any distribution can be induced from Equation 1 by construing a
particular cost function – which is precisely the problem posed by IRL – and
so one might assume any new parameters would be redundant. However, by
imposing structure on our cost and introducing new parameters, we aim to reach
a more natural parametrization. Such a description can critically ease the
optimization search and can encode interpretable information about what is
being learned. This core idea is integral to the algorithm presented in the
next section.

The Boltzmann Rationality model seems very natural, as it coincides with
generic theoretical derivations and many physical systems. In practice
however, the mismatch between this model and actual human behavior bottlenecks
human-robot interaction algorithms, particularly in the case of collaboration.
The problem to which this paper contributes is finding better models of human
behavior.

##  3 Theory

In this section, we introduce mathematical formulation of state dependent
rationality. The derivations begin at the most general form, and we state any
simplifying assumptions along the way.

To better model human behavior, we will introduce new parameters to describe
systematic suboptimality, expanding on the one scalar value $\beta$. First, we
will impose additional structure on our cost, just as we did under Boltzmann
Rationality. We will now consider several human agents, assuming they are all
optimizing the same cost function. In an experiment with a clearly
communicated objective where humans are subjected to the same task, this is a
fair assumption. Then, we will introduce the possibility that each human has
varying suboptimality over states.

###  3.1 Why States?

It is not obvious that suboptimality should be formulated as a function over
states. Often a notion of proficiency of a ”maneuver” is preferred, while in
other cases, a dynamic notion of attention, knowledge, and mood state of a
human is needed. The space of trajectories can be more expressive than states,
but computations over it are usually intractable. In different instantiations
of the reinforcement learning framework, the relevant and rich information
could be in the action space. In practice this is rarely the case, but a
meaningful space to work in is the policy space, as is used for systematic
suboptimality in BPD [8].

While these methods should all be investigated, formulating suboptimality over
states however seems initially the most promising. In many cases states
themselves are the best available representation of when the environment may
become noisy or unfamiliar. They also contain implicit information about
trajectories and maneuvers, as certain states are only reached in specific
traversals towards goals. Particularly, we can make use of the iterative
nature of our optimizers to incorporate the implicit sequential nature of
state data into our training algorithms. Notably, the alignment of state-based
formulation with our data structure, cost formulation, and iterative training
method means we can leverage many of the same conventional computational
tricks to make optimization tractable.

###  3.2 Forward Model

Let’s see how our trajectory distribution looks when we vary $\beta$ over
states $s$. We will parametrize this function just as we did the cost, as a
vector of feature weights $\theta_{\beta}$. It may be useful to set or learn a
separate featurization than the one used for rewards, but we will assume the
features are descriptive enough to compute both suboptimality and rewards.
Equation 3 becomes

$\displaystyle P(\xi\mid\theta,\beta)$
$\displaystyle=\frac{1}{Z}e^{\sum_{s\in\xi}\beta(s)\theta^{T}\phi(s)}$
$\displaystyle P(\xi\mid\theta_{R},\theta_{\beta})$
$\displaystyle=\frac{1}{Z}e^{-\sum_{s\in\xi}\theta_{\beta}^{T}\phi_{\beta}(s)\theta_{R}^{T}\phi_{R}(s)}$
$\displaystyle P(\xi\mid\theta_{R},\theta_{\beta})$
$\displaystyle=\frac{1}{Z}e^{-\sum_{s\in\xi}\theta_{\beta}^{T}\phi(s)\theta_{R}^{T}\phi(s)}$
$\displaystyle P(\xi\mid\theta_{R},\theta_{\beta})$
$\displaystyle=\frac{1}{Z}e^{-\sum_{s\in\xi}\theta_{\beta}^{T}\phi(s)\phi(s)^{T}\theta_{R}}$
$\displaystyle P(\xi\mid\theta_{R},\theta_{\beta})$
$\displaystyle=\frac{1}{Z}e^{-\theta_{\beta}^{T}\left(\sum_{s\in\xi}\phi(s)\phi(s)^{T}\right)\theta_{R}}$

where

$Z=\sum_{\bar{\xi}\in\Xi}e^{-\theta_{\beta}^{T}\left(\sum_{s\in\bar{\xi}}\phi(s)\phi(s)^{T}\right)\theta_{R}}$
(4)

To simplify notation, we define $\Phi_{\xi}=\sum_{s\in\xi}\phi(s)\phi(s)^{T}$
to be the ”feature counts” matrix of trajectory $\xi$. So our trajectory
distribution is

$P(\xi\mid\theta_{R},\theta_{\beta})=\frac{1}{Z}e^{-\theta_{\beta}^{T}\Phi_{\xi}\theta_{R}}$
(5)

###  3.3 Inverse Model

Now we invert the reinforcement learning problem, which can be computed using
Bayesian inference. Let’s consider several humans with the same reward model
but varying state-dependent proficiency. Given rollouts $\xi^{i}_{j}$ from
human $i$ and run $j$, we compute

$\displaystyle P(\theta_{R},\\{\theta^{i}_{\beta}\\}\mid\\{\Xi^{i}\\})$
$\displaystyle=\frac{P(\\{\Xi^{i}\\}\mid\theta_{R},\\{\theta^{i}_{\beta}\\})P(\theta_{R},\\{\theta^{i}_{\beta}\\})}{P(\\{\Xi^{i}\\})}$
(6)

We assume the rollouts are independent of each other up to our human
parameters. Note, this assumption could be broken if the humans gain
proficiency during data collection (between or during trials). Mathematically,
this means

$\displaystyle
P(\\{\Xi^{i}\\}\mid\theta_{R},\\{\theta^{i}_{\beta}\\})=\prod_{i}P(\Xi^{i}\mid\theta_{R},\theta^{i}_{\beta})=\prod_{i}\prod_{j}P(\xi^{i}_{j}\mid\theta_{R},\theta^{i}_{\beta})$

We also assume the human’s parameters are independent of each other, except
for reward which they share.

$\displaystyle P(\theta_{R},\\{\theta^{i}_{\beta}\\})$
$\displaystyle=P(\theta_{R})\prod_{i}P(\theta^{i}_{\beta})$

By plugging these derived equations in, we will now attempt to solve for the
parameters $\theta^{*}$ (where $\theta$ refers to $\theta_{R}$ and
$\\{\theta_{\beta}^{i}\\}$) that achieve the maximum likelihood.

$\displaystyle\theta^{*}$
$\displaystyle=\arg\max_{\theta}\,\,P(\theta\mid\\{\Xi^{i}\\})$
$\displaystyle=\arg\max_{\theta}\,\,\log P(\theta\mid\\{\Xi^{i}\\})$
$\displaystyle=\arg\max_{\theta}\,\,\log P(\\{\Xi^{i}\\}\mid\theta)+\log
P(\theta)-\log P(\\{\Xi^{i}\\})$ $\displaystyle=\arg\max_{\theta}\,\,\log
P(\\{\Xi^{i}\\}\mid\theta)+\log P(\theta)$
$\displaystyle=\arg\max_{\theta}\,\,\sum_{i}\sum_{j}\log
P(\xi^{i}_{j}\mid\theta_{R},\theta^{i}_{\beta})+\log P(\theta)$
$\displaystyle=\arg\min_{\theta}\,\,\sum_{i}\sum_{j}\left(\theta_{R}^{T}\Phi_{\xi^{i}_{j}}\theta_{\beta}^{i}+\log
Z(\theta)\right)-\log P(\theta)$ (7)

This optimization problem is theoretically analyzed in A.

##  4 Experiments

To test the applicability of our theory, there are several experiments we can
try. Unfortunately we failed to run these experiments to produce results, but
the rough experimental design for these is elaborated here.

###  4.1 Environments

We run a simple GridWorld environment to tractably test our math. For the
OverCooked setting, we use the environment provided in [6] and human data
gathered from Mechanical Turk [4] labeled with human IDs. We chunk our data by
layout and human ID, with a total of 8 different tasks, 88 humans, and 8
(sometimes less) rollouts from each agent of 397 timesteps each. Note this
data is all gathered in the collaborative setting, but we are learning only
one agent’s reward and suboptimality parameters, and absorbing the other into
the environment.

###  4.2 Parameter recovery

First, we want to ensure our math is internally consistent, and that our
parameter spaces are meaningful. For this purpose, we create agents operating
directly using our model, and see if we can recover their internal parameters
just from rollout observation. In the GridWorld setting, this seems works
pretty well (see Figure 1), though a compilation of rigorous results has not
been collected.

Figure 1: Posterior Belief of $\theta$ parameters from trajectories generated
by $\theta_{R}=[0,1]$ and $\theta_{\beta}=[1,100]$. For this arbitrary
discrete space, we are able to accurately recover our parameters.

###  4.3 Shared Goals

We can construct another experiment to test our assumption of shared goals. We
assumed that by communicating to each play-tester the same objective, their
actions would align with models that share state feature reward weights. For
any task in which we ground our training from trials on the same objective, we
can validate this assumption experimentally. Using live subjects, we can ask
the humans to label sampled trajectories from different human runs and ask
them to choose the trajectory that best achieves the goal of the task, and see
if there exists a statistical difference between the human labelers.

With offline human data (such as our MTurk dataset), we can at least
qualitatively (if not statistically) compare the learned reward parameters of
trajectories by different humans learned with BR to see how closely they align
compared to those learned from other tasks and random weights.

###  4.4 Generalization of Learned Rewards

We wish to test how well our learned $\theta_{R}$ generalizes, to show it
meaningfully corresponds to the reward model and is independent from human
suboptimality.

We set up an experiment on each of our representative tasks (from Overcooked
and GridWorld). We fit a BR rational model to the data in aggregate and a BSDR
model with each human identity labeled. Then, we optimize the learned
$\theta_{R}$ to find $\pi^{*}_{\theta_{R}}$. By measure the performance of
these learned policies with the original task rewards (known to the
experimenters, but not our algorithms), we can see which achieves better
performance. The cost function that can be opimized to achieve a higher true
reward should be more closely aligned to the true reward.

We can also compare the rewards predicted by our learned parameters evaluated
over some sample human trajectories to the true values.

Note that it seems BSDR has a competitive advantage in this experiment since
it gets an extra dimension of data: human ID’s. Without imposing the structure
of different agents with shared cost, our BSDR parameters are redundant and
have no preference for one corresponding to reward and the other to
suboptimality. Thus BSDR cannot generalize from training on only 1 human. To
equalize the playing field, we can train BR with different $\beta_{i}$ for
each human, and evaluate their performance individually.

####  4.4.1 Action Prediction

In this experiment, we directly compare the action distributions predicted by
our human models. Under our list of tasks, we again train BR, BPD, and BSDR.
We then compare their cross entropy prediction performance of human data to
each other. We can also plot results for a self-play policy, and a random
policy.

###  4.5 Generalization of Learned Suboptimality for Goal Inference

After learning from demonstrator’s performing the same tasks, we should be
able to use our human suboptimality models to better infer new goals and
predict their actions.

###  4.6 Goal Inference

In this experiment, we test the ability to predict a goal from a partial
trajectory, which should be easier if we have a pre-existing accurate model of
their suboptimality.

We use our calculated values of $\beta^{i}$ and $\theta_{\beta}^{i}$ from the
last experiment. We restrict the reward parameters to be one of a small set.
In this case, we will use GridWorld with certain coordinates as goals. Then,
we consider unseen trajectories by the same human used for training the
$\beta$ parameters. Using only the initial portion of the trajectory, we
numerically compute with Bayesian inference the likelihood of each possible
goal. We plot the likelihood of the true goal computed from BR and BSDR
averaged over 8 different trajectories from 20 different humans. We repeat the
trial for the first 25%, 50%, 75%, and 100% of the trajectory. And we again
repeat for different GridWorld environments.

Unfortunately, we did not yet collect human data on GridWorld tasks.

##  5 Future Steps

  * 1.

State-dependent suboptimality may be useful in certain in environments with
relevant information encoded in the state, while irrelevant in others.
Learning a latent representation of states, actions, or maneuvers may be the
best bet. [2] cited meaningful state features as challenging to learn for
unexplored environments, but [15] successfully learned abstract state-
dependent action representations to surpass demonstrator performance. Once we
have a meaningful rich space we can learn suboptimality as a function of that
representation.

  * 2.

One could also investigate the use of more complex suboptimality models such
as different ways of interpolating between optimality and irrationality.

  * 3.

To best gauge human models, we need to work with real human data. Diverse
tasks, environments, and agents are important for experiments aiming to
understand human models.

  * 4.

The question is left of how to take advantage of these human models for
collaboration. Knowledge of $\beta(s)$ could be used to estimate conditional
uncertainty to assess risk in plans. For assistance tasks where the human’s
objective is unknown, an assistance task could be to steer the state towards
the human’s expertise. In a game context, given the human model, one could
explicitly compute the best response strategy to it.

  * 5.

We can experiment using different features for our reward and suboptimality
models.

## Acknowledgements

Anca Dragan and Cassidy Lailaw, the instructors of CS 287H at Berkeley, both
authored background papers, brainstormed on this idea, gave feedback, and
supported me in writing this paper. For them I am very grateful.

##  Appendix A Further Analysis of the Max Likelihood Problem

The following is an incomplete attempt at analytically analyzing the
optimization problem from the MLE in Equation 7. It can be safely ignored.

Note that the $Z(\theta)$ term is omitted from the following computations,
which is an error.

$\displaystyle\theta^{*}$
$\displaystyle=\arg\max_{\theta}\sum_{i}\sum_{j}\left(\theta_{\beta}^{T}\left(\sum_{s\in\xi^{i}_{j}}\phi(s)\phi(s)^{T}\right)\theta_{R}\right)+\log
P(\theta_{R})+\sum_{i}\log P(\theta^{i}_{\beta})$

This quantity can be computed efficiently by flattening $\Xi^{i}$ into a
vector of states visited, and performing tensor multiplications with the
$\theta^{i}$ vector. An einsum summation could work as well. Notice that the
terms involving $\theta^{i}$ can be separated out. So

$\displaystyle\theta^{*}$
$\displaystyle=\arg\min_{\theta^{i}_{\beta}}\left(\sum_{s\in\Xi^{i}}\theta_{R}^{T}\phi(s)\phi(s)^{T}\right)\theta_{\beta}^{i}+\log
P(\theta^{i}_{\beta})$

We can summarize our data with the frequency rates at which human $i$ visits
state $s$: $\rho^{i}(s)$. Then Given a discretized search space of $\theta$
values we can simply search over them all.

To find this numerically for large spaces, let’s compute the gradient.

Let’s define the matrix
$\Phi^{i}=\sum_{j}\sum_{s\in\xi^{i}_{j}}\phi(s)\phi(s)^{T}$.

Let’s assume a prior that is uniform over $\theta$ with unit length, and 0
otherwise. Then our maximum likelihood problem becomes

$\displaystyle\theta^{*}$
$\displaystyle=\arg\min_{\|\theta\|=1}\theta_{R}^{T}\sum_{i}\Phi^{i}\theta_{\beta}^{i}$

Then we see that $\theta_{R}$ that minimizes the quantity is just the unit
vector in the opposite direction from $\sum_{i}\Phi^{i}\theta_{\beta}^{i}$.

$\theta_{R}^{*}=-\frac{\sum_{i}\Phi^{i}\theta_{\beta}^{i}}{\|\sum_{i}\Phi^{i}\theta_{\beta}^{i}\|}$

Then with that fixed, we can go into the minimization over
$\\{\theta_{\beta}^{i}\\}$. We are left with

$\displaystyle\theta^{*}$
$\displaystyle=\arg\max_{\theta}P(\theta_{R},\\{\theta^{i}_{\beta}\\}\mid\\{\Xi^{i}\\})$
$\displaystyle=\arg\min_{\|\theta\|=1}-\frac{1}{\|\sum_{i}\Phi^{i}\theta_{\beta}^{i}\|}(\sum_{i}\Phi^{i}\theta_{\beta}^{i})^{T}(\sum_{i}\Phi^{i}\theta_{\beta}^{i})$
$\displaystyle=\arg\min_{\|\theta\|=1}-\frac{1}{\|\sum_{i}\Phi^{i}\theta_{\beta}^{i}\|}\|\sum_{i}\Phi^{i}\theta_{\beta}^{i}\|^{2}$
$\displaystyle=\arg\max_{\|\theta\|=1}\|\sum_{i}\Phi^{i}\theta_{\beta}^{i}\|$
$\displaystyle=\arg\max_{\|\theta\|=1}\|\sum_{i}\Phi^{i}\theta_{\beta}^{i}\|^{2}$

This is not an obvious expression to maximize, but we can get some insight to
the theoretical analysis. Let’s calculate the derivative of the unconstrained
objective

$\displaystyle\frac{\partial}{\partial\theta_{\beta}^{i}}\|\sum_{i}\Phi^{i}\theta_{\beta}^{i}\|^{2}$
$\displaystyle=\frac{\partial}{\partial\theta_{\beta}^{i}}(\sum_{i}\Phi^{i}\theta_{\beta}^{i})^{T}(\sum_{i^{\prime}}\Phi^{i^{\prime}}\theta_{\beta}^{i^{\prime}})$
$\displaystyle=2(\frac{\partial}{\partial\theta_{\beta}^{i}}\sum_{i}\Phi^{i}\theta_{\beta}^{i})^{T}(\sum_{i^{\prime}}\Phi^{i^{\prime}}\theta_{\beta}^{i^{\prime}})$
$\displaystyle=2{\Phi^{i}}^{T}(\sum_{i^{\prime}}\Phi^{i^{\prime}}\theta_{\beta}^{i^{\prime}})$

The critical points appear when the derivatives are all 0. This occurs at a
minimum when $\sum_{i}\Phi^{i}\theta_{\beta}^{i}=0$. Other critical points
appear whenever $\sum_{i}\Phi^{i}\theta_{\beta}^{i}$ is in the null space of
every ${\Phi^{i}}^{T}$. Given our weight constraint, we can use Lagrange
multipliers to solve the constrained optimization problem. Let
$f(\theta)=\|\sum_{i}\Phi^{i}\theta_{\beta}^{i}\|^{2}$ and
$g(\theta)=\|\theta_{\beta}^{i}\|-1$. Then our likelihood is maximized when
$L(\theta,\lambda)=f(\theta)+\lambda g(\theta)$ is at a stationary point.

$\displaystyle\forall_{i}\,\,\,\,\,\frac{\partial}{\partial\theta_{\beta}^{i}}L(\theta,\lambda)=0\,\,\,$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\frac{\partial}{\partial\lambda}L(\theta,\lambda)=0$
$\displaystyle\forall_{i}\,\,\,\,\,0$
$\displaystyle=2{\Phi^{i}}^{T}(\sum_{i^{\prime}}\Phi^{i^{\prime}}\theta_{\beta}^{i^{\prime}})+2\lambda\theta_{\beta}^{i}\
$ $\displaystyle 1$ $\displaystyle=\|\theta_{\beta}^{i}\|$

These computations can in practice be rather large.

  * 1.

With a large state dimension $D$, $\Phi^{i}$ will have $D^{2}$ entries. It may
be easier to keep the expression decomposed into
$\sum_{s\in\Xi^{i}}\theta_{R}^{T}\phi(s)\phi(s)^{T}$.

  * 2.

For limited data per human, it may be more efficient to use the sum over
states.

  * 3.

In the case of large data per human, and a small state space (or in rare
cases, a small number of states visited per human), it can be efficient to
transform the sum over states into a product of a state feature tensor by
$\rho^{i}(s)$, a vector representing the frequency at which a human visits
state $s$. We can construct an operational quantity like that used in
MaxEntIRL: feature counts $\tilde{\phi}^{i}$.

## References

  * AWHS [16] Henny Admoni, Thomas Weng, Bradley Hayes, and Brian Scassellati.  Robot nonverbal behavior improves task performance in difficult collaborations.  In 2016 11th ACM/IEEE International Conference on Human-Robot Interaction (HRI), pages 51–58, Christchurch, New Zealand, March 2016. IEEE.
  * BSE+ [22] Mark Beliaev, Andy Shih, Stefano Ermon, Dorsa Sadigh, and Ramtin Pedarsani.  Imitation Learning by Estimating Expertise of Demonstrators, June 2022.  arXiv:2202.01288 [cs].
  * CPA [21] Alejandro Chacón, Pere Ponsa, and Cecilio Angulo.  Cognitive Interaction Analysis in Human–Robot Collaboration Using an Assembly Task.  Electronics, 10(11):1317, May 2021.
  * CSH+ [20] Micah Carroll, Rohin Shah, Mark K. Ho, Thomas L. Griffiths, Sanjit A. Seshia, Pieter Abbeel, and Anca Dragan.  On the Utility of Learning about Humans for Human-AI Coordination, January 2020.  arXiv:1910.05789 [cs, stat].
  * DS [79] A. P. Dawid and A. M. Skene.  Maximum Likelihood Estimation of Observer Error-Rates Using the EM Algorithm.  Applied Statistics, 28(1):20, 1979.
  * FHZ+ [21] Matthew C. Fontaine, Ya-Chuan Hsu, Yulun Zhang, Bryon Tjanaka, and Stefanos Nikolaidis.  On the Importance of Environments in Human-Robot Coordination, June 2021.  arXiv:2106.10853 [cs].
  * GYS+ [22] Xiaofeng Gao, Luyao Yuan, Tianmin Shu, Hongjing Lu, and Song-Chun Zhu.  Show Me What You Can Do: Capability Calibration on Reachable Workspace for Human-Robot Collaboration.  IEEE Robotics and Automation Letters, 7(2):2644–2651, April 2022\.  arXiv:2103.04077 [cs].
  * LD [22] Cassidy Laidlaw and Anca Dragan.  The Boltzmann Policy Distribution: Accounting for Systematic Suboptimality in Human Models, April 2022.  arXiv:2204.10759 [cs].
  * LR [20] Yee Yeng Liau and Kwangyeol Ryu.  Task Allocation in Human-Robot Collaboration (HRC) Based on Task Characteristics and Agent Capability for Mold Assembly.  Procedia Manufacturing, 51:179–186, 2020.
  * LR [22] Yee Yeng Liau and Kwangyeol Ryu.  Genetic algorithm-based task allocation in multiple modes of human–robot collaboration systems with two cobots.  The International Journal of Advanced Manufacturing Technology, 119(11-12):7291–7309, April 2022.
  * LZL+ [23] Shufei Li, Pai Zheng, Sichao Liu, Zuoxu Wang, Xi Vincent Wang, Lianyu Zheng, and Lihui Wang.  Proactive human–robot collaboration: Mutual-cognitive, predictable, and self-organising perspectives.  Robotics and Computer-Integrated Manufacturing, 81:102510, June 2023\.
  * RGBL+ [17] S. Robla-Gomez, Victor M. Becerra, J. R. Llata, E. Gonzalez-Sarabia, C. Torre-Ferrero, and J. Perez-Oria.  Working Together: A Review on Safe Human-Robot Collaboration in Industrial Environments.  IEEE Access, 5:26754–26773, 2017.
  * SXCH [20] Harold Soh, Yaqi Xie, Min Chen, and David Hsu.  Multi-Task Trust Transfer for Human-Robot Interaction.  The International Journal of Robotics Research, 39(2-3):233–249, March 2020.  arXiv:1807.01866 [cs].
  * WBBP [10] Peter Welinder, Steve Branson, Serge Belongie, and Pietro Perona.  The multidimensional wisdom of crowds.  pages 2424–2432, 01 2010.
  * YLN [21] Mengjiao Yang, Sergey Levine, and Ofir Nachum.  TRAIL: Near-Optimal Imitation Learning with Suboptimal Data, October 2021.  arXiv:2110.14770 [cs].
  * ZMBD [08] Brian D. Ziebart, Andrew Maas, J. Andrew Bagnell, and Anind K. Dey.  Maximum entropy inverse reinforcement learning.  In Proceedings of the 23rd National Conference on Artificial Intelligence - Volume 3, AAAI’08, page 1433–1438. AAAI Press, 2008.