# Identifying Occluded Agents in Dynamic Games
with Noise-Corrupted Observations

Tianyu Qiu1 and David Fridovich-Keil2 1Tianyu Qiu is with the Department of
Automation in the School of Electronic Information and Electrical Engineering,
Shanghai Jiao Tong University, Shanghai 200240, China.
<EMAIL_ADDRESS>Fridovich-Keil is with the Aerospace Engineering
department, UT Austin. This research is supported by the National Science
Foundation, under grant 2211548<EMAIL_ADDRESS>
###### Abstract

To provide safe and efficient services, robots must rely on observations from
sensors (lidar, camera, etc.) to have a clear knowledge of the environment. In
multi-agent scenarios, robots must further reason about the intrinsic
motivation underlying the behavior of other agents in order to make inferences
about their future behavior. Occlusions, which often occur in robot operating
scenarios, make the decision-making of robots even more challenging. In
scenarios without occlusions, dynamic game theory provides a solid theoretical
framework for predicting the behavior of agents with different objectives
interacting with each other over time. Prior work proposed an inverse dynamic
game method to recover the game model that best explains observed behavior.
However, an apparent shortcoming is that it does not account for agents that
may be occluded. Neglecting these agents may result in risky navigation
decisions. To address this problem, we propose a novel inverse dynamic game
technique to infer the behavior of occluded, unobserved agents that best
explains the observation of visible agents’ behavior, and simultaneously to
predict the agents’ future behavior based on the recovered game model. We
demonstrate our method in several simulated scenarios. Results reveal that our
method robustly estimates agents’ objectives and predicts trajectories for
both visible and occluded agents from a short sequence of noise corrupted
trajectory observation of only the visible agents.

##  I Introduction

Robots depend on sensor observations to be alert to the presence of static and
dynamic obstacles. Yet, sensors are fundamentally limited due to occlusion or
sensing range. In practice, humans can use their prior experience to make
inferences about occluded agents and avoid potential risks. For example, in
Figure 1, the blue pedestrian is running across the road while the green and
red vehicle are driving along the road in the same direction. However, the
blue pedestrian is occluded by the red vehicle, and thus is invisible to the
green vehicle. The green vehicle will likely collide with the blue pedestrian
if the driver maintains its current speed. However, he notices that the red
vehicle beside him brakes. Thus, he infers that someone is moving in the
occluded area, forcing the red vehicle’s deceleration. Consequently, he can
actively brake to avoid the potential collision before the occluded agent
finally comes into view. This example reveals the fact that the behaviors of
occluded agents can be inferred from those of visible agents since both are
affected by their interaction. The green vehicle is then capable of acting
based on the knowledge of the invisible agents.

Figure 1: Traffic scenario, in which the green vehicle’s view is occluded by
the red vehicle so it is not able to observe the occluded blue pedestrian. Our
method utilizes only observations of visible agents’ behavior, solves an
inverse dynamic game to recover a game model that best explains the observed
behavior, and indirectly identifies the behavior of occluded agents. In this
traffic scenario, the driver in the green vehicle notices the deceleration of
the red vehicle next to him. Our method infers that the blue pedestrian is
running across the road within the occluded area. As a consequence, the driver
in the green vehicle can actively brake to avoid a potential collision.

Prior works [1, 2] have investigated multi-agent interaction problems with
dynamic game techniques. Fridovich-Keil et al.[2] proposed an iterative linear
quadratic game algorithm to compute optimal trajectories in multi-agent
noncooperative scenarios. Based on this work, Peters et al.[1] solved the
inverse problem to learn the game model from noisy observations, estimated
agents’ trajectories from the recovered model, and achieved high-quality
estimation performance in occlusion-free scenarios.

Inspired by these works [1, 2], we propose a novel technique that identifies
the unknown parameters of each agent’s cost function and estimates the
trajectories of occluded agents in a Nash game based on realistic sensor
measurements of only visible agents. We evaluate our method in various
scenarios at different levels of observation noise. Results reveal that our
method is noise robust and provides accurate trajectory prediction for both
visible and invisible agents.

##  II Related Work

###  II-A Social Occlusion Inference

Many works have explored methods to make inferences for occluded areas based
on observation of visible agents’ behavior in social settings. These works
utilize the fundamental fact that one’s behavior is affected by the
surrounding environment to a great extent. Representative works often apply
computer vision techniques [3, 4] and occupancy grid map generating techniques
[5, 6] to make inferences about the occluded area. Hara et al. [3] proposed a
CNN that identifies the existence of a person in blind spots with artificial
occlusions in a volleyball game dataset. Subsequent work [4] further extended
to traffic scenarios and proposed a spatio-temporal D-CNN to predict whether a
vehicle would come into view soon based on the tracking of visible pedestrians
from first-person view video input. Afolabi et al.[5] proposed a mapping
framework that incorporates people as sensors to generate an occupancy grid
map for the occluded area. Itkina et al.[6] further trained a conditional
variational autoencoder to generate an occupancy grid map from the observation
of driver trajectories.

However, despite these encouraging achievements in social occlusion inference,
these works only solve the problem of whether the occluded area is occupied,
or whether an agent is coming out from the occluded area. They cannot infer
the actual behavior or intent of agents in occluded areas, which prevents
robots in the scene from making optimal decisions. Other techniques must be
utilized to model the behavior of agents in occluded areas more precisely from
the observation of only visible agents’ social behavior.

###  II-B Planning-based Agent Behavior Modeling

Before we move towards the ultimate goal of precise inference about the
occluded agents’ behavior, we first introduce several planning-based methods
of agent behavior modeling which are closely related to our work. Planning-
based methods make the assumption that humans rationally make decisions about
their trajectories while walking or driving; that is, their trajectories must
explicitly or implicitly follow certain rules and these rules can be applied
to predict their future behaviors. This assumption casts the behavior
inference problem as one of (inverse) optimal control, and the goal becomes to
reconstruct such rules from the observation of agents’ behavior.

####  II-B1 Single agent behavior modeling

With the fundamental investigation of inverse optimal control (IOC) in [7],
researchers extensively applied IOC techniques in modelling single agent
behaviors[8, 9, 10], where an explicit optimal control model was proposed and
parameters were identified to best fit the observations of human behaviors.
The advance of inverse reinforcement learning (IRL) introduced in [11] also
gave rise to applications[12, 13, 14, 15] towards single agent behavior
modeling. Representative work [12] described human behavior based on a Markov
decision process with the principle of maximum entropy. Results were utilized
to improve a robot’s ability to navigate in the presence of pedestrians.

####  II-B2 Multi-agent behavior modeling

Aside from the interaction between agents and the environment, each agent’s
behavior affects that of other agents. A typical framework is multi-agent IRL
generalized from IRL for single agent, under which human social behaviors are
investigated[16, 17]. Subsequent works further make applications in robot
social navigation tasks[18, 19] to improve a robot’s ability to navigate in
the presence of pedestrians, and autonomous driving scenarios[20, 21] to
guarantee safe driving performance.

###  II-C Dynamic Games and Inverse Dynamic Games

Our work builds upon dynamic game theory, which generalizes single agent
optimal control cases to multi-agent scenarios and provides a solid framework
to model multi-agent behaviors. Fridovich-Keil et al.[2] proposed an forward
iterative linear quadratic game to depict interaction between agents in
traffic scenarios and computed feedback Nash strategies. To model human
behavior more accurately, more recent works [1], [22, 23, 24] optimized the
model by solving inverse dynamic games. Le Cleac’h et al.[22] iteratively
solved the inverse dynamic game problem to update a Bayesian estimate of
agents’ cost function parameters by recasting it in a recursive parameter-
estimation framework. Peters et al.[1, 23] jointly optimized player objectives
and trajectory estimates by coupling them through Nash equilibrium constraints
based on noisy, partial state observations.

We emphasize that these works share the same limitation, which is that
observations of all agents are required. To deal with scenarios where agents
are occluded, inspired by [1], we propose a novel inverse dynamic game
technique to identify the unknown weighting parameters in each agent’s cost
function, and simultaneously estimate both visible and invisible agents’
trajectories that best explain the observations of visible agents’
trajectories.

##  III Discrete Time Open-loop Nash Games

TABLE I: Main Symbols and Notations ${J^{i}}$ | Cost function (Objective) of the ${i^{\text{th}}}$ agent
---|---
${M\in\mathbb{N}^{*}}$ | Number of agents in the game
${\theta=(\theta^{1},\cdots,\theta^{M})}$ | Tuple of weighting parameters of all agents
${T\in\mathbb{N}^{*}}$ | Time horizon of the game
${x_{t}^{i}\in\mathbb{R}^{n}}$ | State of the ${i^{\text{th}}}$ agent at time ${t}$
${\textbf{x}^{i}:=(x_{1}^{i},\cdots,x_{T}^{i})}$ | Trajectory of the ${i^{\text{th}}}$ agent
${\textbf{x}_{t}:=(x_{t}^{1},\cdots,x_{t}^{M})}$ | State of all agents at time ${t}$
${\textbf{x}:=(\textbf{x}^{1},\cdots,\textbf{x}^{M})}$ | Tuple of trajectories of all agents
${u_{t}^{i}\in\mathbb{R}^{m}}$ | Control input of the ${i^{\text{th}}}$ agent at time ${t}$
${\textbf{u}^{i}:=(u_{1}^{i},\cdots,u_{T}^{i})}$ | Control sequence of the ${i^{\text{th}}}$ agent
${\textbf{u}_{t}:=(u_{t}^{1},\cdots,u_{t}^{M})}$ | Control input of all agents at time ${t}$
${\textbf{u}:=(\textbf{u}^{1},\cdots,\textbf{u}^{M})}$ | Tuple of control sequences of all agents
${\mathcal{L}^{i}}$ | Lagrangian of the ${i^{\text{th}}}$ agent
${\lambda_{t}^{i}\in\mathbb{R}^{n}}$ | Lagrange multiplier of the ${i^{\text{th}}}$ agent at time ${t}$
${\boldsymbol{\lambda}^{i}:=(\lambda_{1}^{i},\cdots,\lambda_{T}^{i})}$ | Sequence of ${\lambda_{t}^{i}}$ of the ${i^{\text{th}}}$ agent
${\boldsymbol{\lambda}_{t}:=(\lambda_{t}^{1},\cdots,\lambda_{t}^{M})}$ | Tuple of ${\lambda_{t}^{i}}$ of all agents at time ${t}$
${\boldsymbol{\lambda}:=(\boldsymbol{\lambda}^{1},\cdots,\boldsymbol{\lambda}^{M})}$ | Tuple of sequences of ${\lambda_{t}^{i}}$ of all agents
${\mathcal{V}}$ | Set of visible agents in the game
${\mathcal{O}}$ | Set of occluded agents in the game

