printacmref=false ifaamas [AAMAS ’21]Proc. of the 20th International
Conference on Autonomous Agents and Multiagent Systems (AAMAS 2021)May 3–7,
2021London, UKU. Endriss, A. Nowé, F. Dignum, A. Lomuscio (eds.) 2021 2021 ???
Department of Computer Science University of Oxford Department of Computer
Science University of Oxford Department of Computer Science University of
Oxford Department of Computer Science University of Oxford

# Semi-On-Policy Training for Sample Efficient Multi-Agent Policy Gradients

Bozhidar Vasilev<EMAIL_ADDRESS>,  Tarun Gupta
<EMAIL_ADDRESS>,  Bei Peng<EMAIL_ADDRESS>and  Shimon Whiteson
<EMAIL_ADDRESS>
###### Abstract.

Policy gradient methods are an attractive approach to multi-agent
reinforcement learning problems due to their convergence properties and
robustness in partially observable scenarios. However, there is a significant
performance gap between state-of-the-art policy gradient and value-based
methods on the popular StarCraft Multi-Agent Challenge (SMAC) benchmark. In
this paper, we introduce semi-on-policy (SOP) training as an effective and
computationally efficient way to address the sample inefficiency of on-policy
policy gradient methods. We enhance two state-of-the-art policy gradient
algorithms with SOP training, demonstrating significant performance
improvements. Furthermore, we show that our methods perform as well or better
than state-of-the-art value-based methods on a variety of SMAC tasks.

###### Key words and phrases:

Reinforcement Learning; Policy Gradient; Multi-Agent; SMAC

##  1\. Introduction

Many real-world learning problems, such as autonomous vehicle coordination
(Cao et al., 2012) and network packet delivery (Ye et al., 2015), can be
naturally modelled as cooperative multi-agent reinforcement learning (MARL)
problems. Reinforcement learning (RL) algorithms designed for single-agent
settings are typically unsuitable for multi-agent scenarios, because the joint
action space grows exponentially with the number of agents. Furthermore, even
when the joint action space is not prohibitively large, communication
constraints may make it impossible to apply single-agent algorithms directly.

This often makes it necessary to learn decentralised policies, where the
agents make decisions independently, based only on their own local
observations. However, if the agents learn independently, this makes the
environment non-stationary from the viewpoint of any agent, so learning might
not converge. To address this issue, one often adopts the framework of
centralised training with decentralised execution (CTDE) (Oliehoek et al.,
2008; Kraemer and Banerjee, 2016), assuming a simulated or controlled setting.
The CTDE framework allows the algorithm to use any extra information available
during training. However, at test time each agent must make decisions based
only on its private observations.

The StarCraft Multi-Agent Challenge (SMAC) Samvelyan et al. (2019) is a
challenging MARL benchmark in which agents must learn complex coordination
behaviours in an environment with significant partial observability. In recent
years SMAC has become a popular environment for evaluating cooperative MARL
algorithms Rashid et al. (2020a); Mahajan et al. (2019); Wang et al. (2020);
Zhou et al. (2020). Despite the progress in multi-agent policy gradient (MAPG)
methods, most of them still struggle to compete with state-of-the-art value-
based methods on SMAC. Rashid et al. (2020b) show that COMA (Foerster et al.,
2018), a state-of-the-art on-policy MAPG algorithm, severely underperforms
even simple value-based methods. Papoudakis et al. (2020) go as far as using
$10$ times more data for MAPG methods when comparing against value-based ones.
They show that state-of-the-art on-policy MAPG methods, such as CentralV
Foerster et al. (2018) and COMA, suffer from poor performance even on simple
SMAC scenarios.

While off-policy algorithms are generally more sample efficient than on-policy
ones, state-of-the-art off-policy MAPG methods, such as MADDPG (Lowe et al.,
2017) and Multi-Agent Soft Actor-Critic (MASAC) (Haarnoja et al., 2018; Iqbal
and Sha, 2019), also perform poorly on SMAC (Rashid et al., 2020a). The
recently proposed stochastic decomposed policy gradient (DOP) Wang et al.
(2020) algorithm mixes off-policy and on-policy training, by adapting the tree
backup technique (Precup et al., 2000; Munos et al., 2016a) to the multi-agent
setting and taking advantage of a factored critic representation. The authors
claim their method is the first MAPG algorithm that outperforms state-of-the-
art value-based methods on SMAC. However, they only present results for a
small number of SMAC maps.

In this paper, we improve COMA by fixing a representation issue with the
centralised critic, which produces inconsistent value estimates for some joint
actions due to its network architecture. We further improve COMA by changing
its training scheme to use the entire batch of data, rather than minibatches.
We find that our version, COMA with a consistent critic (COMA-CC),
consistently outperforms the original COMA on SMAC, and provide ablation
studies to evaluate the effect of each change.

Furthermore, we introduce semi-on-policy (SOP) training as an effective and
computationally efficient way to address the sample inefficiency of on-policy
policy gradient algorithms. Key to our method is using recent off-policy data
for training, based on an eligibility criterion. We apply SOP training to
CentralV and COMA-CC, and demonstrate significant performance improvements.
Our methods achieve similar or better performance and sample efficiency,
compared to state-of-the-art value-based methods, on a variety of challenging
SMAC tasks. Additionally, our best performing method often outperforms other
state-of-the-art MAPG methods, including DOP.

##  2\. Related Work

##### Value-based methods

A natural value-based MARL method for learning decentralised policies is
independent $Q$-learning Tan (1993), where each agent learns an individual
action-value function and acts based on it; this is typically done using DQN
Mnih et al. (2015); Tampuu et al. (2017). To exploit CTDE, Sunehag et al.
(2018) propose value decomposition networks (VDN), where a centralised but
factored action-value function is learned, which is the sum of the individual
agent utilities. The QMIX algorithm Rashid et al. (2020b) expands this idea by
learning a factored joint action-value function that is only constrained to be
monotonic with respect to each agent’s utility. Recent works, like Weighted-
QMIX Rashid et al. (2020a) and Qatten Yang et al. (2020), aim to address the
limitations imposed by the monotonicity constraint of QMIX, while still using
a monotonic factorisation. On the other hand, works like QTRAN Son et al.
(2019) address the monotonic restriction of QMIX by learning a VDN-factored
joint action-value function along with an unrestricted centralised state-value
function, and QTRAN++ Son et al. (2020) address the gap between the
theoretical and empirical results of QTRAN.

##### Policy gradient methods

The most successful policy gradient methods have been actor-critic methods,
where both a policy and a value function are learned during training, and the
estimates of the value function are used to train the policy Sutton and Barto
(1998). If each agent learns independently, we get the independent actor-
critic (IAC) algorithm, which is an adaptation of the advantage actor-critic
(A2C) Mnih et al. (2016) algorithm to the multi-agent setting. The CentralV
algorithm Foerster et al. (2018) is another adaptation of A2C, but using a
centralised critic that makes use of global state information for learning. To
address the issue of multi-agent credit assignment, the COMA algorithm
Foerster et al. (2018) uses a counterfactual baseline, which marginalises out
an agent’s action when computing the advantage of their action. Zhou et al.
(2020) address the same problem with their technique of learning implicit
credit assignment (LICA) by formulating a centralised critic as a
hypernetwork. To improve sample efficiency, the stochastic decomposed policy
gradients algorithm (DOP) Wang et al. (2020) uses a mixture of on-policy and
off-policy data by adapting the tree backup technique (Precup et al., 2000;
Munos et al., 2016a) to the multi-agent setting and taking advantage of a
factored critic representation. Iqbal and Sha (2019) adapt the soft actor-
critic (SAC) Haarnoja et al. (2018) algorithm to the multi-agent setting and
combine it with an attention mechanism and a centralised critic. The MADDPG
Lowe et al. (2017) algorithm is an adaptation of the DDPG algorithm Lillicrap
et al. (2019) to the multi-agent setting; while originally the algorithm was
used to learn deterministic policies in continuous domains using off-policy
data, it has since been adapted to discrete domains as well Rashid et al.
(2020a); Iqbal and Sha (2019).

##### SMAC

The StarCraft Multi-Agent Challenge (SMAC) (Samvelyan et al., 2019) is a
popular cooperative multi-agent reinforcement learning benchmark. The
challenges require agents to learn complex cooperative behaviours, such as
kiting and focus fire. For a long time, value-based methods have significantly
outperformed policy gradient methods on SMAC. Some recent works, such as LICA
Zhou et al. (2020) and DOP Wang et al. (2020) have shown that actor-critic
methods can sometimes outperform value-based methods on SMAC. However, they
only focus on a limited set of maps, and their algorithms are much more
complex compared to previous work.

Our work focuses on improving the sample efficiency of on-policy policy
gradient algorithms by directly using off-policy data during training, based
on an eligibility criterion. In principle, this makes it applicable to any
such method. In contrast to methods that correct for off-policy data using
importance sampling or returns corrections Wang et al. (2016); Munos et al.
(2016b); Wang et al. (2020); Precup et al. (2000); Harutyunyan et al. (2016);
Foerster et al. (2017), we simply require the generating policy for all
training data to be closely related to the current behaviour policy. This
removes the associated computational cost of applying corrections to the data.
Furthermore, in contrast to proximal policy optimisation (PPO) Schulman et al.
(2017), which uses a ”surrogate” objective to constrain policy updates over
multiple training iterations on a fixed set of on-policy data, our method
optimises the true objective of the base algorithm and uses a rolling-window
replay buffer that contains recent off-policy data, along with on-policy data.
In contrast to Zhou et al. (2020); Wang et al. (2020), we show that a simple
base algorithm like CentralV suffices for many challenging SMAC tasks, with
our method of training.

##  3\. Background

##### Dec-POMDPs

We consider a fully cooperative multi-agent task, which can be modelled as a
Decentralised Partially Observable Markov Decision Process (Dec-POMDP)
(Oliehoek and Amato, 2016) consisting of a tuple $G=\left\langle
S,U,P,r,Z,O,A,\gamma\right\rangle$. The state of the environment at each time
step is $s\in S$. Each of $n$ agents $a\in A\equiv\left\\{1,\dots,n\right\\}$
chooses an action $u^{a}\in U\equiv\left\\{1,\dots,m\right\\}$, forming the
joint action $\mathbf{u}\in\mathbf{U}:=U^{n}$. The environment dynamics are
modelled by a state transition function $P:S\times\mathbf{U}\times
S\rightarrow[0,1]$, where $P(s^{\prime}|\mathbf{u},s)$ is the probability that
the environment transitions to state $s^{\prime}$ given that the current state
is $s$ and the joint action is $\mathbf{u}$. All agents share a global reward
function $r:S\times\mathbf{U}\rightarrow\mathbb{R}$. We consider partially
observable scenarios, where an observation function $O:S\times A\rightarrow Z$
individually determines each agent’s observation $z^{a}$ at each time step.
Each agent has a private action-observation history $\tau^{a}\in T:=(Z\times
U)^{*}$, on which it conditions a stochastic policy $\pi^{a}:U\times
T\rightarrow[0,1]$. The joint policy $\boldsymbol{\pi}:\mathbf{U}\times
T^{n}\rightarrow[0,1]$ admits a joint state-value function
$V^{\boldsymbol{\pi}}(s_{t}):=\mathbb{E}_{\boldsymbol{\pi}}[R_{t}|s_{t}]$

and action-value function
$Q^{\boldsymbol{\pi}}(s_{t},\mathbf{u}_{t})=\mathbb{E}_{\boldsymbol{\pi}}[R_{t}|s_{t},\mathbf{u_{t}}]$,
where $R_{t}=\sum_{i=0}^{\infty}\gamma^{i}r_{t+i}$ is the discounted return.

##### CentralV and COMA

Both CentralV and COMA (Foerster et al., 2018) learn decentralised policies
$\pi^{a}$ for each agent $a$, which induce a joint policy $\boldsymbol{\pi}$.
During training they maximise the objective
$J=V^{\boldsymbol{\pi}}(S_{1})=\mathbb{E}_{\boldsymbol{\pi}}[R_{1}]$,
following the policy gradient:

$\nabla
J=\mathbb{E}_{\boldsymbol{\pi}}\left[\sum_{a}\nabla\log\pi^{a}(u^{a}|\tau^{a})A_{a}^{\boldsymbol{\pi}}(s,\mathbf{u})\right],$
(1)

where $A_{a}^{\boldsymbol{\pi}}(s,\mathbf{u})$ is an advantage function,
which, intuitively, estimates how much better or worse a joint action
$\mathbf{u}$ is compared to the agents’ expected performance under joint
policy $\boldsymbol{\pi}$. As in most cooperative MARL, the agents share
policy parameters (Rashid et al., 2020b; Foerster et al., 2018; Samvelyan et
al., 2019).

