# Scalable Multi-Agent Reinforcement Learning through
Intelligent Information Aggregation

Siddharth Nayak 1, Kenneth Choi 1, Wenqi Ding 1,
Sydney Dolan 1, Karthik Gopalakrishnan 2, Hamsa Balakrishnan 1
1 Massachusetts Institute of Technology
2 Stanford University
{sidnayak, kenchoi, wenqi22, sydneyd<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

We consider the problem of multi-agent navigation and collision avoidance when
observations are limited to the local neighborhood of each agent. We propose
InforMARL, a novel architecture for multi-agent reinforcement learning (MARL)
which uses local information intelligently to compute paths for all the agents
in a decentralized manner. Specifically, InforMARL aggregates information
about the local neighborhood of agents for both the actor and the critic using
a graph neural network and can be used in conjunction with any standard MARL
algorithm. We show that (1) in training, InforMARL has better sample
efficiency and performance than baseline approaches, despite using less
information, and (2) in testing, it scales well to environments with arbitrary
numbers of agents and obstacles.

##  1 Introduction

Reinforcement Learning (RL) has seen wide-ranging successes recently in high-
dimensional robot control DDPG ; RLVisuomotor , solving physics-based control
problems contControlRL , playing Go alphaGo , Chess alphaZero and Atari video
games DQN_Atari ; Atari2 , etc. However, challenges remain in many real-world
applications in which the tasks cannot be handled by a single agent, e.g.,
multi-player games, search-and-rescue drone missions, etc. MARLSurvey . In
such cases, multiple agents may need to work together and share information in
order to accomplish the task MARLseminal . Naïve extensions of single-agent RL
algorithms to multi-agent settings do not work well because of the non-
stationarity in the environment, i.e. the actions of one agent affect the
actions of others IQL1 ; IQL2 . Furthermore, tasks may require cooperation
among the agents. Classical approaches to optimal planning may (1) be
computationally intractable, especially for real-time applications, and (2) be
unable to account for complex interactions and shared objectives between
multiple agents. The ability of RL to learn by trial-and-error makes it well-
suited for problems in which optimization-based methods are not effective. In
particular, multi-agent reinforcement learning (MARL) approaches may be
suitable in these situations due to their fast run-times and superior
performance, and their ability to model shared goals between agents using
appropriate reward structures.

In this paper, we focus on the multi-agent navigation and collision avoidance
problem, where there are $N$ agents navigating to their respective goals in a
2D environment with static and/or dynamic obstacles. The agents are assumed to
be collaborative, i.e. the rewards are shared across all agents. We assume
that an agent can only sense the presence of obstacles or other agents within
a certain limited radius $r$. The overarching objective is for all the agents
to reach their goals in the shortest time possible, while avoiding collisions
with other agents and obstacles. This problem setting is quite general and
arises in many contexts, _e.g_., search-and-rescue robot teams search_rescue ,
environmental monitoring env_monitoring , and drone delivery systems
drone_delivery1 ; drone_delivery2 ; drone_delivery3 .

MARL-based techniques have achieved significant successes in recent times,
_e.g_., DeepMind’s AlphaStar surpassing professional level players in
StarCraft II AlphaStar , OpenAI Five defeating the world-champion in Dota II
Dota2 , etc. The performance of many of these MARL algorithms depends on the
amount of information included in the state given as input to the neural
networks MAPPO . In many practical multi-agent scenarios, each agent aims to
share as little information as possible to accomplish the task at hand. This
structure naturally arises in many multi-agent navigation settings, where
agents may have a desired end goal but do not want to share their information
due to communication constraints or proprietary concerns satelliteInfoSharing1
; satelliteInfoSharing2 . These scenarios result in a decentralized structure,
as agents only have locally available information about the overall system’s
state. In this paper, we focus on the question: “Can we train scalable multi-
agent reinforcement learning policies that use limited local information about
the environment to perform collision-free navigation effectively?"

We propose InforMARL, an approach for solving the multi-agent navigation
problem using graph-reinforcement learning. The main features of our approach
are that it: (1) uses a graph representation of the navigation environment
which enables local information-sharing across the edges of the graph; (2)
transfers well to different numbers of agents; and (3) achieves better sample
complexity in training compared to prior approaches by aggregating relevant
local information from neighbors in the underlying graph.

More broadly, our work (1) demonstrates that graphs provide a valuable
abstraction for multi-agent navigation environments; (2) highlights that more
information (i.e., global information as states) may not necessarily improve
performance, and can, in fact, overwhelm the RL agent networks and lead to
increased sample complexity; and (3) shows how graph architectures can
identify the most valuable information for navigation from local observations
to improve performance and scalability.

The rest of this paper is organized as follows. In Section 2, we briefly
summarize recent developments in MARL under different settings. In Section 3,
we describe InforMARL, our method for multi-agent navigation in settings with
limited observability and communication. In Section 4, we show the
effectiveness of InforMARL in learning policies for navigation by comparing
its performance to other baseline algorithms. We also perform experiments on
transferring InforMARL to scenarios with a different number of entities in the
environment. Finally, in Section 5, we discuss the main conclusions of this
paper and describe promising directions for future work.

##  2 Related Work

In this section, we discuss prior approaches to scaling MARL algorithms, as
well as some recent work on information-sharing between agents in multi-agent
systems.

###  2.1 Scaling MARL

Research on scaling MARL algorithms has broadly followed two main themes: (1)
decentralized execution, and (2) transferring learning between scenarios.

Decentralized MARL: A popular approach to improve scaling is the centralized-
training-decentralized-execution (CTDE) framework. Typically, CTDE frameworks
utilize actor-critic methods ActorCriticSeminal , where the training step uses
a centralized critic that incorporates global information from all actors.
However, during execution, the agents use their own actor networks to select
their actions in a decentralized manner to improve scalability. MADDPG MADDPG
builds upon DDPG DDPG by learning a centralized critic that is provided the
joint state and actions of all agents. MATD3 MATD3 uses a double-centralized
critic model, which helps reduce the over-estimation bias. In MAPPO , the
authors show the effectiveness of PPO in several standard multi-agent
environments. In value-based methods, VDN VDN decomposes a centralized value
function to a sum of individual agent-specific functions. Q-Mix QMIX improves
upon this by imposing a monotonicity requirement on agents’ individual value
functions and using a learnable mixing of the individual functions instead of
only summing them as done in VDN. These MARL algorithms perform well in the
navigation environment when they have access to information about the
positions of all entities in the environment. But, as we demonstrate in this
work, they fail to learn when that information is restricted to just local
neighborhoods around the agents.

Transferability in MARL: It is desirable to have MARL formulations where the
number of agents/entities in the environment doesn’t hinder the performance of
the model. Most of the previous works using MARL for the navigation task
require concatenation of observations of other entities in the environment to
be able to learn meaningful policies. As the neural network size depends on
the state input dimensions used while training, the learned policy fails to
work in scenarios with a different number of entities in the environment. With
the recent success of graph neural networks, many recent works have focused on
leveraging the inherent graph structure present in multi-agent problems and
tackling the issue of transferability. Zhou et al. GNN_IL create a
neighborhood graph according to how close the entities are to each other and
use imitation learning to imitate a greedy behavior for the target-coverage
problem targetCoverage1 ; targetCoverage2 . Similarly, Khan et al. GPG use
graph neural networks with vanilla policy gradient policyGradient for the
formation flying task formationFlying1 ; formationFlying2 . They show that
dynamic graphs have worse performance than static graphs due to the large
number of possible graphs the model has to learn. At the beginning of each
episode, their model determines the connectivity of the agents with each other
and uses the same graph over the whole episode. This works well for the
formation flying task as the graph structure does not change much if the
agents fly in the same formation over an episode.

DGN DGN combines a graph convolutional neural network GCN architecture with
multi-head attention AttentionIsAllYouNeed for a variety of multi-agent
environments, including 2D coverage and tracking. The DGN architecture assumes
that each agent communicates with its three closest neighbors. Communicating
with the three closest agents leads to a graph that is always connected.
However, the three-agent connectivity assumption is highly restrictive in real
life applications, as communication is generally restricted by mutual
separation due to hardware constraints.

We use a distance-based agent-entity graph similar to Entity-Message Passing
(EMP) EMP . However, the EMP architecture assumes that the agents have access
to the positions of all entities in the environment at the beginning of the
episode, which is not possible in the navigation setting where there are
occluded static or dynamic obstacles. By contrast, our model does not require
information on the positions of all entities at the beginning of the episode.

###  2.2 Information Sharing for MARL

For various cooperative tasks where explicit coordination is required to solve
the task, enabling communication to share information across the agents helps
with the performance. In CommNet CommNet , the authors introduce a model to
learn a differentiable communication protocol between multiple agents.
However, this communication does not explicitly model the interactions between
the agents and rather used an averaged state of all the neighboring agents.
VAIN VAIN improves upon CommNet by using an exponential kernel-based attention
to choose which specific messages from other agents to attend to. Similarly,
ATOC ATOC and TarMAC TarMAC use an attention mechanism for communication among
agents but without any restrictions on which agent can communicate with which
others, leading to a centralized setting during execution. GAXNET GAXNet also
uses an attention mechanism AttentionIsAllYouNeed , but additionally allows
for the exchanging of attention weights with other agents to reduce the
attention mismatch between them. While they show that their model works well
for the navigation task, their formulation requires the maximum number of
agents in the environment to be fixed before training. Our work also uses an
attention mechanism for inter-agent communication. However, we do not have any
message passing between objects in the environment (also referred to as
’entities’) and agents; instead, the agents themselves do all the computation
with local information of the entities’ states.

##  3 Description of InforMARL

Our MARL framework for navigation, InforMARL, consists of four modules, as
shown in Figure 1. We describe each in detail below.

Figure 1: Overview of InforMARL. (i) Environment: The agents are depicted by
green circles, the goals are depicted by red rectangles, and the unknown
obstacles are depicted by gray circles. $x^{(i)}_{agg}$ represents the
aggregated information from the neighborhood, which is the output of the GNN.
A graph is created by connecting entities within the sensing-radius of the
agents. (ii) Information Aggregation: Each agent’s observation is concatenated
with $x^{(i)}_{\mathrm{agg}}$. The inter-agent edges are bidirectional, while
the edges between agents and non-agent entities are unidirectional. (iii)
Graph Information Aggregation: The aggregated vector from all the agents is
averaged to get $X_{\mathrm{agg}}$. (iv) Actor-Critic: The concatenated vector
$[o^{(i)},x^{(i)}_{\mathrm{agg}}]$ is fed into the actor network to get the
action, and $X_{\mathrm{agg}}$ is fed into the critic network to get the
state-action values. The code for InforMARL is available at:
https://github.com/nsidn98/InforMARL

###  3.1 Environment

Every object in the environment (also known as an ‘entity’) is assumed to be
either an agent, an obstacle, or a goal. We define an agent-entity graph with
respect to agent $i$ at each time-step $t$, as
$g^{(i)}_{t}\in\mathcal{G}:(\mathcal{V},\mathcal{E})$, where each node
$v\in\mathcal{V}$ is an entity in the environment. The variable
entity_type(j)$\in\\{\texttt{agent},\texttt{obstacle},\texttt{goal}\\}$
determines the type of entity at node $j$. There exists an edge
$e\in\mathcal{E}$ between an agent and an entity if they are within a ‘sensing
radius’ $\rho$ of each other. The agent-agent edges are bidirectional, whereas
the agent-non-agent edges are unidirectional (i.e. messages can only be passed
from the non-agent entity to the agent). In other words, a unidirectional edge
is equivalent to the agent sensing a nearby entity’s state, while a
bidirectional edge is equivalent to a communication channel between agents.
This structure is similar to the agent-entity graph defined in EMP , but
without the assumption that all agents have access to the positions of all
entities at the beginning of each episode. Note that our formulation also
supports cases where disconnected sub-graphs are formed due to the positioning
of the entities in the environment.

The corresponding Decentralized Partially Observable Markov Decision Process
(Dec-POMDP) MDP ; POMDP1 ; POMDP2 is characterized by the tuple $\langle
N,S,O,\mathcal{A},\mathcal{G},P,R,\gamma\rangle$, where $N$ is the number of
agents, $s\in S=\mathbb{R}^{N\times D}$ is the state space of the environment
and $D$ is the dimension of the state, and
$o^{(i)}=O(s^{(i)})\in\mathbb{R}^{d}$ is the local observation for agent $i$,
where $d\leq D$ is the observation dimension. $a^{(i)}\in\mathcal{A}$ is the
action space for agent $i$ and the joint action for all $N$ agents is given by
$A=(a^{(1)},\cdots,a^{(N)})$. Specifically, $a^{(i)}$ is a one-hot vector of
size equal to the number of possible actions. $g^{(i)}\in\mathcal{G}(s;i)$ is
the graph network formed by the entities in the environment with respect to
agent $i$. $P(s^{\prime}|s,A)$ is the transition probability from $s$ to
$s^{\prime}$ given the joint action $A$. $R(s,A)$ is the joint reward
function. $\gamma\in[0,1)$ is the discount factor. The MARL training process
seeks to find an optimal policy,
$\Pi=\left(\pi^{(1)},\cdots,\pi^{(N)}\right)$, where each agent uses a policy
$\pi_{\theta}^{(i)}\left(a^{(i)}|o^{(i)},g^{(i)}\right)$ parameterized by
$\theta$ to determine its action $a^{(i)}$ from its local observation
$o^{(i)}$ and the graph network $g^{(i)}$ that it is a part of, while
optimizing the discounted accumulated reward
$J(\theta)=\operatorname{\mathbb{E}}_{A_{t},s_{t}}\left[\sum\limits_{t}\gamma^{t}R\left(s_{t},A_{t}\right)\right]$.

Agent $i$’s local observation $o^{(i)}$ consists of its position and velocity
in a global frame of reference and the relative position of the agent’s goal
with respect to its position. Each node $j$ on the graph $g^{(i)}$ has node
features
$x_{j}=[p^{j}_{i},v^{j}_{i},p^{\mathrm{goal},j}_{i},\texttt{entity\\_type(j)}]$
where $p^{j}_{i},v^{j}_{i},p^{\mathrm{goal},j}_{i}$ are the _relative_
position, velocity, and position of the goal of the entity at node $j$ with
respect to agent $i$, respectively. If node $j$ corresponds to a
(static/dynamic) obstacle or a goal, we set $p^{\mathrm{goal},j}_{i}\equiv
p^{j}_{i}$. Each edge $e_{ij}$ has an associated edge feature given by the
Euclidean distance between the entities $i$ and $j$. For processing the
entity_type categorical variable, we experimented with both using an embedding
layer word2vec and using one-hot encoding. No significant performance
advantage with one method over the other was found, so we chose to use the
embedding layer. Further analysis revealed that the learned embedding vectors
for entity_type were equidistant from each other when visualized in 2D. This
was to be expected because each of the entity types $(agent,obstacle,goal)$
are equally distinct from another. Future work could include further
refinement of entity_type, _e.g_., adversarial vs. cooperative obstacles,
static vs. dynamic obstacles, etc.

We adopt a similar reward function as used in multi-agent particle environment
(MAPE)111https://github.com/openai/multiagent-particle-envs MAPE where each
agent $i$ gets a reward:
$r_{t}^{(i)}=r_{\mathrm{dist},t}^{(i)}+r_{\mathrm{coll},t}^{(i)}+r_{\mathrm{goal},t}^{(i)}$,
where $r_{\mathrm{dist},t}^{(i)}$ is the negative of the Euclidean distance to
the goal, $r_{\mathrm{coll},t}^{(i)}=-5$ if it collides with any other entity
and zero otherwise, and $r_{\mathrm{goal},t}^{(i)}=+5$ if the agent has
reached the goal and zero otherwise. The joint reward function is defined as
$R(s_{t},A_{t})=\sum_{i=1}^{N}r_{t}^{(i)}$, which encourages cooperation among
all agents.

###  3.2 Information Aggregation

To infer information about the local neighborhood around each agent $i$, we
use a graph neural network (GNN) with a message passing framework
message_passing_GNN . Specifically, we use Unified Message Passing Model
(UniMP) graphTransformerConv , a variant of a graph transformer
AttentionIsAllYouNeed ; graphTransformer where each layer update is as
follows:

$x_{i}^{\prime}=W_{1}\cdot
x_{i}+\sum\limits_{j\in\mathcal{N}(i)}\alpha_{i,j}W_{2}\cdot x_{j},$ (1)

where $x_{k}$ are the node features in the graph, $\mathcal{N}(i)$ is the set
of nodes which are connected to node $i$, $W_{k}$ are learnable weight
matrices and the attention coefficients $\alpha_{i,j}$ are computed via multi-
head dot product attention:

$\alpha_{i,j}=\mathrm{softmax}\left(\dfrac{\left(W_{3}\cdot
x_{i}\right)^{T}\left(W_{4}\cdot x_{j}+W_{5}\cdot
e_{ij}\right)}{\sqrt{\mathrm{c}}}\right),$ (2)

where $W_{k}$ are learnable weight matrices, $e_{ij}$ are edge features for
the edge connecting nodes $i$ and $j$, and $c$ is the output dimension for
that specific layer.

The attention mechanism allows the agents to selectively prioritize messages
coming from their neighbors according to their importance. We use multiple
layers of this message-passing so that information can be propagated between
agents that are higher-order neighbors with each other. For each agent $i$,
this module aggregates information from the neighboring nodes in the graph
into a fixed-sized vector $x^{(i)}_{\mathrm{agg}}$. The concatenated vector
$[o^{(i)},x^{(i)}_{\mathrm{agg}}]$ is given as the input to the actor network.
This architecture allows InforMARL to dynamically adapt to a changing number
of entities in the environment while remaining invariant to the permutation of
the observed entities.

###  3.3 Graph Information Aggregation

While training the model in the CTDE setting, the critic generally gets the
state-action pairs of all individual agents in the environment as a
concatenated vector. To make the training transferable to a variable number of
agents and to aid with curriculum learning curriculum_learning ;
transfer_learning1 ; transfer_learning2 , we replace this concatenation with a
graph information aggregation module. This module is similar to the
‘Information Aggregation’ one in which a GNN aggregates information from the
agent’s neighbors. A global mean pooling operator is applied to aggregate the
updated node features in the graph:

$X_{\mathrm{agg}}=\frac{1}{N}\sum\limits_{i=1}^{N}x^{(i)}_{\mathrm{agg}}$ (3)

Note that $X_{\mathrm{agg}}$ is a vector of fixed size independent of the
number of agents, which is not the case when concatenating the state action-
pairs of all individual agents. This vector is then given as input to the
critic network.

###  3.4 Actor-Critic Networks

The actor and critic networks can be either a multi-layer perceptron (MLP) or
a recurrent neural network (RNN) rnnForPOMDP , using either LSTMs LSTM or GRUs
GRU . Our proposed information aggregation method can be used in conjunction
with any standard MARL algorithm (_e.g_., MADDPG MADDPG , MATD3 MATD3 , MAPPO
MAPPO , QMIX QMIX , VDN VDN , etc.).

##  4 Experiments

Task description: We evaluate our proposed model on the navigation task by
modifying the MAPE MAPE . In this environment, agents move around in a 2D
space following a double integrator dynamics model double_integrator_model .
The action space for each agent is discretized such that it can control unit
acceleration and deceleration in the $x$\- and $y$\- directions. The
navigation task differs from the target-coverage targetCoverage1 ;
targetCoverage2 problem in that each agent navigates to its own specific goal
rather than any of the available goals. Agents start at random locations at
the beginning of each episode; the corresponding goals are also randomly
distributed. Static obstacles are placed randomly in the environment in each
episode. The overarching objective is for every agent to reach its
corresponding goal without colliding with any other entity.

Implementation specifications: We chose to use MAPPO MAPPO as the base MARL
algorithm for InforMARL because it was found to be the best performing of the
standard MARL baselines in the 3 agent-3 obstacle scenario. We implemented
InforMARL by modifying the official codebase for MAPPO MAPPO in PyTorch. The
code for our experiments can be found at https://github.com/nsidn98/InforMARL.
The codebase links to our baseline implementations can be found in Appendix A.
We used the official implementations for most of the baselines considered in
Section 4.2. We prefix the algorithm name with ‘R’ to denote the recurrent
neural network (RNN) version of the algorithm (_e.g_. MAPPO-RMAPPO, MADDPG-
RMADDPG, etc.). We use the same hyperparameters as used in the experiments for
MAPPO, and do not perform parameter tuning on InforMARL for any of the MAPPO-
based parameters. The hyperparameters used in our experiments can be found in
Appendix B. We train the models with five different random seeds, and compute
the means and standard deviations of the rewards.

Amount of information available to agents: The amount of information available
to each agent determines whether or not the agent can learn a meaningful
policy. Although having more information is generally correlated to better
performance, it does not necessarily scale well with the number of agents.
Prior works MAPPO ; MADDPG have typically used a naïve concatenation of the
states of all agents or entities in the environment fed into a neural network.
Such an approach scales poorly (the network input size is determined by the
number of agents) and does not transfer well to scenarios with a different
number of agents than the training environment. We illustrate the dependence
of the learned policies on the amount of information available to agents by
defining three information modes:

  * •

_Local:_ In the local information mode,
$o^{(i)}_{\mathrm{loc}}=[p^{(i)},v^{(i)},p^{(i)}_{\mathrm{goal}}]$ where
$p^{(i)}$ and $v^{(i)}$ are the position and velocity of agent $i$ in a global
frame, and $p^{(i)}_{\mathrm{goal}}$ is the position of the goal relative to
the agent’s position. InforMARL uses this information mode.

  * •

_Global:_ Here,
$o_{\mathrm{glob}}^{(i)}=[p^{(i)},v^{(i)},p^{(i)}_{\mathrm{goal}},p^{(i)}_{\mathrm{other}}]$,
where $p^{(i)}_{\mathrm{other}}$ comprises of the relative positions of all
the other entities in the environment. The scenarios defined in the MAPE (and
consequently, other approaches that use MAPE) use this type of information
mode unless explicitly stated otherwise.

  * •

_Neighborhood:_ In this information mode, agent $i$ observes
$o_{\mathrm{nbd}}^{(i)}=[p^{(i)},v^{(i)},p^{(i)}_{\mathrm{goal}},p^{(i)}_{\mathrm{other}}]$,
where $p^{(i)}_{\mathrm{other}}$ comprises of the relative positions of all
other entities which are within a distance nbd-dist of the agent. The maximum
number of entities within the neighborhood is denoted max-nbd-entities, and so
the dimension of the observation vector is fixed. If there are fewer than max-
nbd-entities within a distance nbd-dist of the agent, we pad this vector with
zeros.

###  4.1 A Motivating Experiment

Using _local_ or _neighborhood_ information modes is preferable, since they
are transferable to scenarios with a different number of entities in the
environment. By contrast, the _global_ information mode is not transferable to
other scenarios. Figure 2 shows the rewards obtained during training using the
three information modes. We use RMAPPO as the base MARL algorithm for the 3
agent-3 obstacle scenario. We chose RMAPPO because it was found to be the
best-performing MARL algorithm in the global information mode.

Figure 2: RMAPPO with the local, neighborhood (with 1, 3, and 5 max-nbd-
entities), and global information given as states. The plots show the rewards
during training for the 3 agent-3 obstacle navigation scenario. The means and
standard deviations of the rewards over training with five random seeds are
shown. Comparing the global information mode to the others, we see that merely
providing local information and a naïve concatenation of neighborhood
information is not sufficient to learn an optimal policy.

We see in Figure 2 that the policy learned with global information is better
than the policies learned with local or neighborhood information. This is
because it has the information necessary to take optimal actions. In the
‘nbd_5’ scenario, if all entities are within the nbd_dist ball, then
$o^{(i)}_{\mathrm{nbd}}\equiv o^{(i)}_{\mathrm{glob}}$, since there are only
five entities in the environment apart from the agent itself. Although the
‘nbd_5’ scenario is almost similar to the global information mode, the
performances are not similar: the global information mode achieves much higher
rewards. This behavior is because the observation $o_{i,\mathrm{nbd}}$ can
change temporally from having information about an entity when it lies within
a nbd_dist ball of the agent, and then getting padded with zeros when it is
out of the ball. This inconsistency in the amount of information accessible to
the agent at every time step causes a significant performance difference
between the policies learned with global and neighborhood information modes.
This experiment motivates us to find a way in which more information can be
leveraged (similar to the global case), but in a manner that does not suffer
from the performance shortcomings of the neighborhood or local information
modes.

###  4.2 Comparing InforMARL to Other Baselines

Figure 3: Comparison of the training performance of InforMARL with the best-
performing baselines using global and local information. The means and
standard deviations of the rewards over training with five random seeds are
shown. InforMARL significantly outperforms most baseline algorithms. Although
RMAPPO has similar performance, it requires global information. Appendix D
presents a complete comparison to more baselines.

As shown in Section 4.1, using local information modes or naïvely
concatenating neighborhood entity information is not sufficient to learn
optimal policies. In this section, we demonstrate that InforMARL can
effectively learn policies for navigation given local information. We then
compare its performance in this problem setting with prior deep MARL
approaches. Specifically, we consider the following classes of methods as
baselines for comparison.

  1. 1.

Graph Policy Gradient (GPG) GPG : GPG uses a graph convolutional neural
network (GCN) GCN to parameterize policies for agents. The authors use the
policy gradient method policyGradient as the base MARL algorithm. We perform
experiments with both dynamic and static graphs. Note: It was shown in GPG
that using a static graph constructed at the beginning of the episode was
better than using a dynamic graph.

  2. 2.

Graph Convolutional Reinforcement Learning (DGN, DGN+ATOC) ATOC ; DGN :
Similar to GPG, these methods use GCNs to capture interactions between the
agents in the environment. A key difference between InforMARL and these two
methods (GPG and DGN) is that while the latter two approaches consider only
agents in their message-passing graph, InforMARL also includes other (non-
agent) entities in the graph.

  3. 3.

Entity Message Passing (EMP) EMP : Similar to InforMARL, EMP uses an agent-
entity graph. However, in contrast to InforMARL, EMP assumes that agents know
the positions of all entities in the graphs (i.e., global information) at the
beginning of the episode.

  4. 4.

Other standard MARL Algorithms: Finally, we compare InforMARL with standard
MARL algorithms, namely, MADDPG MADDPG , MATD3 MATD3 , QMIX QMIX , VDN VDN and
MAPPO MAPPO . In each case, we also consider the recurrent neural network
versions. We focus on results for the global information modes, since our
experiments found that these methods did not learn well with local
information.

Figure 3 shows the training performance of InforMARL and the best-performing
of the baselines mentioned above. For ease of visualization, we plot only the
four best-performing methods with global and local information, respectively.
A complete comparison plot showing all the baseline methods can be found in
Appendix D. Each line corresponds to the mean and standard deviation over five
random seeds. For each method, we consider scenarios with $N=3$, $N=7$, and
$N=10$ agents in the environment.

Figure 3 illustrates that only using RMAPPO (i.e. the RNN version of MAPPO) or
InforMARL are agents able to learn to navigate and get to their goals.
However, unlike RMAPPO which requires global information, InforMARL is able to
achieve this with just local information. Furthermore, we find that InforMARL
requires a similar number of training steps as RMAPPO, despite having access
to much less information.

Algorithm | Inform- | $N=3$ | $N=7$ | $N=10$
---|---|---|---|---
ation | Reward | $T$ | # col | $S\%$ | Reward | $T$ | # col | $S\%$ | Reward | $T$ | # col | $S\%$
mode
MADDPG MADDPG | Global | $-100.73$ | $0.97$ | $1.60$ | 5 | $-206.07$ | $0.98$ | $6.01$ | 0 | -210.43 | 1.00 | 9.26 | 0
RMADDPG | Global | $10.73$ | $0.75$ | $3.20$ | 23 | $-122.069$ | $0.95$ | $9.75$ | 3 | -127.98 | 1.00 | 14.89 | 0
MATD3 MATD3 | Global | $-90.31$ | $0.98$ | $1.11$ | 5 | $-169.08$ | $0.99$ | $1.98$ | 1 | -173.20 | 1.00 | 3.50 | 0
RMATD3 | Global | $105.49$ | $0.51$ | $3.07$ | 67 | $-128.93$ | $0.96$ | $7.40$ | 6 | -131.72 | 0.99 | 11.14 | 1
QMIX QMIX | Global | $-54.24$ | $0.84$ | $2.06$ | 7 | $-288.81$ | $0.97$ | $17.05$ | 0 | -273.46 | 1.00 | 25.31 | 0
RQMIX | Global | $19.21$ | $0.77$ | $1.42$ | 28 | $-83.41$ | $0.85$ | $11.5$ | 12 | -76.98 | 0.96 | 17.04 | 2
VDN VDN | Global | $18.86$ | $0.67$ | $2.29$ | 27 | $39.87$ | $0.64$ | $8.58$ | 23 | 43.23 | 0.73 | 12.54 | 19
RVDN | Global | $64.04$ | $0.62$ | $1.05$ | 45 | $140.94$ | $0.62$ | $6.60$ | 47 | 157.63 | 0.64 | 10.00 | 43
RMAPPO MAPPO | Global | 173.13 | 0.41 | 1.47 | 96 | 327.39 | 0.44 | 8.64 | 88 | 366.81 | 0.44 | 13.21 | 79
GPG (static) GPG | Local | $-67.03$ | $0.96$ | $1.47$ | 7 | $-180.14$ | $0.99$ | $6.41$ | 1 | -182.57 | 1.00 | 9.44 | 0
GPG (dynamic) | Local | $-46.27$ | $0.87$ | $0.43$ | 8 | $-165.91$ | $1.00$ | $3.20$ | 3 | -173.53 | 1.00 | 4.68 | 0
DGN DGN | Local | $32.94$ | $0.59$ | $2.57$ | 32 | $-232.32$ | $0.97$ | $5.39$ | 0 | -243.45 | 1.00 | 8.62 | 0
DGN + ATOC | Local | $67.70$ | $0.66$ | $1.49$ | 35 | $-189.61$ | $0.97$ | $2.20$ | 0 | -201.01 | 1.00 | 4.06 | 0
EMP EMP | Local | $-83.96$ | $0.98$ | $1.35$ | 6 | $-211.90$ | $0.98$ | $5.91$ | 0 | -209.90 | 1.00 | 9.22 | 0
InforMARL | Local | 205.24 | 0.38 | 1.45 | 100 | 399.01 | 0.37 | 6.97 | 100 | 429.14 | 0.39 | 10.50 | 100
Table 1: Comparison of InforMARL with other baseline methods, for scenarios
with 3, 7, and 10 agents in the environment. The results presented represent
the average of 100 test episodes. The following metrics are compared: (a)
Total reward obtained in an episode by all the agents (higher is better). (b)
Fraction of episode taken by the agents to reach the goal, $T$ (lower is
better). (c) The total number of collisions the agents had in the episode, #
col (lower is better). (d) Percent of episodes in which all agents are able to
get to their goals, $S$% (higher is better). The best-performing methods that
use global information (RMAPPO) and local information (InforMARL) are
highlighted. As noted in Section 4.2, the metrics # col and $S$ should be
considered on balance.

We present the following metrics in Table 1. The results represent an average
over 100 test episodes.

  1. 1.

The total rewards obtained by the agents during an episode using the reward
function defined in Section 3.1. A higher value corresponds to better
performance.

  2. 2.

The fraction of an episode that the agents take on average to get to the goal,
denoted $T$. If an agent does not reach its goal, then $T$ is set to be 1
(lower is better).

  3. 3.

Percent of episodes in which all agents are able to get to their goals,
denoted $S$% (higher is better).

  4. 4.

The total number of collisions (both agent-agent + agent-obstacle) that agents
had in an episode, denoted # col. The lower this metric, the better the
performance of the algorithm.

Although having a smaller number of collisions is better, the policies of some
of the baseline algorithms do not significantly move the agents from their
initial position after training and hence do not get to the goal. This leads
to them having a lower number of collisions. Hence, this metric should be
judged with the success rate in context.

The graph-based methods, namely GPG (static and dynamic), DGN (+ATOC), and EMP
do not learn effectively with local information modes. Although both GPG and
DGN use GCNs, they do not perform as well as InforMARL because they use only
agent-agent and not agent-entity graphs. The lack of information about non-
agent entities means that agents cannot maneuver through the environment to
avoid collisions.

In contrast to the results showed in GPG , our results show that the dynamic
graph version of GPG is slightly better than the static graph one. A possible
reason for this discrepancy is that the fixed formation environment considered
in GPG is more amenable to the use of static graphs than the navigation
environment. EMP, which uses a similar agent-entity graph as ours, fails to
learn because of the strong assumption of having access to the positions of
all entities in the environment at the beginning of the episode. The original
implementation of EMP (and the associated environments) included information
about all the entities other than agents in the observation vector.

###  4.3 Scalability of InforMARL

To evaluate the scalability of InforMARL with minimum information, we perform
experiments by testing the models in scenarios with a different number of
agents from those that they were trained on.

Test Train | $n=3$ | $n=7$ | $n=10$
---|---|---|---
$m=3$ | Reward/$m$ | 68.41 | 58.30 | 57.27
$T$ | 0.38 | 0.38 | 0.39
(# col)/m | 0.48 | 0.61 | 0.64
$S\%$ | 100 | 100 | 99
$m=7$ | Reward/$m$ | 61.16 | 57.00 | 58.32
$T$ | 0.38 | 0.37 | 0.38
(# col)/$m$ | 1.15 | 0.99 | 1.03
$S\%$ | 100 | 100 | 100
$m=10$ | Reward/$m$ | 58.59 | 53.25 | 52.10
$T$ | 0.38 | 0.38 | 0.39
(# col)/$m$ | 1.60 | 1.43 | 1.68
$S\%$ | 100 | 99 | 100
$m=15$ | Reward/$m$ | 53.19 | 46.39 | 48.15
$T$ | 0.39 | 0.39 | 0.40
(# col)/$m$ | 2.24 | 2.314 | 2.20
$S\%$ | 100 | 99 | 99
Table 2: Test performance of InforMARL, when trained on scenarios with $n$
agents and tested on scenarios with $m$ agents in the environment. InforMARL
is able to achieve a success rate of nearly 100% while the fraction of
episodes taken to reach the goals ($T$) is almost constant. Note that we
normalize the reward and the number of collisions by $m$ for ease of
comparison.

Table 2 shows the results of testing InforMARL trained on $n$ agents and
tested on $m$ agents. Each scenario is tested over 100 episodes. The number of
obstacles in the environment is randomly chosen from $(0,10)$ at the beginning
of the episode. Based on the findings presented in Table 1) , we did not test
the scalability of other methods since they did not perform well in the local
information mode, or could not handle varying numbers of entities in the
environment.

We see from Table 2 that InforMARL is able to achieve a success rate of 100%
for almost all the scenarios. Furthermore, with InforMARL, the agents are able
to get to their goals within $T\sim 0.39$ of the episode length in all
scenarios. The number of collisions per agent increases, which is to be
expected as the environment becomes denser. We believe that this can be
remedied by using a stricter penalty (negative reward) for collisions instead
of the $-5$ penalty used in our experiments. Alternatively, a more
sophisticated approach with control barrier functions for satisfying safety
constraints could be used to provide a formal safety guarantee in the MARL
setting MARL_control_barrier .

###  4.4 Effect of Sensing Radius

We investigate how the performance of InforMARL depends on the sensing radius,
namely, how much information (i.e. over what neighborhood of the ego-agent) an
agent has access to. As the radius increases to a large value, the graph
becomes fully connected; as the radius decreases to zero, the neighborhood
information mode converges to the local information mode. As seen in Figure 4,
when learning with a small sensing radius (_e.g_. $\rho=\\{0.1,0.2\\}$), the
agents are not able to achieve the same reward as can be achieved with a
larger sensing radius ($\rho=\\{0.5,1,2,5\\}$)222Measurements for $\rho$, a
distance, are in meters. We also note that there are diminishing returns when
increasing the sensing radius from $\rho=0.5$ to $\rho=5$. Since the agents
far away from the ego-agent have little influence on its immediate decisions,
the extra information obtained by increasing the sensing radius does not
improve performance very much.

Figure 4: Diminishing returns in performance gains from increasing the sensing
radius for InforMARL. The dashed lines are the reward values after saturation
for RMAPPO in the global (in green) and local (in red) information modes,
respectively. They are provided for reference.

Figure 5: Training performance of InforMARL with, and without, the graph
information aggregation module, for a 3-agent scenario. The two variants have
similar sample complexities. However, the critic network with the graph
information aggregation module has fewer parameters than the one without this
module.

###  4.5 Ablation Study for Graph Information Aggregation Module

The graph information aggregation module allows our method to perform transfer
learning to more complex environments. In this section, we compare models with
and without the graph information aggregation module. In the absence of the
graph information aggregation module, the states of individual agents are
concatenated to be given as input to the centralized critic. Figure 5 presents
the results of this ablation study. We find that both models have similar
sample complexities and performance. However, the number of parameters for the
critic network with the graph information aggregation module is much smaller
and independent of the number of entities in the environment.

##  5 Conclusions and Future Work

We introduced InforMARL, a novel architecture that uses GNNs for scalable
multi-agent reinforcement learning. We showed that having just local
observations as states is not enough for standard MARL algorithms to learn
meaningful policies. Along with this, we also showed that albeit naïvely
concatenating state information about all the entities in the environment
helps to learn good policies, they are not transferable to other scenarios
with a different number of entities than what it was trained on. InforMARL is
able to learn transferable policies using standard MARL algorithms using just
local observations and an aggregated neighborhood information vector.
Furthermore, it has better sample complexity than other standard MARL
algorithms that use global observation. Future work will include the
introduction of more complex (potentially adversarial) dynamic obstacles in
the environment and adding a safety guarantee layer for the actions of the
agents to avoid collisions. Additionally, the use of InforMARL for curriculum
learning and transfer learning to different environments is a topic of ongoing
research.

## Acknowledgements

The authors would like to thank the MIT SuperCloud supercloud and the Lincoln
Laboratory Supercomputing Center for providing high performance computing
resources that have contributed to the research results reported within this
paper. The NASA University Leadership initiative (grant #80NSSC20M0163)
provided funds to assist the authors with their research, but this article
solely reflects the opinions and conclusions of its authors and not any NASA
entity. This research was sponsored in part by the United States AFRL and the
United States Air Force Artificial Intelligence Accelerator and was
accomplished under Cooperative Agreement Number FA8750-19-2-1000. The views
and conclusions contained in this document are those of the authors and should
not be interpreted as representing the official policies, either expressed or
implied, of the United States Air Force or the U.S. Government. The U.S.
Government is authorized to reproduce and distribute reprints for Government
purposes notwithstanding any copyright notion herein.

## References

  * [1] Timothy P. Lillicrap, Jonathan J. Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra.  Continuous control with deep reinforcement learning.  In Yoshua Bengio and Yann LeCun, editors, 4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings, 2016.
  * [2] Sergey Levine, Chelsea Finn, Trevor Darrell, and Pieter Abbeel.  End-to-end training of deep visuomotor policies.  CoRR, abs/1504.00702, 2015.
  * [3] Nicolas Heess, Greg Wayne, David Silver, Timothy P. Lillicrap, Yuval Tassa, and Tom Erez.  Learning continuous control policies by stochastic value gradients.  CoRR, abs/1510.09142, 2015.
  * [4] David Silver, Aja Huang, Chris J. Maddison, Arthur Guez, Laurent Sifre, George van den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, and et al.  Mastering the game of go with deep neural networks and tree search.  Nature, 529(7587):484–489, 2016.
  * [5] David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, Timothy P. Lillicrap, Karen Simonyan, and Demis Hassabis.  Mastering chess and shogi by self-play with a general reinforcement learning algorithm.  CoRR, abs/1712.01815, 2017.
  * [6] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin A. Riedmiller.  Playing atari with deep reinforcement learning.  CoRR, abs/1312.5602, 2013.
  * [7] V. Mnih, K. Kavukcuoglu, D. Silver, Andrei A. Rusu, J. Veness, Marc G. Bellemare, A. Graves, Martin A. Riedmiller, A. Fidjeland, Georg Ostrovski, Stig Petersen, Charlie Beattie, A. Sadik, Ioannis Antonoglou, Helen King, D. Kumaran, Daan Wierstra, S. Legg, and D. Hassabis.  Human-level control through deep reinforcement learning.  Nature, 518:529–533, 2015.
  * [8] Sven Gronauer and Klaus Diepold.  Multi-agent deep reinforcement learning: a survey.  Artificial Intelligence Review, 55(2):895–943, 2022.
  * [9] Ming Tan.  Multi-agent reinforcement learning: Independent vs. cooperative agents.  In Proceedings of the tenth international conference on machine learning, pages 330–337, 1993.
  * [10] Ming Tan.  Multi-agent reinforcement learning: Independent vs. cooperative agents.  In In Proceedings of the Tenth International Conference on Machine Learning, pages 330–337. Morgan Kaufmann, 1993.
  * [11] Ardi Tampuu, Tambet Matiisen, Dorian Kodelja, Ilya Kuzovkin, Kristjan Korjus, Juhan Aru, Jaan Aru, and Raul Vicente.  Multiagent cooperation and competition with deep reinforcement learning.  CoRR, abs/1511.08779, 2015.
  * [12] Vijay Kumar, D. Rus, and Sanjiv Singh.  Robot and sensor networks for first responders.  IEEE Pervasive Computing, 3(4):24–33, 2004.
  * [13] Matthew Dunbabin and Lino Marques.  Robots for environmental monitoring: Significant advancements and applications.  IEEE Robotics & Automation Magazine, 19(1):24–39, 2012.
  * [14] Hanlin Zhang, Sixiao Wei, Wei Yu, Erik Blasch, Genshe Chen, Dan Shen, and Khanh Pham.  Scheduling methods for unmanned aerial vehicle based delivery systems.  In 2014 IEEE/AIAA 33rd Digital Avionics Systems Conference (DASC), pages 6C1–1–6C1–9, 2014.
  * [15] Kevin Dorling, Jordan Heinrichs, Geoffrey G. Messier, and Sebastian Magierowski.  Vehicle routing problems for drone delivery.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 47(1):70–85, 2017.
  * [16] Michelle Sing Yee Hii, Patrick Courtney, and Paul G. Royall.  An evaluation of the delivery of medicines using drones.  Drones, 3(3), 2019.
  * [17] Oriol Vinyals, Igor Babuschkin, Wojciech M. Czarnecki, Michaël Mathieu, Andrew Dudzik, Junyoung Chung, David H. Choi, Richard Powell, Timo Ewalds, Petko Georgiev, Junhyuk Oh, Dan Horgan, Manuel Kroiss, Ivo Danihelka, Aja Huang, L. Sifre, Trevor Cai, John P. Agapiou, Max Jaderberg, Alexander Sasha Vezhnevets, Rémi Leblond, Tobias Pohlen, Valentin Dalibard, David Budden, Yury Sulsky, James Molloy, Tom Le Paine, Caglar Gulcehre, Ziyun Wang, Tobias Pfaff, Yuhuai Wu, Roman Ring, Dani Yogatama, Dario Wünsch, Katrina McKinney, Oliver Smith, Tom Schaul, Timothy P. Lillicrap, Koray Kavukcuoglu, Demis Hassabis, Chris Apps, and David Silver.  Grandmaster level in starcraft ii using multi-agent reinforcement learning.  Nature, pages 1–5, 2019.
  * [18] Christopher Berner, Greg Brockman, Brooke Chan, Vicki Cheung, Przemyslaw Debiak, Christy Dennison, David Farhi, Quirin Fischer, Shariq Hashme, Christopher Hesse, Rafal Józefowicz, Scott Gray, Catherine Olsson, Jakub Pachocki, Michael Petrov, Henrique Pondé de Oliveira Pinto, Jonathan Raiman, Tim Salimans, Jeremy Schlatter, Jonas Schneider, Szymon Sidor, Ilya Sutskever, Jie Tang, Filip Wolski, and Susan Zhang.  Dota 2 with large scale deep reinforcement learning.  CoRR, abs/1912.06680, 2019.
  * [19] Chao Yu, Akash Velu, Eugene Vinitsky, Jiaxuan Gao, Yu Wang, Alexandre Bayen, and Yi Wu.  The surprising effectiveness of PPO in cooperative multi-agent games.  In Thirty-sixth Conference on Neural Information Processing Systems Datasets and Benchmarks Track, 2022.
  * [20] James Rendleman and Sarah Mountin.  Responsible SSA Cooperation to Mitigate On-orbit Space Debris Risks, page 851–856.  International Conference on Recent Advances in Space Technologies, 2015\.
  * [21] Charity Weeden, Jonathan Goff, Jessica Noble, Jeremy Schiel, Dana Turse, Colin Doughan, John Carrico, and Rob Patterson.  Comments of global newspace operators, 2019.
  * [22] Vijay Konda and John Tsitsiklis.  Actor-critic algorithms.  In S. Solla, T. Leen, and K. Müller, editors, Advances in Neural Information Processing Systems, volume 12. MIT Press, 1999.
  * [23] Ryan Lowe, Yi Wu, Aviv Tamar, Jean Harb, Pieter Abbeel, and Igor Mordatch.  Multi-agent actor-critic for mixed cooperative-competitive environments.  CoRR, abs/1706.02275, 2017.
  * [24] Johannes Ackermann, Volker Gabler, Takayuki Osa, and Masashi Sugiyama.  Reducing overestimation bias in multi-agent domains using double centralized critics.  CoRR, abs/1910.01465, 2019.
  * [25] Peter Sunehag, Guy Lever, Audrunas Gruslys, Wojciech Marian Czarnecki, Vinícius Flores Zambaldi, Max Jaderberg, Marc Lanctot, Nicolas Sonnerat, Joel Z. Leibo, Karl Tuyls, and Thore Graepel.  Value-decomposition networks for cooperative multi-agent learning.  CoRR, abs/1706.05296, 2017.
  * [26] Tabish Rashid, Mikayel Samvelyan, Christian Schröder de Witt, Gregory Farquhar, Jakob N. Foerster, and Shimon Whiteson.  QMIX: monotonic value function factorisation for deep multi-agent reinforcement learning.  CoRR, abs/1803.11485, 2018.
  * [27] Lifeng Zhou, Vishnu D. Sharma, Qingbiao Li, Amanda Prorok, Alejandro Ribeiro, and Vijay Kumar.  Graph neural networks for decentralized multi-robot submodular action selection.  CoRR, abs/2105.08601, 2021.
  * [28] Pratap Tokekar, Volkan Isler, and Antonio Franchi.  Multi-target visual tracking with aerial robots.  In 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 3067–3072, 2014.
  * [29] Philip Dames, Pratap Tokekar, and Vijay Kumar.  Detecting, localizing, and tracking an unknown number of moving targets using a team of mobile robots.  The International Journal of Robotics Research, 36(13-14):1540–1553, 2017.
  * [30] Arbaaz Khan, Ekaterina Tolstaya, Alejandro Ribeiro, and Vijay Kumar.  Graph policy gradients for large scale robot control, 2019.
  * [31] Richard S Sutton, David McAllester, Satinder Singh, and Yishay Mansour.  Policy gradient methods for reinforcement learning with function approximation.  In S. Solla, T. Leen, and K. Müller, editors, Advances in Neural Information Processing Systems, volume 12. MIT Press, 1999.
  * [32] M. Turpin, Nathan Michael, and V. Kumar.  Trajectory design and control for aggressive formation flight with quadrotors.  Autonomous Robots, 33(1):143 – 156, August 2012.
  * [33] Matthew Turpin, Nathan Michael, and Vijay Kumar.  Decentralized formation control with variable shapes for aerial robots.  In 2012 IEEE International Conference on Robotics and Automation, pages 23–30, 2012.
  * [34] Jiechuan Jiang, Chen Dun, Tiejun Huang, and Zongqing Lu.  Graph convolutional reinforcement learning.  In ICLR, 2020.
  * [35] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  CoRR, abs/1609.02907, 2016.
  * [36] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [37] Akshat Agarwal, Sumit Kumar, and Katia P. Sycara.  Learning transferable cooperative behavior in multi-agent teams.  CoRR, abs/1906.01202, 2019.
  * [38] Sainbayar Sukhbaatar, Arthur Szlam, and Rob Fergus.  Learning multiagent communication with backpropagation.  CoRR, abs/1605.07736, 2016.
  * [39] Yedid Hoshen.  Vain: Attentional multi-agent predictive modeling.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [40] Jiechuan Jiang and Zongqing Lu.  Learning attentional communication for multi-agent cooperation.  CoRR, abs/1805.07733, 2018.
  * [41] Abhishek Das, Théophile Gervet, Joshua Romoff, Dhruv Batra, Devi Parikh, Michael G. Rabbat, and Joelle Pineau.  Tarmac: Targeted multi-agent communication.  CoRR, abs/1810.11187, 2018.
  * [42] Won Joon Yun, Byungju Lim, Soyi Jung, Young-Chai Ko, Jihong Park, Joongheon Kim, and Mehdi Bennis.  Attention-based reinforcement learning for real-time UAV semantic communication.  CoRR, abs/2105.10716, 2021.
  * [43] Martin L. Puterman.  Markov Decision Processes.  Wiley, 1994.
  * [44] Leslie Pack Kaelbling, Michael L. Littman, and Anthony R. Cassandra.  Planning and acting in partially observable stochastic domains.  Artificial Intelligence, 101(1):99–134, 1998.
  * [45] Frans A. Oliehoek and Christopher Amato.  A Concise Introduction to Decentralized POMDPs.  Springer Publishing Company, Incorporated, 1st edition, 2016.
  * [46] Tomás Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeffrey Dean.  Distributed representations of words and phrases and their compositionality.  CoRR, abs/1310.4546, 2013.
  * [47] Ryan Lowe, Yi Wu, Aviv Tamar, Jean Harb, Pieter Abbeel, and Igor Mordatch.  Multi-agent actor-critic for mixed cooperative-competitive environments.  Neural Information Processing Systems (NIPS), 2017.
  * [48] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.  Neural message passing for quantum chemistry.  CoRR, abs/1704.01212, 2017.
  * [49] Yunsheng Shi, Zhengjie Huang, Wenjin Wang, Hui Zhong, Shikun Feng, and Yu Sun.  Masked label prediction: Unified massage passing model for semi-supervised classification.  CoRR, abs/2009.03509, 2020.
  * [50] Vijay Prakash Dwivedi and Xavier Bresson.  A generalization of transformer networks to graphs.  CoRR, abs/2012.09699, 2020.
  * [51] Sanmit Narvekar, Bei Peng, Matteo Leonetti, Jivko Sinapov, Matthew E. Taylor, and Peter Stone.  Curriculum learning for reinforcement learning domains: A framework and survey.  CoRR, abs/2003.04960, 2020.
  * [52] Matthew E. Taylor and Peter Stone.  Transfer learning for reinforcement learning domains: A survey.  Journal of Machine Learning Research, 10(56):1633–1685, 2009.
  * [53] Alessandro Lazaric, Marcello Restelli, and Andrea Bonarini.  Transfer of samples in batch reinforcement learning.  pages 544–551, 01 2008.
  * [54] Matthew J. Hausknecht and Peter Stone.  Deep recurrent q-learning for partially observable mdps.  CoRR, abs/1507.06527, 2015.
  * [55] Sepp Hochreiter and Jürgen Schmidhuber.  Long Short-Term Memory.  Neural Computation, 9(8):1735–1780, 11 1997.
  * [56] KyungHyun Cho, Bart van Merrienboer, Dzmitry Bahdanau, and Yoshua Bengio.  On the properties of neural machine translation: Encoder-decoder approaches.  CoRR, abs/1409.1259, 2014.
  * [57] V.G. Rao and D.S. Bernstein.  Naive control of the double integrator.  IEEE Control Systems Magazine, 21(5):86–97, 2001.
  * [58] Zengyi Qin, Kaiqing Zhang, Yuxiao Chen, Jingkai Chen, and Chuchu Fan.  Learning safe multi-agent control with decentralized neural barrier certificates.  CoRR, abs/2101.05436, 2021.
  * [59] Albert Reuther, Jeremy Kepner, Chansup Byun, Siddharth Samsi, William Arcand, David Bestor, Bill Bergeron, Vijay Gadepally, Michael Houle, Matthew Hubbell, Michael Jones, Anna Klein, Lauren Milechin, Julia Mullen, Andrew Prout, Antonio Rosa, Charles Yee, and Peter Michaleas.  Interactive supercomputing on 40,000 cores for machine learning and data analysis.  In 2018 IEEE High Performance extreme Computing Conference (HPEC), pages 1–6. IEEE, 2018.
  * [60] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Z. Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  CoRR, abs/1912.01703, 2019.
  * [61] Matthias Fey and Jan E. Lenssen.  Fast graph representation learning with PyTorch Geometric.  In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.

##  Appendix A Baseline Implementation Sources

We modified the codebases from the official implementations for the GPG, DGN,
EMP, and MAPPO baselines. We also adapted the codebase for MADDPG, MATD3, VDN
and QMIX for our experiments which was benchmarked on other standard
environments. The links for all these implementations are listed below.

  * •

GPG: https://github.com/arbaazkhan2/gpg_labeled

  * •

DGN: https://github.com/jiechuanjiang/pytorch_DGN

  * •

EMP: https://github.com/sumitsk/marl_transfer

  * •

MAPPO: https://github.com/marlbenchmark/on-policy

  * •

MADDPG, MATD3, QMIX, VDN: https://github.com/marlbenchmark/off-policy

##  Appendix B Hyperparameters

We performed a hyperparameter search for these algorithms by varying the
learning-rates, network size and a few algorithm specific parameters. We
observed that the hyperparameters used in the original implementation gave the
best performance, and used those same values for our experiments as well.

Table 3 lists the hyperparameters specific to the InforMARL implementation.
Here, “entity embedding layer dim" and “entity hidden dim" refer to the
embedding layer input and output dimension respectively which is used to
process the entity-type variable in the graph. “add self loop" refers to
whether a self-loop should be added while constructing the agent-entity graph.
“gnn layer hidden dim" is the output dimension of each layer in the graph
transformer. “num gnn heads" and “num gnn layers" are the number of heads in
the attention layer and number of attention layers used in the graph
transformer. “gnn activation" is the activation function used after each layer
in the GNN module.

Hyperparameters | Value
---|---
entity embedding layer dim | 3
entity hidden dim | 16
num embedding layer | 1
add self loop | False
gnn layer hidden dim | 16
num gnn heads | 3
num gnn layers | 2
gnn activation | ReLU
Table 3: Hyperparameters used in InforMARL

Tables 4, 5, 6 list the hyperparameters common for the InforMARL, MAPPO,
MADDPG, MATD3, QMIX and VDN implementations. For MADDPG, MATD3, MAPPO, QMIX,
VDN and InforMARL, “batch size” refers to the number of environment steps
collected before updating the policy via gradient descent. “mini batch” refers
to the number of mini-batches a batch of data is split into. “gain” refers to
the weight initialization gain of the last network layer for the actor
network. “num envs" refers to the number of parallel roll out threads used to
collect state-transition tuples.

Common Hyperparameters | Value
---|---
recurrent data chunk length | 10
gradient clip norm | 10.0
gae lambda | 0.95
gamma | 0.99
value loss | huber loss
huber delta | 10.0
batch size | num envs $\times$ buffer length $\times$ num agents
mini batch size | batch size / mini-batch
gain | 0.01
network initialization | Orthogonal
optimizer | Adam
optimizer epsilon | 1e-5
weight decay | 0
use reward normalization | True
use feature normalization | True
Table 4: Common Hyperparameters used in MAPPO and InforMARL Common Hyperparameters | Value
---|---
gradient clip norm | 10.0
random episodes | 5
epsilon | 1.0 $\rightarrow$ 0.05
epsilon anneal time | 50000 timesteps
train interval | 1 episode
gamma | 0.99
critic loss | mse loss
buffer size | 5000 episodes
batch size | 32 episodes
optimizer | Adam
optimizer eps | 1e-5
weight decay | 0
network initialization | Orthogonal
use reward normalization | True
use feature normalization | True
Table 5: Common Hyperparameters used in MADDPG, MATD3, QMIX, VDN Common Hyperparameters | Value
---|---
num envs | 128
buffer length | 25
num GRU layers | 1
RNN hidden state dim | 64
fc layer hidden dim | 64
num fc | 2
num fc after | 1
Table 6: Common Hyperparameters used in MAPPO, MADDPG, MATD3, QMIX, VDN and
InforMARL

##  Appendix C Computational Environment

Our models were trained on a server with 40 Intel Xeon Gold 6248 @ 2.50 GHz
processor cores and 2 NVIDIA Volta V100 graphics cards. Our code uses PyTorch
pytorch version 1.11.0, CUDA version 11.3, and PyTorch Geometric Library PyG
version 2.0.4.

##  Appendix D Full Comparison

We showcase the learning curves of all the baseline algorithms, including
MADDPG, RMADDPG, MATD3, QMIX, VDN, GPG (static), and DGN in Figure 6.

Figure 6: A thorough comparison of InforMARL with baselines. The algorithms
that use global information modes are represented with dashed lines and the
algorithms that use local information modes are represented with solid lines.
The means and standard deviations of the rewards over training with five
random seeds are shown.


# Broadband Emission from a Kilonova Ejecta-Pulsar Wind Nebula System: Late-
Time X-ray Afterglow Rebrightening of GRB 170817A

J. Ren1,2, Z. G. Dai3,1
1School of Astronomy and Space Science, Nanjing University, Nanjing 210023,
China
2Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University),
Ministry of Education, Nanjing 210023, China
3Department of Astronomy, School of Physical Sciences, University of Science
and Technology of China, Hefei 230026, Anhui, China
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We study the broadband radiation behavior of a kilonova ejecta-pulsar wind
nebula (PWN) system. In this model, we jointly fit the observations of AT
2017gfo in UV-optical-IR bands and the late-time X-ray afterglow of GRB
170817A. Our work shows that a PWN powered by the remnant neutron star (NS)
post GW170817 event could affect the optical transient AT 2017gfo and re-
brighten the late-time X-ray afterglow of GRB 170817A. The PWN radiation will
regulate the trend of future X-ray observations from a flattening to a steep
decline until some other sources (e.g., a kilonova afterglow) become dominant.
The restricted ranges of the central NS parameters in this work are consistent
with the previous works based on the observations of AT 2017gfo only. In
addition, the new fitting result indicates that the NS wind is highly
magnetized. We point out that the radio and X-ray emission from a kilonova
ejecta-PWN system could be an important electromagnetic feature of binary NS
mergers when a long-lived remnant NS is formed. Therefore, observations of a
kilonova ejecta-PWN system will provide important information to inferring the
nature of a merger remnant.

###### keywords:

Gravitational waves – gamma-ray burst: individual (GRB 170817A) – stars:
neutron – pulsars: general

††pubyear: 2021††pagerange: Broadband Emission from a Kilonova Ejecta-Pulsar
Wind Nebula System: Late-Time X-ray Afterglow Rebrightening of GRB 170817A–4

##  1 Introduction

Compact binary mergers are the main sources of gravitational wave (GW) events
in the frequency range of the Advanced Laser Interferometer Gravitational-wave
Observatory (LIGO) and the Advanced Virgo GW detectors. Among them, the
mergers of binary neutron star (NS) and NS-black hole (BH) draw a lot of
attention since they are also potential sources of electromagnetic radiation
(EM). The first GW signal from a NS–NS merger was detected by the advanced
LIGO and Virgo detectors on 2017 August 17 12:41:04 UT (Abbott et al., 2017b).

The nature of the merger remnant of the GW170817 event has been debated so
far. Based on the observational fact that the prompt EM signal is delayed
about 1.7 seconds after the GW signal as well as the sufficiently heavy
gravitational mass of the merger remnant(Abbott et al., 2017a; Goldstein et
al., 2017; Savchenko et al., 2017; Zhang et al., 2018), a short-lived NS
remnant seems likely (e.g., Metzger et al., 2018). Although the time-lag
between the GW and EM signals point to a short-lived NS remnant, recent works
on numerical simulation and analytical calculations have attributed this
phenomenon to the effect of jet propagation (e.g., Beniamini et al., 2020;
Hamidani et al., 2020; Hamidani & Ioka, 2021; Lyutikov, 2020; Lazzati & Perna,
2019; Lazzati et al., 2020; Ren et al., 2020; Pavan et al., 2021, and
reference therein). Therefore, the observations are not inconsistent with the
view that a long-lived NS was formed after the merger. Yu et al. (2013) first
suggested the possibility of the existence of a remnant-NS-powered kilonova,
namely a mergernova. Works have been done under this picture to comprehend the
behavior of AT 2017gfo (e.g., Matsumoto et al., 2018; Li et al., 2018; Yu et
al., 2018; Ren et al., 2019).

The observations of the GRB 170817A afterglow have lasted about four years
(e.g., Makhathini et al., 2020; Balasubramanian et al., 2021; Hajela et al.,
2021a; , 2021). The popular interpretation of the origin of GRB 170817A
afterglow nowadays is an off-axis observed structured jet (e.g., Ren et al.,
2020). Some other scenarios have been proposed, say, a refreshed jet (Lamb et
al., 2020) or an $e^{\pm}$-wind-injected jet (Li & Dai, 2021). About 3.4 years
after the trigger of the GW170817 event, the Chandra X-ray Observatory
detected an X-ray source located at the position of GRB 170817A/AT 2017gfo
(Hajela et al., 2021a, b). But the unabsorbed X-ray flux derived from the
observed data exceeds the conjectural extension based on the off-axis
structured jet model. Since the flux excess is obvious at the X-ray band but
undetected at the radio band, an additional source instead of the GRB jet
external shock to explain the data is naturally needed. At present, the
kilonova afterglow gets increasing attention. It has been supposed to be the
cause of GRB 170817A rebrightening (Hajela et al., 2021a; Nedora et al.,
2021). However, some other explanations, e.g., the $e^{\pm}$-wind-injected jet
(Geng et al., 2018; Li & Dai, 2021), the fall-back accretion of central BH
(Ishizaki et al., 2021a, b; Metzger & Fernandez, 2021), and maybe the counter-
jet (Li et al., 2019; , 2021) are possible scenarios.

In this paper, we introduce a different scenario in which a pulsar wind nebula
(PWN) is powered by the remnant NS of the GW170817 event to explain the
optical transient AT 2017gfo and the late time rebrightening X-ray afterglow
of GRB 170817A. This paper is organized as follows. In Section 2 we present
our numerical model. In Section 3 we describe the fitting process in detail.
We discuss the significance of the fitting results in Section 4. We summarize
our conclusions in Section 5. Throughout this work, we use the notation
$Q=10^{x}Q_{x}$ in the c.g.s. unit unless noted otherwise, and $D_{L}=40$ Mpc
is the luminosity distance of AT 2017gfo/GRB 170817A.

##  2 The model

Similar to the Crab supernova remnant, a PWN could be formed if a pulsar is
left at the center of a merger remnant. In the previous works of an energy-
injected kilonova, the released magnetic dipole (MD) radiation from the
central NS is considered as an injected thermalization energy source (e.g., Yu
et al., 2013, 2018; Kasen et al., 2015). A PWN is formed by the interaction
between the pulsar wind and ejecta after the merger (e.g., Kotera et al.,
2013; Murase et al., 2018). Based on those works, Ren et al. (2019) considered
a PWN embedded in the merger ejecta to explain the behavior of AT 2017gfo. Wu
et al. (2021) adopted the same picture in their study of GRB 160821B afterglow
with slightly different modeling.

The spin-down energy of the central NS first powers a PWN. Then, nonthermal
photons radiated from the PWN will cross the ejecta shell and some photons are
absorbed to thermalize the ejecta material while the other photons escaping
from the PWN could impact the observations. The observed flux thus consists of
the emission $F_{\nu}^{b}$ from the ejecta and the leaked part $F_{\nu}^{\rm
leak}$ from the PWN, i.e.,

$F_{\nu}^{\rm tot}=F_{\nu}^{b}+F_{\nu}^{\rm leak}.$ (1)

The estimates of $F_{\nu}^{b}$ and $F_{\nu}^{\rm leak}$ are given in Sections
2.1 and 2.2, respectively.

After the ejecta becomes transparent to X-ray photons, an additional
observable signal from the PWN may appear. Both the early kilonova behavior
and the late X-ray band observations may be affected. The kilonova ejecta-PWN
system was proposed in the previous work (Ren et al., 2019). We present some
key formulas and new refinements below.

###  2.1 Emission from the ejecta

The dynamics and emission of the quasi-isotropic ejecta are implemented based
on a simplified radiation transfer model given by Kasen & Bildsten (2010) and
Metzger (2019). The merger ejecta expanding homologously is divided into
$N(\gg 1)$ layers with different expansion velocities $v_{i}$, where
$v_{1}=v_{\rm min}$ and $v_{N}=v_{\rm max}$. The location of the $i$th layer
at time $t$ is $R_{i}=v_{i}t$, and the mass of the $i$th layer is
${m_{i}}=\int_{{R_{i}}}^{{R_{i+1}}}{4\pi{r^{2}}{\rho_{{\rm{ej}}}}(r,t)dr}$
with (Nagakura et al., 2014)

$\rho_{\rm ej}(r,t)={\frac{(\delta-3)M_{\rm ej}}{4\pi
R_{\max}^{3}}}\left[\left({R_{\min}\over
R_{\max}}\right)^{3-\delta}-1\right]^{-1}\left({r\over
R_{\max}}\right)^{-\delta},$ (2)

where $M_{\rm ej}$ is the total mass of the ejecta. Evolution of the initial
energy $E_{i}$ for the $i$th layer can be described by

${dE_{i}\over dt}=(1-e^{-\Delta\tau_{i}})e^{-\tau_{i}}\xi L_{\rm
md}+m_{i}\dot{q}_{\rm r}\eta_{{\rm th}}-{E_{i}\over R_{i}}{dR_{i}\over
dt}-L_{i}\;.$ (3)

Here, in each layer, the first term on the right-hand side describes the
absorption of PWN emission, the second term is the radioactive heating rate,
the third term is the adiabatic cooling rate, and the last term is the
radiation cooling rate, respectively. The details about the parameters are
presented as follows.

  1. 1.

The power of the pulsar wind $L_{\rm md}$ from the NS can be estimated by MD
radiation, i.e.,

$L_{\rm md}(t)=L_{\rm md,0}\left(1+{t\over t_{\rm sd}}\right)^{-\alpha}$ (4)

with

$L_{\rm md,0}={B_{\rm p}^{2}R^{6}\Omega_{0}^{4}\over{6c^{3}}}=9.6\times
10^{42}R_{6}^{6}B_{\rm p,12}^{2}P_{0,-3}^{-4}~{}\rm erg\cdot s^{-1},$ (5)

where $\Omega_{0}$, $R$, $B_{p}$, $P_{0}$ and $c$ are the initial angular
frequency, the radius, the surface polar magnetic field, the initial spin
period of the NS, and the speed of light, respectively. Based on the fitting
results of Yu et al. (2018) and Ren et al. (2019), the spin-down timescale
$t_{\rm sd}$ can be taken to be

$t_{\rm sd}={\frac{5c^{5}}{128GI\epsilon^{2}\Omega_{0}^{4}}}=9.1\times
10^{5}\epsilon_{-4}^{-2}I_{45}^{-1}P_{0,-3}^{4}~{}\rm s,$ (6)

with $\alpha=1$ for the GW-dominated spin-down loss regime, where $G$ is the
gravitational constant, $I$ is the stellar moment of inertia, and $\epsilon$
is the NS ellipticity. $\xi$ describes the fraction of $L_{\rm md}$ that can
be absorbed by the ejecta. In addition, $\tau_{i}$ is the optical depth from
the innermost layer to the $i$th layer and can be described by
$\tau_{i}=\sum_{1}^{i-1}{\Delta{\tau_{i}}}$ with
$\Delta\tau_{i}=\int_{R_{i}}^{R_{i+1}}{\kappa\rho(r)dr}$.

  2. 2.

The radioactive power per unit mass $\dot{q}_{\rm r}$ and the thermalization
efficiency of the radioactive power $\eta_{\rm th}$ can be estimated by
(Korobkin et al., 2012; Barnes et al., 2016; Metzger, 2019)

$\dot{q}_{\rm r}=4\times 10^{18}\left[{1\over 2}-{1\over\pi}{\rm
arctan}\left({{t-t_{0}}\over\sigma}\right)\right]^{1.3}\rm
erg{\cdot}s^{-1}{\cdot}g^{-1}$ (7)

and

$\eta_{\rm th}=0.36\left[\exp\left(-0.56t_{\rm day}\right)+{\ln(1+0.34t_{\rm
day}^{0.74})\over 0.34t_{\rm day}^{0.74}}\right],$ (8)

respectively. Here, $t_{0}=1.3$ s, $\sigma=0.11$ s, and $t_{\rm day}=t/1~{}\rm
day$.

  3. 3.

The luminosity of the $i$th layer $L_{i}$ is estimated by

$L_{i}=\frac{E_{i}}{\max\left\\{t_{\rm d}^{i},t_{\rm lc}^{i}\right\\}},$ (9)

where the diffusion timescale $t_{\rm d}^{i}$ of photons reads

$t_{\rm d}^{i}\simeq\frac{\kappa}{\beta
R_{i}c}{\sum\limits_{j=i}^{N-1}{m_{j}}},$ (10)

where $t_{\rm lc}^{i}=R_{i}/c$ is the light crossing time. Here $\beta\simeq
13.7$ is adopted (Arnett, 1982).

The total bolometric luminosity $L_{\rm bol}$ of the ejecta is estimated by

$L_{\rm bol}={\sum\limits_{i=1}^{N-1}}L_{i}.$ (11)

We assume that a blackbody spectrum of the ejecta is emitted from the
photosphere at $R_{\rm ph}$ and the effective temperature $T_{\rm eff}$ is
described as (Yu et al., 2013; Xiao et al., 2017; Li et al., 2018)

$T_{\rm eff}=\left({L_{\rm bol}\over 4\pi\sigma_{\rm SB}R_{\rm
ph}^{2}}\right)^{1/4},$ (12)

where $\sigma_{\rm SB}$ is the Stephan-Boltzmann constant. The photosphere
radius $R_{\rm ph}$ is estimated by setting $\tau_{\rm ph}=\int_{R_{\rm
ph}}^{R_{\max}}\rho(r)dr=1$ since $\tau_{\rm tot}>1$. If $\tau_{\rm tot}\leq
1$, we fix $R_{\rm ph}$ to $R_{\min}$. The flux density at frequency $\nu$
from the ejecta is given by

$F_{\nu}^{b}={2\pi h\nu^{3}\over c^{2}}{1\over\exp(h\nu/kT_{\rm
eff})-1}{R_{\rm ph}^{2}\over D_{L}^{2}},$ (13)

where $h$ is the Planck constant and $k$ is the Boltzmann constant.

###  2.2 Emission from the PWN

At the interface between the shocked and unshocked pulsar wind (“termination
shock”), electrons and positrons (leptons, hereafter) carried in the cold
pulsar wind are accelerated, and the magnetic field is amplified. The
accelerated leptons and the amplified magnetic field fill the PWN out to the
radius $R_{\rm PWN}$. Assuming $\epsilon_{B}$ to describe the fraction of the
magnetic energy density of the total energy density behind the shock, the
magnetic energy density $U_{B}^{\rm PWN}$ in the PWN can be parameterized as
(Tanaka & Takahara, 2010, 2013; Murase et al., 2016)

$U_{B}^{\rm PWN}={B_{\rm PWN}^{2}\over{8\pi}}={3\over{4\pi}}\epsilon_{B}R_{\rm
PWN}^{-3}(t)\int_{0}^{t}{L_{\rm md}(s)ds}.$ (14)

Here $R_{\rm PWN}\sim R_{\min}$ is taken because the deceleration timescale of
the ejecta is much larger than the scope of our calculation. A broken power-
law is adopted to describe the energy distribution of leptons in the PWN
(Murase et al., 2015),

${{d\dot{n}_{e}}\over{d\gamma_{e}}}\propto\left\\{\begin{array}[]{cc}\gamma_{e}^{-q_{1}},&\gamma_{m}\leq\gamma_{e}<\gamma_{b},\\\
\gamma_{e}^{-q_{2}},&\gamma_{b}\leq\gamma_{e}\leq\gamma_{M},\end{array}\right.$
(15)

where $q_{1}\sim 1-2$ ($q_{2}\sim 2-3$) is the low (high)-energy spectral
index, $\gamma_{b}\sim 10^{4}-10^{8}$ is the characteristic Lorentz factor of
the accelerated leptons in the PWN, and $\gamma_{m}$ ($\gamma_{M}$) is the
minimum (maximum) Lorentz factor of leptons. In this work, we assume
$\gamma_{m}=3$ and $\gamma_{M}=\sqrt{9m_{e}^{2}c^{4}/(8B_{\rm PWN}q_{e}^{3})}$
(Kumar et al., 2012), where $q_{e}$ is the charge of leptons, and $m_{e}$ is
the electron mass.

For the synchrotron emission of the PWN, two break frequencies are related
with the leptons’ property, i.e., the characteristic synchrotron frequency
$\nu_{b}$ corresponding to $\gamma_{b}$, and the synchrotron cooling frequency
$\nu_{c}$ to $\gamma_{c}={6\pi m_{e}c}/({\sigma_{\rm T}B_{\rm PWN}^{2}t})$.
The frequencies are expressed by

$\nu_{i}\approx{3\over{4\pi}}\gamma_{i}^{2}{q_{e}B_{\rm
PWN}\over{m_{e}c}},\quad i=b,c,m,M,$ (16)

where $\sigma_{\rm T}$ is the Thomson cross section (Sari et al., 1998). The
synchrotron emission is described as follows correspondingly. In the fast-
cooling regime ($\nu_{c}<\nu_{b}$), the synchrotron emission flux density
$L_{\nu}$ at frequency $\nu$ can be expressed by (Murase et al., 2016)

$\nu L_{\nu}^{{}^{\rm PWN}}\approx{{\xi L_{\rm
md}}\over{2R_{b}}}\left\\{\begin{array}[]{lcc}({\nu_{c}\over\nu_{b}})^{2-q_{1}\over
2}({\nu_{a}\over\nu_{c}})^{3-q_{1}\over 2}({\nu\over\nu_{a}})^{5\over
2},&\nu_{m}\leq\nu\leq\nu_{a},\\\ ({\nu_{c}\over\nu_{b}})^{2-q_{1}\over
2}({\nu\over\nu_{c}})^{3-q_{1}\over 2},&\nu_{a}\leq\nu\leq\nu_{c},\\\
({\nu\over\nu_{b}})^{2-q_{1}\over 2},&\nu_{c}\leq\nu\leq\nu_{b},\\\
({\nu\over\nu_{b}})^{2-q_{2}\over
2},&\nu_{b}\leq\nu\leq\nu_{M}.\end{array}\right.$ (17)

In the slow-cooling regime ($\nu_{c}>\nu_{b}$),

$\nu L_{\nu}^{{}^{\rm PWN}}\approx{{\xi L_{\rm
md}}\over{2R_{b}}}\left\\{\begin{array}[]{lcc}({\nu_{b}\over\nu_{c}})^{3-q_{2}\over
2}({\nu_{a}\over\nu_{b}})^{3-q_{1}\over 2}({\nu\over\nu_{a}})^{5\over
2},&\nu_{m}\leq\nu\leq\nu_{a},\\\ ({\nu_{b}\over\nu_{c}})^{3-q_{2}\over
2}({\nu\over\nu_{b}})^{3-q_{1}\over 2},&\nu_{a}\leq\nu\leq\nu_{b},\\\
({\nu\over\nu_{c}})^{3-q_{2}\over 2},&\nu_{b}\leq\nu\leq\nu_{c},\\\
({\nu\over\nu_{c}})^{2-q_{2}\over
2},&\nu_{c}\leq\nu\leq\nu_{M}.\end{array}\right.$ (18)

Where $R_{b}\simeq(2-q_{1})^{-1}+(q_{2}-2)^{-1}$, the radiation efficiency
$\xi=\eta\epsilon_{e}$ with $\eta=\min\\{1,(\nu_{b}/\nu_{c})^{(q_{2}-2)/2}\\}$
(Fan & Piran, 2006), and $\epsilon_{e}=1-\epsilon_{B}$ is adopted. The
critical synchrotron self-absorption (SSA) frequency $\nu_{a}$ can be
calculated with $\tau_{a}(\nu_{a})=1$, and the SSA optical depth is estimated
by (e.g. Panaitescu & Kumar, 2004; Murase et al., 2014)

$\tau_{a}(\nu)\simeq\xi_{q_{1}}\frac{q_{e}n_{\rm ext}R_{\rm PWN}}{B_{\rm
PWN}\gamma_{m}^{5}}\left({\nu\over\nu_{m}}\right)^{(-q_{1}/2-2)},$ (19)

where $\xi_{q_{1}}\simeq 5/3$, and $n_{\rm ext}\simeq n_{\rm PWN}+n_{\rm ej}$
is the number density of leptons. Here, $n_{\rm PWN}\simeq L_{\rm md}/(4\pi
R_{\rm PWN}^{2}\gamma_{b}m_{e}c^{3})$ and $n_{\rm ej}\simeq 3M_{\rm ej}/(4\pi
m_{p}R_{\rm PWN}^{3})$. Based on Equations (17) and (18), one can have
$\int_{0}^{+\infty}{{L_{\nu}}d\nu}\approx\eta{\epsilon_{e}}{L_{{\rm{md}}}}$.
In this work, the effects of the inverse Compton scattering process and the
SSA heating are ignored. The observed flux from a PWN can be expressed as

$F_{\nu}^{{\rm{leak}}}=\frac{{L_{\nu}}{e^{-{\tau_{\rm{tot}}}}}}{{4\pi
D_{L}^{2}}},$ (20)

where
$\tau_{\rm{tot}}=\sum_{1}^{N-1}{\Delta{\tau_{i}}}=\int_{R_{\min}}^{R_{\max}}{\kappa\rho(r)dr}$.

##  3 fitting method

Based on the numerical model developed with section 2, we jointly fit the
observations of AT 2017gfo in UV-optical-IR bands and the late-time X-ray
observations of GRB 170817A. We implement the Markov Chain Monte Carlo (MCMC)
techniques by the use of the Python package emcee (Foreman-Mackey et al.,
2013).

The fitting dataset in this work is considered to cover together the AT
2017gfo observations and X-ray observations of GRB 170817A from 300 days to
1234 days after GW170817 trigger. The data of AT 2017gfo have been obtained
after the extinction correction111https://kilonova.space/. In addition, the
X-ray data are chosen from (2021) with free photon index $\Gamma$ at the
$0.3-10$ keV band. It is allowed for the possible presence of additional
components that may affect the photon index of late observations.

The opacity generated from the bound-free scattering of X-ray photons is much
larger than that for UV-optical-IR photons when the PWN photons crossed the
ejecta. We consider this effect in detail below. On the one hand, the gray
opacity $\kappa_{\rm opt}$ (ignoring the differences of opacity in different
bands) of UV-optical-IR photons is taken as a free parameter, and the optimal
result can be obtained by the MCMC method. On the other hand, when calculating
the flux of PWN emission at $0.3-10$ keV band, the opacity of ejecta is fixed
as $\kappa_{\rm X}(E)=\kappa_{0}(\frac{E}{\rm 1keV})^{-\beta}$, where $E$ is
the energy of X-ray photons, $\kappa_{0}=\kappa_{\rm X}(\rm 1keV)=7\times
10^{3}\;{\rm cm}^{2}\cdot{\rm g}^{-1}$, and $\beta=1$ is assumed (Chen et al.,
2021). Taking into account the variable opacity, the observed flux from PWN at
$0.3-10$ keV band is calculated by

$\mathcal{F}_{\rm PWN}=\int_{\rm 0.3keV}^{\rm 10keV}F_{\rm
X}^{\rm{leak}}(E)dE.$ (21)

Here

$F_{\rm X}^{\rm{leak}}(E)=\frac{L_{\rm X}(E)e^{-\tau_{\rm tot,X}(E)}}{4\pi
D_{L}^{2}}$ (22)

is the flux density of leakage after X-ray photons cross the ejecta, where the
total optical depth from the ejecta of X-ray photons with energy $E$ is
substituted as $\tau_{\rm{tot,X}}(E)=\int_{R_{\min}}^{R_{\max}}{\kappa_{\rm
X}(E)\rho(r)dr}$.

Except for the PWN radiation, the external shock emission of the jet have
still affect the observations to the remnant of GW170817 event. However, the
post-break light curves of GRB 170817A afterglow in different jet models are
different based on, for example, the specific jet profile, the consideration
of the lateral expansion, and the selection of parameters (e.g., Ren et al.,
2020; Lamb et al., 2020; Li & Dai, 2021; Hajela et al., 2021a; , 2021). It is
difficult to proceed with the work if all the models at once are taken into
account. Therefore, we choose a phenomenological formula to describe the flux
at $0.3-10$ keV band of jet emissions,

$\mathcal{F}_{\rm afterglow}=120\cdot t^{-p}\;{\rm erg\cdot cm^{-2}\cdot
s^{-1}},$ (23)

where $t$ ($>200$ days) is the observed time, and $p=2.16$ is the spectrum
index of electrons in the GRB jet. Here $p$ is chosen from the fits of the
early spectra of GRB 170817A afterglow. Thus, the observed total flux
$\mathcal{F}_{\rm total}$ at $0.3-10$ keV band should be considered by

$\mathcal{F}_{\rm total}=\mathcal{F}_{\rm PWN}+\mathcal{F}_{\rm afterglow}.$
(24)

In Ren et al. (2019), the value of the characteristic Lorentz factor of the
accelerated leptons in the PWN is fixed to $\gamma_{b}=10^{4}$, and the
fraction of the magnetic energy density to the total energy density is fixed
to be $\epsilon_{B}=10^{-2}$. Both of them are taken to be free parameters in
this work. Overall, the eleven physical parameters in our models,
{$q_{1},q_{2},\gamma_{b},\epsilon_{B},L_{\rm md,0},t_{\rm sd},M_{\rm
ej},\kappa_{\rm opt},v_{\min},v_{\max},\delta$}, are fitted as free
parameters.

##  4 Fitting Result and Discussion

The model parameters at the $1\sigma$ confidence level are given in Table 1,
where the projections of the posterior distribution for the parameters are
presented in Figure 1. As shown in Figure 2, the joint fitting shows a good
result for the AT 2017gfo observations, and the additional X-ray component
formed by the PWN could explain the flattening of the X-ray observations of
GRB 170817A.

###  4.1 Comparison with previous fitting

Comparing the results in this paper and Ren et al. (2019) in Table 1, we find
that the optimal values of most of the parameters change slightly.
Nonetheless, the spin-down timescale $t_{\rm sd}$, remains the most important
parameter, changing from $2.34\times 10^{5}$ s to $3.96\times 10^{5}$ s. This
change is based on the enlargement of the dataset caused by the addition of
X-ray data. Using Equation 6, the newly fitting result of the ellipticity of
the central NS can be obtained, $\epsilon=1.52\times
10^{-4}I_{45}^{-1/2}P_{0,-3}^{2}$, which is slightly less than that of Ren et
al. (2019). Meanwhile, the initial spin-down luminosity $L_{\rm md,0}$ is the
same as before, which means the early-time kilonova observations could well
estimate the spin-down luminosity of central NS. However, the degeneracy
between $t_{\rm sd}$ and $L_{\rm md,0}$ also indicates the possibility of the
other parameter combinations, as shown in Figure 1. But we note it does not
affect the conclusion here.

There is an important change in the parameters of the spectral energy
distribution in the PWN. We have gotten the optimal results of $q_{1}=1.85$,
$\gamma_{b}=1.66\times 10^{7}$, and $q_{2}\sim 3$ in this work. Since, in the
early fast-cooling stage, the frequencies of most of the PWN photons are
located within the range of $\nu_{c}\leq\nu\leq\nu_{b}$, the fitting of
$q_{1}$ is determined by the observations of AT 2017gfo, and its changeless
result is easy to understand (see Equation 17). We show the evolution of
$\nu_{b}$, $\nu_{c}$, $\nu_{a}$, and $\nu_{M}$ in Figure 3 for the MCMC
results. As shown in Figure 3, lines of $\nu_{b}$ and $\nu_{c}$ cross each
other at $\sim 3000$ days after the merger, and the conversion frequency is
located in the X-ray band. One can find that the PWN is still in the fast
cooling regime in $\sim 3000$ days after the merger. It means that the
constraints on $q_{2}$ from observations are very weak. This is why $q_{2}$ is
not well constrained, as shown in Figure 1. Based on Equation 17, the photon
index of PWN photons is $\Gamma=1+q_{1}/2=1.92$, suggests a soft spectrum as
seen in the X-ray observations (, 2021).

###  4.2 The trend of X-ray observations

The X-ray afterglow of GRB 170817A is predicted in this work to evolve from
flattening to a steep decline, as shown in the left panel of Figure 2. In
terms of the trend of X-ray observation, the kilonova ejecta-PWN system model
and the BH fallback accretion model (Ishizaki et al., 2021a, b; Metzger &
Fernandez, 2021) have similar predictions. Nevertheless, there are certain
differences in the slope of decay predicted by the two models. The kilonova
ejecta-PWN system model has a steeper decay index than that in the BH fallback
accretion model, i.e., steeper than $t^{-5/3}$. In addition, the BH fallback
accretion model predicts an approximate blackbody spectrum of the additional
component. However, it is a power-law spectrum for the additional PWN
component. Future observations will test our model.

What we should note is that the influence of one or more extra sources on
X-ray observations is possible. The models of the kilonova afterglow (Hajela
et al., 2021a; Nedora et al., 2021), and even the count-jet (Li et al., 2019;
, 2021) are predicted to have an impact during this stage. The combined effect
could complicate the evolution of future X-ray observations. To distinguish
these components, continuous observations at the radio band is an important
approach (see Section 5).

###  4.3 The radio emission from PWN

We also check the light curves on radio bands emitted by the PWN. The model
curves at radio bands shown in Figure 4 are calculated with the MCMC fitting
optimal values of the parameters. We find that the PWN’s radio emission flux
is below the observations, which verifies the dominance of the GRB afterglow
to the observations. Because the late-time radio emission from the PWN is
significantly dimmer than the emission from external forward shock of the jet,
it is not expected to see an effect of the PWN in future observations.

More signals of binary NS mergers will be detected by GW detectors in near
future. Due to the Doppler effect of the large off-axis viewing angle of the
merger jet, the prompt GRB radiation is hard to detect for most parts of
events. Meanwhile, the afterglows of merger jets are also going to be much
dimmer. As shown by our results, however, if a long-lived NS is formed during
the binary NS merger, the kilonova ejecta-PWN system could power a long-
lasting radio emission. This trait considerably increases the probability of
detecting EM counterparts to GW events.

##  5 Conclusions

In this paper, we have studied the broadband radiation behavior of the
kilonova ejecta-PWN system, and fitted jointly the observations of AT 2017gfo
in UV-optical-IR bands and the late-time X-ray observations of GRB 170817A.
The kilonova ejecta-PWN system could explain the observations of AT 2017gfo
based on two effects, namely the leakage of PWN radiation and its heating
influence on the merger ejecta. An additional component shows a signature in
the late-time X-ray observations of GRB 170817A when the ejecta is optically
thin to X-ray photons from the PWN. We have shown that the trend of future
X-ray observations will be dominated by the PWN emission, from flattening to a
steep decline, until other sources have their impact, e.g., the kilonova
afterglow.

The constrained ranges of the parameters of the central NS in this work, i.e.,
the dipole magnetic field strength and the ellipticity, are almost the same as
those in the previous works (Yu et al., 2018; Ren et al., 2019) and are
consistent with the limits set by Ai et al. (2020). It is suggested that a
long-lived NS still exists in the merger remnant center a few years after the
GW170817 event. Additionally, the new fitting result of the characteristic
Lorentz factor of the accelerated leptons indicates that the NS wind is highly
magnetized. Note that the parameter values of the PWN are in a reasonable
range compared with the PWN observations in our galaxy.

We also showed the light curves of radio emission from the PWN obtained by the
optimal values of the fitting parameters. We found that the flux of radio
emission from the PWN is lower than the observations of GW170817 event
counterpart. This indicates that the radio observations are dominated by the
GRB jet afterglow, and not affected by PWN radiation even in the future. This
leads to the difference between the kilonova ejecta-PWN system model and the
kilonova afterglow model. Besides, the kilonova ejecta-PWN system model and
the fallback accretion model have similar predictions of the trend of X-ray
observation but certain differences in the slope of decay predicted by the two
models. A difference is that the decay index in the kilonova ejecta-PWN system
model is steeper than that (e.g, $-5/3$) in the fallback accretion model. In
addition, an approximate blackbody spectrum of the additional component is
suggested by the fallback accretion model but a power-law spectrum is
predicted in the kilonova ejecta-PWN system model. Future continuous
observations will test our model.

What should be pointed out is that the radio and X-ray emission from the
kilonova ejecta-PWN system model could be an important probe of binary NS
mergers from days to years after the GW signals. There is a lack of jet-
related signals to trigger the detectors when the axis of a merger-produced
jet is far away from the line of sight. Differently, the quasi-isotropic,
long-existing emission of a kilonova ejecta-PWN system gives rise to more
advantage in observations. Its flux is related to the luminosity and total
energy released by the central engine. This means that the radiation behavior
of a kilonova ejecta-PWN system would directly provide the information of the
merger remnant itself, possibly revealing the properties of a newborn NS and
an accompanying young PWN.

## acknowledgments

RJ would like to thank Ken Chen, Da-Bin Lin and Yun-Wei Yu for helpful
discussions. We thank the referee for his scrutiny, kind comments and helpful
suggestions to help improve this article. This work was supported by the
National Key Research and Development Program of China (grant No.
2017YFA0402600), the National SKA Program of China (grant No. 2020SKA0120300),
and the National Natural Science Foundation of China (grant No. 11833003).

## data availability

The data underlying this paper will be shared on reasonable request to the
corresponding authors.

## References

  * Abbott et al. (2017a) Abbott B. P., et al., 2017a, ApJ, 848, L12
  * Abbott et al. (2017b) Abbott B. P., et al., 2017b, ApJ, 848, L13
  * Ai et al. (2020) Ai S., Gao H., Zhang B., 2020, ApJ, 893, 146
  * Arnett (1982) Arnett W. D., 1982, ApJ, 253, 785
  * Balasubramanian et al. (2021) Balasubramanian A., et al., 2021, ApJ, 914, L20
  * Barnes et al. (2016) Barnes J., Kasen D., Wu M.-R., Martínez-Pinedo G., 2016, ApJ, 829, 110
  * Beniamini et al. (2020) Beniamini P., Duran R. B., Petropoulou M., Giannios D., 2020, ApJ, 895, L33
  * Chen et al. (2021) Chen M.-H., Li L.-X., Lin D.-B., Liang E.-W., 2021, ApJ, 919, 59
  * Fan & Piran (2006) Fan Y., Piran T., 2006, MNRAS, 369, 197
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Geng et al. (2018) Geng J.-J., Dai Z.-G., Huang Y.-F., Wu X.-F., Li L.-B., Li B., Meng Y.-Z., 2018, ApJ, 856, L33
  * Goldstein et al. (2017) Goldstein A., et al., 2017, ApJ, 848, L14
  * Hajela et al. (2021a) Hajela A., et al., 2021a, arXiv e-prints, p. arXiv:2104.02070
  * Hajela et al. (2021b) Hajela A., et al., 2021b, GRB Coordinates Network, 29375, 1
  * Hajela et al. (2021c) Hajela A., et al., 2021c, GRB Coordinates Network, 31231, 1
  * Hamidani & Ioka (2021) Hamidani H., Ioka K., 2021, MNRAS, 500, 627
  * Hamidani et al. (2020) Hamidani H., Kiuchi K., Ioka K., 2020, MNRAS, 491, 3192
  * Ishizaki et al. (2021a) Ishizaki W., Ioka K., Kiuchi K., 2021a, ApJ, 916, L13
  * Ishizaki et al. (2021b) Ishizaki W., Kiuchi K., Ioka K., Wanajo S., 2021b, ApJ 922, 185
  * Kasen & Bildsten (2010) Kasen D., Bildsten L., 2010, ApJ, 717, 245
  * Kasen et al. (2015) Kasen D., Fernández R., Metzger B. D., 2015, MNRAS, 450, 1777
  * Korobkin et al. (2012) Korobkin O., Rosswog S., Arcones A., Winteler C., 2012, MNRAS, 426, 1940
  * Kotera et al. (2013) Kotera K., Phinney E. S., Olinto A. V., 2013, MNRAS, 432, 3228
  * Kumar et al. (2012) Kumar P., Hernández R. A., Bošnjak Ž., Barniol Duran R., 2012, MNRAS, 427, L40
  * Lamb et al. (2020) Lamb G. P., Levan A. J., Tanvir N. R., 2020, ApJ, 899, 105
  * Lazzati & Perna (2019) Lazzati D., Perna R., 2019, ApJ, 881, 89
  * Lazzati et al. (2020) Lazzati D., Ciolfi R., Perna R., 2020, ApJ, 898, 59
  * Li & Dai (2021) Li L., Dai Z.-G., 2021,ApJ, 918, 52
  * Li et al. (2018) Li S.-Z., Liu L.-D., Yu Y.-W., Zhang B., 2018, ApJ, 861, L12
  * Li et al. (2019) Li L.-B., Geng J.-J., Huang Y.-F., Li B., 2019, ApJ, 880, 39
  * Lyutikov (2020) Lyutikov M., 2020, MNRAS, 491, 483
  * Makhathini et al. (2020) Makhathini S., et al., 2021, ApJ, 922, 154
  * Matsumoto et al. (2018) Matsumoto T., Ioka K., Kisaka S., Nakar E., 2018, ApJ, 861, 55
  * Metzger (2019) Metzger B. D., 2019, Living Reviews in Relativity, 23, 1
  * Metzger & Fernandez (2021) Metzger B. D., Fernandez R., 2021, ApJ, 916, L3
  * Metzger et al. (2018) Metzger B. D., Thompson T. A., Quataert E., 2018, ApJ, 856, 101
  * Murase et al. (2014) Murase K., Thompson T. A., Ofek E. O., 2014, MNRAS, 440, 2528
  * Murase et al. (2015) Murase K., Kashiyama K., Kiuchi K., Bartos I., 2015, ApJ, 805, 82
  * Murase et al. (2016) Murase K., Kashiyama K., Mészáros P., 2016, MNRAS, 461, 1498
  * Murase et al. (2018) Murase K., et al., 2018, ApJ, 854, 60
  * Nagakura et al. (2014) Nagakura H., Hotokezaka K., Sekiguchi Y., Shibata M., Ioka K., 2014, ApJ, 784, L28
  * Nedora et al. (2021) Nedora V., Radice D., Bernuzzi S., Perego A., Daszuta B., Endrizzi A., Prakash A., Schianchi F., 2021, MNRAS, 506, 5908
  * Panaitescu & Kumar (2004) Panaitescu A., Kumar P., 2004, MNRAS, 353, 511
  * Pavan et al. (2021) Pavan A., Ciolfi R., Kalinani J. V., Mignone A., 2021, MNRAS,
  * Ren et al. (2019) Ren J., Lin D.-B., Zhang L.-L., Li X.-Y., Liu T., Lu R.-J., Wang X.-G., Liang E.-W., 2019, ApJ, 885, 60
  * Ren et al. (2020) Ren J., Lin D.-B., Zhang L.-L., Wang K., Li X.-Y., Wang X.-G., Liang E.-W., 2020, ApJ, 901, L26
  * Sari et al. (1998) Sari R., Piran T., Narayan R., 1998, ApJ, 497, L17
  * Savchenko et al. (2017) Savchenko V., et al., 2017, ApJ, 848, L15
  * Tanaka & Takahara (2010) Tanaka S. J., Takahara F., 2010, ApJ, 715, 1248
  * Tanaka & Takahara (2013) Tanaka S. J., Takahara F., 2013, MNRAS, 429, 2945
  * (51) Troja E., O’Connor B., Ryan G., et al., 2021, MNRAS, 510, 1902
  * Wu et al. (2021) Wu G.-L., Yu Y.-W., Zhu J.-P., 2021, A&A, 654, 124
  * Xiao et al. (2017) Xiao D., Liu L.-D., Dai Z.-G., Wu X.-F., 2017, ApJ, 850, L41
  * Yu et al. (2013) Yu Y.-W., Zhang B., Gao H., 2013, ApJ, 776, L40
  * Yu et al. (2018) Yu Y.-W., Liu L.-D., Dai Z.-G., 2018, ApJ, 861, 114
  * Zhang et al. (2018) Zhang B.-B., et al., 2018, Nature Communications, 9, 447

Table 1: Parameters estimated from the MCMC sampling.

[b]  Parameter Constraint Range1 Pervious ${\rm log}_{10}L_{\rm md,0}$
($\rm{erg\cdot s}^{-1}$) $41.317_{-0.023}^{+0.023}$ [39, 43]
$1.85_{-0.23}^{+0.85}\times 10^{41}$ ${\rm log}_{10}t_{\rm sd}$ (s)
$5.598_{-0.044}^{+0.044}$ [3, 7] $2.34_{-0.94}^{+0.57}\times 10^{5}$
$M_{\rm{ej}}/0.01M_{\odot}$ $3.044_{-0.033}^{+0.039}$ [1, 10]
$3.52_{-0.11}^{+0.06}$ $\kappa_{\rm opt}$ ($\rm{cm}^{2}\cdot\rm{g}^{-1}$)
$1.104_{-0.028}^{+0.031}$ [0.1, 10] $1.69_{-0.09}^{+0.06}$ $v_{\rm{min}}/c$
$0.099_{-0.002}^{+0.002}$ [0.05, 0.15] $0.10_{-0.00}^{+0.01}$ $v_{\rm{max}}/c$
$0.473_{-0.020}^{+0.022}$ [0.3, 0.6] $0.34_{-0.01}^{+0.01}$ $\delta$
$3.834_{-0.098}^{+0.108}$ [1, 5] $2.47_{-0.32}^{+0.02}$ $q_{1}$
$1.845_{-0.013}^{+0.012}$ [1, 2] $1.83_{-0.05}^{+0.02}$ $q_{2}$
$2.947_{-0.056}^{+0.032}$ [2, 3] $2.25_{-0.05}^{+0.10}$ ${\rm
log}_{10}{\gamma_{b}}$ $7.221_{-0.046}^{+0.045}$ [4, 8] $10^{4}$ ${\rm
log}_{10}{\epsilon_{B}}$ $-1.579_{-0.064}^{+0.107}$ [-3,-0.5] $0.01$
$\kappa_{\rm X}({\rm 1keV})$ ($\rm{cm}^{2}\cdot\rm{g}^{-1}$) $7\times 10^{3}$
– – $\alpha$ 1 – 1

  * 1

Priors are uniformly distributed.

Figure 1: Posterior probability density contours for the physical parameters
from the MCMC sampling.

Figure 2: Left Panel: The X-ray light curve of GRB 170817A afterglow with the
addition of PWN emission. Data are taken from the free-$\Gamma$ column,
$0.3-10$ keV flux, in Table 1 of (2021). We also mark the new observation at
$\sim 1575$ days with a star, as given by Hajela et al. (2021c), where the
photon index $\Gamma=1.6$ was fixed but the $1\sigma$ uncertainties were not
reported. We notice the data is consistent with our model. Right Panel: The
multiband light curves of AT 2017gfo fitted by our model. The data are taken
from https://kilonova.space/ and the extinction correction is performed.

Figure 3: The evolution of $\nu_{b}$, $\nu_{c}$, $\nu_{a}$, and $\nu_{M}$ of
the PWN using the MCMC optimal results. The picture shows that the PWN is in
the fast-cooling regime in $\sim 3000$ days after the merger. The transition
of $\nu_{b}$ and $\nu_{c}$ in the X-ray band will lead to softening of the
X-ray spectrum.

Figure 4: Multiband radio light curves emitted from the kilonova ejecta-PWN
system of the GW170817 event remnant are calculated by the optimal parameters
estimated from the MCMC sampling. The same bands are represented by the same
colors correspondingly, where the observational data are described with
circles, and the triangles are the upper limits. The data are chosen from
Makhathini et al. (2020) and Balasubramanian et al. (2021).