A discrete time open-loop Nash game with $M$ agents is characterized by the
state ${x_{t}^{i}\in\mathbb{R}^{n}}$ and control inputs
${u_{t}^{i}\in\mathbb{R}^{m}}$, ${\forall i\in[M]}$, at time step $t$. The
dynamics ${\textbf{x}_{t+1}=f(\textbf{x}_{t},\textbf{u}_{t})}$ describe how
the game evolves with each agent’s control input, at each time step
${t\in[T]}$. ${J^{i}:=\sum_{t=1}^{T}g_{t}^{i}(\textbf{x}_{t},\textbf{u}_{t})}$
for each agent evaluates their cost in the game. The game is thus fully
characterized by the tuple of all agents’ cost functions, the initial
condition ${\textbf{x}_{1}}$, and the dynamics, which is denoted by
${\Gamma:=(\\{J^{i}\\}_{i=1}^{M},\textbf{x}_{1},f)}$. In a Nash game, each
agent aims to minimize his cost function, subject to dynamic feasibility
constraints, i.e.

$\displaystyle\min_{\textbf{x},\textbf{u}^{i}}\quad$ $\displaystyle
J^{i}(\textbf{u};\textbf{x}_{1}),$ $\displaystyle\forall i\in[M],$ (1a) s.t.
$\displaystyle\textbf{x}_{t+1}=f(\textbf{x}_{t},\textbf{u}_{t}),$
$\displaystyle\forall t\in[T-1].$ (1b)

Given that each agent decides his control input rationally at equilibrium, the
following inequalities are satisfied:

$J^{i}(\textbf{u}^{*};\textbf{x}_{1})\leq
J^{i}(\textbf{u}^{i},\textbf{u}^{-i*};\textbf{x}_{1}),\forall i\in[M],$ (2)

and ${\textbf{u}^{*}:=(\textbf{u}^{1*},\cdots,\textbf{u}^{M*})}$ is called a
Nash strategy. Equation (2) reveals the fact that no agents will decrease cost
by unilaterally deviating from Nash strategy ${\textbf{u}^{i*}}$[25].

To make these concepts concrete, we introduce the following running example.
Consider $M=2$ pedestrians walking toward their goal, avoiding each other.
${\textbf{x}_{t}:=(x_{t}^{1},x_{t}^{2})}$ denotes the position of both
pedestrians and they follow single-integrator dynamics at time discretization
$\Delta t$:

$\displaystyle x_{t+1}^{i}$ $\displaystyle=\left\\{\begin{aligned}
&p_{x,t+1}^{i}=p_{x,t}^{i}+v_{x,t}^{i}\Delta t,\\\
&p_{y,t+1}^{i}=p_{y,t}^{i}+v_{y,t}^{i}\Delta t,\end{aligned}\right.$ (3)
$\displaystyle=x_{t}^{i}+u_{t}^{i}\Delta t,\quad t\in[T-1],i\in\\{1,2\\},$
$\displaystyle\textbf{x}_{1}$ is known.

${u_{t}^{i}=[v_{x,t}^{i},v_{y,t}^{i}]^{\top}}$ denotes the velocity of each
pedestrian and ${\textbf{x}_{1}}$ is the initial state of two pedestrians.
Each pedestrian’s objective is the sum of the running cost ${g_{t}^{i}}$ over
time, where ${g_{t}^{i}}$ is the combination of different features weighted by
non-negative weighting parameters ${\theta^{i}}$:

$g_{t}^{i}=\sum_{j=1}^{n}\theta_{j}^{i}g_{j,t}^{i}\left\\{\begin{aligned}
&g_{1,t}^{i}=\|x_{t}^{i}-x_{d}^{i}\|_{2}^{2}\\\
&g_{2,t}^{i}=-\log(\|x_{t}^{i}-x_{t}^{-i}\|_{2}^{2})\\\
&g_{3,t}^{i}=\|u_{t}^{i}\|_{2}^{2}\end{aligned}\right.,$ (4)

where ${\|\cdot\|_{2}^{2}}$ denotes the squared Euclidean norm. The pedestrian
aims to get closer to his destination (${g_{1,t}^{i}}$), and to keep far away
from the other agents (${g_{2,t}^{i}}$) without great energy consumption
(${g_{3,t}^{i}}$). In practice, ${g_{t}^{i}}$ can be readily modified to
accommodate different scenarios.

##  IV Our Approach

In this work, we introduce two roles in the dynamic game: the participants and
the observer. The participants compete with each other in the game, while the
observer is outside the game, observing the interaction between the agents.
All participants are observable to each other. However, only some of the
participants are visible to the observer, while others are occluded, and hence
invisible. No matter whether the agents are visible or not, their interaction
remains the same.

In practice, the observer could be a robot wishing to navigating in human-rich
environments. To do so, he must estimate both visible and invisible agents’
objectives and trajectories from the observation of only visible agents, and
thereby improve his navigation performance.

For clarity, we seek to estimate the value of all weighting parameters
${\theta}$ in the game model as well as the trajectories x of all agents that
maximize the likelihood of a given sequence of state observations.
${\textbf{y}^{\mathcal{V}}:=(\textbf{y}^{i}),i\in\mathcal{V}}$:

$\displaystyle\max_{\theta,\textbf{x},\textbf{u}}\quad$ $\displaystyle
p(\textbf{y}^{\mathcal{V}}|\textbf{x},\textbf{u}),$ (5a) s.t.
$\displaystyle(\textbf{x},\textbf{u})\text{ is an OLNE of }\Gamma(\theta),$
(5b) $\displaystyle(\textbf{x},\textbf{u})\text{ is dynamically feasible under
}f,$ (5c)

where ${\theta}$ is the tuple of weighting parameters over all agents, i.e.
${\theta:=(\theta^{1},\cdots,\theta^{M})}$, and
${p(\textbf{y}^{\mathcal{V}}|\textbf{x},\textbf{u})}$ denotes a likelihood
model based on observation ${\textbf{y}^{\mathcal{V}}}$. Note that this
formulation extends that of [1] for our problem. In particular, we have
constructed a new objective (5a) and kept the constraints (5b), (5c) the same
since the visibility to the observer does not interfere with the participants’
interaction.

To solve (5), we need to first solve the OLNE in the forward dynamic game (2).
Since the constraints are the same, akin to the method in [1], we construct
each player’s Lagrangian with additional Lagrange multipliers
${\lambda_{t}^{i}}$:

$\mathcal{L}^{i}=J^{i}+\sum_{t=1}^{T-1}{\lambda_{t}^{i}}^{\top}\left(x_{t+1}^{i}-f(x_{t}^{i},u_{t}^{i})\right).$

The first-order necessary conditions for the optimal solution are given by the
KKT conditions for each agent:

$\displaystyle\textbf{G}(\textbf{x}^{i},\textbf{u}^{i},\boldsymbol{\lambda}^{i}):=\begin{bmatrix}\nabla_{\textbf{x}^{i}}\mathcal{L}^{i}\\\
\nabla_{\textbf{u}^{i}}\mathcal{L}^{i}\\\
x_{t+1}^{i}-f(x_{t}^{i},u_{t}^{i}),t\in[T-1]\end{bmatrix}=\textbf{0},$ (6)
$\displaystyle\forall i\in[M].$

Thus (5b) and (5c) are replaced by (6) and the inverse dynamic game model for
occluded agents evolves into

$\displaystyle\max_{\theta,\textbf{x},\textbf{u},\boldsymbol{\lambda}}\quad$
$\displaystyle p(\textbf{y}^{\mathcal{V}}|\textbf{x},\textbf{u})$ (7a) s.t.
$\displaystyle\textbf{G}(\textbf{x}^{i},\textbf{u}^{i},\boldsymbol{\lambda}^{i})=\textbf{0},\forall
i\in[M].$ (7b)

The objective of (7) is to estimate the unknown weighting parameters and all
agents’ states and actions from the observation of only visible agents’
trajectories under noise corruption. To deal with the observation noise, we
assume that observations take white Gaussian noise, i.e.
${\textbf{n}_{t}\sim\mathcal{N}(\textbf{0},\Sigma)}$, and the observation of
visible agents
${\textbf{y}_{t}^{\mathcal{V}}:=\textbf{x}_{t}^{\mathcal{V}}+\textbf{n}_{t}}$
in this work. We then solve (7) by substituting the likelihood maximization
problem of (7a) with the negative log-likelihood minimization problem with
objective
${\sum_{t\in[T]}\sum_{i\in\mathcal{V}}\|y_{t}^{i}-x_{t}^{i}\|_{2}^{2}}$.

In the next section, we conduct simulation experiments to evaluate the
estimation performance of our method and how it functions under different
levels of observation noise.

##  V Simulation Experiments

We implement our proposed approach in YALMIP [26], a MATLAB interface for
mathematical programming. We use the open-source COIN-OR IPOPT algorithm [27]
as a low-level solver. We conduct Monte Carlo studies to analyze the
performance of our proposed method in several simulated scenarios.

###  V-A Experiment Setup

To demonstrate the performance and robustness of our method, we perform a
sequence of Monte Carlo studies. For each scenario, we fix the weighting
parameters ${\theta}$ for each agent and find the corresponding OLNE
trajectories. We hide the invisible agents’ trajectories and then corrupt the
visible agents’ trajectories with white Gaussian noise. To evaluate the
performance of our method at different levels of noise corruption, we generate
24 sets of random observation sequences at 21 different levels of noise. For
each of the resulting 504 observation sequences, we run our method to recover
estimates of weights for each agent as well as the trajectory for invisible
agents.

###  V-B Evaluation Metrics

To evaluate the performance of weighting parameter recovery, we first follow
[1] and measure the cosine dissimilarity between the unobserved true weighting
parameters $\theta_{\text{true}}$ and the estimated parameters
${\theta_{\text{est}}}$

$D(\theta_{\text{true}},\theta_{\text{est}}):=1-\frac{1}{M}\sum_{i\in[M]}\frac{{\theta_{\text{true}}^{i}}^{\top}\theta_{\text{est}}^{i}}{\|\theta_{\text{true}}^{i}\|_{2}\|\theta_{\text{est}}^{i}\|_{2}},$
(8)

which denotes the degree of dissimilarity between $\theta_{\text{true}}$ and
$\theta_{\text{est}}$.

To evaluate the performance of trajectory estimation for all agents, we
measure the trajectory estimation error between the unobserved true trajectory
based on true weighting parameters and the estimated trajectory based on
recovered weighting parameters for both visible and invisible agents with the
average displacement error (ADE) metric:

$\displaystyle ADE_{\text{visible}}$
$\displaystyle:=\frac{1}{|\mathcal{V}|\cdot
T}\sum_{i\in\mathcal{V}}\sum_{t\in[T]}\|x_{\text{GT},t}^{i}-x_{\text{est},t}^{i}\|_{2},$
(9) $\displaystyle ADE_{\text{invisible}}$
$\displaystyle:=\frac{1}{|\mathcal{O}|\cdot
T}\sum_{i\in\mathcal{O}}\sum_{t\in[T]}\|x_{\text{GT},t}^{i}-x_{\text{est},t}^{i}\|_{2},$

where ${x_{\text{GT},t}^{i}}$ and ${x_{\text{est},t}^{i}}$ denote the ground
truth state and the estimated position of the $i^{\text{th}}$ agent at time
step $t$, respectively. Note that we measure the trajectory estimation error
separately for visible and invisible agents, rather than measuring the total
error, as in [1]. Results reveal the great difference of estimation
performance between the two categories, and hence the necessity of separation.

###  V-C Simulation Experiments

####  V-C1 Identifying an agent with constant velocity

We evaluate our method in a Monte Carlo study of the running example of a
visible agent avoiding collision with an occluded agent with constant
velocity. This example represents a simplified scenario where the occluded
agent does not react to other agents and keeps his velocity constant. In this
example, the problem degrades from solving a dynamic game problem to solving
the following forward optimal control problem

$\displaystyle\min_{\textbf{x},\textbf{u}^{1}}$ $\displaystyle
J^{1}(\textbf{u};\textbf{x}_{1}),$ $\displaystyle\forall i\in[M],$ (10) s.t.
$\displaystyle\textbf{x}_{t+1}^{1}=x_{t}^{1}+u_{t}^{1}\Delta t,$
$\displaystyle\textbf{x}_{t+1}^{2}=x_{t}^{2}+v\Delta t,$ $\displaystyle\forall
t\in[T-1],$

where $v$ takes an unknown constant value, and the inverse optimal control
problem

$\displaystyle\max_{\theta,\textbf{x},\textbf{u},v}$ $\displaystyle
p(\textbf{y}^{1}|\textbf{x},\textbf{u},v)$ s.t.
$\displaystyle(\textbf{x}^{1},\textbf{u}^{1})\text{ is optimal under
\eqref{optimal_control}},$ $\displaystyle(\textbf{x},\textbf{u},v)\text{ is
dynamically feasible under }f.$

Figure 3a shows the parameter estimation performance at different levels of
observation noise. We present the median and the interquartile range (IQR) for
parameter dissimilarity. Our method accurately recovers the unobserved
weighting parameters with which the model best explains the observed noisy
trajectories of visible agents. The result also reveals that our method is
noise-robust since the error remains small ($\times 10^{-3}$ order of
magnitude) even as the noise increases.

Benefiting from the accurate recovery of weighting parameters, our method also
outputs accurate trajectory estimates for both visible and invisible agents.
Figure 3d displays the trajectory estimation performance at different levels
of observation noise. The median estimation error for the visible agent is
under $0.15\text{\,}\mathrm{m}$ and that for the invisible agent is under
$0.35\text{\,}\mathrm{m}$ at all levels of observation noise.

Note that the trajectory estimation error for visible agents is on the order
of the standard deviation of the observation noise, which means that
observation noise contributes the majority of the error. The narrow IQR also
indicates that our method achieves high precision in trajectory estimation for
visible agents. In contrast, the trajectory estimation for invisible agents
suffers from greater median error and larger IQR. The extra error and
uncertainty is due to the lack of observation. This fact reflects the
importance of observation in robot operating scenarios. Still, it is
interesting to note that the estimation error of the occluded agent’s
trajectory is around twice the error of visible agents’ trajectory, which
means our method is fairly accurate.

####  V-C2 Identifying an agent applying a Nash strategy

(a)

(b)

(c)

Figure 2: Demonstration for the 2-agent running example in Section V-C2 at
noise with standard deviation $\sigma=0.05$ for 24 different observation
sequences, where both agents apply open-loop Nash strategy. 2a: Ground truth
trajectory for both agents. 2b: Observations from the observer’s perspective,
where the trajectories of the red agent are corrupted by noise and the blue
agent is invisible. 2c: Trajectories estimated by solving the recovered
dynamic game based on noisy observations. Despite the large observation noise
in 2b, our method can reliably estimate even the invisible agent’s trajectory
in 2c.

In this example, we evaluate our method for collision avoidance with two
agents, both of whom apply open-loop Nash strategies. Figure 2a shows the
ground truth trajectory of both agents. Figure 2b shows the sensor
measurements from the observer’s perspective: the trajectories of the visible
agent are corrupted by noise, while the invisible agent cannot be observed.
Figure 2c shows the predicted trajectories for both agents based on the game
model recovered by our method.

Figure 3b and Figure 3e display the estimation performance for the weighting
parameters and the trajectory for both agents at different levels of
observation noise. Compared with the example in Section V-C1, the median
parameter estimation dissimilarity and median trajectory estimation error for
both agents are higher, and IQR for both metrics is larger in this example. We
attribute this degradation to this example’s added complexity. It is
reasonable to infer that the estimation performance decreases with the number
of agents.

####  V-C3 Identifying an occluded pedestrian in traffic

Next, we consider a more realistic road crossing scenario, depicted in Figure
1. In this scenario, each agent not only tries to keep a safe distance from
other agents but also seeks to move forward while keeping himself in the lane
according to traffic rules. Figure 1 demonstrates that the blue pedestrian and
the driver in the red vehicle are observable to each other. However, the blue
pedestrian is occluded by the red vehicle and hence invisible to the driver in
the green vehicle. Both the blue pedestrian and the red vehicle apply open-
loop Nash strategies: the red vehicle brakes and the blue pedestrian sightly
turns left to avoid collisions. Although the movement of the blue pedestrian
is invisible to the driver in the green vehicle, he can make the inference
that someone is running across the road and predict that agent’s trajectory
based on the observed deceleration of the red vehicle next to him. He can then
also brake in advance to avoid a potential collision with the blue pedestrian.

Figure 3c and Figure 3f display the estimation performance for the weighting
parameters and the trajectory for both agents at different levels of
observation noise. Our method still achieves accurate and noise-robust
estimation of the unknown weighting parameters as well as the visible agent’s
trajectory. Compared with the simpler example in Section V-C2, a more
complicated game model contributes to a greater estimation error of the
invisible agent’s trajectory. Note that in the traffic scenario, the visible
agent’s motion is highly restricted by the traffic rules, and we only have
access to part of the trajectory observation (first 5 steps), therefore the
estimation error of the visible agent’s trajectory is lower than that of the
two agent example in Section V-C2.

##  VI Conclusion & Future Work

In this work, we have proposed a novel method based on the inverse dynamic
game technique to identify behaviors of invisible agents in occluded areas
from noise-corrupted observations of only visible agents. Our proposed method
recovers unobserved weighting parameters in the game model that best explain
the observed trajectories, and simultaneously computes open-loop Nash
trajectories for both visible and invisible agents. The computed results can
be further utilized in agents’ trajectory estimation and prediction. Numerical
results in simulation experiments show that our method is robust to
observation noise and provides an accurate estimation for both weighting
parameters and invisible agents’ trajectories.

Currently, we evaluate our method in simulation experiments. Future work
should investigate techniques to incorporate our method with first-person view
sensing results from real sensors (lidar, camera, etc.) so that our proposed
method can be utilized in real urban traffic scenarios.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Estimation performance for the our method in different examples at
different levels of observation noise. Figures on the top demonstrate the
parameter estimation dissimilarity and figures on the bottom demonstrate the
trajectory estimation error for each example. 3a and 3d: the example in
Section V-C1, where the invisible agent takes a constant velocity. 3b and 3e:
the example in Section V-C2, where both the visible and invisible agent apply
open loop Nash strategies. 3c and 3f: the example of the traffic scenario in
Section V-C3, where both the visible and invisible agent apply open loop Nash
strategies. In each figure, we plot the median dissimilarity (8), displacement
error (9) and their IQR over 24 samples at each level of observation noise.

## References

  * [1] L. Peters, D. Fridovich-Keil, V. Rubies-Royo, C. J. Tomlin, and C. Stachniss, “Inferring objectives in continuous dynamic games from noise-corrupted partial state observations,” _arXiv preprint arXiv:2106.03611_ , 2021.
  * [2] D. Fridovich-Keil, E. Ratner, L. Peters, A. D. Dragan, and C. J. Tomlin, “Efficient iterative linear-quadratic approximations for nonlinear multi-player general-sum differential games,” in _2020 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2020, pp. 1475–1481.
  * [3] K. Hara, H. Kataoka, M. Inaba, K. Narioka, and Y. Satoh, “Recognizing people in blind spots based on surrounding behavior,” in _Proceedings of the European Conference on Computer Vision (ECCV) Workshops_ , 2018, pp. 0–0.
  * [4] K. Hara, H. Kataoka, M. Inaba, K. Narioka, R. Hotta, and Y. Satoh, “Predicting vehicles appearing from blind spots based on pedestrian behaviors,” in _2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC)_. IEEE, 2020, pp. 1–8.
  * [5] O. Afolabi, K. Driggs-Campbell, R. Dong, M. J. Kochenderfer, and S. S. Sastry, “People as sensors: Imputing maps from human actions,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 2342–2348.
  * [6] M. Itkina, Y.-J. Mun, K. Driggs-Campbell, and M. J. Kochenderfer, “Multi-agent variational occlusion inference using people as sensors,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 4585–4591.
  * [7] R. E. Kalman, “When is a linear control system optimal?” in _Joint Automatic Control Conference_ , no. 1, 1963, pp. 1–15.
  * [8] I. Maroger, O. Stasse, and B. Watier, “Walking human trajectory models and their application to humanoid robot locomotion,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 3465–3472.
  * [9] S. Albrecht, K. Ramirez-Amaro, F. Ruiz-Ugalde, D. Weikersdorfer, M. Leibold, M. Ulbrich, and M. Beetz, “Imitating human reaching motions using physically inspired optimization principles,” in _2011 11th IEEE-RAS International Conference on Humanoid Robots_. IEEE, 2011, pp. 602–607.
  * [10] I. Maroger, O. Stasse, and B. Watier, “Inverse optimal control to model human trajectories during locomotion,” _Computer Methods in Biomechanics and Biomedical Engineering_ , vol. 25, no. 5, pp. 499–511, 2022.
  * [11] A. Y. Ng, S. Russell _et al._ , “Algorithms for inverse reinforcement learning.” in _Icml_ , vol. 1, 2000, p. 2.
  * [12] B. D. Ziebart, N. Ratliff, G. Gallagher, C. Mertz, K. Peterson, J. A. Bagnell, M. Hebert, A. K. Dey, and S. Srinivasa, “Planning-based prediction for pedestrians,” in _2009 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2009, pp. 3931–3936.
  * [13] K. M. Kitani, B. D. Ziebart, J. A. Bagnell, and M. Hebert, “Activity forecasting,” in _Computer Vision–ECCV 2012: 12th European Conference on Computer Vision, Florence, Italy, October 7-13, 2012, Proceedings, Part IV 12_. Springer, 2012, pp. 201–214.
  * [14] M. Monfort, A. Liu, and B. Ziebart, “Intent prediction and trajectory forecasting via predictive inverse linear-quadratic regulation,” in _Twenty-Ninth AAAI Conference on Artificial Intelligence_ , 2015.
  * [15] F. Previtali, A. Bordallo, L. Iocchi, and S. Ramamoorthy, “Predicting future agent motions for dynamic environments,” in _2016 15th IEEE International Conference on Machine Learning and Applications (ICMLA)_. IEEE, 2016, pp. 94–99.
  * [16] P. Henry, C. Vollmer, B. Ferris, and D. Fox, “Learning to navigate through crowded environments,” in _2010 IEEE international conference on robotics and automation_. IEEE, 2010, pp. 981–986.
  * [17] M. Pfeiffer, U. Schwesinger, H. Sommer, E. Galceran, and R. Siegwart, “Predicting actions to act predictably: Cooperative partial motion planning with maximum entropy models,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 2096–2101.
  * [18] B. Kim and J. Pineau, “Socially adaptive path planning in human environments using inverse reinforcement learning,” _International Journal of Social Robotics_ , vol. 8, pp. 51–66, 2016.
  * [19] H. Kretzschmar, M. Spies, C. Sprunk, and W. Burgard, “Socially compliant mobile robot navigation via inverse reinforcement learning,” _The International Journal of Robotics Research_ , vol. 35, no. 11, pp. 1289–1307, 2016\.
  * [20] L. Y. Morales, A. Naoki, Y. Yoshihara, and H. Murase, “Towards predictive driving through blind intersections,” in _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_. IEEE, 2018, pp. 716–722.
  * [21] L. Sun, W. Zhan, C.-Y. Chan, and M. Tomizuka, “Behavior planning of autonomous cars with social perception,” in _2019 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2019, pp. 207–213.
  * [22] S. Le Cleac’h, M. Schwager, and Z. Manchester, “Lucidgames: Online unscented inverse dynamic games for adaptive trajectory prediction and planning,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 3, pp. 5485–5492, 2021\.
  * [23] L. Peters, V. Rubies-Royo, C. J. Tomlin, L. Ferranti, J. Alonso-Mora, C. Stachniss, and D. Fridovich-Keil, “Learning players’ objectives in continuous dynamic games from partial state observations,” _arXiv preprint arXiv:2302.01999_ , 2023.
  * [24] J. Li, C.-Y. Chiu, L. Peters, S. Sojoudi, C. Tomlin, and D. Fridovich-Keil, “Cost inference for feedback dynamic games from noisy partial state observations and incomplete trajectories,” _arXiv preprint arXiv:2301.01398_ , 2023.
  * [25] T. Başar and G. J. Olsder, _Dynamic noncooperative game theory_. SIAM, 1998.
  * [26] J. Lofberg, “Yalmip: A toolbox for modeling and optimization in matlab,” in _2004 IEEE international conference on robotics and automation (IEEE Cat. No. 04CH37508)_. IEEE, 2004, pp. 284–289.
  * [27] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical programming_ , vol. 106, no. 1, pp. 25–57, 2006.


${\bf SolarNeutronDecayProtonsobservedinNovember7,2004}$

Yasushi Muraki1,5,∗111*)correspondence<EMAIL_ADDRESS>Jose F.
Vald$\acute{e}s$-Galicia2, Ernesto Ortiz3, Yutaka Matsubara1, Shoichi
Shibata4, Takashi Sako1,6, Satoshi Masuda1, Munetoshi Tokumaru1, Tatsumi Koi4,
Akitoshi Ooshima4,5, Takasuke Sakai7, Tsuguya Naito8, and Pedro Miranda9
222Full author information is available at the end of the article.