In CentralV, the critic approximates the joint state-value function
$V_{\psi}(s_{t})\approx V^{\boldsymbol{\pi}}(s_{t})$, which is used to
estimate the advantage for all agents:

$\displaystyle A_{a}^{\boldsymbol{\pi}}(s_{t},\mathbf{u}_{t})$
$\displaystyle=r_{t}+\gamma
V^{\boldsymbol{\pi}}(s_{t+1})-V^{\boldsymbol{\pi}}(s_{t})$ (2)
$\displaystyle\approx r_{t}+\gamma V_{\psi}(s_{t+1})-V_{\psi}(s_{t})$ (3)

In COMA, the critic approximates the joint action-value function
$Q_{\psi}(s_{t},\mathbf{u}_{t})\approx
Q^{\boldsymbol{\pi}}(s_{t},\mathbf{u}_{t})$, which is used to estimate a
counterfactual advantage for each agent:

$\displaystyle A^{\boldsymbol{\pi}}_{a}(s_{t},\mathbf{u}_{t})$ (4)
$\displaystyle=Q^{\boldsymbol{\pi}}(s_{t},\mathbf{u}_{t})-\sum_{u^{\prime}\in
U}\pi^{a}(u^{\prime}|\tau^{a})Q^{\boldsymbol{\pi}}(s_{t},\mathbf{u}_{t}^{-a},u^{\prime})$
(5) $\displaystyle\approx Q_{\psi}(s_{t},\mathbf{u}_{t})-\sum_{u^{\prime}\in
U}\pi^{a}(u^{\prime}|\tau^{a})Q_{\psi}(s_{t},\mathbf{u}_{t}^{-a},u^{\prime}),$
(6)

where $Q_{\psi}(s_{t},\mathbf{u}_{t}^{-a},u^{\prime})$ is the action-value
function evaluated by keeping the actions of all agents except $a$ fixed, and
replacing agent $a$’s action with $u^{\prime}$. The quantity:

$\sum_{u^{\prime}\in
U}\pi^{a}(u^{\prime}|\tau^{a})Q^{\boldsymbol{\pi}}(s_{t},\mathbf{u}_{t}^{-a},u^{\prime})$
(7)

is called the counterfactual baseline.

In both CentralV and COMA, the critic $f_{\psi}$ is trained to minimise the
loss $\mathcal{L}(\psi)=(y^{(\lambda)}-f_{\psi}(\cdot_{t}))^{2}$, using a
variant of TD($\lambda$) (Sutton and Barto, 1998). The targets
$y^{(\lambda)}=(1-\lambda)\sum_{n=1}^{\infty}\lambda^{n-1}G_{t}^{(n)}$ are a
mixture of $n$-step returns:

$G_{t}^{(n)}=\sum_{i=0}^{n-1}\gamma^{i}r_{t+i}+\gamma^{n}f_{\psi^{\prime}}(\cdot_{t+n}),$
(8)

which are calculated using bootstrapped values from a target network Mnih et
al. (2015). Periodically, the target network parameters are copied from
$\psi$.

##  4\. Method

In this section, we present our semi-on-policy (SOP) training framework and
our COMA-CC algorithm.

###  4.1. Semi-on-policy training

The key insight behind our method is that, if the policy of the agents only
changes slightly after a step of training, then the data generated by the
previous policy can be reused in the next training step. In particular, we
assume that if an episode $e$ was generated following a policy $\pi^{\prime}$,
and $KL(\pi\;\|\;\pi^{\prime})\leq$ kl_threshold, then we can still use the
episode $e$ for training our policy $\pi$.111The choice of using the KL-
divergence between policies to quantify similarity is somewhat arbitrary, and
other criterions could also be used.

Thus, if the base on-policy algorithm induces gradual policy updates, we can
train it using the most recent episodes, despite the fact that some of them
were generated by a different policy. To take advantage of this, we propose
permissive semi-on-policy training (Algorithm 1). Our results on SMAC show
that this is an effective and reliable way of improving the sample efficiency
of on-policy algorithms. For cases when the policy might diverge significantly
during one update, we propose strict semi-on-policy training (Algorithm 2),
which enforces the KL constraint between policies in the replay buffer.
Because our analysis in Section 7.1 shows that permissive SOP training
(henceforth simply called SOP training) is generally stable on SMAC, we defer
evaluating strict SOP training to future work.

1 Let $b$ be the batch-size.

2 Initialise the current policy $\pi$ and replay buffer buf of size $b$.

3 Fill buf with $b$ episodes.

4 while  _there is still time to train_ do

5 Train $\pi$ using all data in buf.*

6 Discard the least recent episode from buf.

7 Sample a new episode from the environment following the current policy $\pi$
and insert it into buf.

8

// (*) this step depends on how the particular policy gradient algorithm
performs policy updates. For actor-critic algorithms: also train the critic.

Algorithm 1 Permissive semi-on-policy training

1 Let $b$ be the batch-size.

2 Let kl_threshold be a hyperparameter.

3 Initialise the current policy $\pi$, replay buffer buf of size $b$ and
corresponding policy buffer pbuf.

4 // Maintain the invariant that buf[i] was generated by following the policy
pbuf[i].

5 while  _there is still time to train_ do

6 Sample new episodes from the environment following the current policy $\pi$
and insert them into buf until it is full. Fill corresponding locations of
pbuf with $\pi$.

7 Train $\pi$ using all data in buf.*

8 Discard the least recent episode from buf and all episodes generated by a
sufficiently different policy.**

9

10// (*) this step depends on how the particular policy gradient algorithm
performs policy updates. For actor-critic algorithms: also train the critic.

// (**) we say $\pi^{\prime}$ is sufficiently different from $\pi$ if
$KL(\pi\;\|\;\pi^{\prime})>$ kl_treshold.

Algorithm 2 Strict semi-on-policy training

####  4.1.1. Implications on batch size

Typical implementations of value-based MARL algorithms use a batch size of
$32$ episodes on SMAC (Rashid et al., 2020b; Rashid et al., 2020a; Mahajan et
al., 2019). Because value-based methods generally use $Q$-learning (Watkins,
1989) to train their value-function estimator (Rashid et al., 2020b; Rashid et
al., 2020a; Mahajan et al., 2019), they can reuse their old data. On-policy
algorithms, however, must discard their data after training, so they typically
opt for a smaller batch size (Rashid et al., 2020b), a massive increase in the
number of episodes sampled from the environment (Zhou et al., 2020), or both
(Papoudakis et al., 2020). If policy updates are stable across a small number
of training steps, our method of training allows us to perform the same number
of training steps as the currently popular off-policy value-based MARL
methods, without requiring more data.

###  4.2. COMA-CC

In this section we describe the improvements our COMA-CC algorithm makes to
the original COMA.

####  4.2.1. Fixing the critic inconsistency

The COMA critic computes the counterfactual $Q$-function for each agent:

$Q^{a}(s_{t},\mathbf{u}_{t}):=\left\langle
Q(s_{t},u_{t}^{a}=1,\mathbf{u}^{-a}_{t}),\dots,Q(s_{t},u_{t}^{a}=|U|,\mathbf{u}^{-a}_{t})\right\rangle$
(9)

varying the action of agent $a$, while keeping other agents’ actions fixed. In
particular it estimates:

$Q^{a}(s_{t},\mathbf{u}_{t})\approx
Q_{\psi}(s_{t},z_{t}^{a},\mathbf{u}_{t-1},\mathbf{u}_{t}^{-a},a)$ (10)

where $Q_{\psi}$ is a neural network with parameters $\psi$. Here $s_{t}$ is
the central state, $z_{t}^{a}$ is the observation of agent $a$ at time step
$t$, $\mathbf{u}_{t-1}$ is the joint action at time step $(t-1)$,
$\mathbf{u}_{t}^{-a}$ is the joint action at time step $t$ with the action of
agent $a$ masked out, and $a$ is the agent id.

Consider a simple scenario where we have two agents, $a_{1},a_{2}$, and they
can perform a single action, $u$. Let
$Q^{1}=f(s_{t},z_{t}^{1}),Q^{2}=f(s_{t},z_{t}^{2})$ be the (neural network)
estimates of the two agents’ action-value functions. Since, in general,
$z_{t}^{1}\neq z_{t}^{2}$, this means that, in general, $Q^{1}(u,u)\neq
Q^{2}(u,u)$. Thus, the agents’ estimates for the actual joint action taken at
that time step are inconsistent with each other. The observation difference
among agents is one high-level inconsistency that causes this issue, while
another is having the agent id as input. The way that $\mathbf{u}_{t}^{-a}$ is
implemented (as a concatenation of one-hot action-vectors with one of the
actions masked out) has a similar impact.

One of the intuitive issues this causes is that two agents might disagree
about whether the same joint action was good or bad. Thus, when they compute
their respective share of the gradient, one will move the policy parameters
towards making the joint action more likely, while the other will move away
from it. Since the agents share parameters, this can cause convergence issues.

To address this issue, COMA-CC uses an alternative critic that estimates:

$Q(s_{t},\mathbf{u}_{t})\approx
Q_{\phi}(s_{t},z_{t}^{1},\dots,z_{t}^{{n}},\mathbf{u}_{t-1},\mathbf{u}_{t})$
(11)

We can use this critic to compute each component of
$Q^{a}(s_{t},\mathbf{u}_{t})$ by inputting the respective counterfactual joint
actions. Furthermore, this critic always receives the same inputs when
estimating the same joint action, resolving the issue of inconsistent
estimates.

The new critic can compute all necessary counterfactual $Q$-values in a single
forward pass; this is done by building the input for each counterfactual
$Q$-value computation and running the concatenated inputs through the network.
For each counterfactual $Q$-value computation, the COMA critic requires $n$
inputs, one for each agent, and the COMA-CC critic requires $nm$ inputs, one
for each agent and counterfactual joint action; therefore, if the internal
network architectures are similar for both critics, the COMA-CC critic must
perform $O(m)$ times more multiplications than the COMA critic. To ensure
this, the concatenated observations $(z_{t}^{1},\dots,z_{t}^{n})$ must be
compressed before they are used as inputs to the critic; this can be done with
an encoding network. Since the running time of our experiments was essentially
the same for both the COMA and COMA-CC critics, we do not include this
optimisation in our implementation.222A more thorough analysis of the
computational complexity of the COMA-CC critic is included in Appendix B.

####  4.2.2. Critic training

The critic in COMA is trained according to the procedure given in Algorithm 3.
This way of training means that each training iteration uses only batch_size
steps of agent-environment interaction to compute the gradients: for example,
if the batch_size is $8$ episodes, and each episode has length $128$, COMA
will train on minibatches of $8$ steps until it trains on all data from the
episodes. Our ablative study in Section 7.2 indicates that this destabilises
learning; therefore, we train COMA-CC on the whole batch at once. The change
necessary is to the for-loop, and is given in Algorithm 4.

1 Initialise the critic parameters $\phi$, target critic parameters
$\hat{\phi}$ = $\phi$, and policy parameters $\theta$

2 …

3 [Let $T$ be the maximum length of an episode.]

4 // For all agents $a\in A$ simultaneously:

5 Calculate TD($\lambda$) targets $y_{t}^{a}$ using the target critic
$\hat{\phi}$

6 for  _t = T to 1_ do

7 $\Delta Q_{t}^{a}=y_{t}^{a}-Q_{\phi}(s_{t},\mathbf{u}_{t})$

8 $\Delta\phi=\nabla_{\phi}(\Delta Q_{t}^{a})^{2}$ // calculate critic
gradient

9 $\phi=\phi-\alpha\Delta\phi$ // update critic parameters

10 Every $C$ steps reset $\hat{\phi}=\phi$

11

12…

Algorithm 3 COMA Critic Training (Foerster et al., 2018)

1 $\Delta\phi=0$

2 for  _t = T to 1_ do

3 $\Delta Q_{t}^{a}=y_{t}^{a}-Q_{\phi}(s_{t},\mathbf{u}_{t})$

4 $\Delta\phi=\Delta\phi+\nabla_{\phi}(\Delta Q_{t}^{a})^{2}$ // accumulate
critic gradient

5

6$\phi=\phi-\alpha\Delta\phi$ // update critic parameters

7 Every $C$ steps reset $\hat{\phi}=\phi$

Algorithm 4 Alternative COMA Critic Training

##  5\. Experimental setup

