# Collective eXplainable AI: Explaining Cooperative Strategies and Agent
Contribution in Multiagent Reinforcement Learning with Shapley Values

Alexandre Heuillet§,
Corresponding author: A. Heuillet. e-mail: alexandre.heuillet@universite-
paris-saclay.fr Université Paris-Saclay, France  Fabien Couthouis§,
Ubisoft, France  Natalia Díaz-Rodríguez  University of Granada, Spain

###### Abstract

While Explainable Artificial Intelligence (XAI) is increasingly expanding more
areas of application, little has been applied to make deep Reinforcement
Learning (RL) more comprehensible. As RL becomes ubiquitous and used in
critical and general public applications, it is essential to develop methods
that make it better understood and more interpretable. This study proposes a
novel approach to explain cooperative strategies in multiagent RL using
Shapley values, a game theory concept used in XAI that successfully explains
the rationale behind decisions taken by Machine Learning algorithms. Through
testing common assumptions of this technique in two cooperation-centered
socially challenging multi-agent environments environments, this article
argues that Shapley values are a pertinent way to evaluate the contribution of
players in a cooperative multi-agent RL context. To palliate the high overhead
of this method, Shapley values are approximated using Monte Carlo sampling.
Experimental results on Multiagent Particle and Sequential Social Dilemmas
show that Shapley values succeed at estimating the contribution of each agent.
These results could have implications that go beyond games in economics,
(e.g., for non-discriminatory decision making, ethical and responsible AI-
derived decisions or policy making under fairness constraints). They also
expose how Shapley values only give general explanations about a model and
cannot explain a single run, episode nor justify precise actions taken by
agents. Future work should focus on addressing these critical aspects.

###### Index Terms:

Reinforcement Learning, Explainable Artificial Intelligence, Responsible
Artificial Intelligence, Shapley values

§§footnotetext: Equal contribution

##  I Introduction

Over the last few years, Reinforcement Learning (RL) has been a very active
research field. Many RL-related works focused on improving performance and
scaling capabilities by introducing new algorithms and optimizers [1, 2]
whereas very few tackled the issue of explainability in RL. However,
explainability in Machine Learning (ML) and deep learning is increasingly
becoming a pressing issue, as it concerns general public trust, and the
transparency of algorithms now conditions the deployment of and ML and RL, in
industry and daily life. As a consequence, Explainable Artificial Intelligence
(XAI) arises as a new field that strives to bring explainability to every ML
aspect, from linear classifiers [3] and time series predictors [4, 5] to RL
[6, 7]. Even though some works provide explanations for specific situations
[8, 9], RL models still lack a general explainability framework similar to
SHAP [10] or LIME [3]. These methods, although designed in XAI towards generic
ML predictive models, can bring a broader form of explainability to RL.

In this work, inspired by SHAP (SHapley Additive exPlanations) [10], the
possibilities offered by the mathematical framework of Shapley values [11] to
explain RL models were explored, with a focus on multi-agent cooperative
environments, often called Common Games. As a kind of collaborative RL, the
improvement of the learning process when agents interact with each other,
usually yields better results than training each agent in isolation [12].
These are challenging scenarios where all participants must cooperate in order
to achieve a common goal. In the particular case where agents must gather
resources from a common pool without being greedy, Garrett Hardin [13] defined
The Tragedy of the Commons: if an agent uses slightly more resources than it
should, this might be inconsequential. However, if every agent starts
following this logic then the consequences can become dire, with the common
pool being exhausted and no one being able to gather resources anymore. This
is why the Common Games are of special interest to study in this context. The
objective is finding ways to evaluate and understand how agents cooperate and
share resources in multiagent RL, through explainability methods such as SHAP.

As a matter of fact, this work could have implications that go beyond games,
since RL-based systems are increasingly used to solve critical problems. In
particular, studying social dilemmas and explaining the contribution of each
policy [14], agent, or model feature becomes relevant in many societal
problems. For instance, it could provide useful insights in economics (of
social structures), allocating resources or designing resilience programs
(e.g., for climate change, non-discriminatory decision making, ethical and
fair policy design, or to achieve the sustainable development goals).

Our main hypothesis is that Shapley values can be a pertinent way to explain
the contribution of an agent in multi-agent RL cooperative settings, and that
valuable insights, especially for the RL developer audience, can be derived
from this analysis. This article features experiments (see Section V)
conducted on Multiagent Particle [15] and Sequential Social Dilemma [16]
environments, and shows that Shapley values can accurately answer the
following research questions:

  * •

Can Shapley values be used to determine how much each agent contributes to the
global reward? (RQ1). If the answer to RQ1 is yes:

  * •

Does the proposed Monte Carlo based algorithm empirically offer a good
approximation of Shapley values? (RQ2)

  * •

What is the best method to replace an agent missing from the coalition (e.g.,
a random action, am action chosen randomly from another player or the “no
operation” action)? (RQ3)

Our experimental setup is brought forward as well in order to show the
limitations of the Shapley framework in order to make deep RL more
explainable. In fact, this approach cannot explain particular notions of a
multi-agent learning model, such as the contribution of a specific episode or
a specific action taken by an agent at a given point in time of its training.
This is due to the requirements inherent to training multi-agent RL, and the
practical design limitations of the Shapley framework, which are brought
upfront. In particular, as discussed in Section IV, Shapley values only yield
an average metric of each player’s contribution to the overall reward and
thus, to obtain this average contribution metric, one must compare the
cooperation of players during several games (or episodes in RL). However, due
to the frequently stochastic nature of (simulation and real RL) environments,
the non-deterministic behaviour of different agents or non-identical
conditions –non inherent to the agent’s policy at consideration–, it makes
concrete episodes, or concrete actions within an episode, not comparable.

This article presents the following contributions:

  * •

A study of the mathematical notion of Shapley values and how it is able to
provide quantitative explanations about the individual contribution of agents
in a cooperative multi-agent RL environment.

  * •

The application of an XAI global model-agnostic[17] method to explain multi-
agent cooperative RL models using Monte Carlo (MC) approximated Shapley
values.

  * •

A set of experiments that demonstrates the applicability and usefulness of
Shapley values and how they can provide insights that can enable a better
comprehension of emergent behaviours [18] in cooperative settings.

The rest of this article is structured as follows: Section II presents a short
survey on cooperative multi-agent RL and explainable RL, Section III presents
some preliminaries about RL and Shapley values, Section VI discusses the
experimental study setup and the general usage of Shapley values in a multi-
agent RL setting and, finally, Section VII presents conclusions on our carried
out experiments and gives some insights on promising lines of future work.

##  II Related Work

Cooperative multi-agent RL has been studied in different settings, e.g., the
emergence of different behaviors in the context of the Commons Tragedy game
[19]. These studies show, for instance, that certain inequity aversion
improves intertemporal social dilemmas [20]. However, these works analyze the
game from a theoretical point of view and not from the XAI angle. This article
specifically focuses on studying to which degree the most relevant factors (or
those contributing the most) in a black-box deep RL model can be explained.
Notably, this article focuses on pointing out at particular agents, episodes,
agents, actions or policies.

Recent techniques to attain eXplainable Reinforcement Learning (or XRL) [6]
can be categorized in two main families or categories: transparent methods or
post-hoc explainability methods (according to XAI taxonomies in [17]). On the
one hand, transparent algorithms include, by definition, every ML model that
is understandable by itself, such as a decision-tree. On the other hand, post-
hoc explainability includes all methods that craft an explanation of an RL
algorithm after its training, such as LIME [3], BreakDown [21] or SHAP [10].
Other studies [22, 23, 24] try to transpose Shapley values into XAI but, to
the best of our knowledge, none applies SHAP to explain the specific issues of
cooperative multi-agent RL. In fact, most XRL methods are based on transparent
algorithms [6].

Wang et al. [25] developed an approach to solve global reward games in a
multi-agent RL context by making use of Shapley values to distribute the
global reward more efficiently across all agents.Their focus was not on
explainability nor interpretability but performance.

All Shapley-based XAI methods listed above consider model features as
participants of a cooperative multiplayer game. In this article, Shapley
values are applied to a cooperative multi-agent RL context by considering
agents as players instead, an approach closer to the original game theory
method presented by Lloyd Shapley in 1953 [11].

Exact computation methods for both BreakDown and SHAP exist only for linear
regression and tree ensemble models. In more complex models, the dependence of
these methods on the number of samples or number of subsets of predictors $p$
to be used makes the approximated vs exact computation of contributions to be
different and, potentially, point in opposite directions [21]. This indicates
that these generic methods are not the universal response to XAI yet.

This article, despite some remaining issues against attaining a general XRL
framework, shows that in a cooperative multi-agent RL setting, Shapley values
can be used to accurately estimate the contributions of different agents.

##  III Explainability and Shapley Values to explain Multi-Agent RL in
Cooperative Settings

###  III-A Explainability in Cooperative RL

With the rapid growth of RL research and industrial applications (deployed,
e.g., in autonomous systems [26], or robotics [27, 28]), we have witnessed
over the past few years a need for Explainable RL. These needs have rapidly
risen, since being able to understand and justify the decisions of such models
is legally and morally necessary for their broad diffusion.

The subfield of XRL that focuses on multi-agent cooperative games has recently
gained significant attention with emerging concepts such as social
learning[29, 30] in a RL context. While studying social interactions between
entities in cooperative games is originally typical in sociology or economics
[31, 32], some AI researchers realized that they could potentially provide
explanations or improve the efficiency of their models this way.

Perolat et al. [19] sought to conduct new behavioral experiments using RL
agents in video-game-like cooperative environments instead of human subjects.
More accurately, they studied games which put the emphasis on common-pool
resources (CPR). They found that agents learn new emergent behaviors and that
some strategies can arise when some agents are excluded from the CPR. In
addition, they came up with metrics that quantify social outcomes such as
sustainability, equality, peace or efficiency for RL models.

Following the same direction, Jaques et al. [33] proposed a framework to
achieve better coordination and communication among agents by rewarding agents
on the basis of causal influence (i.e., actions leading to big changes in
other agents’ behavior). Their empirical results show that agents that choose
their actions carefully in order to influence others lead to better
coordination and thus, better global performance in socially challenging
settings where cooperation is paramount.

Exploring this aspect further, Ndousse et al. [34] analyzed the behavior of
independent RL agents in multi-agent environments and found out that model-
free agents do not use social learning. Thus, they introduce a model-based
auxiliary loss that allows agents to learn from other well-performing agents
to improve themselves. In addition to outperforming the experts, these agents
were also able to achieve better zero-shot performance than those which did
not rely on social learning when transferred to another task.