###### Abstract

We have found an interesting event registered by the solar neutron telescopes
installed at high mountains in Bolivia (5250 m a.s.l.) and Mexico (4600 m
a.s.l.). The event was observed November 7th of 2004 in association with a
large solar flare of magnitude X2.0. Some features in our registers and in two
satellites (GOES 11 and SOHO) reveal the presence of electrons and protons as
possible products of neutron decay. Solar neutron decay protons (sndp) were
recorded on board ISEE3 satellite in June 3rd, 1982 . On October 19th, 1989,
the ground level detectors installed in Goose Bay and Deep River revealed the
registration of solar neutron decay protons (sndp). Therefore this is the
second example that such an evidence is registered on the Earth′s surface.

${\bf keywords}$: Solar neutron decay protons, Solar flare, Solar Energetic
Particles, Particle acceleration

##  1 Introduction

Gamma rays and neutrons propagate freely in the interplanetary medium when
emitted as secondary products of solar explosion events. Therefore they may
provide information regarding the condition of production site and mechanism
acceleration. Solar neutrons have been observed in space crafts and ground
based detectors (see, e.g. Galicia et al. 2009, Dorman 2010, Kamiya et al
2019, Muraki et al. 2020 and reference there in). In order to push forward
this study, we have international solar neutron telescope (SNT) network at
high mountains in the world. Evenson et al. (1983) reported the discovery of
interplanetary protons by the decay of solar flare neutrons. Shea et al.(1991)
found signals in neutron monitors for the event on 19 October 1989 that may be
interpreted as the detection of relativistic protons that were the decay
products of solar neutrons. Observing protons produced by the neutron decay in
flight, very accurate energy spectra may be obtained. These spectra will be
close to the source spectra.

In this paper we report registers obtained at two Solar Neutron Telescopes
installed at high mountains, combined with the observations of two spacecrafts
(GOES 11 and SOHO). They could be interested as products of the solar neutrons
decay that were produced at the 7 November 2004 X2 solar flare.

The plan of the paper is as follows: in section 2 we give a brief description
of the Solar Neutron Telescopes (SNT) at Bolivia and Mexico; section 3 is
dedicated to a description of the observations in spacecrafts and on the
Earth′s surface. In section 4, we give a plausible interpretation of the
results. Then section 5 is dedicated to convert the fluxes observed on earth′s
surface to fluxes on the top of the atmosphere to compare with the spacecraft
observations. The sites where solar neutrons may have decayed are analyzed in
section 6, to finish with our conclusions in section 7.

##  2 The two mountain detectors located on American continent

Two solar neutron detectors are located at Mt. Sierra Negra (4,600 m, 19.0∘N
97.3∘W) and Mt. Chacaltaya (5,250 m, 16.3∘S 68.1∘W). The detector installed at
Mt. Chacaltaya in Bolivia is comprised of 4 $m^{2}$ plastic scintillator with
40 cm thickness and the detector installed at Mt. Sierra Negra, Mexico is
composed of the 4 $m^{2}$ plastic scintillator with 30 cm thickness. In
addition, four layers of proportional counters are installed underneath the
scintillator for the identification of the arrival direction of charged
particles (Valdes-Galicia et al., 2004). These SNTs were constructed for the
detection of solar neutrons. Both instruments distinguish neutral incidents
from charged particles.

##  3  Main features of the observed events

1\. $\it Signals{\,}recorded{\,}in{\,}the{\,}ground{\,}level{\,}detectors$:
increases of the counting rate were recognized in the 5-minute-value of the
solar neutron telescopes (SNT) located in Bolivia from 15:48 to 16:06 UT, and
in Mexico from 15:51 to 17:09 UT of November 7th, 2004. The statistical
significance of each excesses was 3.7$\sigma$ and 12$\sigma$ respectively.
Both data are shown in Figures 1 and 2. The excesses of signal of the SNTs
were observed almost at the same time, assuming them to be due to the solar
flare of X2.0. Therefore, the increase cannot be explained by very plausibly
originated at fluctuations of the background signal.

2\. $\it
Signals{\,}recorded{\,}in{\,}the{\,}detectors{\,}onboard{\,}the{\,}satellites$:
Figure 3 shows the counting rate of the protons with $>$10 MeV observed by the
GOES 11 satellite. An increase is apparently recognized between 15:50 UT and
16:00 UT (16$\sigma$). An appropriate explanation on this bump is that the
detector onboard the GOES 11 satellite received solar neutron decay protons
(sndp). Figure 4 presents another candidate of neutron decay products detected
by the electron detector COSTEP onboard the SOHO satellite that stays near the
Lagrange point L1 (SOHO-COSTEP-EPHIN website). The excess was recorded by the
COSTEP instrument during 15:58 UT - 16:17 UT. Solar neutron decay protons
(sndp) were also observed by the COSTEP. In the proton channels 4.3-7.8 MeV
and 7.8-25 MeV, the excess of the sndp was clearly observed. However, in the
highest channel 25-35 MeV, the excess was not definitely recognized. On the
other hand, in the electron channels 0.25-0.70 MeV, 0.67-3.0 MeV and 2.64-10.4
MeV, the excess was clearly see. These may be another evidence of neutron
decay products via n$\rightarrow p+e^{-}+\nu_{e}$ process. Therefore neutron
decay protons and electrons were observed not only inside but also outside of
the magnetosphere.

##  4  Interpretations on the event

We make the following hypotheses in order to explain the observations.

(a) The increase of the 5-minute counting rate registered during 15:48-16:07
UT at Chacaltaya was induced by solar neutrons. These neutrons were produced
by the impulsive flare with X2.0. One of the most probable production time was
presumed when the derivative of the GOES X-ray intensity showed the first
maximum at 15:47:00 UT (Figure 5). Figure 5 shows that the flare reached X2.0
intensity via three step increase. The tendency is quite similar to a previous
event observed in April 15, 2001 (Muraki et al. 2008).

(b) The increase of the counting rate observed between 15:51 and 17:09 UT at
Mt. Sierra Negra was produced by solar neutron decay protons (sndp) in the
space between the CME front and the Earth. Here we specify the preceding CME
as CME1. Since two large CMEs were produced: one in November 6th at 00:34 UT
and another in 7th at 16 UT. Possible acceptable area for the sndp by each
detector is pictorially shown in Figure 6. The distance D1 in Figure 6 is
estimated to be $10^{7}$ km. First, the sndp produced near the magnetosphere
arrived, then the sndp produced near the CME1 transported to the Earth. This
is due to the fact that charged particles produced in the interplanetary space
are transported toward the Earth along the inter-planetary magnetic field.
That is one of the reasons why the excess of Mt. Sierra Negra continued for 78
minutes and was slightly delayed with respect to the neutron arrival time. In
comparison with the June 3rd, 1982 event, the sndp of November 7th, 2004 event
was observed rather short period ($\sim$1/10). One of the reasons may due to
the difference of decay space of the sndp by each detector. An explanation
will be given in the next section.

The hypotheses (a) and (b) were introduced for two main reasons. One arises
from the difference of the atmospheric thickness of the neutron path length
between the two observatories. At the flare time, there were 200 g/$cm^{2}$
difference of atmospheric thickness between the two observatories for the
passage of neutrons (Dorman et al., 1999, Tsuchiya et al. 1999, 2001).
Therefore the flux at Mt. Sierra Negra should be 7 times less than that of
Chacaltaya ($\sim$1/7), although both detectors are located at nearly the same
altitude. (Let us remind that the Sun was situated above South America.)
However, both excesses were detected with nearly the same intensity.

The second reason is based on the difference of the event duration. The signal
of Chacaltaya was observed for 18 minutes, while the excess of Sierra Negra
continued for 78 minutes (as shown by the horizontal arrow in Figure 2).
Neutrons with 50 MeV energy are expected to arrive to the Earth 18 minutes
later than the fastest neutron, if neutrons were produced. Neutrons with
energy 50 MeV are possible to detect by the Chacaltaya detector, since the
threshold energy is set at higher than $>$40 MeV. On the other hand, neutrons
with the energy of 30 MeV arrive at the Earth after 22 minutes (1,300 seconds)
later than the highest energy neutron. The threshold energy of the Sierra
Negra SNT was set at $>$30 MeV. However, the increase continued for 78
minutes. Therefore it would be difficult to explain the signals observed at
Mt. Sierra Negra by the direct hit of solar neutrons.

The excess of the counting rate at Chacaltaya was $\sim$100
events/($m^{2}$minute) in the $>$40 MeV channel between 15:48 UT and 16:06 UT.
On the other hand, the S1 channel of Sierra Negra (the channel of charged
particles with the energy higher than $>$30 MeV) showed 133 events /
($m^{2}$·minute) as the excess between 15:51-16:27 UT. The flux was observed
with nearly the same intensity as that of Chacaltaya, and the excess start
time was almost the same. However, the excess duration was completely
different. Therefore, in order to explain both enhancements consistently, we
may introduce the assumption that the excess Chacaltaya was produced by the
direct arrival of solar neutrons, while the excess Sierra Negra was produced
by solar neutron decay protons (sndp).