Figure 1. Comparing MAPG methods with QMIX. Batch size $32$. Easy maps (top):
[2s3z, 3s5z, 1c3s5z], Hard maps (bottom): [bane_vs_bane, 2c_vs_64zg], Super-
hard maps (bottom): [MMM2].

6 plots with titles 2s3z, 3s5z, 1c3s5z, bane_vs_bane, 2c_vs_64zg, MMM2. X-axis
labelled T (mil). Y-axis labelled Median Test Win %. Legend contains CentralV-
SOP, COMA-CC-SOP, QMIX, DOP.

We evaluate all methods on the challenging SMAC benchmark (Samvelyan et al.,
2019), using a batch size of either $8$ or $32$ episodes. In all comparisons,
all methods use the same batch size; we indicate this in the figure notes. We
run $4$ random seeds per map for each method.

For QMIX and DOP we use the author-provided open-source
implementations.333https://github.com/oxwhirl/pymarl and
https://github.com/TonghanWang/DOP. We create fresh implementations of
CentralV and COMA-CC, based on the existing CentralV and COMA
implementations.444https://github.com/oxwhirl/pymarl. The changes we make are
outlined in Section 4. Additionally, we set $\gamma=1$ when computing the
advantage function for CentralV, which we find helps stabilise policy updates;
we include results for the version without this change in Appendix C.2.

All experiments, including baselines, use StarCraft II Version 2.4.10.555At
the time of writing, the latest version of the PyMARL framework provided by
the original SMAC authors uses SC2.4.10. Older versions of the framework used
SC2.4.6.2.69232. The SMAC authors note that results are not always comparable
between versions. We have retested all baselines to ensure fair comparisons.
We use the same hyperparameters, network architectures, and evaluation
methodology as in the original SMAC paper Samvelyan et al. (2019), unless
stated otherwise. We list these in detail in Appendix A. Doing so, we have
attempted to eliminate all possible confounding factors in our COMA-CC
ablation study and in our comparisons of CentralV and COMA-CC. As such, our
results should be interpreted as purely comparative, and we believe that
better absolute results are achievable in practice, with proper tuning.

We focus on $6$ maps in this section: $3$ easy maps `2s3z`, `3s5z`, `1c3s5z`,
$2$ hard maps `bane_vs_bane`, `2c_vs_64zg`, and $1$ super-hard map `MMM2`. The
easy maps were chosen because of the surprisingly low performance of CentralV
and COMA reported previously in the literature Samvelyan et al. (2019);
Papoudakis et al. (2020); Rashid et al. (2020b); Wang et al. (2020).666Some of
these papers use the older version of StarCraft II. However, our experiments
indicate that the low performance is not due to version number alone. The hard
maps include a challenging symmetric and asymmetric scenario, and the super-
hard map tests the limits of our method. For each map, we report the median
test win-rate across all seeds, and we shade the area between the 25th and
75th percentile. For aggregate results, we plot the averaged median win-rate
across all seeds and all $14$ SMAC maps for the given algorithm. Finally, when
we refer to the entire SMAC challenge, we mean the $14$ maps explicitly
mentioned in the original paper Samvelyan et al. (2019); in particular, we do
not consider the additional maps distributed with the package. The algorithms
that use semi-on-policy training are suffixed with ”-SOP”.

Figure 2. Averaged median test win rate for policy gradient algorithms and
QMIX across all 14 SMAC scenarios.

1 plot with title All Maps Median. X-axis labelled T (mil). Y-axis labelled
Averaged Median Test Win %. Legend contains CentralV-SOP, COMA-CC-SOP, QMIX,
DOP.

Figure 3. Comparing CentralV with SOP training and with purely on-policy
training. Batch size $8$.

6 plots with titles 2s3z, 3s5z, 1c3s5z, bane_vs_bane, 2c_vs_64zg, MMM2. X-axis
labelled T (mil). Y-axis labelled Median Test Win %. Legend contains CentralV-
SOP and CentralV.

##  6\. Results

Figure 4. Comparing COMA-CC with SOP training and with purely on-policy
training. Batch size $8$.

6 plots with titles 2s3z, 3s5z, 1c3s5z, bane_vs_bane, 2c_vs_64zg, MMM2. X-axis
labelled T (mil). Y-axis labelled Median Test Win %. Legend contains COMA-CC-
SOP and COMA-CC.

Figure 1 shows that CentralV-SOP and COMA-CC-SOP perform on par with the
state-of-the-art value-based method QMIX on various SMAC maps, and even
outperform it on some maps. In particular, QMIX fails to converge on
`bane_vs_bane`, whereas all MAPG methods manage to solve the scenario.
Additionally, CentralV-SOP achieves the same performance as QMIX on
`2c_vs_64zg`, outperforming all other MAPG methods and converging much faster.
The best performing MAPG method is CentralV-SOP, which outperforms other MAPG
methods on most maps. Furthermore, COMA-CC-SOP significantly outperforms the
original COMA. Results for all 14 SMAC maps are included in Appendix C.

Figure 2 shows the aggregate results on the entire SMAC challenge. QMIX
performs the best on average. However, as previously discussed, it is not the
best on all maps. CentralV-SOP and DOP achieve approximately the same
performance at the end of training, with CentralV-SOP converging much faster.
One significant difference between CentralV and COMA or DOP is that the former
learns a centralised $V$-function, while the latter learn a centralised
$Q$-function, which could be more difficult to learn in the multi-agent
setting due to the exponentially sized joint-action space.

Figures 3 and 4 show the results comparing SOP training to on-policy training
for CentralV and COMA-CC respectively. We see that the algorithms trained with
SOP data consistently outperform the corresponding on-policy versions,
demonstrating better sample efficiency and stable learning.

##  7\. Analysis

###  7.1. SOP Training

In this section we provide more insight into when SOP training works and how
to best utilise it. Figure 5 shows a plot of the maximum KL divergence between
the current policy $\pi$ and any other policy in the buffer for various maps.
To estimate the KL divergence, we use an unbiased estimator Schulman (2020):

$KL(p\;\|\;q)=\mathbb{E}_{x\sim
p}\left[\frac{q(x)}{p(x)}-1-\log\frac{q(x)}{p(x)}\right]$ (12)

As before, we plot the median across seeds and shade the area between the 25th
and 75th percentile. In the map 3s5z we can see that the divergence is
extremely low: less than $4$ in most cases, and going up to at most $10$ for
the 75th percentile. This fits our intuitive description that we can use past
data, as long as the current policy has not diverged too much from the policy
that generated it. However, even divergences as high as $150$ can still yield
stable training, as the results for 2c_vs_64zg indicate. The reason is that
the action space in 2c_vs_64zg is much larger than most other maps, and there
are many actions with a similar effect; thus, while the KL divergence may
rise, the underlying behaviour remains approximately the same.

In practice we find that sudden big changes to the current policy lead to
instability when using SOP training. There are different ways to remedy this:
for example modifying critic training or advantage estimation, lowering the
learning rate, or using strict SOP training (Algorithm 2) with an appropriate
kl_threshold. We hypothesise that, in most cases, permissive SOP training
(Algorithm 1) would be preferable, since, if the policy changes are unstable,
strict SOP training would prune too much of the old data and essentially
degenerate to on-policy training. The KL estimator is also an important factor
to consider when determining appropriate values of kl_threshold. Additionally,
our results from Figure 5 indicate that, if used as a hyperparameter, it
should be tuned on a per-application basis (i.e., per map in the case of
SMAC). However, our results from Section 6 indicate that permissive SOP
training suffices for SMAC.

Figure 5. Maximum estimated KL divergence between the current policy $\pi$ and
any other policy $\pi^{\prime}$ in the buffer. Median over all seeds. BS
indicates batch size.

3 plots with titles 3s5z, 2c_vs_64zg, MMM2. X-axis labelled T (mil). Y-axis
labelled Max KL Divergence. Legend contains CentralV-SOP BS8 and CentralV-SOP
BS32.

###  7.2. COMA Ablations

Figure 6. Averaged median test win rate for COMA-CC ablations across all 14
SMAC scenarios. BS indicates batch size.

1 plot with title All Maps Median. X-axis labelled T (mil). Y-axis labelled
Averaged Median Test Win %. Legend contains COMA-CC BS8, COMA-SOP BS32, COMA-
WholeBatch-SOP BS32, COMA-CC-SOP BS8, COMA-CC-SOP BS8, COMA-CC-SOP BS32, COMA
BS32, COMA BS8.

We now analyse the effect of each change we made to the COMA algorithm. Figure
6 shows the averaged median test win rate across all $14$ SMAC scenarios. We
use the same naming convention as before, with the addition of the COMA-
WholeBatch algorithm, which incorporates new training from Section 4.2.2, but
does not use the new critic proposed in Section 4.2.1.

We find that COMA-CC-SOP outperforms all ablations. Training on the entire
batch at once and using SOP training are crucial to achieving the best results
with COMA. We find that it is important to have both additions, as COMA-SOP
(just doing SOP training) and COMA-CC (just doing whole-batch training),
produce far weaker results compared to the ablations that incorporate both
changes. Fixing the critic representation also yields some improvements,
though they are less significant than we expected, as COMA-WholeBatch-SOP
performs similarly to COMA-CC-SOP. Additionally, we do not find significant
differences between using a batch size of $32$ episodes, compared to using a
batch size of $8$ episodes. Overall, COMA-CC-SOP achieves a more than $200\%$
higher average median test win-rate compared to the original COMA on the SMAC
benchmark. Additionally, it exhibits more stable learning.

##  8\. Conclusion

In this paper we proposed a simple to implement enhancement that yielded
significant improvements to the performance and sample efficiency of both
multi-agent policy gradient algorithms we examined: namely, using semi-on-
policy training. With this enhancement, we established CentralV-SOP as a
strong benchmark for multi-agent policy gradient algorithms on SMAC, showing
that it is competitive with state-of-the-art value-based and policy gradient
methods. The fact that CentralV-SOP often outperforms COMA and DOP, two
methods that address the multi-agent credit assignment problem, indicates that
credit assignment might be an insignificant factor on SMAC.