However, even if these works managed to extract useful information from
studying social interactions between agents, the literature lacks a general
framework that could automatically provide explanations about the level of
performance of each agent and their added value in the cooperative game, such
as SHAP [10] does for the features of a ML model.

###  III-B Shapley Values and RL in Cooperative Settings

Shapley values originate from game theory. They evaluate importance in terms
of contribution of each participant in a cooperative game, in order to help
split a shared payout in a fair way [11]. The concept which is key here is to
be able to form “coalitions” (or subsets) of players in order to measure the
performance of each player in every possible team situation (for instance,
“player A”, “player A and player B” or “player B and player C”).

Formally, a coalitional game $C=(N,v)$ is defined by a set $N$ of players with
$|N|=n$ (the number of players) and a function $v$, that maps a coalition of
players $S$ to a real number, corresponding to the total expected sum of
payouts the members of $S$ can obtain through cooperation:
$v:2^{N}\Rightarrow\mathbb{R}$, where $v(\varnothing)=0$ and $\varnothing$ is
the empty set. Thus, $v$ is denominated the gain function of the considered
game.

The idea is to quantify how much players cooperate in a coalition and how much
profit they from this cooperation [35]. According to the Shapley value
definition [11], the contribution added by player $i$ in a coalition $S$ in a
coalitional game $(N,v)$ is given by Eq. 1:

$\small\phi_{i}(v)=\sum_{S\subseteq
N\setminus\\{i\\}}\frac{|S|!(n-|S|-1)!}{n!}(v(S\cup\\{i\\})-v(S))$ (1)

A more interpretable equivalent formula to express the Shapley value for
player $i$, rewritten with binomial coefficients [36], is:

$\small\phi_{i}(v)=\frac{1}{n}\sum_{S\subseteq
N\setminus\\{i\\}}\binom{n-1}{|S|}^{-1}(v(S\cup\\{i\\})-v(S))$ (2)

In summary, the Shapley value of a feature (or player) is the mean marginal
contribution (the term: $v(S\cup\\{i\\})-v(S)$ in (2)) of all possible
coalitions, or average change in the prediction that the coalition already in
the room receives when the feature value joins them. It satisfies desirable
properties [11, 37] lacking in other XAI techniques: Efficiency (the sum of
the Shapley values of all players equals the shared payout of the grand
coalition), Dummy, or dummy feature [35, 37] (if a feature does not change the
predicted value, –e.g., the RL global reward– regardless of which coalition of
feature values it is added to, then its Shapley value is equal to zero),
Symmetry (the contributions of two feature values should be the same if they
contribute equally to all possible coalitions) and Linearity (the contribution
of a coalition of features should be the sum of the individual contributions
of the features that compose the coalition).

Shapley values must be interpreted as follows: it is the average contribution
of a feature to the prediction in different coalitions. Note that it is not
the difference in prediction when a feature would be removed from the model
[35]. As shown in Eq. 1, for $|N|=n$ players, the exact computation of Shapley
values for a specific participant requires computing the average of $2^{n-1}$
possible coalitions, which is computationally expensive, especially when
considering all players. Therefore, this would mean computing $n(2^{n-1})$
coalitions in total to obtain values for all players. While this value can be
reduced to $2^{n}$ coalitions (for all players, if the algorithm is optimized
to avoid computing the same coalition multiple times), the cost remains
exponential with respect to the number of players. In fact, the number of
agents in RL environments can vary greatly from a few [16, 15] to hundreds
[38] where exact Shapley values are prohibitively expensive to compute.
Furthermore, estimating Shapley values in a stochastic RL environment requires
sampling multiple episodes in order to estimate all marginal contributions,
which further worsens the computational cost.

In terms of computational efficiency, assuming that simulating a game is done
in constant time $O(1)$, finding the exact solution to this problem becomes
difficult, as the number of coalitions exponentially increases as more
features are added. Despite Shapley value computation being an NP-hard problem
[39], Shapley distributes the feature attribution fairly, i.e., allowing
contrastive explanations. For instance, it permits the comparison of a
prediction to another feature subset prediction, or a single data point.

In spite of the broad applicability of Shapley values, they have some
theoretical limitations. In particular, Shapley values are only a way to
obtain an average metric of each player’s contribution to the overall reward
(i.e., payout). In order to obtain this average metric one must compare the
cooperation of players during several games (or episodes in RL). As a
consequence, it obviously cannot explain the contribution of a concrete single
episode to the learning process nor explain one specific action taken by a
player.

##  IV Monte Carlo Approximation of Shapley Values

Since the complexity of computing Shapley values grows exponentially with
respect to the number of players (as discussed in Section III-B), in order to
keep the computation time manageable, contributions can be computed for only a
subset of all possible coalitions. The Shapley value $\phi_{i}(v)$ can be
approximated by Monte Carlo sampling in order to apply it to any type of
classification or regression model f [40] as follows:

$\small\hat{\phi_{i}}(\hat{f})=\frac{1}{M}\sum_{m=1}^{M}(\hat{f}(x^{m}_{+i})-\hat{f}(x^{m}_{-i}))\approx\phi_{i}(\hat{f})$
(3)

where $\hat{f}(x^{m}_{+i})$ is the model prediction (or gain function in game
theory) for input $x$ with a random number of feature values replaced by
feature values from a random data point, except for the respective value of
feature $i$, and $M$ is the number of marginal contributions to estimate in
order to compute the Shapley value for one feature. The value $x^{m}_{-i}$ is
almost identical to the last one, but the value $x_{i}^{m}$ is also taken from
this randomly sampled data point.

In this work, our contribution consists of adapting Eq. 3 to the multi-agent
RL setting, by replacing input features with agent actions. While in Eq. 3,
$\hat{f}$ is the function approximated by a classification or regression model
[40], $\hat{f}$ is instead considered to represent the global reward obtained
by agents from a random subset (or coalition) of all agents on a sample
episode. This leads to the following reformulation of Shapley value:

$\small\hat{\phi_{i}}^{RL}(r)=\frac{1}{M}\sum_{m=1}^{M}(r^{m}_{+i}-r^{m}_{-i})\approx\phi_{i}(r)$
(4)

where $r^{m}_{+i}$ corresponds to the global reward obtained by simulating one
sample episode with a random subset of players where player $i$ is present,
and $r^{m}_{-i}$ is the global reward obtained by simulating one episode with
the same subset than in $r^{m}_{+i}$, except that the current player $i$ has
been removed from the subset.

Three approaches to exclude players from a coalition (i.e., let the absent
player take “substitute” actions) are explored:

  1. a.

Replace: A missing agent’s actions are replaced by those of a randomly chosen
player among the trained agents which are present in the coalition (and
ideally with the same role as the missing one). This is the direct translation
from the standard application of Shapley values in ML [10] (against the
traditional use in game theory where it is often possible to completely remove
a player) since they replace missing feature values by ones randomly selected
among present (non-zero) features.

  2. b.

Random: Letting the absent player act by taking random actions.

  3. c.

NoOp: Replacing the actions of the missing agent by “noop” (no operation),
i.e., letting the agent do nothing, and not move.

The estimation of the Shapley value is repeated for each player. Thus, the
algorithm must roll out $2M$ times per player ($2Mn$ roll outs in total, where
$n$ is the total number of players in the game). At the end of the process,
one Shapley value per agent policy is obtained, indicating each player’s
average contribution to the grand coalition global reward (i.e., the reward
collectively obtained by all agents working simultaneously) on the sampled
episodes (as in Monte Carlo method, Algorithm 1). Hence the complexity of this
method only depends on $M$ ($O(M)$) as illustrated in Table II. The Shapley
value estimation [40] in Eq. 3 allows the model to conclude, for instance, “On
average, the contribution of Player 1 to the team has an impact of +0.6 on the
global reward”. This allows us to quantify and rank how relevant each player
is in terms of cooperation and contribution to the overall common goal.
Algorithm 1 describes the process to estimate Shapley values via Monte Carlo
sampling of coalitions of players:

Algorithm 1 Monte Carlo approximation of Shapley values applied to a multi-
agent RL context with shared payout

0: List: $agents$

0: Integer: $M$ (number of coalition permutations to be used)

0: List: $shapley\\_values$

1: $shapley\\_values\leftarrow empty\\_list()$

2: for $i\leftarrow 1$ to $length(agents)$ do

3: $marginal\\_contributions\leftarrow empty\\_list()$

4: for $m\leftarrow 1$ to $M$ do

5: $coal\\_i\leftarrow sample\\_coalition(agents[i])$

6: $coal\\_no\\_i\leftarrow remove\\_from\\_list(coal\\_i,agents[i])$

7: $r_{+i}\leftarrow rollout(coal\\_i)$

8: $r_{-i}\leftarrow rollout(coal\\_i)$

9: $add\\_to\\_list(marginal\\_contributions,(r_{+i}-r_{-i}))$

10: end for

11: $shapley\\_value\\_i\leftarrow mean(marginal\\_contributions)$

12: $add\\_to\\_list(shapley\\_values,shapley\\_value\\_i)$

13: end for

14: return $shapley\\_values$

In Section II and with the Shapley value estimation in Eq. 3, different
approaches to apply Shapley values as an XAI method were presented. While it
is possible to use the notion of contribution on classification or regression
model predictions, it is also possible to use it on an RL reward r as a final
payout that needs to be explained. The analogy could also be easy to
understand if we were to use a value function–as in [25]–instead of roll outs,
to evaluate each player’s contribution. However, with this generic approach of
using reward as a payout, it becomes intuitive: the more a player (i.e.,
feature in XAI) is important, the more its presence will lead to a higher
reward on average. This contribution ranking scheme–computed by making
coalitions of agents to estimate marginal contributions on the mean final
reward–can be used to rank agents in order of average importance in the team,
as explained above. This is the approach presented in this article.

##  V Experimental Study

###  V-A Context and Hypotheses

A straightforward application of Shapley values to endow an AI model with a
notion of explainability consists of considering the features of a model as
participants of a cooperative game, and the final prediction as the shared
payout [10, 41]. Computing the Shapley value of each feature provides us with
its weight in the final decision.

In this setup, our aim is to explain a deep multi-agent RL model using Shapley
values to determine the contribution of each agent to the group’s global
reward. The participants of the cooperative game will be agents, and the
shared payout the global reward obtained by the agents at the end of an
episode. By using the sampling described in Section IV to estimate the Shapley
value for each agent, it can be expected that each agent would obtain a
Shapley value proportional to its contribution to the collaborative task. It
would then be possible to answer the research questions defined in Section I
(RQ1, RQ2, RQ3).

Experiments were conducted §§§The repository linking to our experiments can be
found here: https://github.com/Fabien-Couthouis/XAI-in-RL using two multi-
agent RL environments and three different RL algorithms:

  * •