We may take into account another effect for the estimation of the flux of the
sndp, since the energy of sndp is expected to be around a few GeV. After
transportation of the sndp in the magnetosphere, some of them may be trapped
by the Earth′s magnetic field, but some of them will penetrate and arrive at
Mt. Sierra Negra. The characteristic energy is called as the cut-off energy or
Sp$ddot{o}$rer limit (rigidity) for low energy protons. The cut-off energy of
Mt. Sierra Negra and Mt. Chacaltaya are estimated as 3$\sim$4 GeV and
11$\sim$12 GeV respectively (Shea and Smart, 2000).

The differential energy spectrum of protons near the cut-off energy has been
measured in space by PAMELA (2009, 2016) and AMS (2000) detectors
independently. The result is shown in Figure 7 by the green triangles for the
cut-off energy of 3-4 GeV. Therefore we will make an expected flux of sndp
near the cut-off energy (the black circles in Figure 7), by multiplying the
PAMELA′s observed differential flux of 3-4 GeV (the green triangles in Figure
7) (Casolino 2007, PAMELA 2009) by the expected neutron energy spectrum of
En-4dEn (the red diamond in Figure 7). The energy spectrum of sndp beyond the
cut-off energy is expected to reflect the neutron energy spectrum. However, in
the low energy region less than the cut-off energy, the energy spectrum of
sndp is predicted to have almost constant value to the neutron spectrum of
En-4dEn. The magnetic latitude of the observatory is estimated as 30∘N. The
difference between the geographical latitude arises from the location of the
magnetic-pole located at northern Canada (79∘N) over$\sim$260∘E line (100∘W
line) in 2001.

(c) The increases observed by the GOES 11 satellite during 15:50-16:00 UT
(Figure 3) was also produced by the neutron decay protons. They were decay
products of high energy neutrons in the energy range between 80 and 400 MeV.
If these neutrons were produced at 15:47 UT, from the observed time, the
parent neutrons had the energy between 80 and 400 MeV.

It may be interesting to know that the flux of neutron decay protons differs
two geostationary satellites. The increase was not observed in the detector
onboard the GOES 10 satellite. As for the longitude of both satellites on the
Earth, the GOES 11 was situated at 114∘W, while the GOES 10 was located at
135∘W. In other words, the GOES 11 satellite was located in the right above
American continent (just over the longitude of Mexico City).

(d) Another bump of electron and proton components around 16 UT was observed
by the COSTEP instrument onboard the SOHO satellite (Figure 4). They were also
produced by the neutron decays in the outer magnetosphere. Now we focus on the
neutron decay electrons. Observation of electrons with energies higher than
2.64 MeV implies that initial neutrons should have a Lorentz factor
$\gamma$=1.75 according to earlier predictions (Koi et al. 1993, Dorman 2010).
It implies that original neutrons must have the kinetic energy higher than 700
MeV. So initial protons must be accelerated beyond 1 GeV at the impulsive
phase of the flare to produce high energy neutrons, En = 700 MeV. The
kinematics of neutron decay protons and electrons is given in Appendix 1.

##  5 Flux Conversion observed by the ground level detectors to the top of
atmosphere

In order to compare each flux, in this section, we convert the flux measured
by the ground level detectors to the top of the atmosphere.

a) ${\it
Deriving{\,}flux{\,}of{\,}solar{\,}neutrons{\,}at{\,}the{\,}top{\,}of{\,}atmosphere{\,}over{\,}Chacaltaya.}$

In this section, let us derive the differential energy spectrum of solar
neutrons at the top of the atmosphere. At 16 UT ( the local time at 12 LT),
the atmospheric thickness over Chacaltaya is estimated as 550 g/$cm^{2}$. The
differential energy spectrum is derived with use of one-minute value of the
counting rate of the channel with E$>$40 MeV of the solar neutron detector. In
order to derive the differential energy spectrum, we introduce the hypothesis
that neutrons were produced impulsively at the solar atmosphere. The flight
time of neutrons depends on its kinetic energy. We fixed the production time
at 15:47:00 UT. Then, the kinetic energy of neutrons is determined from the
flight time (=arrival time $–$15:47:00 UT), so that the differential flux of
solar neutrons is derived as a function of energy. The observed energy
spectrum is given in Figure 8(a).

After we derived the observed energy spectrum of solar neutrons by the
detector, we converted it to the flux at the top of the atmosphere. For this
purpose we used two correction factors, the detection efficiency of solar
neutrons by the SNT (Watanabe 2005) and the neutron attenuation curve in the
atmosphere (Shibata 1994). Combining these two correction factors into one
curve, a correction curve is made. The results is shown in Figure 8(b) as a
function of incoming neutron energy.

Dividing the energy spectrum observed at Mt. Chacaltaya (Figure 8a) by the
correction curve (Figure 8b), finally the differential energy spectrum of
solar neutrons at the top of the atmosphere of Chacaltaya has been derived.
The differential energy spectrum at the top of the atmosphere is presented in
Figure 8(c). The differential flux has been already normalized to the flux per
unit area ($/m^{2}$).

b) $\it
Flux{\,}of{\,}solar{\,}neutron{\,}decay{\,}protons{\,}over{\,}at{\,}Sierra{\,}Negra$

Here we estimate the flux of the sndp over Mt. Sierra Negra. We use the
observed results in space by PAMELA (2007, 2009) and AMS (2000). These results
show that cosmic rays less energy than the cut-off energies have been
observed. The green points of Figure 7 represent such effect. The events
recorded by the Mt. Sierra Negra Solar Neutron Telescope between 15:51 UT and
15:54 UT may correspond to the highest of the solar neutron decay protons
(sndp). Therefore, they may carry information on the flux near the cut-off
energy of protons around Ecut$\sim$3 GeV. Incoming protons with energy about
$\sim$3 GeV make nuclear interactions with air nuclei. As a result, neutral
pions are produced with an energy around $\sim$1 GeV. These neutral pions
immediately decay into two gamma rays and the gamma rays will make the
electromagnetic cascade shower in the air. The tails of the cascade shower
enter into the scintillator and gamma rays involved in the shower are
converted into electrons and positrons. From the intensity of electron and
positron signals, the arrival intensity of protons over Mt. Sierra Negra may
be estimated.

Figure 7 (the black circles) suggests that the intensity of the sndp between
Ep=1-3 GeV may be observed with an equal intensity and the spectrum is
estimated to be flat. Therefore, we choose energy bin width ($\Delta$E) of
3000 MeV (=4000-1000 MeV). Then we are able to estimate the differential flux
of the sndp at 1-4 GeV as (4.4$\pm$3.0) / ($m^{2}$.min.MeV) at the top of the
atmosphere. On the other hand the flux of solar neutrons at the top of the
atmosphere over Mt. Chacaltaya is estimated as 12.8, 1.8, and 0.5
events/($m^{2}$.min.MeV) at En=1, 2 and 3 GeV respectively. If we compare
these numbers, the estimated solar neutron decay proton flux over Mt. Sierra
Negra shows the Chacaltaya neutron flux. Given the assumptions made, we could
say that the agreement is fairly good. Each flux is summarized in Table 1.

We are now preparing end to end simulation based on the GEANT4 code. In this
simulation, The primary energy spectrum measured by the PAMELA detector and
the attenuation of gamma rays in the atmosphere will be taken into account. We
expect that with these simulations the error bars on the flux of sndp will be
reduced. Furthermore if you look carefully to Figure 2, you may notice another
enhancement: the counting rate of the anti-counter registering charged
particles with energy higher than 30 MeV, and the lower detector located under
the scintillator (L1 channel) indicate an increase from 18:00 UT (14 LT in
Bolivia and 12 LT in Mexico). It is possible that protons were further
accelerated into high energies by the shock acceleration mechanism (Tsuneta
and Naito 1998). We will discuss this matter in the next paper, together with
the results obtained by the GEANT4 simulation.

##  6  Production place of solar neutron decay protons

Before we discuss the production place of sndp in the interplanetary space,
let us summarize the general situation of the interplanetary space around 16UT
on November 7th 2004. Around the flare time, a very short gamma ray burst was
observed by the gamma ray burst monitors onboard INTEGRAL and WIND satellite
(INTEGRAL web site, WIND web site). According to INTEGRAL SPIACS data, a very
short gamma ray burst was recorded at 15:49:30 UT. The KONUS detector onboard
the WIND satellite also detected the short gamma ray signal and they detected
2.2 MeV line gamma rays. The 2.2 MeV line gamma rays are emitted when a
neutron is captured by a proton and to form a deuterium. The RHESSI satellite
detected hard X-rays after 16:05 UT. In Figure 9, we present two images of the
flare at 15:36 UT and 15:48 UT. The images were taken by an ultraviolet
telescope onboard SOHO. Within 16 minutes, an arch was emphasized in the image
at 15:48 UT. It is shown by the white artificial arc in the left side of
Figure 9. Particles may possibly be accelerated within this arc.

Now let us describe briefly the characteristics of the CME1. According to the
magnetometer measurement onboard ACE satellite (ACE web site) that stays at
the Lagrange point L1, the maximum field strength of the CME1 was 40 nT. This
first CME1 (on the right side in Figure 6) was produced by the M9.3 flare of
November 6th, not by the X2.0 flare of November 7th. According to the IPS
observation, the CME2 arrived on the Earth around 9UT of November 9th
(Tokumaru 2013). The CME2 produced at 16UT is depicted near the Sun in Figure
6. The CME1 was already expanded in the interplanetary space at 16 UT on
November 7. From the record of the ACE satellite, the diameter of the CME1 may
be estimated. We estimated it as $\ell$1= $10^{7}$ km. The width of the CME1
was calculated by multiplying the time (9,000 sec) by speed (1,150 km/sec).
Then we can estimate the maximum momentum of charged particles that will be
trapped inside the CME1. By putting the numerical values into a simple
equation pc=300 H$\rho$, as H=4$\times$10-4gauss, $\rho$= 0.5 $\times$1012cm,
we get pc=6$\times 10^{10}$ eV=60 GeV. The value is quite high. In other
words, low energy charged particles may be trapped by the “CME barrier”, so
that neutron decay protons produced between the Sun and the CME1 could not
arrive the detector located near the Earth. Therefore, the observed sndp must
be produced in between the front side of CME1 and the Earth. This region is
pictorially shown in Figure 6 as D1. (Drawing is not proportional to the
actual distance. ) The direct line approximation of the decay path may be
guaranteed by the preceding calculations, at least during 60 minutes (Sakai
and Muraki 1993, Sakai et al. 1997).