## References

  * (1)
  * Cao et al. (2012) Yongcan Cao, Wenwu Yu, Wei Ren, and Guanrong Chen. 2012\.  An Overview of Recent Progress in the Study of Distributed Multi-Agent Coordination.  _IEEE Transactions on Industrial Informatics_ 9 (07 2012).  https://doi.org/10.1109/TII.2012.2219061
  * Foerster et al. (2017) Jakob Foerster, Nantas Nardelli, Gregory Farquhar, Triantafyllos Afouras, Philip H. S. Torr, Pushmeet Kohli, and Shimon Whiteson. 2017\.  Stabilising Experience Replay for Deep Multi-Agent Reinforcement Learning. In _Proceedings of the 34th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 70)_ , Doina Precup and Yee Whye Teh (Eds.). PMLR, 1146–1155.  http://proceedings.mlr.press/v70/foerster17b.html
  * Foerster et al. (2018) Jakob N. Foerster, Gregory Farquhar, Triantafyllos Afouras, Nantas Nardelli, and Shimon Whiteson. 2018\.  Counterfactual multi-agent policy gradients. In _Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence_. AAAI Press, 2974–2982.
  * Haarnoja et al. (2018) Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine. 2018.  Soft Actor-Critic: Off-Policy Maximum Entropy Deep Reinforcement Learning with a Stochastic Actor. In _International Conference on Machine Learning_. 1856–1865.
  * Harutyunyan et al. (2016) Anna Harutyunyan, Marc G. Bellemare, Tom Stepleton, and Rémi Munos. 2016.  Q($$\lambda$$) with Off-Policy Corrections.  _CoRR_ abs/1602.04951 (2016).  arXiv:1602.04951
  * Iqbal and Sha (2019) Shariq Iqbal and Fei Sha. 2019.  Actor-Attention-Critic for Multi-Agent Reinforcement Learning. In _Proceedings of the 36th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 97)_ , Kamalika Chaudhuri and Ruslan Salakhutdinov (Eds.). PMLR, Long Beach, California, USA, 2961–2970.
  * Kraemer and Banerjee (2016) Landon Kraemer and Bikramjit Banerjee. 2016.  Multi-agent reinforcement learning as a rehearsal for decentralized planning.  _Neurocomputing_ 190 (2016), 82–94.
  * Lillicrap et al. (2019) Timothy P. Lillicrap, Jonathan J. Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra. 2019.  Continuous control with deep reinforcement learning.  arXiv:1509.02971 [cs.LG]
  * Lowe et al. (2017) Ryan Lowe, Yi Wu, Aviv Tamar, Jean Harb, Pieter Abbeel, and Igor Mordatch. 2017\.  Multi-Agent Actor-Critic for Mixed Cooperative-Competitive Environments.  _Neural Information Processing Systems (NIPS)_ (2017), 6382–6393.
  * Mahajan et al. (2019) Anuj Mahajan, Tabish Rashid, Mikayel Samvelyan, and Shimon Whiteson. 2019.  MAVEN: Multi-Agent Variational Exploration. In _Advances in Neural Information Processing Systems_. 7611–7622.
  * Mnih et al. (2016) Volodymyr Mnih, Adrià Puigdomènech Badia, Mehdi Mirza, Alex Graves, Tim Harley, Timothy P. Lillicrap, David Silver, and Koray Kavukcuoglu. 2016.  Asynchronous Methods for Deep Reinforcement Learning. In _Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48_ (New York, NY, USA) _(ICML’16)_. JMLR.org, 1928–1937.
  * Mnih et al. (2015) Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A. Rusu, Joel Veness, Marc G. Bellemare, Alex Graves, Martin Riedmiller, Andreas K. Fidjeland, Georg Ostrovski, Stig Petersen, Charles Beattie, Amir Sadik, Ioannis Antonoglou, Helen King, Dharshan Kumaran, Daan Wierstra, Shane Legg, and Demis Hassabis. 2015.  Human-level control through deep reinforcement learning.  _Nature_ 518, 7540 (01 Feb 2015), 529–533.  https://doi.org/10.1038/nature14236
  * Munos et al. (2016a) Remi Munos, Tom Stepleton, Anna Harutyunyan, and Marc Bellemare. 2016a.  Safe and Efficient Off-Policy Reinforcement Learning. In _Advances in Neural Information Processing Systems_ , D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (Eds.), Vol. 29. Curran Associates, Inc., 1054–1062.
  * Munos et al. (2016b) Remi Munos, Tom Stepleton, Anna Harutyunyan, and Marc Bellemare. 2016b.  Safe and Efficient Off-Policy Reinforcement Learning. In _Advances in Neural Information Processing Systems_ , D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (Eds.), Vol. 29. Curran Associates, Inc.
  * Oliehoek and Amato (2016) Frans A. Oliehoek and Christopher Amato. 2016.  _A Concise Introduction to Decentralized POMDPs_ (1st ed.).  Springer Publishing Company, Incorporated.
  * Oliehoek et al. (2008) Frans A Oliehoek, Matthijs TJ Spaan, and Nikos Vlassis. 2008\.  Optimal and approximate Q-value functions for decentralized POMDPs.  _Journal of Artificial Intelligence Research_ 32 (2008), 289–353.
  * Papoudakis et al. (2020) Georgios Papoudakis, Filippos Christianos, Lukas Schäfer, and Stefano V. Albrecht. 2020\.  Comparative Evaluation of Multi-Agent Deep Reinforcement Learning Algorithms.  arXiv:2006.07869 [cs.LG]
  * Precup et al. (2000) Doina Precup, Richard S. Sutton, and Satinder P. Singh. 2000\.  Eligibility Traces for Off-Policy Policy Evaluation. In _Proceedings of the Seventeenth International Conference on Machine Learning_ _(ICML ’00)_. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 759–766.
  * Rashid et al. (2020a) Tabish Rashid, Gregory Farquhar, Bei Peng, and Shimon Whiteson. 2020a.  Weighted QMIX: Expanding Monotonic Value Function Factorisation for Deep Multi-Agent Reinforcement Learning. In _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (Eds.), Vol. 33. Curran Associates, Inc., 10199–10210.
  * Rashid et al. (2020b) Tabish Rashid, Mikayel Samvelyan, Christian Schroeder de Witt, Gregory Farquhar, Jakob Foerster, and Shimon Whiteson. 2020b.  Monotonic Value Function Factorisation for Deep Multi-Agent Reinforcement Learning.  _Journal of Machine Learning Research_ 21, 178 (2020), 1–51.
  * Samvelyan et al. (2019) Mikayel Samvelyan, Tabish Rashid, Christian Schroeder de Witt, Gregory Farquhar, Nantas Nardelli, Tim G. J. Rudner, Chia-Man Hung, Philip H. S. Torr, Jakob Foerster, and Shimon Whiteson. 2019.  The StarCraft Multi-Agent Challenge. In _Proceedings of the 18th International Conference on Autonomous Agents and MultiAgent Systems_ (Montreal QC, Canada) _(AAMAS ’19)_. International Foundation for Autonomous Agents and Multiagent Systems, Richland, SC, 2186–2188.
  * Schulman (2020) John Schulman. 2020\.  _Approximating KL Divergence_.  http://joschu.net/blog/kl-approx.html
  * Schulman et al. (2017) John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. 2017.  Proximal Policy Optimization Algorithms.  arXiv:1707.06347 [cs.LG]
  * Son et al. (2020) Kyunghwan Son, Sungsoo Ahn, Roben Delos Reyes, Jinwoo Shin, and Yung Yi. 2020.  QTRAN++: Improved Value Transformation for Cooperative Multi-Agent Reinforcement Learning.  arXiv:2006.12010 [cs.LG]
  * Son et al. (2019) Kyunghwan Son, Daewoo Kim, Wan Ju Kang, David Earl Hostallero, and Yung Yi. 2019.  QTRAN: Learning to Factorize with Transformation for Cooperative Multi-Agent Reinforcement Learning. In _Proceedings of the 36th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 97)_ , Kamalika Chaudhuri and Ruslan Salakhutdinov (Eds.). PMLR, 5887–5896.
  * Sunehag et al. (2018) Peter Sunehag, Guy Lever, Audrunas Gruslys, Wojciech Marian Czarnecki, Vinicius Zambaldi, Max Jaderberg, Marc Lanctot, Nicolas Sonnerat, Joel Z. Leibo, Karl Tuyls, and Thore Graepel. 2018.  Value-Decomposition Networks For Cooperative Multi-Agent Learning Based On Team Reward. In _Proceedings of the 17th International Conference on Autonomous Agents and MultiAgent Systems_ (Stockholm, Sweden) _(AAMAS ’18)_. International Foundation for Autonomous Agents and Multiagent Systems, Richland, SC, 2085–2087.
  * Sutton and Barto (1998) Richard S. Sutton and Andrew G. Barto. 1998.  _Introduction to Reinforcement Learning_ (1st ed.).  MIT Press, Cambridge, MA, USA.
  * Tampuu et al. (2017) Ardi Tampuu, Tambet Matiisen, Dorian Kodelja, Ilya Kuzovkin, Kristjan Korjus, Juhan Aru, Jaan Aru, and Raul Vicente. 2017\.  Multiagent cooperation and competition with deep reinforcement learning.  _PLOS ONE_ 12, 4 (04 2017), 1–15.  https://doi.org/10.1371/journal.pone.0172395
  * Tan (1993) Ming Tan. 1993.  Multi-Agent Reinforcement Learning: Independent vs. Cooperative Agents. In _In Proceedings of the Tenth International Conference on Machine Learning_. Morgan Kaufmann, 330–337.
  * Wang et al. (2020) Yihan Wang, Beining Han, Tonghan Wang, Heng Dong, and Chongjie Zhang. 2020.  Off-Policy Multi-Agent Decomposed Policy Gradients.  arXiv:2007.12322 [cs.LG]
  * Wang et al. (2016) Ziyu Wang, Victor Bapst, Nicolas Heess, Volodymyr Mnih, Rémi Munos, Koray Kavukcuoglu, and Nando de Freitas. 2016.  Sample Efficient Actor-Critic with Experience Replay.  _CoRR_ abs/1611.01224 (2016).  arXiv:1611.01224
  * Watkins (1989) Christopher Watkins. 1989\.  _Learning From Delayed Rewards_.  Ph.D. Dissertation. King’s College, Cambridge.
  * Yang et al. (2020) Yaodong Yang, Jianye Hao, Ben Liao, Kun Shao, Guangyong Chen, Wulong Liu, and Hongyao Tang. 2020.  Qatten: A General Framework for Cooperative Multiagent Reinforcement Learning.  arXiv:2002.03939 [cs.MA]
  * Ye et al. (2015) Dayong Ye, Minjie Zhang, and Yun Yang. 2015.  A multi-agent framework for packet routing in wireless sensor networks.  _Sensors (Basel, Switzerland)_ 15, 5 (28 Apr 2015), 10026–10047.  https://doi.org/10.3390/s150510026 25928063[pmid].
  * Zhou et al. (2020) Meng Zhou, Ziyu Liu, Pengwei Sui, Yixuan Li, and Yuk Ying Chung. 2020.  Learning Implicit Credit Assignment for Cooperative Multi-Agent Reinforcement Learning.  arXiv:2007.02529 [cs.LG]

##  Appendix A Experimental setup

All environment configurations for SMAC are kept at their default values. We
use the same network architecture and hyperparameters as in Samvelyan et al.
(2019), unless stated otherwise. The configuration is the same for both
CentralV and COMA-CC. We present it here for completeness.

###  A.1. Network architecture

The policy network is a recurrent network which contains a GRU of hidden
dimension $64$ with a fully-connected linear layer before and after. To speed
up learning, parameters are shared across agents. The critic architecture for
both algorithms is a fully-connected feedforward neural network, with $2$
hidden layers of $128$ units each, and a final layer of $1$ unit; as discussed
in Section 4.2.1, the original COMA critic has a final layer of $|U|$ units.
On each training iteration, we first update the critic, and then the actor,
performing a single gradient descent step in both cases; the original
implementations of these algorithms used minibatch gradient descent to update
their critics, as in Algorithm 3 from Section 4.2.2. Both networks apply a
ReLU non-linearity between all but the final layer.

###  A.2. Hyperparameters

For both CentralV and COMA-CC (including ablations) we use the following
hyperparameters: for TD($\lambda$) we set $\lambda=0.8$ for all maps, with
discount factor $\gamma=0.99$. Both algorithms update the target critic every
$200$ training iterations. The actor and critic are optimised separately using
an RMSProp optimiser; in both cases we set the optimiser parameters as
follows: learning rate $=0.005$, $\alpha=0.99$, and $\varepsilon=1$e-$5$. To
enforce exploration, both algorithms use an $\epsilon$-floor scheme (as in
Foerster et al. (2018)) for the agents’ policies, with $\epsilon$ linearly
annealed from $0.5$ to $0.01$ over 100k time steps. As in Foerster et al.
(2018); Samvelyan et al. (2019), at test time, the algorithms perform greedy
action selection, choosing the actions with highest probability.

##  Appendix B Computational complexity of the COMA-CC critic

Multiplying a $n\times k$ matrix with a $k\times m$ matrix requires $O(n\times
k\times m)$ multiplications. The COMA critic is a neural network, which takes
$i_{1}$ inputs, has two hidden layers of size $h$, and produces $m$ outputs.
The COMA-CC critic takes $i_{2}$ inputs, has two hidden layers of size $h$,
and produces a single output. Thus, on a single input, the COMA critic
requires $O(hi_{1})+O(h^{2})+O(hm)$ multiplications to compute its outputs and
the COMA-CC critic requires $O(hi_{2})+O(h^{2})+O(h)$ multiplications.

Note that $h$ is typically a fixed constant. Furthermore, note that
$m=O(\min\left\\{i_{1},i_{2}\right\\})$, since the joint action is included in
the input. Thus, we get that the COMA critic requires $O(i_{1})$
multiplications and the COMA-CC critic requires $O(i_{2})$ multiplications.
For each baseline computation, the COMA critic takes $n$ inputs, one for each
agent, so requires $O(n\times i_{1})$ multiplications; the COMA-CC critic
takes $nm$ inputs, one for each agent and counterfactual joint action, so
requires $O(n\times m\times i_{2})$ multiplications. Therefore, assuming
$i_{1}=\Theta(i_{2})$, the COMA-CC critic network performs $O(m)$ times as
many multiplications as the COMA critic.

##  Appendix C Full results

In this section, we present our full set of results.

###  C.1. Actor-critic algorithms vs QMIX

The per-map comparisons of all algorithms we examine are presented in Figure
9.

###  C.2. CentralV-default

Here we show the results for using SOP training on CentralV without using the
$\gamma=1$ trick when computing advantages. We refer to this version as
CentralV-default. The results are shown in Figure 7. SOP training still
manages to improve the performance of the algorithm. The increased training
instability can be explained by unstable policy updates, as shown in Figure 8.

Figure 7. Comparing CentralV-default with SOP training and with purely on-
policy training. Batch size $8$.