Multiagent Particle (Predator Prey scenario) [15]: A light environment with
continuous observations and a discrete action space. Predator-Prey scenario
was evaluated with three predators and a single prey. In this scenario, the
prey is faster and aims to avoid being caught by predators. The prey is
positively rewarded when escaping predators, and negatively rewarded when
caught, and vice versa for predators. Both types of agents are rewarded
negatively when trying to overpass the screen boundaries or hitting an
obstacle. As it was done by the authors of Multiagent Particle [15], the prey
was trained using the DDPG [42] algorithm §§§DDPG implementation from
repository https://github.com/openai/maddpg/ whereas predators were trained
using MADDPG [15]§§§MADDPG implementation: https://github.com/openai/maddpg/.
For this scenario, 5 different models with the same hyperparameters were
trained in order to obtain meaningful values despite the stochastic nature of
MADDPG, as suggested in [43] and detailed in the Appendix.

  * •

Sequential Social Dilemmas (Harvest scenario) [16]: Another light environment
with continuous observations and discrete actions that proposes scenarios
emphasizing social interactions and cooperation between agents. In the Harvest
scenario, agents must cooperate to harvest the maximum number of apples, while
being careful to not “kill” apple trees by collecting all apples they contain,
as this would prevent this tree from spawning apples further. As it was done
by the authors of Sequential Social Dilemmas open source implementation of the
environment [44], all agents were trained using the Asynchronous Advantage
Actor-Critic (A3C) [1]§§§A3C implementation: https://github.com/ray-
project/ray/tree/master/rllib algorithm. The best reported hyperparameters and
training protocol showcased in [44, 33] were used.

Two experiments were conducted on the Predators-Prey scenario: 1) Using
default settings provided by the authors of the environment [15] (to verify
RQ1, RQ2 and RQ3), and 2) with different speeds for each predator (to further
confirm RQ1 and RQ2). Speeds used on each experiment are presented in Table I.
For Harvest, three experiments were also conducted: first, the Shapley values
of a simple model trained with 6 agents using the default settings suggested
by [44, 33] were computed to further investigate RQ1 and, especially,
determine the minimal number of agents required for this task. Then, Shapley
values were recomputed using the same model but modified according to
information extracted from the Shapley analysis obtained during the first
experiment, in order to confirm the validity of such information. In addition,
social outcome metrics from [19] are reported to have a more fine-grained view
of the payout notion (instead of merely a global reward). The following
metrics were implemented:

  * •

Efficiency (Eq. 5): Measures the total sum of all rewards obtained by all
agents.

  * •

Equality (Eq. 6): Measures the statistical dispersion intended to represent
inequality (Gini coefficient [45]).

  * •

Sustainability (Eq. 7): Defined as the average time $t_{s}\in t$ at which
rewards are collected.

Considering $N$ independent agents, let $\\{r_{t}^{i}|t=1,...,T\\}$ be the
sequence of rewards obtained by the i-th agent over an episode of duration $T$
timesteps. Its return is given by $R^{i}=\sum_{t=1}^{T}r^{i}_{t}$. Thus, the
equations describing the social metrics are as follows:

$\small\text{Efficiency U}=\mathbb{E}[\frac{\sum_{i=1}^{N}R^{i}}{T}]$ (5)
$\small\text{Equality
E}=1-\frac{\sum_{i=1}^{N}\sum_{j=1}^{N}|R^{i}-R^{j}|}{2N\sum_{i=1}^{N}R^{i}}$
(6) $\small\text{Sustainability
S}=\mathbb{E}[\frac{1}{N}\sum_{i=1}^{N}t^{i}],\text{where
}t^{i}=\mathbb{E}[t|r_{t}^{i}>0]$ (7)

Finally, additional experiments were executed on different Harvest checkpoints
to explore how agents cooperate. Each experiment will be detailed in the
subsections below.

Table I: Speed settings for each agent used in Experiment 1 and 2 (RQ1, RQ2,
RQ3), conducted on the Predator-Prey setting of Multiagent Particle [15].

| Prey | Pred. 1 (slow) | Pred. 2 (medium) | Pred. 3 (fast)
---|---|---|---|---
Speed in exp. #1 | 1.3 | 1.0 | 1.0 | 1.0
Speed in exp. #2 | 1.3 | 0.2 | 0.8 | 2.0

Table II: Table reporting Shapley value computation times for the Harvest [44] experiments (5 agents) run on a 6-core AMD Ryzen 5 5600X CPU, using Algorithm 1. Computation time grows proportionally to $M$. M | 100 | 500 | 1000
---|---|---|---
Computation time |  $\approx 1$ hour |  $\approx 5$ hour |  $\approx 10$ hour

###  V-B Experiment 1: Agents with Identical Settings in Multiagent Particle

####  V-B1 Environment Settings

The goal of this experiment is to answer RQ1 and RQ2 (Section I), i.e.,
whether the Shapley values of predator agents correlate with the number of
times they catch the prey, and whether these are a close approximation to the
exact Shapley values.

First, when training a model of three predators using default settings on the
Predator-Prey scenario, one could think that each predator should provide a
similar contribution, as predators do not have significant differences on
their speed, action space or training method. However, in contrast to this
assumption, statistics in Figure 1 show that the performance of each predator
agent (i.e., the number of times each predator catches the prey) varies
significantly: Predator 3 has a larger contribution than Predator 1, which has
a higher contribution than Predator 0. In fact, the trained MADDPG model
developed a strategy on which Predator 3 and Predator 2 perform better than
Predator 0 at catching the prey. This can be explained by the fact that MADDPG
provides the same reward to all agents, instead of only rewarding the highest-
contributing agent. This can bias the training, as explained in [25]. On the
contrary, rewarding only the agent who contributed the most does not highlight
or recognize team strategies where the contribution of every agent was
critical to the shared payout.

Figure 1: Predator-Prey environment: Agents are endowed with the same speed.
Left plot shows predator agents mean performance comparison out of 10,000
sample episodes over 5 different models (run 2,000 sample episodes each) while
right plot presents the Monte Carlo estimation of Shapley values obtained for
each predator agent (M=1,000; “random” player exclusion method; mean over 5
models).

####  V-B2 Shapley Values Analysis

Figure 2 shows Shapley values computed for each agent on each of the five
models, over 1,000 sample episodes per model. This can be observed in a more
convenient way in Figure 1, which shows Shapley values computed for each agent
on a single model. As hypothesized, the decreasing order of agents’ Shapley
values is the following: Predator 2, Predator 1 and Predator 0. Thus, this
first experiment supports RQ1, since Shapley values are able to correctly map
contributions to agents usefulness in a cooperative multi-agent setting.

####  V-B3 Comparison of Approximated Shapley Values with the Exact Shapley
Values

In this subsection, the Monte Carlo approximation of Shapley values is
compared with the exact computation of Shapley values to verify RQ2. For that,
Eq. 2 is used to compute Shapley values: marginal contributions are estimated
as the mean global reward obtained by the coalition of players over a high
number of episodes–1,000 episodes in this experiment–using the player
exclusion mechanisms in Section IV. Indeed, a large number of samples is
needed for each coalition due to the stochastic nature of the environment,
which leads to high variance in results. As depicted in Fig. 2, Shapley values
estimated by Monte Carlo sampling are very close to the real Shapley values.
The difference of only 5% on average for all agents and models therefore
supports RQ2. In addition, MC approach is simpler to implement and with a
complexity that does not grow exponentially with the number of agents in the
RL environment.

Figure 2: Predator-Prey environment, same agents’ speeds. Comparison of MC
approximation of Shapley values with the exact Shapley values. Each point
relates to one of the 5 different models run.

###  V-C Experiment 2: Introducing Variations in Agents Speeds in Multiagent
Particle

####  V-C1 Environment Settings

This experiment introduces variations in the agent settings with the objective
of disturbing the actual distribution of contributions between agents obtained
in Experiment 1 and see if, as claimed in RQ1, this change will be reflected
in the computed Shapley Values (i.e., the aim is to ensure that Shapley values
correlate with the observed contributions). Thus, the speed of each predator
will differ from the default one and so, it can be expected that the faster
agent will catch the prey more often and contribute the most to the global
reward. In this manner, our goal is to obtain a clear hierarchy between agents
that correlates with Shapley Values so they reflect the agents’ observed
behavior coherently. The following speeds are arbitrarily set for each agent:
Predator 0 (slow): 0.2, Predator 1 (medium): 0.8, Predator 2 (fast): 2.0.
Statistics in Figure 3 show the performance of each predator agent in terms of
the number of times each predator catches the prey. As expected, the faster
agent (Predator 3) presents a higher contribution than Predator 1, which
exhibits a higher contribution than Predator 0 (the slowest). More precisely,
the ranking in speed is reflected in the ranking of contributions , which is
the ideal setting to test if the assumption in RQ1 is valid. In addition, the
real Shapley values for these settings were also computed in order to further
verify RQ2.

Figure 3: Predator-Prey environment, variable agents’ speeds. The left plot
shows predator agents’ performance comparison for 10,000 sample episodes and 5
different models (2,000 sample episodes for each of the five trained models)
while the right plot presents the Monte Carlo estimation of Shapley values
(M=1,000) for each predator agent (averaged over 5 model runs) with the
“random” player exclusion option.

####  V-C2 Shapley Values Analysis

Figure 5 shows the Shapley values computed for each agent on each of the five
models, over 1,000 sample episodes per model. This can be observed in a more
convenient way in Figure 3, which shows the Shapley values computed for each
agent on a single model. Following RQ1, the order of agents’ Shapley values is
the following: Predator 2, Predator 1, and finally Predator 0. These Shapley
values accurately correlate with the number of times each of the agents caught
the prey (see Figure 3). In addition, as noticeable in Fig. 4, when making a
single agent’s speed vary, its Shapley value grows proportionally: the faster
the agent is, the higher its Shapley value is. As expected, it makes sense,
since it can more easily catch the prey, hence contributing more to the
overall payout. Thus, this experiment further supports RQ1, since Shapley
values accurately correlate to the exact (observed) distribution of
contributions.

Figure 4: Predator-Prey environment, variable agents’ speeds. Monte Carlo
approximation of Shapley values (M=1,000) obtained by a single predator agent
with respect to its speed.

####  V-C3 Comparison of Approximated Shapley Values with Exact Shapley Values