We examine “another barrier” for charged particles. It is the wall produced by
the magnetosphere. The length of the magnetic-sheath $\ell$3 is estimated as
$\ell$3= 2.7$\times$$R_{earth}$ = 17,000 km = 1.7$\times$109cm. The field
strength is estimated as between 20 nT $–$ 40 nT. Therefore, we choose 30 nT.
Again, putting these values into the equation of pc=300 H$\rho$, we will get
pc= 15 MeV. This time, the threshold energy is quite low, however protons with
the energy less than 15 MeV produced between the CME1 and the bow shock cannot
penetrate inside the magnetosphere. In present case, the initial neutrons
observed by the GOES detector had an energy between 80 $-$ 400 MeV. Therefore
the decay products (sndp) between Ep=80 MeV and 400 MeV can penetrate the
magneto-sheath. (For the reference, we provide the length of $\ell$2 and
$\ell$4 as l2=$10^{7}$ km =0.066 AU, and $\ell$4= 56,700 km respectively. The
flight length from the front edge of the CME1 to the mountain detector D1 is
estimated as D1= 1$\times 10^{7}$km + 67,000 km$\approx$1$\times$107 km.)

In the events observed in 1982 June 3, the sndp were registered for almost 12
hours (Evenson et al. 1983, 1990), while in the present event the excess was
observed only during 75 minutes. The difference may be due to the existence of
the preceding CME (CME1). In the 1982 June 3rd flare, there was no preceding
CME (Solar Geophysical Data website). However in the event of 1989 October
19th, the sndp were observed for 30 minutes and the signals of sndp were
masked by the strong accelerated proton beam produced by the X13 flare (Shea
et al. 1991).

##  7 Conclusive Remarks

We have here presented evidence corresponding to registers in two spacecrafts
and two ground based detectors. Taken together, observations admit an
interpretation that would be consistent with the observation of solar neutron
decay protons. Solar neutron decay protons were first reported in 1981 by Paul
Evenson, Peter Meyer, Roger Pyle. Ruffolo and other researchers discussed
solar neutron decay electrons and protons afterwards (Ruffolo 1991,
Dr$\ddot{o}$ge et al. 1996). However, we do not know any report of this kind
of events after that of 19, October 1989. To the best of our knowledge, there
are only three early reports on the detection of neutron decay protons and
electrons onboard satellites and one report of neutron decay protons by ground
level detectors in October 19th, 1989. The event of October 1989, was also
discussed (Koi et al. 1993). So, present event may be the second case where
neutron decay protons were registered by a ground-based detector. Further
study is necessary to understand this phenomenon deeply. However, we wanted to
publish this quick report to call community attention in the search of neutron
decay protons and electrons.

##  8 Competing interests

The authors declare that they have no competing interests.

##  9 Author contributions

YM, JFVG, YMa, SS, TSako, and TS have constructed the Mt. Sierra Negra Solar
Neutron Telescope, while YMa, TSako, YM and PM constructed the Chacaltaya
solar neutron detector. YM, JFVG, EO and SS have made the data analysis. SM
and MT provided the solar image and solar wind data respectively. SS prepared
Appendix 1. All members jointed discussions.

##  10 Acknowledgments

The authors acknowledge INAOE personnel and authorities for their continued
support in providing the services needed to keep the Mt. Sierra Negra SNT
functioning well. We also acknowledge UNAM-PA‘IIT partial support through
grant IN-104115. The authors express the acknowledgment to the staffs of
physics department of UMSA, who have kept the Solar Neutron Detector located
at Mt. Chacaltaya under good condition for 28 years. We express sincere thanks
to the satellite teams; particularly SOHO (LASCO, COSTEP, and EIT instruments)
and ACE satellites, who provided valuable data available from the web site.
This work has been carried out, being based on the international solar neutron
telescope (SONTEL) network data reserved in the storage system of ISEE, Nagoya
University.

##  11 Author details

1) Institute for Space Earth Environmental Research, Nagoya University, Nagoya
464-8601, Japan, 2) Instituto de Geofisica, UNAM, 04510, Mexico D. F., Mexico,
3) Escuela Nacional de Ciencias de la Tierra, UNAM. Ciudad de Mexico. 04510.
Mexico, 4) Engineering Science Laboratory, Chubu University, Kasugai, Aichi
487-0027, Japan, 5) Astronomy Observatory, Chubu University, Kasugai, Aichi
487-0027, Japan, 6) Institute for Cosmic Ray Research, The University of
Tokyo, Kashiwa, Chiba 277-8582, Japan, 7) Physical Science laboratory, College
of Industrial Technologies, Nihon University, Narashino, Chiba 275-0006,
Japan, 8) Department of Information Science, Yamanashi Gakuin University, Kofu
400-8575, Yamanashi, Japan, and 9) Department of Physics, UMSA, La Paz,
Bolivia

##  12 References

$\ast$) ACE web site;
http:// www.srl.caltech.edu/ACE/ASC/DATA/level3/mag/ACECpec.cgi?LATEST=1
$\ast$) Alex A. (2018) Script Lecture PHY 432, Physics with muons: from atomic
physics to Solid state physics (psi ch.) (Access is available from website)
$\ast$) AMS collaboration; Alearaz et al. (2000) Protons in near earth orbits.
Physics Letters B 472:215-226 https://doi.org/10.1016/S0370-2693(99)01427-6
$\ast$) Burman R.L. and Smith E.S (1989) LA-11502-MS, Parametrization of Pion
production and Reaction cross-section at LAMPF energies (Access is available
from website.)
$\ast$) Casolino M et al. (2007) Two Years of Flight of the Pamela Experiment:
Results and Perspectives. J. Phys. Soc. Jpn. Suppl. A 78:35-40, M. Casolino et
al. (2007) Observation of primary, trapped and quasi trapped particles with
PAMELA experiment. Proceed. 30th ICRC (Merida) 1:709-712, Also in the
rapporteur talk by E.O. Flückiger (2007) Ground Level Events and Terrestrial
Effects. Proceed. 30th ICRC (Merida) 6: 239-253
$\ast$) Dorman LI, Valdés-Galicia JF, Dorman VI (1999) Numerical simulation
and analytical description of solar neutron transport in the Earth’s
atmosphere, Journal Geophysical Research 104 (A10): 22417-22426.
Doi.org/10.1029/1999jA900182
$\ast$) Dorman LI, Dorman IV, and Valdés-Galicia JF (1997) Simulation of solar
neutron scattering and attenuation in the Earth’s atmosphere for different
initial zenith angles. Proceed. 25th ICRC (Durban) 1:25-28
$\ast$) Dorman L (2010) Solar Neutrons and Related Phenomena (Springer). pp
374-377doi: 10.1007/978-90-481-3737-4
$\ast$) Droge W, Ruffolo D, and Klecker B (1996) Observation of electrons from
the decay of solar flare neutrons. ApJ 464: L87-L90
$\ast$) Evenson P, Meyer P, and Roger Pyle K (1983) Protons from the decay of
solar flare neutrons. ApJ 274:875-882 $\ast$) Evenson P, Kroger R, Meyer P,
and Reames D (1990) Solar Neutron Decay Proton Observations in Cycle 21. ApJ
Suppl. 73:273-277
$\ast$) INTEGRAL web site:
https://www.isdc.unige.ch/integral/science/grb$\sharp$ACS
$\ast$) Kamiya K, Koga K. Matsumoto H., Masuda S, Muraki Y., Tajima H, and
Shibata S (2019) Solar Neutrons observed from September 4 to 10 2017 by SEDA-
FIB, PoS (ICRC2019) 1150.
$\ast$) Koi T. et al.(1993) Prediction of electrons as decay products of solar
neutrons. Proceeding of 23rd ICRC (Calgary) 3:151-154
$\ast$) Muraki Y., Matsubara Y., Masuda S., Sakakibara S., Sako T., Watanabe
K., Bütikofer R., Flückiger E.O., Chilingarian A., Hovsepyan G., Kakimoto F.,
Terasawa T., Tsunesada Y., Tokuno H., Velarde A., Evenson P., Poirier J., and
Sakai T. (2008) Detection of high-energy solar neutrons and protons by ground
level detectors on April 15, 2001. Astroparticle Physics 29:229-242
Doi:10.1016/j.astropartphys.2007.12.007
$\ast$) Muraki Y. Valdés-Galicia JF, González LX, Kamiya K, Katayose Y, Koga
K, Matsumoto H, Masuda S, Matsubara Y, Nagai Y, Ohnishi M, Ozawa S, Sako T,
Shibata, S, Takita M, Tanaka Y, Tsuchiya H, Watanabe K, and Zhang JL (2020)
Possible detection of solar gamma-rays by ground-level detectors in solar
flares on 2011 March 7. Pub. Astron. Soc. Japan 72:1-17 Doi:
10.1093/pasj/psz141
$ast$) PAMELA collaboration; N. De Simone et al. (2009) Comparison of models
and measurements of protons of trapped and secondary origin with PAMELA
experiment. Proceed. 31st ICRC (Lodz) icrc0795
Adriani et al. (2016) PAMELA’s measurements of geomagnetic cutoff variations
during the 14 December 2006 storm.
Space Weather 14 210-220 doi:10.1002/2016SW001364
$\ast$) Ruffolo D (1991) Interplanetary transport of decay protons from solar
flare neutrons. Astrophysical Journal 388:688-698 Doi: 10.1086/170756
$\ast$) Sakai T. and Muraki Y. (1993) Solar neutron Decay Proton. Proceed of
23rd ICRC (Calgary) 3:147S. Website:
articles.ads.abc.harvard.edu/pdf/1993ICRC….3..1475
$\ast$) Sakai T, Kato M, and Muraki Y (1997) Propagation of Solar neutron
decay protons near the Earth. Geomag. Geoelectr. 49:1105-1113
$\ast$) Shea MA, Smart DF, Wilson MD, and Flückiger EO (1991) Possible ground-
level measurement of solar neutron decay protons during the 19 October 1989
solar cosmic ray event. Geophys. Res. Lett. 18 : 829-832.
Doi.org/10.1029/90GL02668
$\ast$) Shea MA and Smart DF (2000) Cosmic Ray Implications for Human Health.
Space Science Review 93:187-205
$\ast$) Shibata S (1994) Propagation of solar neutrons through the atmosphere
of the Earth. Journal of Geophysical Research 99: 6651-6665
$\ast$) SOHO CME catalog website: https://cdawweb.gsfc.nasa.gov
$\ast$)SOHO-COSTEP-EPHIN : http://www2.physik.unikiel.de/SOHO/phpeph/EPHIN.htm
$\ast$) Solar Geophysical Data web site:
https://www.ngdc.noaa.gov/stp
/space-weather/online-publications/stp sgd/1982/sgd8212c.pdf
$\ast$) Tsuchiya H.et al. (1999) Detection efficiency of new solar neutron
detector. Proceed. 26th ICRC (Salt Lake City) 7: 363
$\ast$) Tsuchiya H et al. (2001) Arrival of solar neutrons from large zenith
angle. Proceed. 27th ICRC (Hamburg) 8: 3056-3059
$\ast$) Tsuneta S. and Naito T (1998) Fermi acceleration at the fast shock in
a solar flare and the inclusive loop-top hard X-ray source. APJL 67:495
$\ast$) Tokumaru, M (2013) Three-dimensional exploration of the solar wind
using observations of interplanetary scintillation, Proceedings of the Japan
Academy Ser. B, 89(2):67-79 doi:10.2183/pjab.89.67.
$\ast$) Valdés-Galicia JF, Muraki Y., Sako T, Musalem O, Huertado A, Gonzalez
X, Matsubara Y, Watanabe K, Hirano N, Tateiwa N, Shibata S, and Sakai T (2004)
An improved solar neutron telescope installed at a very high altitude in
México. NIMA 535: 656-664 doi.org/10.1016/j.nima.2004.06.148
$\ast$) Valdés-Galicia JF, Muraki Y, Watanabe K, Matsubara Y, Sako T, Gonzalez
X, Musalem O, Huertado A (2009) Solar neutron events as a tool to study
particle acceleration at the Sun, Advances in Space Research, 43: 565-572.
https://doi.org/10.1016/j.asr2008.09.023
$\ast$) Watanabe K. (2005) Proceeding of the Cosmic-ray Research Section of
Nagoya University 46(2), pp1-249. Solar neutron Events associated with Large
Solar Flares in Solar Cycle 23 (PhD thesis in English)
$\ast$) WIND web site: https://asd.gsfc.nasa.gov