6 plots with titles 2s3z, 3s5z, 1c3s5z, bane_vs_bane, 2c_vs_64zg, MMM2. X-axis
labelled T (mil). Y-axis labelled Median Test Win %. Legend contains CentralV-
default On-Policy Training and CentralV-default SOP Training.

Figure 8. Comparing KL divergence of CentralV and CentralV-default with SOP
training. Batch size $8$.

6 plots with titles 2s3z, 3s5z, 1c3s5z, bane_vs_bane, 2c_vs_64zg, MMM2. X-axis
labelled T (mil). Y-axis labelled Max KL Divergence. Legend contains CentralV-
SOP and CentralV.

Figure 9. Comparing MAPG methods with QMIX on the entire SMAC challenge. Batch
size $32$.

Series of plots comparing CentralV-SOP, COMA-CC-SOP, QMIX, and DOP.


# The Stellar-mass Function of Long Gamma-Ray Burst Host Galaxies

Guang-Xuan Lan  Purple Mountain Observatory, Chinese Academy of Sciences,
Nanjing 210023, China<EMAIL_ADDRESS><EMAIL_ADDRESS>School of Astronomy
and Space Science, University of Science and Technology of China, Hefei
230026, China  Guangxi Key Laboratory for Relativistic Astrophysics, Nanning
530004, China  Department of Astronomy, School of Physics, Peking University,
Beijing 100871, China  Jun-Jie Wei  Purple Mountain Observatory, Chinese
Academy of Sciences, Nanjing 210023, China<EMAIL_ADDRESS><EMAIL_ADDRESS>School of Astronomy and Space Science, University of Science and Technology of
China, Hefei 230026, China  Guangxi Key Laboratory for Relativistic
Astrophysics, Nanning 530004, China  Ye Li  Purple Mountain Observatory,
Chinese Academy of Sciences, Nanjing 210023, China<EMAIL_ADDRESS><EMAIL_ADDRESS>Hou-Dun Zeng  Purple Mountain Observatory, Chinese Academy of
Sciences, Nanjing 210023, China<EMAIL_ADDRESS><EMAIL_ADDRESS>Key
Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory,
Chinese Academy of Sciences Nanjing 210034, China<EMAIL_ADDRESS>Xue-Feng Wu
Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210023,
China<EMAIL_ADDRESS><EMAIL_ADDRESS>School of Astronomy and Space Science,
University of Science and Technology of China, Hefei 230026, China

###### Abstract

Long gamma-ray bursts (GRBs) have been discussed as a potential tool to probe
the cosmic star formation rate (SFR) for a long time. Some studies found an
enhancement in the GRB rate relative to the galaxy-inferred SFR at high
redshifts, which indicates that GRBs may not be good tracers of star
formation. However, in these studies, the GRB rate measured at any redshift is
an average over all galaxies at that epoch. A deep understanding of the
connection between GRB production and environment also needs to characterize
the population of GRB host galaxies directly. Based on a complete sample of
GRB hosts, we constrain the stellar-mass function (SMF) of GRB hosts, and
examine redshift evolution in the GRB host population. Our results confirm
that a strong redshift evolution in energy (with an evolution index of
$\delta=2.47^{+0.73}_{-0.89}$) or in density ($\delta=1.82^{+0.22}_{-0.59}$)
is needed in order to account for the observations. The GRB host SMF can be
well described by the Schechter function with a power-law index
$\xi\approx-1.10$ and a break mass $M_{b,0}\approx 4.9\times 10^{10}$ ${\rm
M}_{\odot}$, independent of the assumed evolutionary effects. This is the
first formulation of the GRB host SMF. The observed discrepancy between the
GRB rate and the galaxy-inferred SFR may also be explained by an evolving SMF.

Gamma-ray bursts (629) — Galaxies (573) — Star formation (1569) — Stellar mass
functions (1612)

##  1 Introduction

As the extremely powerful explosions of massive stars (Woosley & Bloom, 2006),
long-duration gamma-ray bursts (GRBs) are detectable out to very high
redshifts (Salvaterra et al., 2009b; Tanvir et al., 2009; Cucchiara et al.,
2011). The redshift distribution of long GRBs and the properties of their host
galaxies can therefore be used to probe the evolution of the star formation
rate (SFR) density with cosmic time (e.g., Totani 1997; Wijers et al. 1998;
Lamb & Reichart 2000; Porciani & Madau 2001). Relevant studies in this field
have been carried out for more than two decades, but it remains an open
question whether long GRBs can be good tracers of cosmic star formation. By
studying the GRB redshift distribution, there is a general agreement about the
fact that the comoving GRB rate does not purely follow the cosmic SFR, but
exhibits an enhancement of GRB detections at high redshifts ($z\gtrsim 2-3$)
with respect to the expectation from the star formation history (e.g., Daigne
et al. 2006; Guetta & Piran 2007; Le & Dermer 2007; Kistler et al. 2008; Li
2008; Qin et al. 2010; Wanderman & Piran 2010; Virgili et al. 2011; Lu et al.
2012; Robertson & Ellis 2012; Wang 2013; Lien et al. 2014; Wei et al. 2014;
Perley et al. 2015, 2016c; Vergani et al. 2015; Japelj et al. 2016; Wei & Wu
2017; Palmerio et al. 2019; Lan et al. 2019, 2021; Ghirlanda & Salvaterra
2022). Several possible mechanisms have been proposed to reconcile the
discrepancy between the GRB rate and the SFR, such as cosmic metallicity
evolution (Langer & Norman, 2006; Li, 2008; Salvaterra et al., 2012; Lloyd-
Ronning et al., 2019; Ghirlanda & Salvaterra, 2022), an evolving initial mass
function of stars (Xu & Wei, 2009; Wang & Dai, 2011), and an evolving
luminosity function of GRBs (e.g., Lloyd-Ronning et al. 2002; Firmani et al.
2004; Salvaterra & Chincarini 2007; Salvaterra et al. 2009a, 2012; Cao et al.
2011; Tan & Wang 2015; Pescalli et al. 2016; Lan et al. 2019, 2021).

Among the possible mechanisms, cosmic metallicity evolution appears to be a
compelling one. The progenitors of long GRBs are expected to be preferentially
localized in low-metallicity environments, where the loss of mass and angular
momentum are lessened, allowing a relativistic jet to form (MacFadyen &
Woosley, 1999; Hirschi et al., 2005; Yoon & Langer, 2005; Woosley & Heger,
2006; Yoon et al., 2006). And indeed it has been shown that long GRBs tend to
occur in lower metallicity (as well as lower stellar-mass) galaxies, with
metallicity $Z\leq 0.3Z_{\odot}$ (see Perley et al. 2016a and references
therein). However, some works suggested that GRBs could occur in high-
metallicity host environments (e.g., Levesque et al. 2010; Elliott et al.
2012, 2013; Savaglio et al. 2012; Hao & Yuan 2013). It is important to note
that previous studies of the GRB redshift distribution are necessarily
imprecise. Specifically, their derived GRB rate at any redshift is an average
over all galaxies at that epoch (a diverse population covering orders of
magnitude in metallicities, SFRs, etc.) and cannot disclose which types of
galaxies contribute most to the GRB rate. To thoroughly understand the link
between GRB production and environment, one also is required to characterize
the population of GRB host galaxies directly (Perley et al., 2016c).

One of the difficulties in attempting to connect GRB observations and their
environments to particular progenitor systems is obtaining a complete sample
of GRB host galaxies. Recently, Perley et al. (2016b, c) carried out a new
survey of the GRB host galaxy population. Their survey provides the first host
galaxy sample that is unbiased and sufficiently large to statistically examine
redshift evolution in the host galaxy population in detail. In this work, we
make use of the complete sample presented in Perley et al. (2016b, c) to
quantify the stellar-mass function (SMF) of GRB host galaxies. Meanwhile, we
use the observed redshift distribution of the complete sample to probe and
constrain the evolution of the GRB host population in redshift.

The structure of the paper is as follows: In Section 2, we introduce the GRB
host sample used for our analysis. We describe our analysis method in Section
3, while the model results are presented in Section 4. Finally, in Section 5
we draw our conclusions. Throughout this work we adopt a flat $\Lambda$CDM
model with cosmological parameters $H_{0}=70$ km s-1 Mpc-1, $\Omega_{\rm
m}=0.3$, and $\Omega_{\Lambda}=0.7$.

##  2 The Sample

The redshifts of about 40% of all GRBs detected by the Swift satellite have
been measured (Lan et al., 2021). Though this is a significant improvement
compared to the pre-Swift era, from the point of view of characterizing the
GRB host population or redshift distribution, the sample is still far from
being considered unbiased and complete. The Swift GRB Host Galaxy Legacy
Survey (SHOALS) built a new unbiased sample of long GRB host galaxies, which
is composed of 119 GRBs at $0.03<z<6.29$ drawn from the Swift GRB catalog
using a series of selection criteria and requiring the Swift/Burst Alert
Telescope (BAT) fluence to exceed $S>10^{-6}$ erg $\rm cm^{-2}$ (Perley et
al., 2016b). Their optimized observability cuts combined with host galaxy
follow-up enabled redshift measurements for 110 of these bursts (representing
a completeness of 92%). Subsequently, Perley et al. (2016c) presented rest-
frame near-infrared (NIR) luminosities and stellar masses for this uniform
sample of 119 GRB host galaxies through the Spitzer Infrared Array Camera
(IRAC) observations. As the first large, highly complete, and multiband survey
of an unbiased GRB host galaxy sample, SHOALS provides unprecedented insight
into galaxy evolution and cosmic history (Perley et al., 2016b, c).

In this work, we use the SHOALS sample to investigate the stellar-mass
distribution of GRB hosts across cosmic history and its possible implications
for interpreting the observed GRB redshift distribution. We calculate the
isotropic-equivalent energies $E_{\rm iso}$ in the 1–$10^{4}$ keV rest-frame
energy range for these 110 GRBs having redshifts. Note that there is a low-
luminosity/energy burst (GRB 060218) with $E_{\rm iso}\leq 10^{49}$ erg. As
the low-luminosity/energy bursts have been claimed to be a distinct population
(Soderberg et al., 2004; Cobb et al., 2006; Chapman et al., 2007; Liang et
al., 2007), GRB 060218 is discarded in our analysis. The remaining sample
includes 109 GRBs, whose stellar mass–redshift distribution is presented in
the left panel of Figure 1. It is also noteworthy that 30 out of these 109 GRB
host galaxies only have stellar-mass upper limits (Perley et al., 2016c). In
principle, the galaxies with stellar-mass upper limits should be removed in
the determination of the stellar-mass distribution of GRB hosts. However, the
direct abandonment of these hosts would introduce a new bias on shaping the
GRB redshift distribution. In the following, we will simply treat the upper
mass limits as the stellar-mass measurements for these 30 GRB hosts.

Figure 1: Left panel: the stellar mass–redshift distribution of 109 GRB host
galaxies in the SHOALS sample. The inverted triangles represent the hosts with
stellar-mass upper limits. The dotted line represents the stellar-mass
threshold adopted in our calculations (corresponding to an apparent-magnitude
limit $m_{\rm obs}=25.5$ mag; see text and Equations (9) and (10)). Right
panel: stellar mass as a function of the absolute magnitude at a wavelength of
$3.6/(1+z)$ $\mu$m for 79 GRB host galaxies with measured stellar mass and
redshift. The best-fitting curve is parameterized in Equation (10).

##  3 Analysis Method

In order to infer the model-free parameters, a maximum likelihood method
firstly proposed by Marshall et al. (1983) is adopted. The likelihood function
$\mathcal{L}$ is defined as (e.g., Chiang & Mukherjee 1998; Narumoto & Totani
2006; Ajello et al. 2009, 2012; Abdo et al. 2010; Zeng et al. 2014, 2016; Lan
et al. 2019, 2021)

$\mathcal{L}=\exp\left(-N_{\rm exp}\right)\prod^{N_{\rm
obs}}_{i=1}\Phi(E_{{\rm iso},i},\;M_{\ast,i},\;z_{i},\;t_{i})\;,$ (1)

where $N_{\rm obs}$ is the number of the observed sample, $N_{\rm exp}$ is the
expected number of GRB detections, and $\Phi(E_{\rm iso},\;M_{\ast},\;z,\;t)$
is the observed rate of GRBs per unit time at redshift $z\sim z+\mathrm{d}z$
with energy $E_{\rm iso}\sim E_{\rm iso}+\mathrm{d}E_{\rm iso}$ and stellar
mass $M_{\ast}\sim M_{\ast}+\mathrm{d}M_{\ast}$, which can be expressed as