In this section, the same experiment as described in Subsection V-B3 was
conducted to further verify RQ2. The Monte Carlo approximation of Shapley
values was compared with the exact or complete computation of Shapley values.
Here again, Fig. 5 showcases that Shapley values estimated by Monte Carlo
sampling (with $M=1,000$) are very close to the real ones with an average
difference of 8% between the approximated values and the real ones for
Predator 1 and Predator 2, when computed with the same replacement method.
However, this percentage reaches 53% when considering only Predator 0, because
its Shapley values are very close to 0 with a small standard deviation.
Therefore, a small difference in value leads to a large difference in
percentage. Therefore, approximated Shapley values are close to the real ones
while being simpler to implement and with a complexity that does not grow with
the number of agents in the RL environment, supporting RQ2 again.

Figure 5: Predator-Prey environment, variable agents’ speeds. Comparison of
the Monte Carlo approximation of Shapley values with the real Shapley values.
Each point represents one of the five model runs..

###  V-D Experiment 3: Explaining Agent Contribution in the Harvest
Environment

####  V-D1 Environment Settings

In this experiment, the Monte Carlo Shapley Value computation method is
applied to another multi-agent environment (i.e., Harvest). The goal here is
to test this method in a more complex use case where it could prove useful to
explain and extract insightful information from a trained model in which all
agents share the same settings. This data is leveraged in order to attempt to
further verify RQ1 (presented in Section I).

Default settings of Harvest were used: 6 agents trying to collect apples on a
pre-configured map. At first glance, this seems quite a large number of agents
for a map this small (39*15, with 159 apples initially). Thus, one can
hypothesize that some agents are superfluous, not contributing much to the
global reward, and may even prevent other agents from elaborating effective
strategies together, obstructing them in their movements.

####  V-D2 Shapley Values Analysis

First, Figure 6 clearly highlights the fact that Agent 5 does not seem to
bring much added value to the team (i.e., its Shapley value is close to 0)
while all other agents seem to contribute a near equal amount to the global
reward. In fact, while watching the agents play, it appears that Agent 5 is
left unaccounted for, wandering on the map, not harvesting any apple and
sometimes randomly hitting the map border, which grants it a negative reward.
This may indicate that the default setting with 6 agents is unnecessary or
unproductive: training only 5 agents could be enough to provide the same level
of performance, will require less hardware and be less time consuming.
Moreover, the fact that the first five agents contribute equally may prove
that the A3C [1] model found a satisfying solution to distribute tasks among
agents, with the exception of Agent 5 who is left unaccounted for, maybe, and
reasonably, because none of its actions could help increase the global reward.
In conclusion,this synergy between agents, observed both in the obtained
global reward and the nearly equal partition of Shapley values among agents,
shows that agents are actually cooperating with each other.

Figure 6: Harvest environment: The left plot shows the Monte Carlo estimation
of Shapley values obtained for each agent (M=1,000, “noop” action selection
method). The right plot displays Shapley values computed over agents with the
same settings but without Agent 5 present. Figure 7: Harvest environment.
Monte Carlo estimation of Shapley values (M=1,000) for each agent using each
of the three agent substitution methods. The same settings were used for all
agents.

####  V-D3 Following The Insight Given by Shapley Values

In Subsection V-D2, the Shapley values analysis indicated that Agent 5 does
not contribute at all to solve the Harvest environment. Thus, it was decided
to re-run the Shapley value computation on the same model; this time with
Agent 5 removed (i.e., completely deactivated and not appearing in the
environment map) to check, whether, as suspected, the global reward and
agents’ contribution remained the same. Figure 6 shows that Shapley values
distribution when Agent 5 is deactivated stays nearly identical to the
previous setting where Agent 5 is included in the game. Some very minor
variation in Shapley values can be attributed to stochasticity, since each
estimation of Shapley values does not yield the exact same results every time.
In addition, the global reward also remains stable at around 450 reward units.
This means that removing Agent 5 did not have any negative effect on the game
and corroborates our hypothesis that its participation was not productive. In
conclusion, a valuable insight was successfully derived from the analysis of
Shapley values, strengthening the validity of the hypothesis in RQ1.

####  V-D4 Social Outcome Metrics: Analysis of the Social Behavior Between
Agents

This section presents an analysis of how each social outcome metric
(introduced in [19]§§§The Peace metric, also presented in [19], was not
included because it relies on a specific mechanism (i.e., a time-out period
during which a tagged agent cannot harvest apples anymore) that was deemed as
of limited utility by the authors, since agents would learn very quickly not
to use it. (whose definitions are given in Subsection V-A) can be explained
with our approach. The goal here is to further explore RQ1 and see whether
agents leverage social learning and cooperate with each other in practice.
Note that the social and sustainability metrics (as originally defined in
[19]) require credit assignment to be computed in a per-agent basis, while
this is not the case for Shapley values, which only needs to be given a global
reward.

Figure 8: Harvest environment: Evolution of the Shapley values and social
metrics from [19] over different training episodes with A3C model. From top to
bottom are displayed the Shapley values (noop action selection method; MC
estimation with M=500), the mean of those Shapley values of all agents, the
efficiency, equality, and sustainability metrics.

For these experiments, over 100 training episodes were run, the three social
outcome metrics presented above were computed, as well as the Shapley values
(using Monte Carlo approximation and M=500), at different steps of the
training (i.e., every 1,000 episodes), in order to analyze the evolution of
those values during training. Results are presented in Figure 8. As the
Efficiency metric is nothing else than the mean of the Shapley values (in
expectation), it can be observed with no surprise that the mean of the Shapley
values of all agents has the same evolution than the efficiency, and directly
correlates with it. The slight differences, such as the peak at episode 4,000,
are due to the stochastic nature of the environment, the choice of $M$ in
Monte Carlo approximation and number of episodes used to compute the metrics.

From previous experiment it can be concluded that the mean of Shapley values
of all agents is a metric that explains the agent’s afficiency and does not
assume any credit assignment among agents. As stated in Equation 4, a shared
global reward is enough. Equality evolves in the same manner (with a peak at
episode 4,000 and a drop at episode 5,000). This metric decreases a bit at
episode 2,000: this fall is not captured by the efficiency nor the mean of
Shapley Values. However, while looking at the Shapley values of each agent in
Figure 8, it is clear that agent 3 obtains a lower reward than other agents.
Using the Shapley values of all agents, an explanation can thus be that the
decrease in equality at episode 2,000 is caused by agent 3, which contributes
little to the global reward (in comparison with other agents). It could not
have been possible to explain this behavior using the equality metric in a
shared global reward context and herein the value of applying Shapley analysis
in this context. As Shapley values are computed using the global reward of the
total team of agents, no link can be established between the sustainability
metric and Shapley values. Especially because sustainability is defined as the
average time (time-step in our use cases) at which rewards are collected, and
this is independent from the value of these rewards.

This experiment showed that Shapley values can effectively capture both
metrics; especially 1) if agents get high rewards (Efficiency) and 2) if
rewards are shared equally among all agents (Equality). However, Shapley
values cannot tell if agents are obtaining their reward continuously (i.e.,
displaying a sustainable behaviour). Contrarily to the social outcome metrics,
Shapley values can be computed even though the reward is globally shared among
agents. This is a huge advantage when the environment does not allow precise
credit assignment.

###  V-E Choosing a Player Exclusion Method

In this subsection, results gathered using the different player exclusion
options (defined in Section IV) during the three experiments above are
analyzed in order to answer RQ3.

When looking at Figures 2, 5 and 7, it becomes clear that the three player
exclusion mechanisms lead to Shapley values that, when considered
individually, are coherent for the excluded agent with respect to the others.
However, in the referential of a single agent, the standard deviation between
their respective values is important. In particular, there is a significant
gap between noop and the other two methods. Figure 7 shows there is an average
gap of 73.1 reward units between noop and random_player_action, while there is
only an average gap of 20.5 reward units between random and
random_player_action). This can be explained by the fact that randomly moving
agents disturb the game significantly more than immobilized ones, as they can
get negative rewards by hitting the map borders in Multiagent Particle or
killing trees in Harvest (i.e., harvesting all apples that are contiguous).
Thus, when using random or replace strategy, a majority of coalitions are
“parasited” by these negative rewards that contribute towards lowering the
global reward and lead to an overall lower Shapley value than the noop method
(as observed in Figures 2, 5 and 7). Therefore, in that context, noop action
selection seems to be the most faithful method to get Shapley values assessing
the agents’ true contributions closely, and free from random and unwanted
negative rewards.

##  VI Discussion

This article demonstrated the usefulness of Shapley values and their Monte
Carlo approximation for explaining RL models in cooperative settings. These
values provide a form of explanation, i.e., continuous values that are
understandable by researchers and developers, since they represent a portion
of the reward value of the agents team, partitioned according to each agent’s
contribution. They could also provide explanations for the general public that
may perceive them as an intrinsic “value” of each agent, making them
accountable for the effectiveness of the system. Moreover, Shapley values
could be a good way to detect biases in the training of an RL model, since
they require analyzing the individual behavior of each agent and this could
highlight disparities between their different strategies and abilities.

Concerning the player exclusion method to replace missing agents from a
coalition, noop (no-operation) action seems to be the most neutral, and
interaction-free method when the environment offers this possibility, since
methods using a substitution mechanism mandated by random-selection of actions
are prone to get high negative rewards and interfere in the game. Social
interaction between agents was also explored and this investigation showed
that Shapley values are able to effectively capture both efficiency and
equality metrics, while they are still able to be computed even though the
reward is globally shared between agents. This is a huge advantage when the
environment does not enable fair individual-level credit assignment. In
consequence, it can be asserted that Shapley Values are an effective way to
explain the contributions of RL agents, and, to some extent, the relationships
between them.

However, our approach is limited to multi-agent cooperative RL and, in its
current form, cannot be applied to competitive and single-agent models. In
addition, it cannot be used to explain an agent’s actions, their
sustainability in time, nor explain a specific episode of interest, as it only
provides an average metric for the contribution of each agent in a cooperative
game, with the total of Shapley values corresponding to the mean global reward
of the grand coalition (i.e., the one containing all agents). Thus, it must be
considered as a way to get a first ranking of contributions of agents in a
model. Finally, while the Monte Carlo method to estimate Shapley values (see
Section IV) is more efficient than computing the exact Shapley values, it
still remains time consuming. Future work should seek to keep accurate value
estimation of SHAP values while accelerating their computational
approximation.

##  VII Conclusion and Future Work

The three research questions were positively answered, with experiments
conducted in two socially challenging multi-agent RL environments (Harvest
from Sequential Social Dilemmas [44, 16] and Particle Multiagent [15]) and two
different RL algorithms (MADDPG [15] and A3C [1]). Experiments showed thatthe
computation of Shapley values could be a potential breakthrough elucidating
understanding towards attaining multi-agent XRL environments. They can
efficiently assess the contribution of agents to the global reward in
cooperative settings. They also provide insightful information about the
agents’ behaviour and their social interactions.