##  13 Appendix 1

Appendix 1 —kinematics of neutron decay electrons —

Let us define the parameters that are used for the estimation of the highest
energy of electrons from neutrons. Here $m_{n}$, $m_{p}$, and $m_{e}$
represents the mass of neutron, proton, and electron respectively. $E_{p}^{*}$
and $E_{e}^{*}$ are the energy of proton and electron in the center of
momentum system (in the neutron rest frame). In this calculation, we take the
neutrino mass zero. Now we study an extreme case; proton and neutrino are
emitted backward, while electron is emitted forward. Then the electron energy
in the center of momentum system is expressed by

$E_{e}^{*}=((m_{n}c^{2}-m_{p}c^{2})^{2}+me^{2}c^{4})/(2(m_{n}c^{2}-m_{p}c^{2})).\qquad(1)$

We know $(m_{n}c^{2}-m_{p}c^{2})=1.29\,{\rm MeV}$ and $m_{e}c^{2}=0.511\,{\rm
MeV}$. Putting these values in Equation (1), we get $E_{e}^{*}=0.746\,{\rm
MeV}$.

From a relation, $(P_{e}^{*}c)^{2}=E_{e}^{*2}-m^{2}c^{4}$, and
$P_{e}^{*}c=0.544\,{\rm MeV}$.

The Lorentz transform is expressed as

$E_{e}=\gamma(E_{e}^{*}+\gamma\beta P_{e}^{*}c\cos{{\theta}^{*}}),\qquad(2)$

where $\gamma$ represents the Lorentz factor of incident neutron.
So the electron energy in the laboratory frame Ee can be expressed by

$E_{e}=\gamma(0.746+\gamma\beta\times 0.544)\,{\rm MeV}\ \ {\rm in\ case}\
\cos{{\theta}^{*}}=1.\qquad(3)$

The threshold energy of COSTEP detector onboard SOHO satellite is 2.64 MeV.
Therefore solving the equation, we will get the minimum Lorentz factor of
solar neutrons($\gamma$) by

$2.64=\gamma(0.746+\gamma\beta\times 0.554).\qquad(4)$

When we put $\gamma=1.75$ and $\beta=0.821$ in Eq.(4), we will get
$E_{e}\approx 2.67\,{\rm MeV}$.

Let us provide another extreme case, i.e., the proton is emitted forward and
the electron and neutrino are emitted backward. The maximum energy of the
proton in the center of momentum system is expressed by

$E_{p}^{*}=(m_{n}^{2}c^{4}+m_{p}^{2}c^{4}-m_{e}^{2}c^{4})/2m_{n}c^{2}.\qquad(5)$

Since $E_{p}^{*}=m_{p}c^{2}$ and $P_{p}^{*}\approx 0$, in Eq.(2), we put
$E_{p}^{*}$ instead of $E_{e}^{*}$,
then we get

$E_{p}=\gamma m_{p}c^{2}=(m_{p}/m_{n})E_{n}.\qquad(6)$

Therefore the proton emitted the forward direction takes almost all the energy
(0.999) of the neutron.

Figure 1:  (left). The 5-minute-value of the counting rate of the solar
neutron detector located at Mt. Chacaltaya in Bolivia (5,250m a.s.l.). From
top to bottom, it presents the counting rate of charged particles with the
deposit energy higher than $>$40, $>$80, $>$120, and $>$160 MeV respectively.
The threshold energy of the detector is calibrated by using the deposit energy
of the minimum ionizing particles like muons. The arrow corresponds to the
time of the enhanced counting rate, around 16 UT and the vertical dotted line
represents the GOES flare start time. The horizontal presents the local time.
Figure 2:  (right). The 5-minute-value of the counting rate of the solar
neutron telescope located at Mt. Sierra Negra in Mexico (4,600m a.s.l.). From
top to bottom, the picture corresponds to the counting rate of all charged
particles, charged particles with energy higher than 30 MeV, neutral particles
with energy higher than 30 MeV, and the proportional counter located
underneath the scintillator. We call it L1 channel (Layer 1). The L1 channel
is triggered with the signal of S1 ($>$30MeV). The L1 channel is used for the
identification of neutrons and gamma rays. The vertical dotted line presents
the GOES flare start time (UT) and the horizontal arrow indicates the time of
the enhanced counting rate.  Figure 3: (top). The time variation of the proton
intensity with energy higher than 10 MeV measured by the proton counter
onboard the GOES 11 satellite. An enhancement can be recognized 15:50-16:00 UT
(indicated by an arrow). The enhancement corresponds to the impulsive phase of
the solar flare X2.0 and may be produced by the neutron decay protons.  Figure
4: (bottom) The time variation of the electron intensity with the energy
2.64-10.4 MeV measured by the COSTEP detector onboard the SOHO satellite. The
enhanced time is indicated by the arrow. If the enhancement was produced by
solar neutron decay electrons, the initial neutron should have the Lorentz
factor $\gamma$=1.75 according to Koi et al. (1993).  Figure 5: The one-minute
time profile of GOES X-ray intensity from 15:20 UT to 17:10 UT November 7,
2004. The red line and the blue line correspond to the X-rays with the wave
length 1-8${\AA}$ and 0.5-4${AA}$ respectively. The green line represents the
derivative of the short band of X-rays. We assume that neutrons were produced
instantaneously at 15:47:00 UT, when there is a change of slope of the
emission.  Figure 6: The interplanetary plasma scenario around 16 UT of
November 7th, 2004 is depicted. had not yet arrived near the Earth at 16 UT,
situated about 1${\times}10^{7}$ km (${\ell}$2) away from the (${\ell}$3) of
the Earth. The front of the magnetosphere was extended at about 9${R_{e}arth}$
(${\ell}$4). The green star presents the position of the GOES satellite and
the brown mark corresponds to Mt. Sierra Negra observatory. D1 represents the
distance from the CME front to the mountain laboratory, and D2
(=2.2${\times}10^{4}$ km) depicts the distance from the bow of the magneto
pause to the GOES satellite Figure 7: The expected solar neutron decay proton
(sndp) flux at Mt. Sierra Negra (the black circles). The red diamonds
represent a possible production spectrum of solar neutrons for the energy
spectrum of En-4dEn. The green triangle presents the observed proton flux by
the PAMELA detector near the cut-off energy of 3-4 GeV. Beyond the cut-off
energy, the energy spectrum of the sndp represents the neutron spectrum
Figure 8: (a) (left) the original energy spectrum observed at Mt. Chacaltaya,
(b) (center) the correction factor for the energy spectrum to the observed
spectrum. (c) (right) The estimated energy spectrum at the top of the
atmosphere over Mt. Chacaltaya in units of flux/(cm2·min. MeV). Figure 9:
Solar images taken by the EIT telescope onboard SOHO satellite at 15:36 UT
(left) and 15:48 UT (right). The UV wave length is 19.5 nm. We assume that
neutrons were instantaneously produced at 17:47 UT. If you compare the two
images (a) and (b), you will notice a plasma arcade in the image of 15:48 UT.
The arc is shown by a white arc in the image of 15:36 UT and possibly it had
an important role for the particle acceleration. Figure 10: The sky projection
maps of the solar wind speed derived from the observations of the
interplanetary scintillation device (IPS) using the 327-MHz multi-station
system of Nagoya University (Tokumaru, 2013). The data were accumulated during
the periods (left side) between 22 UT of November 7, and 7 UT of November 8,
2004, and (right side) between 22 UT of November 8 and 7 UT of November 9,
2004 respectively. Four IPS antennas observed IPS for the same radio sources
simultaneously. The solar wind speed was determined by detecting the time lag
between the variations of the intensity at the separated stations. The center
of the map corresponds to the location of the Sun, and circles indicate the
relative positions of lines-of-sight for radio sources. The radius of the
circle represents the solar wind speed. Red circles denote high speed ($>$1000
km/s) plasma flow. The solar offset distances of high-speed data for November
8-9 (right panel) are larger than those for November 7-8 (left panel). This
fact suggests that the expanding CME was approaching near the Earth. T he
shock in association with the X2.0 flare arrived on Earth at 09 UT in November
9th Figure 11:  The image of the CME-CME interaction taken by the SOHO LASCO
telescope at 16:54 UT on November 7th, 2004. The central white circle
represents the solar surface. Figure 12:  The CME emission log from November
6th 00 UT to November 8th 00 UT observed by the SOHO coronagraph. Two high
speed CMEs are recognized in the plot in association with the two large solar
flares at November 6th (at 00: 34UT) and November 7th (at 16 UT). Before the
emission of the CME at 16 UT, another CME was emitted around 12 UT of November
7th with rather slow speed. As a result, they will collide at 17 UT. The
collision time is indicated by the vertical line in the figure. The CME
emission log is shown together with the GOES X-ray data as reference. The
numbers in the GOES X-ray data indicate the flare position on the solar
surface. The data are available from SOHO home page. Figure 13: The flux of
X-rays observed by the RHESSI satellite. At 18:15 UT a small bump can be
recognized.


# What Does DALL-E 2 Know About Radiology?

Lisa C. Adams
Department of Radiology
Stanford University School of Medicine
Stanford, CA, USA
<EMAIL_ADDRESS>
Felix Busch11footnotemark: 1
Department of Radiology
Charité – Universitätsmedizin Berlin
Berlin, Germany
<EMAIL_ADDRESS>
Daniel Truhn
Department of Diagnostic and Interventional Radiology
University Hospital Aachen
Aachen, Germany
<EMAIL_ADDRESS>
Marcus R. Makowski
Department of Diagnostic and Interventional Radiology
School of Medicine and Klinikum Rechts der Isar
Technical University of Munich
Munich, Germany
<EMAIL_ADDRESS>
Hugo JWL. Aerts
Artificial Intelligence in Medicine (AIM) Program
Mass General Brigham, Harvard Medical School
Boston, MA, USA

Departments of Radiation Oncology and Radiology
Dana-Farber Cancer Institute and Brigham and Women’s Hospital
Boston, MA, USA

Radiology and Nuclear Medicine,
CARIM & GROW, Maastricht University
Maastricht, the Netherlands
<EMAIL_ADDRESS>
Keno K. Bressem
Department of Radiology
Charité – Universitätsmedizin Berlin
Berlin, Germany

Berlin Institute of Health at Charité – Universitätsmedizin Berlin
Berlin, Germany
<EMAIL_ADDRESS>

contributed equally

###### Abstract