$\displaystyle\Phi\left(E_{\rm iso},\;M_{\ast},\;z,\;t\right)$
$\displaystyle=$
$\displaystyle\frac{\mathrm{d}^{4}N}{\mathrm{d}t\mathrm{d}z\mathrm{d}M_{\ast}\mathrm{d}E_{\rm
iso}}=\frac{\mathrm{d}^{4}N}{\mathrm{d}t\mathrm{d}V\mathrm{d}M_{\ast}\mathrm{d}E_{\rm
iso}}\times\frac{\mathrm{d}V}{\mathrm{d}z}$ (2) $\displaystyle=$
$\displaystyle\frac{\Delta\Omega}{4\pi}\frac{\psi(z)}{1+z}\varphi(M_{\ast},z)\phi(E_{\rm
iso})\times\frac{\mathrm{d}V}{\mathrm{d}z}\;,$

where $\Delta\Omega=1.33$ sr is the $Swift$/BAT field of view, $\psi(z)$ is
the comoving formation rate of GRBs (in units of $\rm Mpc^{-3}$ $\rm
yr^{-1}$), $(1+z)^{-1}$ accounts for the cosmological time dilation,
$\phi(E_{\rm iso})$ is the normalized isotropic-equivalent energy function of
GRBs, and $\varphi(M_{\ast},z)$ is the normalized SMF of GRB hosts, which
evolves with redshift dependents on the assumed model (see more below). ${\rm
d}V(z)/{\rm d}z=4\pi cD_{L}^{2}(z)/[H_{0}(1+z)^{2}\sqrt{\Omega_{\rm
m}(1+z)^{3}+\Omega_{\Lambda}}]$ is the comoving volume element in a flat
$\Lambda$CDM model, where $D_{L}(z)$ is the luminosity distance.

Since in the collapsar model the births of GRBs just mean the deaths of short-
lived massive stars, the GRB formation rate $\psi(z)$ should in principle be
related to the cosmic SFR $\psi_{\star}(z)$, i.e.,

$\psi(z)=\eta\psi_{\star}(z)\;,$ (3)

where $\eta$ is the GRB formation efficiency in units of ${\rm
M}_{\odot}^{-1}$. The SFR (in units of ${\rm M}_{\odot}$ $\rm yr^{-1}$ $\rm
Mpc^{-3}$) can be parameterized as (Hopkins & Beacom, 2006; Li, 2008)

$\psi_{\star}(z)=\frac{0.0157+0.118z}{1+(z/3.23)^{4.66}}\;.$ (4)

For the GRB energy function $\phi(E_{\rm iso})$, we employ a broken power-law
form