Nonetheless, numerous issues remain to be explained in future work. Different
interpretations of Shapley values to further explain deep RL issues must be
explored to increase the levels of explanation granularity. Robustness and
reproducibility remain a critical issue for XRL (and XAI in a more general
sense), and other statistical methods could prove very useful for that
purpose, as presented in [46, 47, 48] (e.g., Winsorised or trimmed
estimators). Moreover, Shapley values could also be combined with a robust
model selection measure (such as the Lorenz Zonoids [22]). Besides, Shapley
values or other additive and non-additive methods could be used not only to
explain the roles taken by agents when learning a policy to achieve a
collaborative task, but also to detect defects in agents while training, or in
the fed data. Furthermore, the dynamic nature of RL (vs. the static settings
of most ML models where only a single data point needs to be explained) could
be taken into account in order to create a novel approach that evaluates the
contributions of agents through time (e.g. during evaluation time). Here,
“temporal” Shapley values could be approximated with a model as in [25].
However, one of the main advantages of SHAP being a post-hoc XAI method (i.e.,
being agnostic to the RL algorithm) would be lost, as the Shapley prediction
model would be dependent on the policy learning model used. Finally, a
different contribution ranking scheme than the one presented in Section IV
could also be proposed (i.e., accounting for more complex objective metrics to
better highlight the order of importance in the team). For instance, each set
of observations per agent could be ranked in order of quality or average
didactic importance in order to assess the learning agents more fairly, with
respect to the quality of the data they were exposed to.