Generative models such as DALL-E 2 could represent a promising future tool for
image generation, augmentation, and manipulation for artificial intelligence
research in radiology provided that these models have sufficient medical
domain knowledge. Here we show that DALL-E 2 has learned relevant
representations of X-ray images with promising capabilities in terms of zero-
shot text-to-image generation of new images, continuation of an image beyond
its original boundaries, or removal of elements, while pathology generation or
CT, MRI, and ultrasound images are still limited. The use of generative models
for augmenting and generating radiological data thus seems feasible, even if
further fine-tuning and adaptation of these models to the respective domain is
required beforehand.

_Keywords_ Deep Learning $\cdot$ DALL-E 2 $\cdot$ Diffusion Models $\cdot$
Radiology

##  1 The potential impact of generative models in radiology

DALL-E 2 is a novel deep learning model for text-to-image generation, first
introduced by OpenAI in April 2022 (Ramesh et al., 2022). The model has
recently gained widespread public interest due to its ability to create
photorealistic images solely from short written inputs (Kather et al.,
2022)(Conwell and Ullman, 2022)(Marcus et al., 2022). Trained on billions of
text-image pairs extracted from the internet, DALL-E 2 learned a wide range of
representations, which it can recombine to create novel images, exceeding the
variability it has seen in the training data, even from implausible prompts
(e.g., “a painting of a pelvic x-ray by Claude Monet”)(Ramesh et al.,
2022)(Nichol et al., 2021). Given these powerful generative capabilities of
DALL-E 2 it raises the question if these can be transferred to the medical
domain to create or augment data, as it is often sparse and hard to come by.
Given its size, it is to be assumed, that the training data for DALL-E 2 also
contained radiological images and the model might have learned about
composition and structure of X-rays and maybe even cross-sectional imaging and
ultrasound.
Therefore, this comment systematically explores the radiological domain
knowledge embedded in DALL-E 2, by creating and manipulating radiological
images of X-rays, CT, MRI and ultrasound.

##  2 Generating radiological images from short descriptive texts

To explore the general radiological domain knowledge embedded in DALL-E 2, we
first tested the how well radiological images can be created from short text
descriptive texts using the wording “An x-ray of” and one-word description of
the anatomical area. For the head, chest, shoulder, abdomen, pelvis, hand,
knee, and ankle, we let DALL-E 2 synthesize four artificial X-ray images,
which are displayed in Figure 1. The coloring of the images and the general
structure of the bones appeared realistic, and the overall anatomical
propositions were correct, indicating the presence fundamental concepts of
X-ray anatomy. However, upon close inspection, it was noted that the
trabecular structure of the bone appeared random and did not follow the course
of mechanical stress as in real X-rays. Sometimes smaller bones were either
missing, like the fibula in multiple X-rays of the knee, or fused into a
single instance like several carpal or tarsal bones. Still, the overall
structure of the depicted anatomic area was recognizable. In rare cases,
additional bones and joints were generated, e.g., an extra index finger in a
hand X-ray. The model had the most difficulties generating joints correctly as
they appeared to be either fused or the joint surfaces were distorted.
Compared to X-rays, the generation of cross-sectional images from CT or MRI
was very limited (See Figure 2). While the requested anatomical structure was
conceivable in the images, the overall impression was primarily nonsensical.
Instead of a cross-sectional slice, a photo of an entire CT or MRI examination
printed on radiographic film was often displayed. For ultrasound, the
requested anatomical structures were not recognizable, instead all images
resembled obstetric images. Still, the images showed concepts of
CT/MRI/ultrasound images, indicating that representations of these modalities
are present in DALL-E 2.

##  3 Changing areas in radiological images while maintaining the anatomical
structure

To investigate how detailed the radiological knowledge is, we tested how well
DALL-E 2 can reconstruct missing areas in a radiological image (inpainting).
Therefore, we chose radiographs of the pelvis, ankle, chest, shoulder, knee,
wrist, and thoracic spine and erased specific areas before providing the
residuals to DALL-E 2 for restoration (see Figure 3). The accompanying prompts
were identical to the text-based imaging described above. DALL-E 2 provided
realistic image replacements that were nearly indistinguishable from the
original radiographs for the pelvis, thorax, and thoracic spine. However, the
results were not as convincing when a joint was included in the image area.
For the ankle and wrist, the number of tarsal bones and the structure varied
greatly and deviated from realistic representations. For the shoulder, DALL-E
2 failed to reconstruct the glenoid cavity and articular surface of the
humerus. In one case, a foreign body was incorporated into the shoulder that
remotely resembled a prosthesis. For knee reconstruction, the model omitted
the patella but retained the bicondylar structure of the femur.

##  4 Extending radiological images beyond image borders

As the previous two experiments showed, that DALL-E 2 possesses fundamentals
of standard radiological images, we now aimed to explore how well anatomical
knowledge is embedded in the model. For this we randomly selected X-rays of
the abdomen, chest, pelvis, knee, different spine areas, wrist bones, and hand
let DALL-E 2 extend these images beyond their borders, as this task requires
the model to understand the location as well as anatomical proportions. The
accompanying prompt was chosen depending on the anatomical regions being
generated. Again, the style of the augmented images represented a realistic
depiction of X-ray (see Figure 4). It was possible to generate a complete
whole-body X-ray using just an image of the knee as a starting point. However,
the further the distance between the original image and the generated area,
the less detailed the images became. Anatomical proportions, such as the
length of the femoral bone or size of the lungs, remained realistic, still
finer details such as the number of lumbar vertebrates were inconsistent. The
model generally performed best when creating anterior and posterior views,
while lateral views were more challenging and produced poorer results.

##  5 Generation of pathological images

The generation of pathological images, e.g., showcasing fractures,
intracranial bleedings, or tumors was limited. We tested this with the
generation of bone fractures on X-ray, which mostly led to distorted images of
normal radiographs, similar to those shown in Figure 1. Furthermore, DALL-E 2
has a filter that prevents the generation of harmful content, preventing us
from testing most pathologies, as words like “hemorrhage” triggered this
filter.

##  6 Conclusion and Outlook

We could show that DALL-E 2 can generate X-ray images that are close to
authentic X-ray images in style and anatomical proportions. From this, we
conclude that relevant representations for X-rays were learned by DALL-E 2
during training. However, it seems only representations for radiographs
without pathologies are present or admitted by the filter, as generating
pathological images, e.g., of fractures and tumors, was limited. Beyond,
generative capabilities for cross-sectional imaging as well as ultrasound were
poor.
Access to data is critical in deep learning and generally, the larger the
dataset, the better the performance of the models trained upon it. However,
especially in radiology, there exists not unified large database from which
such dataset could be created but instead data is split between multiple
institutions and privacy concerns additionally prevent merging of this data.
Here, synthetic data brings the promise to solve these issues, allowing the
creation much larger datasets than those currently available substantially
accelerating the development of new radiological deep learning tools (Jordon
et al., 2020)(Jordon et al., 2022). Previous approaches on data generation
using generative adversarial networks (GANs) were challenged by instabilities
during training and the inability to capture sufficient diverytiy of the
training data (Kodali et al., 2017). In contrast novel diffusion based
generative models, of which DALL-E 2 is a prominent representative, solved
these issues, beating GANs on image generation (Dhariwal and Nichol, 2021).
Furthermore, the amount of control and flexibility diffusion models can
provide during data generation, makes them exciting for radiological data
augmentation, as images may be altered, and pathologies may be added or
removed just by text prompts. However, our experiments showed, that DALL-E 2
was not yet able to create convincing pathological images, which is likely
explained by a lack of those images in the training data. Radiological images
containing a wide variety of pathologies are more likely to be shared on
professional sites, might be under restricted access or in special data
formats such as DICOM. Altogether probably leading to only a small number of
those images included on the training data.
Therefore, more work is needed to explore diffusion models designed to create
and modify medical data. Although DALL-E 2 is not openly available, other
models of similar architecture, such as Stable Diffusion, are available to the
public(Ho and Salimans, 2022)(Saharia et al., 2022). Future research should
focus on fine-tuning these models to medical data and incorporate medical
terminology to create powerful models for data generation and augmentation in
radiological research.

## References

  * Ramesh et al. [2022] Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen.  Hierarchical text-conditional image generation with clip latents.  _arXiv preprint arXiv:2204.06125_ , 2022.
  * Kather et al. [2022] Jakob Nikolas Kather, Narmin Ghaffari Laleh, Sebastian Foersch, and Daniel Truhn.  Medical domain knowledge in domain-agnostic generative ai.  _NPJ digital medicine_ , 5(1):1–5, 2022.
  * Conwell and Ullman [2022] Colin Conwell and Tomer Ullman.  Testing relational understanding in text-guided image generation.  _arXiv preprint arXiv:2208.00005_ , 2022.
  * Marcus et al. [2022] Gary Marcus, Ernest Davis, and Scott Aaronson.  A very preliminary analysis of dall-e 2.  _arXiv preprint arXiv:2204.13807_ , 2022.
  * Nichol et al. [2021] Alex Nichol, Prafulla Dhariwal, Aditya Ramesh, Pranav Shyam, Pamela Mishkin, Bob McGrew, Ilya Sutskever, and Mark Chen.  Glide: Towards photorealistic image generation and editing with text-guided diffusion models.  _arXiv preprint arXiv:2112.10741_ , 2021.
  * Jordon et al. [2020] James Jordon, Alan Wilson, and Mihaela van der Schaar.  Synthetic data: Opening the data floodgates to enable faster, more directed development of machine learning methods.  _arXiv preprint arXiv:2012.04580_ , 2020.
  * Jordon et al. [2022] James Jordon, Lukasz Szpruch, Florimond Houssiau, Mirko Bottarelli, Giovanni Cherubin, Carsten Maple, Samuel N Cohen, and Adrian Weller.  Synthetic data–what, why and how?  _arXiv preprint arXiv:2205.03257_ , 2022.
  * Kodali et al. [2017] Naveen Kodali, Jacob Abernethy, James Hays, and Zsolt Kira.  On convergence and stability of gans.  _arXiv preprint arXiv:1705.07215_ , 2017.
  * Dhariwal and Nichol [2021] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  _Advances in Neural Information Processing Systems_ , 34:8780–8794, 2021.
  * Ho and Salimans [2022] Jonathan Ho and Tim Salimans.  Classifier-free diffusion guidance.  _arXiv preprint arXiv:2207.12598_ , 2022.
  * Saharia et al. [2022] Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily Denton, Seyed Kamyar Seyed Ghasemipour, Burcu Karagol Ayan, S Sara Mahdavi, Rapha Gontijo Lopes, et al.  Photorealistic text-to-image diffusion models with deep language understanding.  _arXiv preprint arXiv:2205.11487_ , 2022.

##  7 Figures

Figure 1: Sample generated anatomical structures in X-ray from short text
descriptions using DALL-E 2. Figure 2: Sample text-to-image generated
anatomical structures in CT, MRI, and ultrasound using DALL-E 2. Figure 3:
Reconstructed areas of different anatomical locations in X-rays using DALL-E
2. The yellow-bordered regions of the original images were erased before
providing the remnant image for reconstruction. Figure 4: Extending X-ray
images of different anatomical regions beyond their borders using DALL-E 2.
The original X-ray is shown in the yellow boxes, with the remainder of the
images being generated by DALL-E 2.