$\phi(E_{\rm iso})=\frac{A}{\ln(10)E_{\rm
iso}}\left\\{\begin{array}[]{l}{\left(\frac{{{E_{\rm
iso}}}}{{{E_{b}}}}\right)^{a}};\,\,{E_{\rm iso}}\leq{E_{b}}\\\
{\left(\frac{{{E_{\rm iso}}}}{{{E_{b}}}}\right)^{b}};\,\,{E_{\rm
iso}}>{E_{b}}\;,\end{array}\right.$ (5)

where $A$ is a normalization factor, and $a$ and $b$ are the power-law indices
before and after the break energy $E_{b}$.

Empirically, the galaxy SMF can be well described by a Schechter form (Drory &
Alvarez, 2008). Therefore, we use the same Schechter expression for the
stellar-mass distribution of GRB host galaxies:

$\varphi(M_{\ast})=\varphi_{0}\left(\frac{M_{\ast}}{M_{b,0}}\right)^{\xi}\exp\left(-\frac{M_{\ast}}{M_{b,0}}\right)\frac{1}{M_{b,0}}\;,$
(6)

where $\varphi_{0}$ is a normalization factor, $\xi$ is the power-law index,
and $M_{b,0}$ is the break mass.

Considering the fluence threshold that is used to define the SHOALS sample
(i.e., $S_{\rm lim}=10^{-6}$ erg $\rm cm^{-2}$ in the 15–150 keV), the
expected number of GRB detections should be

$\displaystyle N_{\rm exp}=$ $\displaystyle\frac{\Delta\Omega
T}{4\pi}\int^{z_{\rm
max}}_{0}\frac{\psi(z)}{1+z}\frac{\mathrm{d}V(z)}{\mathrm{d}z}\mathrm{d}z\int^{M_{\ast,\rm
max}}_{M_{\ast,\rm lim}(z)}\varphi(M_{\ast},z)\mathrm{d}M_{\ast}$ (7)
$\displaystyle\times\int^{E_{\rm max}}_{{\rm max}[E_{\rm min},E_{\rm
lim}(z)]}\phi(E_{\rm iso})\mathrm{d}E_{\rm iso}\;,$

where $T\sim 7$ yr is the observational period of Swift that covers the SHOALS
sample and $M_{\ast,\rm lim}(z)$ represents the stellar-mass threshold of
Spitzer-IRAC at redshift $z$. Since $z<10$ and $M_{\ast}<10^{12}$ ${\rm
M}_{\odot}$ for the current GRB host population, we adopt a maximum redshift
$z_{\rm max}=10$ and a maximum stellar mass $M_{\ast,\rm max}=10^{12}$ ${\rm
M}_{\odot}$. The energy function is assumed to extend between minimum and
maximum energies $E_{\rm min}=10^{51}$ erg and $E_{\rm max}=10^{56}$ erg. The
energy threshold appearing in Equation (7) can be computed by

$E_{\rm lim}(z)=\frac{4\pi D_{L}^{2}(z)S_{\rm
lim}}{1+z}\frac{\int^{10^{4}/(1+z)\;{\rm keV}}_{1/(1+z)\;{\rm
keV}}EN(E)\mathrm{d}E}{\int^{150\;{\rm keV}}_{15\;{\rm
keV}}EN(E)\mathrm{d}E}\;,$ (8)

where $N(E)$ is the GRB photon spectrum. To describe the burst spectrum, we
employ a typical Band function with low- and high-energy spectral indices $-1$
and $-2.3$, respectively (Band et al., 1993; Preece et al., 2000; Kaneko et
al., 2006). For a given isotropic-equivalent energy $E_{\rm iso}$, the peak
energy of the spectrum $E_{p}$ is estimated through the empirical
$E_{p}$–$E_{\rm iso}$ correlation (Amati et al., 2002; Nava et al., 2012),
i.e., $\log\left[E_{p}(1+z)\right]=-29.60+0.61\log E_{\rm iso}$.

For each GRB host galaxy (with measured redshift) from the SHOALS sample, the
absolute magnitude at a wavelength of $3.6/(1+z)$ $\mu$m can be calculated by

$M_{3.6/(1+z)}=m_{\rm obs}-\mu(z)+2.5\log(1+z)\;,$ (9)

where $m_{\rm obs}$ is the observed IRAC aperture magnitude and $\mu(z)$ is
the distance modulus at $z$. It has been suggested that the quantity
$M_{3.6/(1+z)}$ alone can be used as a reasonable stellar-mass proxy (Perley
et al., 2016c). As shown in the right panel of Figure 1, there is a strong
correlation between the absolute magnitudes $M_{3.6/(1+z)}$ and the
logarithmic stellar masses $\log M_{\ast}$ of GRB host galaxies. For
simplicity, here we use a second-order polynomial to model the $\log
M_{\ast}$–$M_{3.6/(1+z)}$ relation

$\log M_{\ast}=p_{0}+p_{1}M_{3.6/(1+z)}+p_{2}M_{3.6/(1+z)}^{2}\;,$ (10)

where the best-fitting parameters are $p_{0}=3.59_{-1.13}^{+1.13}$ ,
$p_{1}=-0.11_{-0.11}^{+0.11}$, and $p_{2}=0.0087_{-0.0025}^{+0.0025}$ at 68%
confidence level. The observed magnitude $m_{\rm obs}$ can then be translated
to stellar-masses $M_{\ast}$ at different $z$ via Equations (9) and (10). We
show that the stellar-mass threshold $M_{\ast,\rm lim}$ appearing in Equation
(7) may be approximated well using a magnitude limit of Spitzer-IRAC, $m_{\rm
obs,lim}=25.5$ mag, according to Equations (9) and (10) (see the dotted line
in the left panel of Figure 1).

Figure 2: Distributions of redshift (top panel), isotropic-equivalent energy
(middle panel), and host galaxy mass (bottom panel) of 109 GRBs in the SHOALS
sample (data points, with the number of detections in each corresponding bin
indicated by a solid point with Possion error bars). Different curves show the
expected distributions from different best-fitting models: no evolution model
(purple dashed lines), energy evolution model (blue dotted-dashed lines),
density evolution model (green dotted lines), and mass evolution model (black
solid lines).

Considering that (i) GRBs may not be good tracers of star formation, (ii) the
GRB energy function may evolve with redshift, and (iii) there is a strong
evolution in the galaxy SMF (Fontana et al., 2006; Drory & Alvarez, 2008;
Duncan et al., 2014), here we explore four different evolution scenarios. In
the first no evolution model, the GRB formation rate strictly follows the
cosmic SFR, and neither the GRB energy function nor the GRB host SMF evolves
with redshift. In the energy evolution model, while the GRB rate is still
proportional to the SFR and the GRB host SMF is still taken to be nonevolving,
the break energy in the GRB energy function increases with redshift as
$E_{b,z}=E_{b,0}(1+z)^{\delta}$. In the density evolution model, the GRB rate
follows the SFR in conjunction with an extra evolving factor, i.e.,
$\psi(z)=\eta\psi_{\star}(z)(1+z)^{\delta}$, and both the GRB energy function
and the GRB host SMF are taken to be nonevolving. In the last mass evolution
model, while the GRB rate is still fully consistent with the SFR and the break
energy in the GRB energy function is still a constant, the GRB host SMF
evolves with redshift as the following relations (Drory & Alvarez, 2008):

$\displaystyle\varphi(M_{\ast},z)$ $\displaystyle=$
$\displaystyle\varphi_{*}(z)\left(\frac{M_{\ast}}{M_{b}(z)}\right)^{\xi}\exp\left(-\frac{M_{\ast}}{M_{b}(z)}\right)\frac{1}{M_{b}(z)}$
$\displaystyle\varphi_{*}(z)$ $\displaystyle=$
$\displaystyle\varphi_{0}(1+z)^{\delta}\;,$ (11) $\displaystyle\log M_{b}(z)$
$\displaystyle=$ $\displaystyle\log M_{b,0}+\gamma\ln(1+z)$

where $\delta$ and $\gamma$ are two additional evolution parameters. When
$\delta\rightarrow 0$ and $\gamma\rightarrow 0$, Equation (3) reduces to
Equation (6).

We optimize the model-free parameters by maximizing the likelihood function
(Equation (1)). Since there are several parameters in our models, the Markov
Chain Monte Carlo (MCMC) technology is applied to give multidimensional
parameter constraints from the observational data. Here we employ the MCMC
code from GWMCMC,111https://github.com/grinsted/gwmcmc which is an
implementation of the Goodman & Weare 2010 Affine invariant ensemble MCMC
sampler (Goodman & Weare, 2010; Foreman-Mackey et al., 2013).

##  4 Model Results

Table 1: Best-fitting parameters in different models with the sample of 109
GRB host galaxies

Model |  | $\log\eta$
---
$a$ | $b$ | $\log E_{b}$ | $\xi$ | $\log M_{b,0}$ | Evolution parameters | $\ln\mathcal{L}$ | AIC
|  | (${\rm M}_{\odot}^{-1}$)
---
|  | (erg) |  | $({\rm M}_{\odot})$ |  |  |
No evolution | $-7.82^{+0.38}_{-0.19}$ | $-0.06^{+0.01}_{-0.22}$ | $-1.23^{+0.18}_{-1.45}$ | $53.27^{+0.70}_{-0.02}$ | $-0.96^{+0.08}_{-0.08}$ | $10.59^{+0.14}_{-0.07}$ | $\cdot\cdot\cdot$ | -115.17 | 242.34
Energy evolution | $-7.52^{+0.58}_{-0.76}$ | $-0.10^{+0.03}_{-0.47}$ | $-1.07^{+0.39}_{-0.55}$ | $51.91^{+1.68}_{-0.27}$ | $-1.09^{+0.21}_{-0.10}$ | $10.68^{+0.23}_{-0.20}$ | $\delta=2.47^{+0.73}_{-0.89}$ | -106.56 | 227.12
Density evolution | $-7.70^{+0.97}_{-0.29}$ | $-0.46^{+0.34}_{-0.04}$ | $-1.30^{+0.55}_{-0.94}$ | $53.52^{+0.55}_{-0.67}$ | $-1.10^{+0.06}_{-0.14}$ | $10.69^{+0.23}_{-0.08}$ | $\delta=1.82^{+0.22}_{-0.59}$ | -104.08 | 222.16
Mass evolution | $-7.49^{+0.46}_{-0.59}$ | $-0.39^{+0.18}_{-0.12}$ | $-1.23^{+0.17}_{-1.17}$ | $53.51^{+0.54}_{-0.25}$ | $-1.14^{+0.10}_{-0.06}$ | $11.24^{+0.25}_{-0.43}$ |  $\delta=1.96^{+0.44}_{-0.47}$, $\gamma=-0.38^{+0.28}_{-0.19}$ | -103.35 | 222.70

_Note._

Errors show the 68% containment regions around the best-fitting values.

Using the above analysis method, we can infer the values of each model’s free
parameters, including the GRB formation efficiency $\eta$, the GRB energy
function, the host galaxy SMF, and the evolution parameters. The best-fitting
parameter values together with their $1\sigma$ uncertainties for different
models are listed in Table 1. The last two columns report the natural-
logarithmic-likelihood value $\ln\mathcal{L}$ and the Akaike information
criterion (AIC) score, respectively. The AIC score of each fitted model is
evaluated as (Akaike, 1974; Liddle, 2007)

${\rm AIC}=-2\ln\mathcal{L}+2f\;,$ (12)

where $f$ is the number of model-free parameters. If there are two or more
models for the same data, $\mathcal{M}_{\rm 1},\mathcal{M}_{\rm
2},...,\mathcal{M}_{\rm N}$, and they have been separately fitted, the one
with the smallest AIC score is the one most preferred by this criterion. With
the${\rm AIC}_{i}$ characterizing model $\mathcal{M}_{i}$, the unnormalized
confidence that this model is true is the Akaike weight $\exp(-{\rm
AIC}_{i}/2)$. The relative probability of $\mathcal{M}_{i}$ being the correct
model is then

$P(\mathcal{M}_{i})=\frac{\exp(-{\rm AIC}_{i}/2)}{\exp(-{{\rm
AIC_{1}}}/2)+\cdot\cdot\cdot+\exp(-{\rm AIC_{N}}/2)}\;.$ (13)

###  4.1 No evolution model

In the no evolution scenario, long GRBs are assumed to trace the cosmic star
formation and both the GRB energy function and the GRB host SMF are taken to
be nonevolving. Figure 2 shows the observed distributions of redshift $z$,
isotropic-equivalent energy $E_{\rm iso}$, and host galaxy mass $M_{\ast}$ of
109 GRBs in the SHOALS sample. The results of our fitting from the no
evolution model are indicated with purple dashed lines. Although the fits of
the $E_{\rm iso}$ and $M_{\ast}$ distributions look good, the expectation from
this model does not provide a good representation of the observed $z$
distribution of the SHOALS sample. In particular, the expected $z$
distribution is found to peak at a lower redshift than observed, and the rate
of GRBs at high-$z$ is underestimated. According to the AIC model selection
criterion, the no evolution model can be safely excluded as having a
probability of only $\sim 10^{-5}$ of being correct compared to the other
three models (more fully described below).

###  4.2 Energy evolution model

Evolution in the GRB energy function can lead to an enhancement of the
high-$z$ GRB detection, providing a feasible way to reconcile model results
with the observations. In this model, we consider the possibility that the
break energy is an increasing function of the reshift, which is
$E_{b,z}=E_{b,0}(1+z)^{\delta}$. Note that here the GRB rate is still
proportional to the SFR and the GRB host SMF is assumed to be nonevolving. We
find that a strong energy evolution with $\delta=2.47^{+0.73}_{-0.89}$ is
required to reproduce the observed distributions (blue dotted-dashed lines in
Figure 2). Using the AIC model selection criterion, we find that among the
evolution models, this one is somewhat disfavoured statistically with a
relative probability of $\sim 5\%$.

###  4.3 Density evolution model

An increase in the GRB formation rate with redshift can also enhance the
number of GRB detections at high-$z$. Here, we assume that the GRB rate
follows the SFR in conjunction with an extra evolving factor characterized by
$(1+z)^{\delta}$, i.e., $\psi(z)=\eta\psi_{\star}(z)(1+z)^{\delta}$. Note that
in this model both the GRB energy function and the GRB host SMF are taken to
be nonevolving. We find that a strong density evolution with
$\delta=1.82^{+0.22}_{-0.59}$ reproduces the observed distributions (green
dotted lines in Figure 2) quite well. The AIC shows that the density evolution
model is statistically preferred with a relative probability of $\sim 54\%$.

###  4.4 Mass evolution model

Here we investigate whether an evolving SMF of GRB hosts can reconcile the
discrepancy between the GRB rate and the SFR. In this model, while the GRB
formation rate still follows the cosmic SFR and the GRB energy function also
does not evolve with redshift, the Schechter parameters of the GRB host SMF
are assumed to be evolving with redshift as
$\varphi_{*}(z)=\varphi_{0}(1+z)^{\delta}$ and $\log M_{b}(z)=\log
M_{b,0}+\gamma\ln(1+z)$. The fitting results from the mass evolution model are
shown as the black solid lines in Figure 2, which agree with the observed
distributions very well. That is, a strong redshift evolution in the GRB host
SMF (with evolution indices of $\delta=1.96^{+0.44}_{-0.47}$ and
$\gamma=-0.38^{+0.28}_{-0.19}$) is required to reproduce the observations
well. We find that $\delta=0$ is ruled out at $\sim 4.2\sigma$ level and
$\gamma=0$ is ruled out at $\sim 1.4\sigma$ level, confirming the significant
dependence of Schechter parameters on redshift. Using the galaxy stellar mass
data spanning $0<z<5$ from the FORS Deep Field, Drory & Alvarez (2008) also
found a strong evolution in mass, in good agreement with our results from the
GRB host galaxies. According to the AIC, the mass evolution model is slightly
disfavoured compared to the density evolution model, but the differences are
statistically insignificant ($\sim 41\%$ for the former versus $\sim 54\%$ for
the latter).

It should be underlined that the $(1+z)^{\delta}$ term in the GRB host SMF
(see Equation (3)) can be referred to as the GRB formation efficiency $\eta$
instead of the normalization $\varphi_{*}(z)$. In other words, this term can
be taken out from the mass integral in Equation (3) as it does not depend on
the stellar mass. However, in addition to this evolutionary term, a modest
redshift evolution in the break mass of the SMF (with an evolution index of
$\gamma=-0.38^{+0.28}_{-0.19}$) is identified by the observations, serving as
an important characteristic for distinguishing between the density evolution
model and the mass evolution model.

Table 2: Best-fitting parameters in Different Models Using the Refined Sample
of 109 GRB hosts (Including 79 hosts with known stellar masses, 26 with mock
stellar masses, and 4 with stellar-mass upper limits)

Model |  | $\log\eta$
---
$a$ | $b$ | $\log E_{b}$ | $\xi$ | $\log M_{b,0}$ | Evolution parameters | $\ln\mathcal{L}$ | AIC
|  | (${\rm M}_{\odot}^{-1}$)
---
|  | (erg) |  | $({\rm M}_{\odot})$ |  |  |
No evolution | $-7.63^{+0.48}_{-0.27}$ | $-0.14^{+0.11}_{-0.11}$ | $-1.38^{+0.29}_{-0.62}$ | $53.44^{+0.24}_{-0.20}$ | $-1.02^{+0.09}_{-0.08}$ | $10.67^{+0.13}_{-0.12}$ | $\cdot\cdot\cdot$ | -114.46 | 240.92
Energy evolution | $-7.30^{+1.01}_{-0.44}$ | $-0.36^{+0.30}_{-0.20}$ | $-1.06^{+0.29}_{-0.30}$ | $52.12^{+1.03}_{-0.62}$ | $-1.16^{+0.12}_{-0.12}$ | $10.78^{+0.30}_{-0.14}$ | $\delta=2.72^{+0.80}_{-0.63}$ | -105.06 | 224.12
Density evolution | $-7.00^{+0.88}_{-0.64}$ | $-0.45^{+0.15}_{-0.10}$ | $-2.26^{+0.99}_{-1.31}$ | $54.03^{+0.22}_{-0.56}$ | $-1.22^{+0.12}_{-0.11}$ | $10.79^{+0.27}_{-0.16}$ | $\delta=1.71^{+0.48}_{-0.33}$ | -102.15 | 218.30
Mass evolution | $-6.86^{+0.57}_{-0.92}$ | $-0.24^{+0.14}_{-0.15}$ | $-1.47^{+0.46}_{-0.71}$ | $53.48^{+0.49}_{-0.32}$ | $-1.22^{+0.11}_{-0.12}$ | $11.38^{+0.43}_{-0.49}$ |  $\delta=2.03^{+0.67}_{-0.67}$, $\gamma=-0.45^{+0.36}_{-0.37}$ | -100.63 | 217.26

Figure 3: Same as Figure 2, but for the refined sample of 109 GRB hosts
(including 79 hosts with known stellar masses, 26 with mock stellar masses,
and 4 with stellar-mass upper limits). Model results are shown as in Table 2.

Note that there are 30 GRB hosts with stellar-mass upper limits in our
analysis sample, and we simply treated these upper mass limits as the stellar-
mass measurements. To examine whether these upper mass limits can reliably be
included for our analysis, we perform Monte Carlo simulations on their stellar
masses. For each GRB host with measured $z$ and stellar-mass upper limit
$M_{\ast,\rm upper}$, we mock the stellar mass $M_{\ast,\rm mock}$ according
to the mass probability distribution of 79 hosts having measured stellar
masses. With the measured $z$, we infer the stellar-mass threshold
$M_{\ast,\rm lim}$ using Equations (9) and (10). If $M_{\ast,\rm
lim}(z)<M_{\ast,\rm mock}<M_{\ast,\rm upper}$, we pick up this simulated event
as a triggered GRB host with mass detection. Since there are four GRB hosts
(GRBs 060522, 060927, 081029, and 091208B) with $M_{\ast,\rm upper}\leq
M_{\ast,\rm lim}(z)$, their upper mass limits are still treated as the
stellar-mass measurements. While the stellar masses of the remaining 26 are
generated through Monte Carlo simulations, the simulation process is repeated
1000 times to ensure the final constraints are unbiased.

Following the same analysis method described in Section 3, we also use the
refined sample of 109 GRB hosts (including 79 hosts with measured $M_{\ast}$,
26 with mock $M_{\ast,\rm mock}$, and 4 with $M_{\ast,\rm upper}$) to
construct the host galaxy SMF in different evolution models. The best-fitting
results are displayed in Table 2 and Figure 3. We find that the expectations
from our evolution models give a good description of the observed
distributions of the refined sample, whereas the no evolution model provides a
poor match for the observed redshift distribution. According to the AIC, we
can safely discard the no evolution model as having a probability of $\sim
10^{-6}$ of being correct compared to the other three models. This further
verifies the reliability of our analysis method, strengthening our
conclusions. That is, the inclusion of 30 GRB hosts with stellar-mass upper
limits would not change the main conclusions of our paper.

##  5 Summary

Whether long GRBs are biased tracers of the star formation activity remains an
open question. Studies of the GRB redshift distribution seem to indicate that
the comoving GRB rate density exhibits similar behavior to the star formation
history at low redshifts ($z\lesssim 2$), but shows an enhancement of GRBs at
high redshifts ($z\gtrsim 2-3$) compared to what would be expected if the GRB
rate followed the cosmic SFR strictly (e.g., Daigne et al. 2006; Kistler et
al. 2008; Perley et al. 2015). However, in the technique of star formation
models, the GRB rate measured at any epoch is an average over all galaxies at
that time (a diverse galaxy population spanning orders of magnitude in
metallicities, SFRs, etc.). Thus, it is difficult to inspect which types of
galaxies contribute most to the GRB formation rate. To get a better
understanding of the connection between GRB production and environment, one
has to characterize the population of GRB host galaxies directly (Perley et
al., 2016c).

Based on the complete GRB host sample presented in Perley et al. (2016b, c),
we construct the SMF of GRB host galaxies as well as the isotropic-equivalent
energy function and redshift distribution of GRBs in the frameworks of
different evolution models. The maximum likelihood method is applied to
perform a joint analysis for the distributions of redshift, isotropic-
equivalent energy, and host galaxy mass of the GRB sample, and the MCMC
algorithm is adopted to provide the best-fitting parameters in each model. Our
results confirm that GRBs must have experienced some kind of redshift
evolution, being more luminous or more population in the past than present
day. In order to match the observed distributions, the typical GRB energy
should increase to $(1+z)^{2.47^{+0.73}_{-0.89}}$ or the GRB rate density
should rise to $(1+z)^{1.82^{+0.22}_{-0.59}}$ on top of the known cosmic
evolution of the SFR. We find that the GRB host SMF can be well described by
the Schechter function with a power-law index $\xi\approx-1.10$ and a break
mass $M_{b,0}\approx 4.9\times 10^{10}$ ${\rm M}_{\odot}$, independent of the
assumed evolutionary effects.

However, since the SMF of all galaxies is suggested to be evolving with
redshift, the intrinsic form of the stellar-mass distribution of GRB host
galaxies may also have a redshift dependence. We also explore models in which
the GRB host SMF evolves with redshift. In this case, we find that the
redshift distribution of GRBs with measured redshift, energy, and stellar mass
can be successfully fitted by the star formation history, which involves with
an evolving SMF of GRB hosts. In order to reproduce the observed redshift
distribution well, the normalization coefficient of the SMF should increase to
$(1+z)^{1.96^{+0.44}_{-0.47}}$ and the break mass should evolve with redshift
according to $\log
M_{b}(z)=11.24^{+0.25}_{-0.43}-0.38^{+0.28}_{-0.19}\ln(1+z)$. Our results
suggest that long GRBs may be biased tracers of star formation and the
evolving SMF may also explain the discrepancy between the GRB rate and the
SFR.

The current sample size of GRB hosts is still too small to exhibit an obvious
redshift evolution in the stellar mass distribution. Ongoing large-area
surveys will be able to increase the host galaxy sample size. With more
abundant observational information in the future, the properties of the
population of GRB host galaxies and their possible evolution with redshift
will be better characterized, and determinations on the SMF of GRB hosts, as
discussed in this work, will be considerably improved.

We thank the anonymous referee for providing constructive comments that have
led to a significant improvement in the presentation of the material in this
paper. This work is partially supported by the National Natural Science
Foundation of China (grant Nos. 11725314, 12041306, and 12103089), the Key
Research Program of Frontier Sciences (grant No. ZDBS-LY-7014) of Chinese
Academy of Sciences, the Major Science and Technology Project of Qinghai
Province (2019-ZJ-A10), the China Manned Space Project (CMS-CSST-2021-B11),
the Natural Science Foundation of Jiangsu Province (grant Nos. BK20221562 and
BK20211000), and the Guangxi Key Laboratory for Relativistic Astrophysics.

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 720, 435
  * Ajello et al. (2009) Ajello, M., Costamante, L., Sambruna, R. M., et al. 2009, ApJ, 699, 603
  * Ajello et al. (2012) Ajello, M., Shaw, M. S., Romani, R. W., et al. 2012, ApJ, 751, 108
  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
  * Amati et al. (2002) Amati, L., Frontera, F., Tavani, M., et al. 2002, A&A, 390, 81
  * Band et al. (1993) Band, D., Matteson, J., Ford, L., et al. 1993, ApJ, 413, 281
  * Cao et al. (2011) Cao, X.-F., Yu, Y.-W., Cheng, K. S., & Zheng, X.-P. 2011, MNRAS, 416, 2174
  * Chapman et al. (2007) Chapman, R., Tanvir, N. R., Priddey, R. S., & Levan, A. J. 2007, MNRAS, 382, L21
  * Chiang & Mukherjee (1998) Chiang, J., & Mukherjee, R. 1998, ApJ, 496, 752
  * Cobb et al. (2006) Cobb, B. E., Bailyn, C. D., van Dokkum, P. G., & Natarajan, P. 2006, ApJ, 645, L113
  * Cucchiara et al. (2011) Cucchiara, A., Levan, A. J., Fox, D. B., et al. 2011, ApJ, 736, 7
  * Daigne et al. (2006) Daigne, F., Olive, K. A., Silk, J., Stoehr, F., & Vangioni, E. 2006, ApJ, 647, 773
  * Drory & Alvarez (2008) Drory, N., & Alvarez, M. 2008, ApJ, 680, 41
  * Duncan et al. (2014) Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960
  * Elliott et al. (2012) Elliott, J., Greiner, J., Khochfar, S., et al. 2012, A&A, 539, A113
  * Elliott et al. (2013) Elliott, J., Krühler, T., Greiner, J., et al. 2013, A&A, 556, A23
  * Firmani et al. (2004) Firmani, C., Avila-Reese, V., Ghisellini, G., & Tutukov, A. V. 2004, ApJ, 611, 1033
  * Fontana et al. (2006) Fontana, A., Salimbeni, S., Grazian, A., et al. 2006, A&A, 459, 745
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Ghirlanda & Salvaterra (2022) Ghirlanda, G., & Salvaterra, R. 2022, arXiv e-prints, arXiv:2206.06390
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Guetta & Piran (2007) Guetta, D., & Piran, T. 2007, J. Cosmology Astropart. Phys, 2007, 003
  * Hao & Yuan (2013) Hao, J.-M., & Yuan, Y.-F. 2013, ApJ, 772, 42
  * Hirschi et al. (2005) Hirschi, R., Meynet, G., & Maeder, A. 2005, A&A, 443, 581
  * Hopkins & Beacom (2006) Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142
  * Japelj et al. (2016) Japelj, J., Vergani, S. D., Salvaterra, R., et al. 2016, A&A, 590, A129
  * Kaneko et al. (2006) Kaneko, Y., Preece, R. D., Briggs, M. S., et al. 2006, ApJS, 166, 298
  * Kistler et al. (2008) Kistler, M. D., Yüksel, H., Beacom, J. F., & Stanek, K. Z. 2008, ApJ, 673, L119
  * Lamb & Reichart (2000) Lamb, D. Q., & Reichart, D. E. 2000, ApJ, 536, 1
  * Lan et al. (2021) Lan, G.-X., Wei, J.-J., Zeng, H.-D., Li, Y., & Wu, X.-F. 2021, MNRAS, 508, 52
  * Lan et al. (2019) Lan, G.-X., Zeng, H.-D., Wei, J.-J., & Wu, X.-F. 2019, MNRAS, 488, 4607
  * Langer & Norman (2006) Langer, N., & Norman, C. A. 2006, ApJ, 638, L63
  * Le & Dermer (2007) Le, T., & Dermer, C. D. 2007, ApJ, 661, 394
  * Levesque et al. (2010) Levesque, E. M., Kewley, L. J., Graham, J. F., & Fruchter, A. S. 2010, ApJ, 712, L26
  * Li (2008) Li, L.-X. 2008, MNRAS, 388, 1487
  * Liang et al. (2007) Liang, E., Zhang, B., Virgili, F., & Dai, Z. G. 2007, ApJ, 662, 1111
  * Liddle (2007) Liddle, A. R. 2007, MNRAS, 377, L74
  * Lien et al. (2014) Lien, A., Sakamoto, T., Gehrels, N., et al. 2014, ApJ, 783, 24
  * Lloyd-Ronning et al. (2019) Lloyd-Ronning, N. M., Aykutalp, A., & Johnson, J. L. 2019, MNRAS, 488, 5823
  * Lloyd-Ronning et al. (2002) Lloyd-Ronning, N. M., Fryer, C. L., & Ramirez-Ruiz, E. 2002, ApJ, 574, 554
  * Lu et al. (2012) Lu, R.-J., Wei, J.-J., Qin, S.-F., & Liang, E.-W. 2012, ApJ, 745, 168
  * MacFadyen & Woosley (1999) MacFadyen, A. I., & Woosley, S. E. 1999, ApJ, 524, 262
  * Marshall et al. (1983) Marshall, H. L., Tananbaum, H., Avni, Y., & Zamorani, G. 1983, ApJ, 269, 35
  * Narumoto & Totani (2006) Narumoto, T., & Totani, T. 2006, ApJ, 643, 81
  * Nava et al. (2012) Nava, L., Salvaterra, R., Ghirlanda, G., et al. 2012, MNRAS, 421, 1256
  * Palmerio et al. (2019) Palmerio, J. T., Vergani, S. D., Salvaterra, R., et al. 2019, A&A, 623, A26
  * Perley et al. (2016a) Perley, D. A., Niino, Y., Tanvir, N. R., Vergani, S. D., & Fynbo, J. P. U. 2016a, Space Sci. Rev., 202, 111
  * Perley et al. (2015) Perley, D. A., Perley, R. A., Hjorth, J., et al. 2015, ApJ, 801, 102
  * Perley et al. (2016b) Perley, D. A., Krühler, T., Schulze, S., et al. 2016b, ApJ, 817, 7
  * Perley et al. (2016c) Perley, D. A., Tanvir, N. R., Hjorth, J., et al. 2016c, ApJ, 817, 8
  * Pescalli et al. (2016) Pescalli, A., Ghirlanda, G., Salvaterra, R., et al. 2016, A&A, 587, A40
  * Porciani & Madau (2001) Porciani, C., & Madau, P. 2001, ApJ, 548, 522
  * Preece et al. (2000) Preece, R. D., Briggs, M. S., Mallozzi, R. S., et al. 2000, ApJS, 126, 19
  * Qin et al. (2010) Qin, S.-F., Liang, E.-W., Lu, R.-J., Wei, J.-Y., & Zhang, S.-N. 2010, MNRAS, 406, 558
  * Robertson & Ellis (2012) Robertson, B. E., & Ellis, R. S. 2012, ApJ, 744, 95
  * Salvaterra & Chincarini (2007) Salvaterra, R., & Chincarini, G. 2007, ApJ, 656, L49
  * Salvaterra et al. (2009a) Salvaterra, R., Guidorzi, C., Campana, S., Chincarini, G., & Tagliaferri, G. 2009a, MNRAS, 396, 299
  * Salvaterra et al. (2009b) Salvaterra, R., Della Valle, M., Campana, S., et al. 2009b, Nature, 461, 1258
  * Salvaterra et al. (2012) Salvaterra, R., Campana, S., Vergani, S. D., et al. 2012, ApJ, 749, 68
  * Savaglio et al. (2012) Savaglio, S., Rau, A., Greiner, J., et al. 2012, MNRAS, 420, 627
  * Soderberg et al. (2004) Soderberg, A. M., Kulkarni, S. R., Berger, E., et al. 2004, Nature, 430, 648
  * Tan & Wang (2015) Tan, W.-W., & Wang, F. Y. 2015, MNRAS, 454, 1785
  * Tanvir et al. (2009) Tanvir, N. R., Fox, D. B., Levan, A. J., et al. 2009, Nature, 461, 1254
  * Totani (1997) Totani, T. 1997, ApJ, 486, L71
  * Vergani et al. (2015) Vergani, S. D., Salvaterra, R., Japelj, J., et al. 2015, A&A, 581, A102
  * Virgili et al. (2011) Virgili, F. J., Zhang, B., Nagamine, K., & Choi, J.-H. 2011, MNRAS, 417, 3025
  * Wanderman & Piran (2010) Wanderman, D., & Piran, T. 2010, MNRAS, 406, 1944
  * Wang (2013) Wang, F. Y. 2013, A&A, 556, A90
  * Wang & Dai (2011) Wang, F. Y., & Dai, Z. G. 2011, ApJ, 727, L34
  * Wei & Wu (2017) Wei, J.-J., & Wu, X.-F. 2017, International Journal of Modern Physics D, 26, 1730002
  * Wei et al. (2014) Wei, J.-J., Wu, X.-F., Melia, F., Wei, D.-M., & Feng, L.-L. 2014, MNRAS, 439, 3329
  * Wijers et al. (1998) Wijers, R. A. M. J., Bloom, J. S., Bagla, J. S., & Natarajan, P. 1998, MNRAS, 294, L13
  * Woosley & Bloom (2006) Woosley, S. E., & Bloom, J. S. 2006, ARA&A, 44, 507
  * Woosley & Heger (2006) Woosley, S. E., & Heger, A. 2006, ApJ, 637, 914
  * Xu & Wei (2009) Xu, C.-y., & Wei, D.-m. 2009, Chinese Astron. Astrophys., 33, 151
  * Yoon & Langer (2005) Yoon, S. C., & Langer, N. 2005, A&A, 443, 643
  * Yoon et al. (2006) Yoon, S. C., Langer, N., & Norman, C. 2006, A&A, 460, 199
  * Zeng et al. (2016) Zeng, H., Melia, F., & Zhang, L. 2016, MNRAS, 462, 3094
  * Zeng et al. (2014) Zeng, H., Yan, D., & Zhang, L. 2014, MNRAS, 441, 1760