## References

  * [1] V. Mnih, A. P. Badia, M. Mirza, A. Graves, T. Lillicrap, T. Harley, D. Silver, and K. Kavukcuoglu, “Asynchronous Methods for Deep Reinforcement Learning,” in _International Conference on Machine Learning_. PMLR, 2016, pp. 1928–1937.
  * [2] L. Espeholt, H. Soyer, R. Munos, K. Simonyan, V. Mnih, T. Ward, Y. Doron, V. Firoiu, T. Harley, I. Dunning _et al._ , “Impala: Scalable Distributed Deep-RL with Importance Weighted Actor-Learner Architectures,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1407–1416.
  * [3] M. T. Ribeiro, S. Singh, and C. Guestrin, “‘Why Should I Trust You?’ Explaining the Predictions of any Classifier,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2016, pp. 1135–1144.
  * [4] S. El-Sappagh, J. M. Alonso, S. R. Islam, A. M. Sultan, and K. S. Kwak, “A Multilayer Multimodal Detection and Prediction Model based on Explainable Artificial Intelligence for Alzheimer’s Disease,” _Scientific Reports_ , vol. 11, no. 1, pp. 1–26, 2021.
  * [5] U. Schlegel, H. Arnout, M. El-Assady, D. Oelke, and D. A. Keim, “Towards A Rigorous Evaluation of XAI Methods On Time Series,” in _IEEE/CVF International Conference on Computer Vision Workshop (ICCVW)_ , 2019, pp. 4197–4201.
  * [6] A. Heuillet, F. Couthouis, and N. Díaz-Rodríguez, “Explainability in Deep Reinforcement Learning,” _Knowledge-Based Systems_ , vol. 214, p. 106685, 2021. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0950705120308145
  * [7] E. Puiutta and E. M. Veith, “Explainable Reinforcement Learning: A Survey,” in _International Cross-Domain Conference for Machine Learning and Knowledge Extraction_. Springer, 2020, pp. 77–95.
  * [8] P. Madumal, T. Miller, L. Sonenberg, and F. Vetere, “Explainable Reinforcement Learning through a Causal Lens,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 03, 2020, pp. 2493–2500.
  * [9] S. Greydanus, A. Koul, J. Dodge, and A. Fern, “Visualizing and Understanding Atari Agents,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1792–1801.
  * [10] S. M. Lundberg and S.-I. Lee, “A Unified Approach to Interpreting Model Predictions,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , 2017, pp. 4768–4777.
  * [11] L. Shapley, “A Value for N-Person Games,” _Contributions to the Theory of Games_ , no. 28, pp. 307–317, 1953.
  * [12] A. Andres, E. Villar-Rodriguez, A. D. Martinez, and J. Del Ser, “Collaborative exploration and reinforcement learning between heterogeneously skilled agents in environments with sparse rewards,” in _2021 International Joint Conference on Neural Networks (IJCNN)_ , 2021, pp. 1–10.
  * [13] G. Hardin, “The Tragedy of the Commons,” _Journal of Natural Resources Policy Research_ , vol. 1, no. 3, pp. 243–253, 2009.
  * [14] M. Chica, J. M. Hernández, and J. Bulchand-Gidumal, “A Collective Risk Dilemma for Tourism Restrictions under the COVID-19 Context,” _Scientific Reports_ , vol. 11, no. 1, pp. 1–12, 2021.
  * [15] R. Lowe, Y. Wu, A. Tamar, J. Harb, P. Abbeel, and I. Mordatch, “Multi-Agent Actor-Critic for Mixed Cooperative-Competitive Environments,” _Neural Information Processing Systems (NIPS)_ , 2017. [Online]. Available: https://arxiv.org/pdf/1706.02275.pdf
  * [16] J. Z. Leibo, V. Zambaldi, M. Lanctot, J. Marecki, and T. Graepel, “Multi-Agent Reinforcement Learning in Sequential Social Dilemmas,” in _Proceedings of the 16th Conference on Autonomous Agents and MultiAgent Systems_ , ser. AAMAS ’17. International Foundation for Autonomous Agents and Multiagent Systems, 2017, p. 464–473.
  * [17] A. B. Arrieta, N. Díaz-Rodríguez, J. Del Ser, A. Bennetot, S. Tabik, A. Barbado, S. García, S. Gil-López, D. Molina, R. Benjamins _et al._ , “Explainable Artificial Intelligence (XAI): Concepts, Taxonomies, Opportunities and Challenges Toward Responsible AI,” _Information Fusion_ , vol. 58, pp. 82–115, 2020.
  * [18] K. K. Ndousse, D. Eck, S. Levine, and N. Jaques, “Emergent social learning via multi-agent reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 7991–8004.
  * [19] J. Perolat, J. Z. Leibo, V. Zambaldi, C. Beattie, K. Tuyls, and T. Graepel, “A Multi-Agent Reinforcement Learning Model of Common-Pool Resource Appropriation,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , 2017, pp. 3646–3655.
  * [20] E. Hughes, J. Leibo, M. Phillips, K. Tuyls, E. Duenez-Guzman, A. Castaneda, I. Dunning, T. Zhu, K. McKee, R. Koster _et al._ , “Inequity Aversion Improves Cooperation in Intertemporal Social Dilemmas,” in _Advances in Neural Information Processing Systems 31_ , vol. 31. Neural Information Processing Systems Foundation, Inc., 2018, pp. 1–11.
  * [21] M. Staniak and P. Biecek, “Explanations of Model Predictions with live and breakDown Packages,” _The R Journal_ , vol. 10, no. 2, pp. 395–409, 2018\. [Online]. Available: 10.32614/RJ-2018-072
  * [22] P. Giudici and E. Raffinetti, “Shapley-Lorenz Decompositions in eXplainable Artificial Intelligence,” _SSRN Electronic Journal_ , 01 2020.
  * [23] D. Shim, Z. Mai, J. Jeong, S. Sanner, H. Kim, and J. Jang, “Online Class-Incremental Continual Learning with Adversarial Shapley Value,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 11, 2021, pp. 9630–9638.
  * [24] M. Sundararajan and A. Najmi, “The Many Shapley Values for Model Explanation,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 9269–9278.
  * [25] J. Wang, Y. Zhang, T.-K. Kim, and Y. Gu, “Shapley Q-value: A Local Reward Approach to Solve Global Reward Games,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 05, 2020, pp. 7285–7292.
  * [26] B. R. Kiran, I. Sobh, V. Talpaert, P. Mannion, A. A. Al Sallab, S. Yogamani, and P. Pérez, “Deep Reinforcement Learning for Autonomous Driving: A Survey,” _IEEE Transactions on Intelligent Transportation Systems_ , 2021\.
  * [27] H. Nguyen and H. La, “Review of Deep Reinforcement Learning for Robot Manipulation,” in _Third IEEE International Conference on Robotic Computing (IRC)_ , 2019, pp. 590–595.
  * [28] T. Lesort, V. Lomonaco, A. Stoian, D. Maltoni, D. Filliat, and N. Díaz-Rodríguez, “Continual Learning for Robotics: Definition, Framework, Learning Strategies, Opportunities and Challenges,” _Information Fusion_ , vol. 58, pp. 52–68, 2020.
  * [29] D. Lee, N. Jaques, J. C. Kew, D. Eck, D. Schuurmans, and A. Faust, “Joint attention for multi-agent coordination and social learning,” _CoRR_ , vol. abs/2104.07750, 2021. [Online]. Available: https://arxiv.org/abs/2104.07750
  * [30] K. Ndousse, D. Eck, S. Levine, and N. Jaques, “Multi-agent social reinforcement learning improves generalization,” _arXiv e-prints_ , pp. arXiv–2010, 2020.
  * [31] J. Duffy and J. Ochs, “Cooperative Behavior and the Frequency of Social Interaction,” _Games and Economic Behavior_ , vol. 66, no. 2, pp. 785 – 812, 2009, special Section In Honor of David Gale. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0899825608001395
  * [32] A. M. Colman, “Cooperation, Psychological Game Theory, and Limitations of Rationality in Social Interaction,” _The Behavioral and Brain Sciences_ , vol. 26(2), p. 139–198, 2003.
  * [33] N. Jaques, A. Lazaridou, E. Hughes, C. Gulcehre, P. Ortega, D. Strouse, J. Z. Leibo, and N. De Freitas, “Social Influence as Intrinsic Motivation for multi-agent deep reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 3040–3049.
  * [34] K. Ndousse, D. Eck, S. Levine, and N. Jaques, “Learning Social Learning,” in _NeurIPS Workshop on Cooperative AI_ , 2020. [Online]. Available: https://arxiv.org/abs/2010.00581
  * [35] C. Molnar, _Interpretable Machine Learning_ , 2019, https://christophm.github.io/interpretable-ml-book/.
  * [36] _The Shapley Value: Essays in Honor of Lloyd S. Shapley_. Cambridge University Press, 1988. [Online]. Available: http://www.library.fa.ru/files/Roth2.pdf
  * [37] E. Friedman and H. Moulin, “Three Methods to Share Joint Costs or Surplus,” _Journal of Economic Theory_ , vol. 87, no. 2, pp. 275 – 312, 1999. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0022053199925346
  * [38] T. Chu, S. Qu, and J. Wang, “Large-Scale Multi-Agent Reinforcement Learning using Image-Based State Representation,” in _IEEE 55th Conference on Decision and Control (CDC)_ , 2016, pp. 7592–7597.
  * [39] U. Faigle and W. Kern, _The Shapley Value for Cooperative Games under Precedence Constraints_ , ser. Memorandum. University of Twente, Faculty of Mathematical Sciences, 1992, no. 1025.
  * [40] E. Štrumbelj and I. Kononenko, “Explaining Prediction Models and Individual Predictions with Feature Contributions,” _Knowledge and Information Systems_ , vol. 41, no. 3, pp. 647–665, 2014. [Online]. Available: https://moodle.telekom.ftn.uns.ac.rs/pluginfile.php/13342/mod_folder/content/0/Feature%20importance%20paper.pdf?forcedownload=1
  * [41] A. Tallón-Ballesteros and C. Chen, “Explainable AI: Using Shapley Value to Explain Complex Anomaly Detection ML-Based Systems,” _Machine Learning and Artificial Intelligence: Proceedings of MLIS 2020_ , vol. 332, p. 152, 2020.
  * [42] T. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous Control with Deep Reinforcement Learning,” _CoRR_ , vol. abs/1509.02971, 2016.
  * [43] P. Henderson, R. Islam, P. Bachman, J. Pineau, D. Precup, and D. Meger, “Deep Reinforcement Learning that Matters,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 32, no. 1, 2018.
  * [44] E. Vinitsky, N. Jaques, J. Leibo, A. Castenada, and E. Hughes, “An Open Source Implementation of Sequential Social Dilemma Games,” https://github.com/eugenevinitsky/sequential_social_dilemma_games/issues/182, 2019, gitHub repository.
  * [45] C. Gini, “Variabilità e Mutabilità,” _Reprinted in Memorie di metodologica statistica (Ed. Pizetti E_ , 1912.
  * [46] P. J. Huber, _Robust statistics_. John Wiley & Sons, 2004, vol. 523.
  * [47] R. A. Maronna, R. D. Martin, V. J. Yohai, and M. Salibián-Barrera, _Robust Statistics: Theory and Methods (with R)_. John Wiley & Sons, 2019.
  * [48] F. R. Hampel, E. M. Ronchetti, P. J. Rousseeuw, and W. A. Stahel, _Robust Statistics: The Approach based on Influence Functions_. John Wiley & Sons, 2011, vol. 196.

## Acknowledgements

We thank Frédéric Herbreteau, Adrien Bennetot and Léo Heidelberger for their
help and support. N. Díaz-Rodríguez is currently supported by the Spanish
Government Juan de la Cierva Incorporación contract (IJC2019-039152-I).

##  VIII Supplementary Material

Supplementary material can be accessed here: https://bit.ly/3xG7ZXy


# Uncovering the behavior of Hf2Te2P and the candidate Dirac metal Zr2Te2P

K. -W. Chen,1,2 S. Das,1,2 D. Rhodes,1,2 S. Memaran,1,2 T. Besara,1 T.
Siegrist,1,3 E. Manousakis,1,2 L. Balicas,1 and R. E. Baumbach1 1National High
Magnetic Field Laboratory, Florida State University  2Department of Physics,
Florida State University  3Department of Chemical and Biomedical Engineering,
Florida State University

###### Abstract

Results are reported for single crystal specimens of Hf2Te2P and compared to
its structural analogue Zr2Te2P, which was recently proposed to be a potential
reservoir for Dirac physics.[Ji  _et al._ , 2015] Both materials are produced
using the iodine vapor phase transport method and the resulting crystals are
exfoliable. The bulk electrical transport and thermodynamic properties
indicate Fermi liquid behavior at low temperature for both compounds. Quantum
oscillations are observed in magnetization measurements for fields applied
parallel but not perpendicular to the $c$-axis, suggesting that the Fermi
surfaces are quasi-two dimensional. Frequencies are determined from quantum
oscillations for several parts of the Fermi surfaces. Lifshitz-Kosevich fits
to the temperature dependent amplitudes of the oscillations reveal small
effective masses, with a particularly small value $m^{*}$ $=$ 0.046$m_{0}$ for
the $\alpha$ branch of Zr2Te2P. Electronic structure calculations are in good
agreement with quantum oscillation results and illustrate the effect of a
stronger spin-orbit interaction going from Zr to Hf. These results suggest
that by using appropriate tuning parameters this class of materials may deepen
the pool of novel Dirac phenomena.

##  I Introduction

Efforts in condensed matter physics during the past century have focused on
states of matter that are characterized by their symmetry. This encompasses
phenomena such as magnetism, superconductivity, charge density wave order,
nematicity, etc. and was foundational for development of society transforming
modern technologies: e.g., those based on rare earth magnets and
superconductors. Another type of ordering has recently gained prominence,
where the electronic state is characterized by its topology instead of a
broken symmetry. The first examples of this were the quantum Hall effect and
fractional quantum Hall effect (e.g., in GaAs/AlGaAs) von Klitzing (1986);
Stormer (1999) and during the past decade this field has rapidly expanded to
include new families of topological materials, particularly those with three
dimensional crystalline structures. These three dimensional Dirac materials
host topologically protected linearly dispersing (“photon-like”) electronic
bands both in the bulk and as surface states. Wehling  _et al._ (2014) There
are various origins for this behavior including strong spin-orbit coupling
(e.g., topological insulators Bi2Te3 Zhang  _et al._ (2009); Chen  _et al._
(2009); Hsieh  _et al._ (2009) and Bi2Se3 Zhang  _et al._ (2009); Xia  _et
al._ (2009); Hsieh  _et al._ (2009)), crystalline symmetries (e.g., Dirac
semimetals Cd3As2 Liu  _et al._ (2014a); Wang  _et al._ (2013) and Na3Bi Wang
_et al._ (2012); Liu  _et al._ (2014b)), combinations of these effects (Weyl
semimetals (Nb,Ta)(P,As) Weng  _et al._ (2015); Lv  _et al._ (2015); Xu  _et
al._ (2015); Shekhar  _et al._ (2015)), and strong electronic correlations
(Kondo insulators SmB6 Dzero  _et al._ (2010); Xu  _et al._ (2013); Jiang  _et
al._ (2013)).

Given their novel physics and potential practical applications,Wehling  _et
al._ (2014); Liang  _et al._ (2015); Shekhar  _et al._ (2015) there is great
interest in diversifying the pool of 3D Dirac materials. Taking inspiration
from earlier successes in condensed matter physics, this will be accomplished
largely through tuning of known systems (and new materials as they are
uncovered) using control parameters such as magnetic fields, applied pressure,
chemical substitution, nanostructuring, and device fabrication. This route not
only gives insight into parent compounds, but also has a high likelihood of
uncovering transformative phenomena. This approach has already proven
successful for some topological systems: in Cd3As2 both applied pressure He
_et al._ (2015) and point contact experiments Aggarwal  _et al._ (2016)
stabilize superconductivity and may provide a route towards studies of
Majorana fermions.

Zr2Te2P was recently introduced as a possible strongly topological metal with
multiple Dirac cones: i.e., a Dirac metal when exfoliated into a few atomic
layers.Ji  _et al._ (2015) Of particular interest is that this compound
crystallizes in the well-known tetradamite structure and is easily exfoliated,
opening the way for device development. In order to explore the effect of a
strengthened spin-orbit interaction in this system, we undertook to synthesize
related compounds with larger Z elements. Here we report the synthesis and
bulk characterization of the isoelectronic analogue Hf2Te2P, which forms in
the same structure, is easily exfoliated, and supports similar bulk electrical
transport and thermodynamic behavior. Quantum oscillations are observed in
both compounds starting near magnetic fields of $H$ = 3 T, opening the
possibility of experimentally comparing their Fermi surfaces and
characterizing the Dirac behavior through the Berry’s phase. We additionally
present electronic structure calculations which reveal the effect of a
strengthening spin-orbit interaction going from Zr to Hf.

##  II Methods

Single crystals of Zr2Te2P and Hf2Te2P were grown from elements as previously
described for the Zr version.Tschulik  _et al._ (2009) Stoichiometric
polycrystalline precursor material was first prepared by reacting the raw
elements at 1000 oC for 24 hours. The resulting powders (1.5 g and 1.3 g for
Zr2Te2P and Hf2Te2P, respectively) were subsequently sealed under vacuum with
46 mg of iodine in quartz tubes with 18 mm diameter, 2 mm thickness, and 10 cm
length. The tubes were placed in a resistive tube furnace with a temperature
gradient with hot and cold zone temperatures of 900 and 800 oC, respectively.
The ampoules were held under these conditions for three weeks and then
quenched to room temperature. Large hexagonal single crystals formed near the
cold zone. For Zr2Te2P the linear dimensiona are as large as 10 $\times$ 10
$\times$ 0.4 mm. Hf2Te2P crystals were slightly smaller, with dimensions 3
$\times$ 3 $\times$ 1 mm (Fig. 1).

X-ray diffraction measurements were collected using a Scintag PAD-V
$\Theta$:2$\Theta$ diffractometer using Cu K-alpha radiation. The sample
consisted of crushed crystals spread randomly on a glass slide. Due to the
tendency of the crystals to separate as plates in the $ab$-plane, completely
random orientation was not achieved, and the XRD pattern shows signs (001)
texture. The pattern was fit with the WINPREP software Stahl  using the known
parameters of the isostructural compound Zr2Te2P as starting parameters.

Flakes of Zr2Te2P and Hf2Te2P were exfoliated from bulk crystals using the
micromechanical cleavage technique with adhesive scotch-tape and transferred
onto a $p$-doped silicon wafer substrate covered with a 270 nm layer of SiO2.
Atomically-thin crystals were identified under an optical microscope using
optical contrast. Atomic Force Microscopy (AFM) imaging was performed using an
Asylum Research MFP-3D AFM system. Images were taken about 2 hours after
exfoliation, during which the samples were exposed to air.

Zero magnetic field electrical resistivity $\rho$ was measured using the 3He
option in Quantum Design Physical Properties Measurement System for
temperatures 400 mK $<$ $T$ $<$ 300 K. Several individual crystals were
measured for each concentration, which revealed a high degree of batch
uniformity. Heat capacity $C$ measurements were performed using the same
apparatus for temperatures 400 mK $<$ $T$ $<$ 20 K. Magnetization $M(T,H)$
measurements were undertaken for single crystals for temperatures $T$ $=$ 1.8
- 300 K under an applied magnetic field of $H$ $=$ 10 kOe applied parallel
$\parallel$ and perpendicular $\perp$ to the $c$-axis using a Quantum Design
Magnetic Property Measurement System. Magnetic susceptibility $\chi$ is
defined as the ratio $M/H$.

Electronic structure calculations were performed using the Vienna _ab-initio_
simulation package Shishkin  _et al._ (2007); Fuchs  _et al._ (2007); Shishkin
and Kresse (2007, 2006) (VASP) within the generalized gradient approximation
(GGA). The contribution from the spin-orbit coupling is included in the
calculations. The Perdew-Burke-Ernzerhof (PBE) exchange correlation functional
Perdew  _et al._ (1996) and the projected augmented wave (PAW) methodology
Blöchl (1994) were used to describe the core electrons. The 4 _s_ , 4 _p_ , 5
_s_ , and 4 _d_ electrons for Zr, the 5 _s_ , 5 _p_ , 6 _s_ , and 5 _d_
electrons for Hf, the 4 _d_ , 5 _s_ and 5 _p_ electrons for Te, and the 3 _s_
and 3 _p_ electrons for P were treated as valence electrons in all
calculations. The energy cut off for the plane-wave basis was chosen to be 400
meV. A total of 240 bands and a $k-$point mesh of $20\times 20\times 4$ were
used for the self-consistent ground state calculations. A total of 100
$k-$points were chosen between each pair of special $k-$points in the
Brillouin-zone for the band-structure calculations. The Fermi surfaces were
generated with a $k-$point mesh of $20\times 20\times 4$, using the
eigenvalues obtained from VASP and were visualized using the XCrysden
software. Kokalj (2003)

##  III Results and Discussion

Results from powder X-ray diffraction measurement of crushed Hf2Te2P crystals
are shown in Fig. 1, along with a pattern fit. Parameters of the isostructural
compound Zr2Te2P, which crystallizes in the space group R$\bar{3}$m (#166),
were used as starting parameters. The fit yielded lattice parameters $a$ $=$
3.7946(2) Å and $c$ $=$ 29.140(1) Å, and displayed no traces of impurity
phases. The lattice parameters are smaller than the Zr2Te2P parameters ($a$
$=$ 3.8119(3) Å and $c$ $=$ 29.189(3) Å),Tschulik  _et al._ (2009) indicative
of the smaller ionic radius of Hf compared to Zr. Teatum  _et al._ (1960)

Atomic force microscopy (AFM) results for exfoliated flakes of Zr2Te2P and
Hf2Te2P are shown in Figs. 1d,e. Edge step profiles reveal that the crystals
are easily reduced to single- or several-unit cell flakes, where the edge step
can be the same size as a single $c$-axis unit cell. Samples that are left
exposed to air are seen to change over time, revealing that these specimens
are air sensitive.

Figure 1: (a) Unit cell of the $M$2Te2P ($M$ $=$ Zr and Hf) compounds. (b)
Single crystal specimens of Hf2Te2P. (c) Powder x-ray diffraction pattern for
Hf2Te2P. The solid red line is a patternf fit to the data to determine lattice
parameters. (d) Atomic Force Microscopy (AFM) image collected from an
atomically thin flake of Zr2Te2P on a SiO2 substrate. Inset: Height Profile
(along the red line shown in the image) indicating a thickness of
approximately 6nm equivalent to two unit cells. (e) AFM image for Hf2Te2P.
Inset: Height profile for two-unit cell flake.

Electrical resistance normalized to its room temperature value
$R/R_{\rm{300K}}$ as a function of temperature $T$ data are shown in Fig. 2a
for Zr2Te2P and Hf2Te2P, where the electrical current was applied in the basal
$ab$-plane. Metallic behavior is observed for both materials, with room
temperature resistivities $\rho$300K $\approx$ 60 $-$ 80 $\mu$$\Omega$cm,
depending on uncertainty of the geometric factor. $R/R_{\rm{300K}}$ decreases
linearly with decreasing temperature down to $T$ $\approx$ 50 K, where the
resistance begins to saturate in a Fermi-liquid like manner. For both
compounds, the resistivity saturates near 3 $-$ 5 $\mu$$\Omega$cm, giving a
residual resistivity ratio $RRR$ = $R_{\rm{300K}}$/$R_{\rm{0}}$ $\approx$ 20\.
This value for $RRR$ is not particularly large, and might be taken to indicate
appreciable disorder scattering. Alternatively, it is possible that as-cast
crystals have large $RRR$ but rapidly begin to degrade after being exposed to
atmosphere. As shown below, quantum oscillations are easily observed in
magnetization measurements, which is unexpected if there is large disorder
scattering.

The temperature dependences of the heat capacity $C$ divided by $T$ are shown
in Figs. 2b,c for Zr2Te2P and Hf2Te2P. Fermi liquid behavior is seen for both
compounds, where the expression $C/T$ $=$ $\gamma$ $+$ $\beta$$T^{2}$
describes the data for 500 mK $<$ $T$ $<$ 10 K. Fits to the data give $\gamma$
$=$ 5.4 mJ/mol-K2 and 4.7 mJ/mol-K2 and $\beta$ $=$ 0.62 and 0.71 mJ/mol-K4
for Zr2Te2P and Hf2Te2P, respectively. That the electronic coefficients of the
heat capacity $\gamma$ are similar for both compounds is consistent with
electronic structure calculations presented below, where the total density of
states at the Fermi energy remains roughly constant upon going from Zr to Hf.
From the expression $\beta$ $=$ $r$1944($T/\theta_{\rm{D}}$)J/mol-K, where $r$
is the number of atoms per formula unit, the Debye temperatures
$\theta_{\rm{D}}$ are calculated to be 250 K and 239 K for the Zr and Hf
compounds, respectively. This result is consistent with the increased atomic
mass of Hf in comparison to Zr.

Figure 2: (a)Room temperature normalized electrical resistance
$R/R_{\rm{300K}}$ for $M$2Te2P ($M$ $=$ Zr and Hf). (b) Heat capacity $C$
divided by temperature $T$ vs. $T^{2}$ for $M$2Te2P ($M$ $=$ Zr and Hf)
showing linear behavior as expected for a Fermi liquid. (c) $C/T$ vs $T$ for
$M_{2}$Te2P ($M$ = Zr and Hf.)

Fig. 3 displays magnetization $M$ as a function of the field $H$ applied
either parallel $\parallel$ or perpendicular $\perp$ to the $c$-axis and for
both compounds. $M(H)$ is diamagnetic for $H$ $\perp$ $c$, whereas curves for
$H$ $\parallel$ $c$ initially are paramagnetic with increasing $H$, but
eventually take on negative slopes, and even become negative for Hf2Te2P. The
initial positive slopes can be attributed to a small fraction of paramagnetic
impurities. In the insets, the magnetic susceptibilities $\chi$ $=$ $M/H$ for
Zr2Te2P and Hf2Te2P are shown, where diamagnetic behavior with little
anisotropy is observed at large temperatures. Weak upturns are seen at low
$T$, which may be due to Curie tails from paramagnetic impurities. Quantum
oscillations are observed starting from low fields for both compounds when $H$
$\parallel$ $c$. The very low field for onset of oscillations indicates the
high sample quality. QOs are not seen when $H$ $\perp$ $c$, giving strong
evidence that the Fermi surface is two dimensional. While this is expected
from the band structure calculations presented below, further angle resolved
measurements will be useful to illustrate the degree of two dimensionality.

Figure 3: (a) Magnetization $M$ vs. magnetic field $H$ for temperature $T$ =
1.8 K for magnetic field $H$ applied parallel $\parallel$ and perpendicular
$\perp$ to the $c$-axis for Zr2Te2P. (b) Magnetic susceptibility $\chi$ $=$
$M/H$ vs. $T$ for $H$ = 10 kOe applied $\parallel$ and $\perp$ to the $c$-axis
for Zr2Te2P. (c) $M$ vs. $H$ for 1.8 K for $H$ $\parallel$ and $\perp$ to the
$c$-axis for Hf2Te2P (d) $\chi$ vs. $T$ for $H$ = 10 kOe applied $\parallel$
and $\perp$ to the $c$-axis for Hf2Te2P.

The quantum oscillations in the magnetization $\Delta$$M$ are isolated by
subtracting a 3rd order polynomial background as shown in Fig. 4 for different
temperatures $T$ between 1.8 K and 30 K over the field range $H$ $=$ 3 T $-$ 7
T. Fast Fourier transforms (FFT) for $\Delta$$M$ at various temperatures were
performed using a Hanning window to extract the oscillation frequencies (Figs.
4b,e). The FFTs yield frequencies $F_{\alpha}$ $=$ 11.4 T and $F_{\beta}$ $=$
89.9 T for Zr2Te2P and $F_{\gamma}$ $=$ 106.8 T and $F_{\eta}$ $=$ 146.6 T for
Hf2Te2P, as summarized in Table 1. The temperature dependences of the FFT
amplitude are shown in Figs. 4c,f. The solid lines are fits to the data
according to the thermal damping factor $R_{\rm{T}}$ of the Lifshitz-Kosevich
formula using the expression,

$R_{\rm{T}}=\frac{2\pi^{2}k_{\rm{B}}Tm^{*}/\hbar
eH}{sinh(2\pi^{2}k_{\rm{B}}Tm^{*}/\hbar eH)}$ (1)

where $m^{*}$ is the effective cyclotron mass. The fits yield small effective
masses as collected in Table 1. Particularly noteworthy is the very small
$m^{*}$ $=$ 0.046$m_{0}$ for the $\alpha$ branch of Zr2Te2P at
$F_{\rm{\alpha}}$ $=$ 11.4 T, where $m_{0}$ is the free electron mass.

Table 1: Summary of results from fast Fourier transform (FFT) of quantum oscillations for Zr2Te2P and Hf2Te2P. Frequencies $F$(T) and effective masses $m^{*}$($m_{0}$) are listed. |  $F_{\rm{\alpha}}$(T), $m_{\alpha}^{*}$($m_{0}$) |  $F_{\rm{\beta}}$, $m_{\beta}^{*}$ |  $F_{\rm{\gamma}}$, $m_{\gamma}^{*}$ |  $F_{\rm{\eta}}$, $m_{\eta}^{*}$
---|---|---|---|---
Zr | 11.4, 0.046 | 89.9, 0.22 | — | —
Hf | — | — | 106.8, 0.27 | 146.6, 0.22
Figure 4: (a) Background subtracted magnetization $\Delta$$M$ for Zr2Te2P
showing de Hass-van Alphen (dHvA) oscillations for different temperatures $T$
between 1.8 K and 30 K in a field range $H$ $=$ 3 T $-$ 7 T. $\Delta$$M$ is
obtained by subtracting a 3rd order polynomial background from the raw data.
(b) Fast Fourier transforms (FFT ) for $\Delta$$M$ at various $T$s for
Zr2Te2P. (c) Amplitude of the peaks observed in FFT spectra as a function of
the temperature for Zr2Te2P. The solid lines are fits to the data using the
Lifshitz-Kosevich formula. (d) $\Delta$$M$ for Hf2Te2P. (e) FFTs of
$\Delta$$M$ at various $T$s (using a Hanning window) for Hf2Te2P. (f) $T$
dependence of the FFT amplitude for $\Delta$$M$ for Hf2Te2P.

The calculated Fermi surfaces and band structures of Zr2Te2P and Hf2Te2P are
shown in Fig. 5 and Fig. 7, respectively, revealing that they are quite
similar. Fig. 6 compares the effect of the spin-orbit interaction on the
pockets near the $\Gamma$ point obtained in GGA with and without spin-orbit
coupling for the Zr (top row) and Hf (bottow row) based 221 compound. Notice
that in going from GGA-spin (left column) to GGA-spin+SOC (right column) the
size of the orbits are not significantly affected for the Zr based compound.
In the case of the Hf based compound, in the absence of SOC the smaller orbit
is estimated to be significantly (about a factor of 2) larger than the
smallest observed value of the oscillation frequency ($\sim$ 100 T). However,
as can be seen from the bottom row of Fig. 6 the effect of SOC is to reduce
the lowest orbit by a significant factor which brings the value of the
smallest frequency ($\sim$ 80 T) in the observed range. Rough estimates of the
size of the three lowest frequencies for the Zr based compound are $\sim$ 30
T, 60 $-$ 70 T, 110 $-$ 120 T. For the case of the Hf based compound the three
lowest frequencies are roughly 70 $-$ 80 T, 110 $-$ 120 T, and $\sim$ 200 T.
Fig. 7 presents the band structure for the Zr and Hf based compounds obtained
with the inclusion of the SOC. They are very similar, where there are small
differences in the size of the energy separation between the bands near Fermi
level at the $\Gamma$ point. This is also true at the K point for those bands
closest to the Fermi level.

Figure 5: (color-online) Comparison between the Fermi surface of the Zr (a) Hf
(b) based 221 compounds.

Figure 6: (color-online) Comparison between the Fermi pockets near the
$\Gamma$ point obtained in GGA with (right column) and without (left column)
spin-orbit coupling for the Zr (top row) and Hf (bottow row) based 221
compound.

Figure 7: (color-online)The electronic band structure of (a) the Zr-based and
(b) Hf-based compound obtained with the inclusion spin-orbit interaction.

##  IV Conclusion

We have used iodine vapor phase transport to grow single crystals of the
compounds Zr2Te2P and Hf2Te2P. Measurements of the bulk electrical transport
and thermodynamic properties indicate Fermi liquid behavior at low $T$. Both
compounds readily exhibit quantum oscillations in magnetization measurements,
from which we have determined the frequencies associated with several parts of
the Fermi surfaces. Lifshitz-Kosevich fits to the temperature dependent
amplitudes of the oscillations reveals small effective masses, with a
particularly small value for the $\alpha$ branch of Zr2Te2P. Electronic
structure calculations are presented that are in good agreement with these
results. Calculations also indicate that measurements in higher magnetic
fields will uncover larger Fermi surface sections with higher oscillation
frequencies. Fits to the quantum oscillation data also provide preliminarly
evidence that some bands have a Berry’s phase. This will be presented in
detail in a separate manuscript.

These results make both Zr2Te2P and Hf2Te2P intriguing candidates for further
study as potential hosts for Dirac physics. Of particular interest is to
determine whether the calculated Dirac points that are away from the Fermi
energy have an impact on the physical properties, or if the Fermi energy can
be tuned to be closer to them. These materials are also particularly well
suited to spectroscopic investigations such as ARPES, given the ease with
which they are cleaved. This property also opens the possibility of device
development. Finally, we point out that this family of materials likely can be
expanded to include a variety of chemical analogues: e.g., $M$2Te2As ($M$ $=$
Ti, Zr, and Hf) are likely to form.

##  V Acknowledgements

This work was performed at the National High Magnetic Field Laboratory
(NHMFL), which is supported by National Science Foundation Cooperative
Agreement No. DMR-1157490, the State of Florida and the DOE. A portion of this
work was supported by the NHMFL User Collaboration Grant Program (UCGP). LB,
DR, and SM are supported by DOE-BES through award DE-SC0002613. TB and TS are
supported by DOE-BES through award DE-SC0008832.

## References

  * Ji  _et al._ (2015) H. Ji, I. Pletikosi, Q. D. Gibson, G. Sahasrabudhe, T. Valla, and C. R. J., arXiv:1512.00528 (2015).
  * von Klitzing (1986) K. von Klitzing, Rev. Mod. Phys. 58, 519 (1986).
  * Stormer (1999) H. L. Stormer, Rev. Mod. Phys. 71, 875 (1999).
  * Wehling  _et al._ (2014) T. Wehling, A. Black-Schaffer, and A. Balatsky, Advances in Physics 63, 1 (2014).
  * Zhang  _et al._ (2009) H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nat Phys 5, 438 (2009).
  * Chen  _et al._ (2009) Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, Science 325, 178 (2009).
  * Hsieh  _et al._ (2009) D. Hsieh, Y. Xia, D. Qian, L. Wray, J. H. Dil, F. Meier, J. Osterwalder, L. Patthey, J. G. Checkelsky, N. P. Ong, A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nature 460, 1101 (2009).
  * Xia  _et al._ (2009) Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nat Phys 5, 398 (2009).
  * Liu  _et al._ (2014a) Z. K. Liu, J. Jiang, B. Zhou, Z. J. Wang, Y. Zhang, H. M. Weng, D. Prabhakaran, S.-K. Mo, H. Peng, P. Dudin, T. Kim, M. Hoesch, Z. Fang, X. Dai, Z. X. Shen, D. L. Feng, Z. Hussain, and Y. L. Chen, Nat Mater 13, 677 (2014a).
  * Wang  _et al._ (2013) Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Phys. Rev. B 88, 125427 (2013).
  * Wang  _et al._ (2012) Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, and Z. Fang, Phys. Rev. B 85, 195320 (2012).
  * Liu  _et al._ (2014b) Z. K. Liu, B. Zhou, Y. Zhang, Z. J. Wang, H. M. Weng, D. Prabhakaran, S.-K. Mo, Z. X. Shen, Z. Fang, X. Dai, Z. Hussain, and Y. L. Chen, Science 343, 864 (2014b).
  * Weng  _et al._ (2015) H. Weng, C. Fang, Z. Fang, B. A. Bernevig, and X. Dai, Phys. Rev. X 5, 011029 (2015).
  * Lv  _et al._ (2015) B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian, and H. Ding, Phys. Rev. X 5, 031013 (2015).
  * Xu  _et al._ (2015) S.-Y. Xu, I. Belopolski, D. S. Sanchez, C. Zhang, G. Chang, C. Guo, G. Bian, Z. Yuan, H. Lu, T.-R. Chang, P. P. Shibayev, M. L. Prokopovych, N. Alidoust, H. Zheng, C.-C. Lee, S.-M. Huang, R. Sankar, F. Chou, C.-H. Hsu, H.-T. Jeng, A. Bansil, T. Neupert, V. N. Strocov, H. Lin, S. Jia, and M. Z. Hasan, Science Advances 1, 1501092 (2015).
  * Shekhar  _et al._ (2015) C. Shekhar, A. K. Nayak, Y. Sun, M. Schmidt, M. Nicklas, I. Leermakers, U. Zeitler, Y. Skourski, J. Wosnitza, Z. Liu, Y. Chen, W. Schnelle, H. Borrmann, Y. Grin, C. Felser, and B. Yan, Nat Phys 11, 645 (2015).
  * Dzero  _et al._ (2010) M. Dzero, K. Sun, V. Galitski, and P. Coleman, Phys. Rev. Lett. 104, 106408 (2010).
  * Xu  _et al._ (2013) N. Xu, X. Shi, P. K. Biswas, C. E. Matt, R. S. Dhaka, Y. Huang, N. C. Plumb, M. Radović, J. H. Dil, E. Pomjakushina, K. Conder, A. Amato, Z. Salman, D. M. Paul, J. Mesot, H. Ding, and M. Shi, Phys. Rev. B 88, 121102 (2013).
  * Jiang  _et al._ (2013) J. Jiang, S. Li, T. Zhang, Z. Sun, F. Chen, Z. Ye, M. Xu, Q. Ge, S. Tan, X. Niu, M. Xia, B. Xie, Y. Li, X. Chen, H. Wen, and D. Feng, Nat Commun 4, 3010 (2013).
  * Liang  _et al._ (2015) T. Liang, Q. Gibson, M. N. Ali, M. Liu, R. J. Cava, and N. P. Ong, Nature Materials 14, 280 (2015).
  * He  _et al._ (2015) L. P. He, Y. T. Jia, S. J. Zhang, X. C. Hong, C. Q. Jin, and S. Y. Li, arXiv:1502.02509 (2015).
  * Aggarwal  _et al._ (2016) L. Aggarwal, A. Gaurav, G. S. Thakur, Z. Haque, A. K. Ganguli, and G. Sheet, Nat Mater 15, 32 (2016).
  * Tschulik  _et al._ (2009) K. Tschulik, M. Ruck, M. Binnewies, E. Milke, S. Hoffmann, W. Schnelle, B. P. T. Fokwa, M. Gillesen, and S. P., Eur. J. Inorg. Chem 2009, 3102 (2009).
  * (24) K. Stahl, “Winprep,” Lyngby, Denmark.
  * Shishkin  _et al._ (2007) M. Shishkin, M. Marsman, and G. Kresse, Phys. Rev. Lett. 99, 246403 (2007).
  * Fuchs  _et al._ (2007) F. Fuchs, J. Furthmüller, F. Bechstedt, M. Shishkin, and G. Kresse, Phys. Rev. B 76, 115109 (2007).
  * Shishkin and Kresse (2007) M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (2007).
  * Shishkin and Kresse (2006) M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (2006).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, , and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kokalj (2003) A. Kokalj, Comp. Mater. Sci. 28, 155 (2003), code available from http://www.xcrysden.org.
  * Teatum  _et al._ (1960) E. Teatum, J. Waber, and K. A. Gschneidner, “Compilation of calculated data useful in predicting metallurgical behavior of the elements in binary alloy systems,” Los Alamos Scientific Laboratory of the University of California (1960).