# Toward Policy Explanations for Multi-Agent Reinforcement Learning

Kayla Boggess1 Sarit Kraus2&Lu Feng1
1University of Virginia
2Bar-Ilan University
{kjb5we<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Advances in multi-agent reinforcement learning (MARL) enable sequential
decision making for a range of exciting multi-agent applications such as
cooperative AI and autonomous driving. Explaining agent decisions is crucial
for improving system transparency, increasing user satisfaction, and
facilitating human-agent collaboration. However, existing works on explainable
reinforcement learning mostly focus on the single-agent setting and are not
suitable for addressing challenges posed by multi-agent environments. We
present novel methods to generate two types of policy explanations for MARL:
(i) policy summarization about the agent cooperation and task sequence, and
(ii) language explanations to answer queries about agent behavior.
Experimental results on three MARL domains demonstrate the scalability of our
methods. A user study shows that the generated explanations significantly
improve user performance and increase subjective ratings on metrics such as
user satisfaction.

##  1 Introduction

Recent years have witnessed a growing body of research in multi-agent
reinforcement learning (MARL), enabling sequential decision making for a range
of exciting multi-agent applications such as cooperative AI Dafoe et al.
(2020) and autonomous driving Kiran et al. (2021). Generating explanations
about agent decisions is crucial for improving system transparency, increasing
user satisfaction, and facilitating human-agent collaboration Kraus et al.
(2020); Chakraborti et al. (2020). However, existing works on explainable
reinforcement learning (RL) mostly focus on the single-agent setting Wells and
Bednarz (2021); Heuillet et al. (2021); Puiutta and Veith (2020).

Generating explanations for MARL agents that interact with each other in a
common environment poses significant challenges. The combinatorial nature of
MARL (i.e., the joint state/action space grows exponentially with the number
of agents) leads to scalability issues. Explanations should provide adequate
information about agent behavior, including the interaction (e.g.,
cooperation) among multiple agents, for user understanding. Furthermore,
explanations should avoid redundant information that may overwhelm or confuse
users, thus decreasing user satisfaction and trust.

To tackle these challenges, we develop novel methods to generate two types of
policy explanations for MARL: (i) policy summarization, and (ii) query-based
language explanations. Our methods rely on first building an abstract
representation of MARL policy as a multi-agent Markov decision process (MMDP),
which can be obtained by abstracting samples observed during the MARL policy
evaluation. The proposed method generates a summarization about the most
probable sequence of agent behavior under a given MARL policy, by finding the
most probable path through the MMDP abstraction and extracting information
about agent cooperation and task sequence. The generated policy summarizations
can help users to have a global view of agent decisions and support human-
agent collaboration (e.g., users may adjust their workflow based on agents’
task sequence).

Additionally, we developed methods for generating language explanations to
answer three types of queries about agent behavior, including “When do
[agents] do [actions]?”, “Why don’t [agents] do [actions] in [states]?”, “What
do [agents] do in [conditions]?” Such explanations can enable users to
understand specific agent decisions, debug faulty agent behavior, and refine
user mental models. Our work is inspired by the method proposed in Hayes and
Shah (2017), which computes a minimal Boolean logic expression covering states
satisfying the query criteria, and converts the Boolean expression to
explanations via language templates. However, we find that a naive adaptation
of this method for MARL generates explanations with redundant information and
has limited scalability. Further, the generated explanations do not
necessarily capture agent cooperation, which is imperative for explaining
MARL. We proposed improved methods that address these limitations by
leveraging MARL domain knowledge (e.g., agent cooperation requirements) to
filter relevant agent states and actions.

We applied a prototype implementation of the proposed methods to three
benchmark MARL domains: (i) multi-robot search and rescue, (ii) multi-robot
warehouse, and (iii) level-based foraging Papoudakis et al. (2021).
Experimental results demonstrate that the proposed methods can generate policy
summarizations and query-based explanations for large MARL environments with
up to 19 agents.

Finally, we conducted a user study to evaluate the quality of generated
explanations. We measured user performance on correctly answering questions
based on explanations, to test user understanding of agent behavior. We also
collected user subjective ratings on _explanation goodness metrics_ Hoffman et
al. (2018). The results show that the generated explanations significantly
improve user performance and increase subjective ratings on various metrics
including user satisfaction.

##  2 Related Work

Explaining agent decision making has recently emerged as a focus area within
the explainable AI paradigm. Chakraborti et al. (2020) provides a survey about
this emerging landscape of explainable decision making. Kraus et al. (2020)
proposes Explainable Decisions in Multi-Agent Environments (xMASE) as a new
research direction, emphasizing many challenges of generating multi-agent
explanations, such as accounting for agent interactions and user satisfaction.

Explainable RL has been attracting increasing interest, as shown in several
recent surveys Wells and Bednarz (2021); Heuillet et al. (2021); Puiutta and
Veith (2020). In particular, Wells and Bednarz (2021) points out the lack of
user studies as a major limitation across existing works. Moreover, current
approaches mostly focus on the single-agent setting, while generating
explanations for MARL has received scant attention so far. Kazhdan et al.
(2020) extracts MARL model as abstract argumentations, which do not explicitly
consider agent cooperation on the same tasks as in our work.

The proposed methods are inspired by several prior works. Topin and Veloso
(2019) develops Abstracted Policy Graphs (i.e., Markov chains of abstract
states) for explaining single-agent RL. Our work differs in that we adopt MMDP
to represent multi-agent policy abstractions. Amir and Amir (2018) summarizes
agent behavior by extracting trajectories from agent simulations and
visualized them as videos. Inspired by this idea, we created GIF animations as
the baseline for evaluating policy summarizations in the user study. But we
extracted the sequence of agent actions differently, focusing on the agent
cooperation. Hayes and Shah (2017) generates RL policy descriptions to answer
queries. We adapted this method to multi-agent environments and proposed
significant improvements, which will be described in Section 4.

##  3 Policy Abstraction and Summarization

We describe how to build an abstract representation of MARL policy in Section
3.1, and present the method for generating policy summarization in Section
3.2.

###  3.1 Policy Abstraction

In the context of MARL, a group of $N$ agents interact with each other in a
common environment and make decisions influenced by the joint states of all
agents. Agent decisions can be captured by a joint policy
$\pi:{\mathcal{X}}\to\Delta({\mathcal{A}})$, which is a function mapping the
set of joint states ${\mathcal{X}}=\\{(x^{1},\dots,x^{N})\\}$ to a
probabilistic distribution over the set of joint actions
${\mathcal{A}}=\\{(a^{1},\dots,a^{N})\\}$, where $x^{i}$ (resp. $a^{i}$)
denotes the state (resp. action) of agent $i$. Once a policy is trained,
agents can act upon it in any given state. But there is a lack of a global
view of the entire policy. Further, the size of the policy grows exponentially
with the number of agents and state variables. To address these issues, we
propose to build an abstract representation of the policy as the basis for
generating explanations about agent behavior.

We use the multi-agent Markov decision process (MMDP) framework to represent
MARL policy abstraction. Formally, an MMDP is a tuple
$({\mathcal{S}},{\mathcal{A}},{\mathcal{T}})$, where
${\mathcal{S}}=\\{(s^{1},\dots,s^{N})\\}$ is the joint (abstract) state space,
${\mathcal{A}}$ is the joint action space, and ${\mathcal{T}}$ is the
transition function. Let ${\mathcal{F}}$ be a set of Boolean predicates
indicating features of the MARL domain. We denote by $f(x^{i})=1$ if an agent
state $x^{i}$ satisfies a feature predicate $f\in{\mathcal{F}}$. An abstract
state $s^{i}$ is then given by the satisfaction of all feature predicates
$f\in{\mathcal{F}}$, where each bit of the binary encoding of
$s^{i}\in\mathbb{N}$ corresponds to the satisfaction of a predicate
$f(x^{i})$. Thus, the choice of features affects the abstraction level and
should include adequate information for explanations. In this work, we assume
that users specify a set of feature predicates for a given MARL domain.

Once an MARL policy is trained, we build an MMDP during the policy evaluation
stage. For each sample $(\mathbf{x},\mathbf{a},\mathbf{x}^{\prime})$,
determine an MMDP transition
$\mathbf{s}\xrightarrow{\mathbf{a}}\mathbf{s}^{\prime}$ by finding the
abstract state $\mathbf{s}$ (resp. $\mathbf{s^{\prime}}$) corresponding to
$\mathbf{x}$ (resp. $\mathbf{x^{\prime}}$). When policy evaluation terminates
(e.g., converging to the expected reward), compute the transition probability
${\mathcal{T}}(\mathbf{s},\mathbf{a},\mathbf{s}^{\prime})$ via frequency
counting.

Properties. The resulting MMDP is a _sound_ abstraction of the MARL policy
because, by construction, every MMDP transition with non-zero probability
corresponds to at least one sampled policy decision. The state space size
$|{\mathcal{S}}|$ is bounded by ${\mathcal{O}}({2^{|{\mathcal{F}}|}}^{N})$,
depending on the number of agents $N$ and feature predicates
$|{\mathcal{F}}|$. In practice, a trained MARL policy may only induce a small
set of reachable states.

Figure 1: Example MARL domain of multi-robot search and rescue.

###### Example 1

Figure 1(a) shows an example MARL domain where three robotic agents cooperate
to complete search and rescue tasks. Rescuing the victim requires the
cooperation of an unmanned aerial vehicle (UAV) and an unmanned ground vehicle
(UGV). Any agent can fight the fire, which is blocked by the wall and
obstacle. Removing the obstacle requires the cooperation of two UGVs. Given a
trained MARL policy, we build an MMDP abstraction with 6 feature predicates
indicating whether each task is detected or completed (e.g.,
$\mathsf{victim\\_detect}$, $\mathsf{victim\\_complete}$). An agent can only
detect a task in a neighboring grid (e.g., UAV detects the victim in Figure
1(a)). The resulting MMDP has 63 (reachable) states and 577 transitions.

###  3.2 Policy Summarization

A policy abstraction containing hundreds of states and transitions is too
complex for humans to understand. An alternative way of communicating agent
behavior is to show execution traces; however, a lengthy trace may be
burdensome for users to review. To overcome these limitations, we develop a
method to generate policy summarization, illustrating the agent cooperation
and task sequence for the most probable sequence of agent behavior under a
given MARL policy.

Algorithm 1 Generating Policy Summarization

Input: policy abstraction
${\mathcal{M}}=({\mathcal{S}},{\mathcal{A}},{\mathcal{T}})$, task completion
predicates ${\mathcal{F}}_{c}\subseteq{\mathcal{F}}$
Output: policy summarization ${\mathcal{Z}}$

1:${\mathcal{Z}}\leftarrow\\{\\}$

2:Compute the most probable path $\rho$ through ${\mathcal{M}}$

3:for $0\leq t\leq|\rho|$ do

4: $y\leftarrow$ new array

5: for $1\leq i\leq N$ do

6: $y[i]\leftarrow$ {}

7: for $f\in{\mathcal{F}}_{c}$ do

8: if agent state $s^{i}_{t}$ in the path $\rho$ satisfies $f$ then

9: insert $f$ to $y[i]$

10: insert non-empty array $y$ to ${\mathcal{Z}}$

11:return ${\mathcal{Z}}$

Algorithm 1 shows the proposed method, which takes the input of a policy
abstraction ${\mathcal{M}}$ and a set of predicates ${\mathcal{F}}_{c}$
representing the completion of tasks (subgoals) in a given MARL domain. The
first step is to compute the most probable path
$\rho=\mathbf{s}_{0}\xrightarrow{\mathbf{a}_{0}}\mathbf{s}_{1}\xrightarrow{\mathbf{a}_{1}}\cdots$
from the initial state to a goal state in the MMDP ${\mathcal{M}}$, which
represents the most probable sequence of agent decisions under the policy.
This problem can be solved by converting the MMDP to a directed weighted graph
with edge weight
$e(\mathbf{s},\mathbf{a},\mathbf{s^{\prime}})=-\log{\mathcal{T}}(\mathbf{s},\mathbf{a},\mathbf{s^{\prime}})$
for each transition, and then applying the Dijkstra’s algorithm Dijkstra et
al. (1959) to find the shortest path.

Next, the algorithm loops through every joint state $\mathbf{s}_{t}$ in the
path $\rho$ to extract the agent cooperation and task sequence. At each step
$t$, the algorithm checks if an agent state $s^{i}_{t}$ satisfies any task
completion predicate $f\in{\mathcal{F}}_{c}$ and inserts completed task $f$
into the array element $y[i]$ (line 4-9). An agent only satisfies a task
completion predicate at step $t$ when it finishes the task and receives a
reward. We assume that if a task is completed via the cooperation of multiple
agents, they must satisfy the task predicate $f$ at the same step $t$ and each
receive a portion of the reward. Thus, the agent cooperation is represented as
multiple elements of the array $y$ sharing the same task. Only non-empty
arrays containing completed tasks are inserted into the summarization
${\mathcal{Z}}$. When the algorithm terminates, the generated summarization is
visualized as a chart, with each column corresponding to a non-empty $y$-array
and each row representing an agent’s task sequence.

Properties. The generated policy summarization ${\mathcal{Z}}$ is _sound_ ,
because it is derived from the most probable path of a _sound_ policy
abstraction (see Section 3.1). The complexity of computing the most probable
path is bounded by ${\mathcal{O}}(|{\mathcal{S}}|^{2})$, following the
complexity of the Dijkstra’s algorithm and depending on the MMDP state space
size. The rest of Algorithm 1 is bounded by ${\mathcal{O}}(|\rho|\cdot
N\cdot|{\mathcal{F}}_{c}|)$, depending on the path length and the number of
agents and tasks.

###### Example 2

We apply Algorithm 1 using the policy abstraction and task predicates from
Example 1. There are 8 states in the most probable path from the initial state
(i.e., all agents starting in the green grid) to a goal state (i.e., all tasks
have been completed). Figure 1(b) visualizes the generated summarization, with
column names (i.e., T1, T2, T3) indicating the sequence of task completions:
UGV2 and UAV cooperate to rescue the victim; next, UGV1 and UGV2 cooperate to
remove the obstacle; and lastly, UAV fights the fire.

##  4 Query-Based Explanations

While policy summarization provides a global view of the agent behavior under
a MARL policy, users may also query about specific agent decisions. In this
section, we develop methods for generating language explanations to answer the
following three types of queries:

  * •

_“When do [agents] do [actions]?”_ for identifying conditions for action(s) of
a single or multiple agent(s) .

  * •

_“Why don’t [agents] do [actions] in [states]?”_ for understanding differences
in expected and observed behaviors of a single or multiple agent(s).

  * •

_“What do [agents] do in [predicates]?”_ for revealing agent behavior under
specific conditions described by the given predicates.

Our work is inspired by a method developed in Hayes and Shah (2017) to
generate query-based explanations for single-agent RL. In the following, we
propose new methods to tackle limitations posed by adapting this baseline
method to multi-agent environments.

###  4.1 Explanations for When Query

Algorithm 2 presents both the baseline and proposed methods for answering
“When do agents $G_{q}$ do actions $A_{q}$?”, where $G_{q}$ and $A_{q}$ are
sets of agents and actions, respectively. The text in blue highlights changes
about relevancy filters (RF) for the proposed method (called WithRF) compared
to the baseline (called NoRF).

WithRF starts the algorithm (line 1-5) by identifying relevant agents $G$,
features $F$, and action sets $A$ based on domain knowledge (e.g., agent
cooperation requirements). For example, consider a query “When does UAV rescue
the victim?”. The domain knowledge is that rescuing the victim requires the
cooperation of a UAV and a UGV (Example 1). Thus, the relevant agent set $G$
is {$\mathsf{UVA,UGV\\_1,UGV\\_2}$}. The relevant feature set $F$ is
{$\mathsf{victim\\_detect,victim\\_complete}$}, while predicates about the
fire and obstacle are irrelevant. The relevant action sets $A$ is an array
with each element representing one possible set of agent actions required for
cooperation:
$[\\{\mathsf{UAV\\_rescue},\mathsf{UGV_{1}\\_rescue}\\},\\{\mathsf{UAV\\_rescue},\mathsf{UGV_{2}\\_rescue}\\}]$,
which can be generated based on the aforementioned domain knowledge about
agent cooperation requirements.

Both NoRF and WithRF loop through all the joint states
$\mathbf{s}\in{\mathcal{S}}$ of the policy abstraction MMDP and check all the
_enabled_ (i.e., with non-zero transition probability) joint actions
$\mathbf{a}$ in state $\mathbf{s}$. In line 9 of Algorithm 2, NoRF checks if
$\mathbf{a}$ is _compatible_ with $A_{q}$; that is, every agent action $a\in
A_{q}$ is contained in the joint action $\mathbf{a}=(a^{1},\dots,a^{N})$. By
contrast, WithRF checks if $\mathbf{a}$ is compatible with at least one set of
relevant actions contained in the array $A$. Following the previous example,
NoRF checks if $\mathbf{a}$ contains $\mathsf{UAV\\_rescue}$, while WithRF
checks if $\mathbf{a}$ contains
$\\{\mathsf{UAV\\_rescue},\mathsf{UGV_{1}\\_rescue}\\}$ _or_
$\\{\mathsf{UAV\\_rescue},\mathsf{UGV_{2}\\_rescue}\\}$. Since each element of
$A$ is a super-set of $A_{q}$, the WithRF check is more restrictive.

Algorithm 2 Generating Query-Based Explanations

Input: policy abstraction $({\mathcal{S}},{\mathcal{A}},{\mathcal{T}})$, query
“when do agents $G_{q}$ do actions $A_{q}$?”
Output: explanations ${\mathcal{E}}$

1:$G\leftarrow\\{\\}$; $F\leftarrow\\{\\}$; $A\leftarrow[\\{\\}]$

2:for all agent action $a^{i}\in A_{q}$ do

3: insert all relevant agents of $a^{i}$ to $G$

4: insert all relevant features of $a^{i}$ to $F$

5: insert all relevant action sets of $a^{i}$ to $A$

6:$V\leftarrow\\{\\}$; $\bar{V}\leftarrow\\{\\}$

7:for all joint state $\mathbf{s}\in{\mathcal{S}}$ do

8: for all joint action $\mathbf{a}$ enabled in $\mathbf{s}$ do

9: if $\mathbf{a}$ is compatible with $A_{q}$ [replace $A_{q}$ with $A$] then

10: insert $\mathbf{s}$ to $V$

11: else

12: insert $\mathbf{s}$ to $\bar{V}$

13:$B_{1}\leftarrow$ States2Boolean($V$); $B_{0}\leftarrow$
States2Boolean($\bar{V}$)

14:$\phi\leftarrow$ Quine-McCluskey(ones=$B_{1}$, zeros=$B_{0}$)

15:translate $\phi$ to explanations ${\mathcal{E}}$ via language templates

16:return ${\mathcal{E}}$

17:function State2Boolean($W$)

18: $B\leftarrow\\{\\}$

19: for all $\mathbf{s}=(s^{1},\dots,s^{N})\in W$ do

20: for $1\leq i\leq N$ [replace with $i\in G$] do

21: $C\leftarrow$ feature predicate valuations of state $s^{i}$

22: for $f\in{\mathcal{F}}$ [replace with $f\in F$] do

23: insert $C(f)$ to $B$

24: return $B$

If a state $\mathbf{s}$ has at least one enabled action $\mathbf{a}$ passing
the aforementioned checks, $\mathbf{s}$ is inserted to the target states set
$V$; and to the non-target states set $\bar{V}$ otherwise. The intuition is
that the generated explanations should describe target states satisfying the
query criteria and exclude conditions of non-target states. WithRF poses
further restrictions that target states need to satisfy criteria captured by
relevant actions $A$, such as agent cooperation requirements. Thus,
explanations generated by WithRF can provide information about agent
cooperation, which may be missed by NoRF explanations.

Next, the algorithm converts the states set $V$ (resp. $\bar{V}$) to a list of
Boolean formulas $B_{1}$ (resp. $B_{0}$) via the function described in line
17-24. Given a joint state $\mathbf{s}=(s^{1},\dots,s^{N})$, NoRF finds
valuations of every feature predicates $f\in{\mathcal{F}}$ for all agent state
$s^{i}$ and insert them to the list $B$. By contrast, WithRF only inserts to
$B$ the valuations of relevant features $f\in F$ in relevant agent states
$s^{i}$ for all $i\in G$. Following the previous example, WithRF only
considers Boolean formulas related to relevant feature predicates
{$\mathsf{victim\\_detect,victim\\_complete}$}, filtering out features related
to the fire and obstacle.

Lastly, the algorithm supplies Boolean formulas $B_{1}$ and $B_{0}$ to the
Quine-McCluskey algorithm Quine (1952) and obtains a minimized Boolean
formula, which can be translated into language explanations following Hayes
and Shah (2017). The runtime of Quine-McCluskey grows exponentially with the
number of variables. Thus, WithRF is more efficient than NoRF, due to the
decreased number of Boolean variables. Moreover, filtering out irrelevant
agents and features helps WithRF to prevent redundant information in the
generated explanations.

Properties. Following the Quine-McCluskey, the complexity of NoRF is bounded
by
${\mathcal{O}}\big{(}3^{N\cdot|{\mathcal{F}}|}/\ln(N\cdot|{\mathcal{F}}|)\big{)}$.
The complexity of WithRF is reduced to
${\mathcal{O}}\big{(}3^{|G|\cdot|F|}/\ln(|G|\cdot|F|)\big{)}$.

Query |  Explanations generated by NoRF (baseline) |  Explanations generated by WithRF (proposed)
---|---|---
When does UAV rescue the victim? |  UAV rescues the victim when UAV detects the victim _and_ UGV1 does not detect the fire, _or_ UAV detects the victim _and_ UGV1 does not detect the obstacle. |  UAV rescues the victim when UAV detects the victim _and_ UGV1 detects the victim, _or_ UAV detects the victim _and_ UGV2 detects the victim.
Why don’t UGV1 and UGV2 remove the obstacle in this state? |  UGV1 and UGV2 don’t remove the obstacle in this state because UGV1 does not detect the obstacle. |  UGV1 and UGV2 don’t remove the obstacle in this state because UGV1 does not detect the obstacle _and_ UGV2 does not detect the obstacle.
What does UAV do when it detects the victim? |  UAV can rescue the victim, move, or wait when it detects the victim. |  UAV is most likely to rescue the victim when it detects the victim.
Table 1: Examples of query-based explanations

###### Example 3

Table 1 (first row) shows the explanations generated by NoRF and WithRF for a
when query. The NoRF explanation contains redundant information about the fire
and obstacle that are irrelevant to the query. The WithRF explanation
completely captures the required agent cooperation for the query, which is
missed by the NoRF explanation.

###  4.2 Explanations for Why Not Query

The query “Why don’t agents $G_{q}$ do actions $A_{q}$ in the joint State
$\mathbf{s}_{q}$?” can be answered by modifying Algorithm 2 as follows. In
line 10, adding $\mathbf{s}$ to $\bar{V}$ instead of $V$. Remove line 11-12
and add a new line for inserting the query state $\mathbf{s}_{q}$ to $V$. The
modified algorithm (see Appendix A) generates an explanation describing the
differences between the observed behavior in the target query state
$\mathbf{s}_{q}$ and the expected behavior of states with actions compatible
with the query actions $A_{q}$ (NoRF) or relevant action sets $A$ (WithRF).
The complexity of the modified algorithm follows Algorithm 2.

###### Example 4

Table 1 (second row) shows the explanations generated by NoRF and WithRF for a
why not query about the behavior of two agents UGV1 and UGV2 in the state
shown in Figure 1(a). The WithRF explanation captures the required agent
cooperation for removing the obstacle, while the NoRF explanation fails to
provide such information.

###  4.3 Explanations for What Query

To answer the query “What do agents $G_{q}$ do when satisfying predicates
$F_{q}$?”, we first identify all the satisfying joint states
$\mathbf{s}=(s^{1},\dots,s^{N})$; that is, for all $i\in G_{q}$, agent state
$s^{i}$ satisfies predicates $F_{q}$. The baseline NoRF method is to generate
a list of all possible enabled actions for agents $G_{q}$ in these states. The
proposed WithRF method improves the baseline by filtering agent actions that
are relevant to predicates $F_{q}$ and finding the most likely relevant agent
actions from the list via frequency counting. Since the proposed methods111See
Appendix A for the pseudos code.  do not need to call the Quine-McCluskey, the
complexity of both NoRF and WithRF are only bounded by
${\mathcal{O}}(|G_{q}|\cdot|{\mathcal{S}}|\cdot|{\mathcal{A}}|)$, depending on
the number of query agents, joint state space, and joint action space of the
policy abstraction.

###### Example 5

Table 1 (third row) shows the explanations generated by NoRF and WithRF for a
what query. The WithRF explanation is more concise than the NoRF explanation
and only contains relevant action for the query predicate.

##  5 Computational Experiments

Case Study |  |  | MMDP |  |  | Summarization |  |  | When Query |  |  | Why Not Query |  |  | What Query
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  |  | Path |  | Chart |  |  | NoRF |  | WithRF |  |  | NoRF |  | WithRF |  |  | NoRF |  | WithRF
Domain | $N$ |  |  | $|{\mathcal{S}}|$ | $|{\mathcal{T}}|$ |  |  | $|\rho|$ |  | $|{\mathcal{Z}}|$ |  |  | $|{\mathcal{E}}|$ | Time (ms) |  | $|{\mathcal{E}}|$ | Time (ms) |  |  | $|{\mathcal{E}}|$ | Time (ms) |  | $|{\mathcal{E}}|$ | Time (ms) |  |  | $|{\mathcal{E}}|$ | Time (ms) |  | $|{\mathcal{E}}|$ | Time (ms)
| 3 |  |  | 63 | 577 |  |  | 8 |  | 3x2 |  |  | 4 | 72.4 |  | 4 | 1.2 |  |  | 1 | 86.8 |  | 2 | 0.8 |  |  | 3 | 1.7 |  | 1 | 1.3
SR | 4 |  |  | 732 | 5,048 |  |  | 23 |  | 4x4 |  |  | - | timeout |  | 6 | 31.4 |  |  | - | timeout |  | 3 | 14.9 |  |  | 3 | 15.7 |  | 1 | 14.5
| 5 |  |  | 839 | 4,985 |  |  | 24 |  | 5x4 |  |  | - | timeout |  | 10 | 73.6 |  |  | - | timeout |  | 4 | 23.4 |  |  | 6 | 24.6 |  | 1 | 19.8
| 2 |  |  | 8 | 62 |  |  | 5 |  | 2x2 |  |  | 4 | 2.6 |  | 4 | 2.7 |  |  | 4 | 2.0 |  | 2 | 0.4 |  |  | 3 | 0.1 |  | 1 | 0.1
RWARE | 4 |  |  | 16 | 387 |  |  | 5 |  | 4x1 |  |  | 12 | 3,321.0 |  | 5 | 82.7 |  |  | 6 | 23.0 |  | 3 | 1.3 |  |  | 3 | 0.1 |  | 1 | 0.1
| 19 |  |  | 114 | 1,500 |  |  | 15 |  | 19x2 |  |  | - | timeout |  | 7 | 3,890.6 |  |  | - | timeout |  | 4 | 56.5 |  |  | 3 | 21.7 |  | 1 | 20.1
| 2 |  |  | 5 | 13 |  |  | 4 |  | 2x2 |  |  | 2 | 0.7 |  | 2 | 0.7 |  |  | 3 | 0.4 |  | 2 | 0.3 |  |  | 2 | 0.1 |  | 1 | 0.1
LBF | 4 |  |  | 15 | 355 |  |  | 5 |  | 4x1 |  |  | 10 | 3,598.5 |  | 2 | 1.4 |  |  | 8 | 3,695.0 |  | 2 | 1.6 |  |  | 3 | 0.6 |  | 1 | 0.6
| 9 |  |  | 482 | 5,841 |  |  | 13 |  | 9x2 |  |  | - | timeout |  | 2 | 200.2 |  |  | - | timeout |  | 2 | 101.3 |  |  | 2 | 24.9 |  | 1 | 19.8

Table 2: Experimental results on three MARL domains (timeout set as one hour).

We implemented and applied the proposed methods to three MARL domains. The
first domain is multi-robot search and rescue (SR) similar to Example 1. The
second and third domains are benchmarks taken from Papoudakis et al. (2021).
Multi-robot warehouse (RWARE) considers multiple robotic agents cooperatively
delivering requested items. Level-based foraging (LBF) considers a mixed
cooperative-competitive game where agents must navigate a grid world to
collect randomly scattered food. Our implementation used the Shared Experience
Actor-Critic Christianos et al. (2020) for MARL policy training and
evaluation. All models were trained and evaluated to 10,000 steps, or until
converging to the expected reward, whichever occurred first. The experiments
were run on a laptop with a 1.4 GHz Quad-Core Intel i5 processor and 8 GB RAM.

Table 2 shows experimental results. For each domain, we report the number of
agents $N$ and the number of states $|{\mathcal{S}}|$ and transitions
$|{\mathcal{T}}|$ of generated policy abstraction MMDP. It is unsurprising
that the MMDP size grows exponentially with the number of agents. We report
the most probable path length $|\rho|$ and the chart size $|{\mathcal{Z}}|$ of
generated policy summarizations, which are more compact and easier to
interpret than complex MMDP abstractions. All summarizations were generated
within 1 second (thus not shown in the table). Additionally, we compare NoRF
and WithRF methods in terms of the number of clauses in the generated query-
based explanations and runtime. The results show that WithRF is more succinct
in general and has better scalability than NoRF. In particular, NoRF failed to
generate explanations for “when” and “why not” queries within an hour for
large cases with more than 4 agents, while WithRF generated explanations for
all cases within seconds. Both methods generate explanations for “what”
queries efficiently, thanks to the lower complexity than other queries (see
Section 4).

In summary, experimental results demonstrate that the proposed methods can
generate policy summarizations and query-based explanations for large MARL
domains (e.g., RWARE with 19 agents, which is the largest number of possible
agents in the provided environments).

##  6 User Study

We conducted a user study 222This study was approved by University of Virginia
Institutional Review Boards IRB-SBS #4701. Study details (e.g., questionnaire,
interface) are included in Appendix B.  via Qualtrics to evaluate the quality
of generated explanations. We describe the study design in Section 6.1 and
analyze results in Section 6.2.

###  6.1 Study Design

Participants. We recruited 116 eligible participants (i.e., fluent English
speakers over the age of 18) through university mailing lists. 62.1% of
participants self-identified as male, 37.1% as female, and 0.8% preferred not
to say. The age distribution is 76(18-24), 31(25-34), 7(35-49), 2(50-64).
Participants were instructed to answer multiple-choice questions about agent
behavior for multi-robot search and rescue tasks. They were incentivized with
bonus payments to answer questions correctly based on the provided
explanations. To ensure data quality, attention checks were injected during
the study.

Independent variables. We employed a within-subject study design with the
explanation generation methods as independent variables. Participants were
asked to complete two trials for evaluating policy summarizations. They were
presented with charts generated by Algorithm 1 in one trial, and GIF
animations illustrating the most probable sequence of agent behavior (i.e.,
visualization of the most probable path in the policy abstraction) in the
other trial. For each trial, there were two questions about agent behavior in
various environments (i.e., 3$\times$6 and 6$\times$6 grid world). Questions
used in the two trials are different but had similar difficulty. All
participants were presented with the same set of four randomly generated
questions for summarization trials. To counterbalance the ordering confound
effect, they were randomly assigned to answer the first two questions based on
either charts or GIF, and the other two questions based on the remaining
method. Additionally, participants were asked to complete two trials for
evaluating query-based explanations generated by NoRF and WithRF methods, with
6 questions (2 environments $\times$ 3 query types) in each trial.
Participants answered the same set of 12 randomly generated questions for
query-based trials, and were randomly assigned to different groups similarly
to summarization trials.

Dependent measures. We measured _user performance_ by counting the number of
correctly answered questions in each trial. In addition, at the end of each
trial, participants were asked to rate in a 5-point Likert scale (1 - strongly
disagree, 5 - strongly agree) about _explanation goodness_ metrics (i.e.,
understanding, satisfaction, detail, completeness, actionability, reliability,
trust) Hoffman et al. (2018).

Hypotheses. We make the following hypotheses in this study.

  * •

H1: Chart-based summarizations lead to better user performance than GIF-based.

  * •

H2: Chart-based summarizations yield higher user ratings on explanation
goodness metrics than GIF-based.

  * •

H3: Query-based explanations generated by WithRF lead to better user
performance than those by NoRF.

  * •

H4: Query-based explanations generated by WithRF yield higher user ratings on
explanation goodness metrics than those by NoRF.

###  6.2 Results Analysis

Figure 2: Mean and SD of participant ratings about policy summarizations (“*”
indicates statistically significant difference).

We used a paired t-test to evaluate hypotheses H1 and H3, and used the
Wilcoxon Signed-rank test to evaluate hypotheses H2 and H4. We set the
significant level as $\alpha=0.05$.

Evaluating policy summarizations. Participants answered more questions
correctly with chart-based summarizations (M=1.8 out of 2, SD=0.6) than GIF-
based (M=0.9 out of 2, SD=0.4), with statistical significance ($t$(462)=-15.8,
$p\leq$0.01, $d$=1.5). _Thus, the data supports H1._

Figure 2 shows average participant ratings about summarizations. Chart-based
summarizations yield higher ratings on the perceived completeness than GIF-
based with statistical significance ($W$=371.5, $Z$=-2.4, $p\leq$0.02,
$r$=-0.2). But no significant difference was found regarding the other
metrics. _Thus, the data partially supports H2._

Evaluating query-based explanations. Participants answered more questions
correctly with explanations generated by WithRF (M=5.2 out of 6, SD=1.7) than
NoRF (M=2.3 out of 6, SD=1.0), with statistical significance ($t$(1390)=-21.1,
$p\leq$0.01,$d$=2.0). _Thus, the data supports H3._

Figure 3 shows that participants gave higher average ratings to WithRF
explanations than NoRF explanations. The Wilcoxon test found significant
differences on all metrics: understanding ($W$=319.5, $Z$=-4.9, $p\leq$0.01,
$r$=-0.3), satisfaction ($W$=266.0, $Z$=-7.0, $p\leq$0.01, $r$=-0.5), detail
($W$=484.0, $Z$=-3.7, $p\leq$0.01, $r$=-0.2), completeness ($W$=494.5,
$Z$=-6.4, $p\leq$0.01, $r$=-0.4), actionability ($W$=167.0, $Z$=-6.9,
$p\leq$0.01, $r$=-0.5), reliability ($W$=382.5, $Z$=-3.6, $p\leq$0.01,
$r$=-0.2), and trust ($W$=217.0, $Z$=-3.4, $p\leq$0.01, $r$=-0.2). _Thus, the
data supports H4._

Figure 3: Mean and SD of participant ratings about query-based explanations
(“*” indicates statistically significant difference).

Discussion. In summary, the data supports all hypotheses, while H2 is only
partially supported because the statistical test found no significant
differences between chart-based and GIF-based summarizations on most metrics.
However, Figure 2 shows that participants rated chart-based summarizations
close to 4 (agree) on all metrics, and above GIF-based ratings on all metrics
except understanding, reliability, and trust. This may be because users showed
a strong preference toward the moving nature of GIF animations and the
visualized effects of agents completing tasks. But watching a GIF can be more
time-consuming and less informative than a quick glance at the chart. This is
supported by the results that participants were able to answer more questions
correctly with chart-based summarizations, and they rated this method
significantly higher on completeness (i.e., providing needed information).
Meanwhile, query-based explanations generated by the proposed WithRF method
led to significantly better user performance and higher user ratings on all
metrics, because users prefer succinct WithRF explanations with adequate
information about agent behavior and cooperation. By contrast, NoRF
explanations do not necessarily provide essential information about agent
cooperation for correctly answering questions, and may contain redundant
information that decreases user satisfaction.

##  7 Conclusion

In this work, we developed methods to generate policy summarizations and
query-based explanations for MARL. Experimental results on three MARL domains
demonstrate the scalability of our methods. Evaluation via a user study shows
that our generated MARL policy explanations can improve user understanding
about agent behavior and enable them to answer more questions correctly, while
maintaining very positive ratings on explanation goodness metrics.

The proposed methods are independent from the inner working of MARL methods,
only relying on policy abstractions that can be built via observing samples
during the policy evaluation. Although we only used one MARL method in our
experiments, the proposed methods can be used to generate policy explanations
for different MARL methods. As part of the future work, we plan to apply the
proposed methods to a wide range of MARL methods and domains.

## Acknowledgments

This work was supported in part by U.S. National Science Foundation grant
CCF-1942836, U.S. Office of Naval Research grant N00014-18-1-2829, Israel
Science Foundation grant 1958/20, EU Project TAILOR grant 952215, and by the
Data Science Institute at Bar-Ilan University. Any opinions, findings, and
conclusions or recommendations expressed in this material are those of the
author(s) and do not necessarily reflect the views of the grant sponsors.

## References

  * Amir and Amir [2018] Dan Amir and Ofra Amir.  Highlights: Summarizing agent behavior to people.  In Proceedings of the 17th International Conference on Autonomous Agents and Multi Agent Systems, pages 1168–1176, 2018.
  * Chakraborti et al. [2020] Tathagata Chakraborti, Sarath Sreedharan, and Subbarao Kambhampati.  The emerging landscape of explainable automated planning & decision making.  In IJCAI, pages 4803–4811, 2020.
  * Christianos et al. [2020] Filippos Christianos, Lukas Schäfer, and Stefano V Albrecht.  Shared experience actor-critic for multi-agent reinforcement learning.  In Thirty-fourth Conference on Neural Information Processing Systems, pages 10707–10717. Curran Associates Inc, 2020.
  * Dafoe et al. [2020] Allan Dafoe, Edward Hughes, Yoram Bachrach, Tantum Collins, Kevin R McKee, Joel Z Leibo, Kate Larson, and Thore Graepel.  Open problems in cooperative ai.  arXiv preprint arXiv:2012.08630, 2020.
  * Dijkstra et al. [1959] Edsger W Dijkstra, others, et al.  A note on two problems in connexion with graphs.  Numerische mathematik, 1(1):269–271, 1959.
  * Hayes and Shah [2017] Bradley Hayes and Julie A Shah.  Improving robot controller transparency through autonomous policy explanation.  In 2017 12th ACM/IEEE International Conference on Human-Robot Interaction (HRI, pages 303–312. IEEE, 2017.
  * Heuillet et al. [2021] Alexandre Heuillet, Fabien Couthouis, and Natalia Díaz-Rodríguez.  Explainability in deep reinforcement learning.  Knowledge-Based Systems, 214:106685, 2021.
  * Hoffman et al. [2018] Robert R Hoffman, Shane T Mueller, Gary Klein, and Jordan Litman.  Metrics for explainable ai: Challenges and prospects.  arXiv preprint arXiv:1812.04608, 2018.
  * Kazhdan et al. [2020] Dmitry Kazhdan, Zohreh Shams, and Pietro Liò.  Marleme: A multi-agent reinforcement learning model extraction library.  In 2020 International Joint Conference on Neural Networks (IJCNN), pages 1–8. IEEE, 2020.
  * Kiran et al. [2021] B Ravi Kiran, Ibrahim Sobh, Victor Talpaert, Patrick Mannion, Ahmad A Al Sallab, Senthil Yogamani, and Patrick Pérez.  Deep reinforcement learning for autonomous driving: A survey.  IEEE Transactions on Intelligent Transportation Systems, 2021.
  * Kraus et al. [2020] Sarit Kraus, Amos Azaria, Jelena Fiosina, Maike Greve, Noam Hazon, Lutz Kolbe, Tim-Benjamin Lembcke, Jorg P Muller, Soren Schleibaum, and Mark Vollrath.  Ai for explaining decisions in multi-agent environments.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 13534–13538, 2020.
  * Papoudakis et al. [2021] Georgios Papoudakis, Filippos Christianos, Lukas Schäfer, and Stefano V Albrecht.  Benchmarking multi-agent deep reinforcement learning algorithms in cooperative tasks.  In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 1), 2021.
  * Puiutta and Veith [2020] Erika Puiutta and Eric MSP Veith.  Explainable reinforcement learning: A survey.  In International Cross-Domain Conference for Machine Learning and Knowledge Extraction, pages 77–95. Springer, 2020.
  * Quine [1952] Willard V Quine.  The problem of simplifying truth functions.  The American mathematical monthly, 59(8):521–531, 1952.
  * Topin and Veloso [2019] Nicholay Topin and Manuela Veloso.  Generation of policy-level explanations for reinforcement learning.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 2514–2521, 2019.
  * Wells and Bednarz [2021] Lindsay Wells and Tomasz Bednarz.  Explainable ai and reinforcement learning—a systematic review of current approaches and trends.  Frontiers in artificial intelligence, 4:48, 2021.

##  Appendix A Algorithms

The text in blue highlight changes about relevancy filters (RF) for the
proposed WithRF method compared to the baseline NoRF method.

Algorithm 3 Answering “why not” query

Input: policy abstraction
${\mathcal{M}}=({\mathcal{S}},{\mathcal{A}},{\mathcal{T}})$, query “why don’t
agents $G_{q}$ do actions $A_{q}$ in the joint State $\mathbf{s}_{q}$”
Output: language explanations ${\mathcal{E}}$

1:$G\leftarrow\\{\\}$; $F\leftarrow\\{\\}$; $A\leftarrow[\\{\\}]$

2:for all agent action $a^{i}\in A_{q}$ do

3: insert all relevant agents of $a^{i}$ to $G$

4: insert all relevant features of $a^{i}$ to $F$

5: insert all relevant action sets of $a^{i}$ to $A$

6:$V\leftarrow\\{\\}$; $\bar{V}\leftarrow\\{\\}$

7:insert $\mathbf{s}_{q}$ to $V$

8:for all joint state $\mathbf{s}\in{\mathcal{S}}$ do

9: for all joint action $\mathbf{a}$ enabled in $\mathbf{s}$ do

10: if $\mathbf{a}$ is compatible with $A_{q}$ [replace $A_{q}$ with $A$] then

11: insert $\mathbf{s}$ to $\bar{V}$

12:$B_{1}\leftarrow$ States2Boolean($V$); $B_{0}\leftarrow$
States2Boolean($\bar{V}$)

13:$\phi\leftarrow$ Quine-McCluskey(ones=$B_{1}$, zeros=$B_{0}$)

14:translate $\phi$ to explanations ${\mathcal{E}}$ via language templates

15:return ${\mathcal{E}}$

16:function State2Boolean($W$)

17: $B\leftarrow\\{\\}$

18: for all $\mathbf{s}=(s^{1},\dots,s^{N})\in W$ do

19: for $1\leq i\leq N$ [replace with $i\in G$] do

20: $C\leftarrow$ feature predicate valuations of state $s^{i}$

21: for $f\in{\mathcal{F}}$ [replace with $f\in F$] do

22: insert $C(f)$ to $B$

23: return $B$

Algorithm 4 Answering “what” query

Input: policy abstraction
${\mathcal{M}}=({\mathcal{S}},{\mathcal{A}},{\mathcal{T}})$, query “what do
agents $G_{q}$ do when satisfying predicates $F_{q}$”
Output: language explanations ${\mathcal{E}}$

1:$A_{q}\leftarrow\\{\\}$

2:$\alpha\leftarrow$ find all relevant actions of predicates $F_{q}$

3:for all joint state $\mathbf{s}=(s^{1},\dots,s^{N})\in{\mathcal{S}}$ do

4: if $s^{i}$ satisfies $F_{q}$ for all $i\in G_{q}$ then

5: for all $\mathbf{a}=(a^{1},\dots,a^{N})$ enabled in $\mathbf{s}$ do

6: for all $i\in G_{q}$ do

7: insert $a^{i}$ to $A_{q}$ [only if $a^{i}\in\alpha$]

8:$A_{q}\leftarrow$ the most frequent action for each agent in $A_{q}$

9:generate explanations ${\mathcal{E}}$ with $A_{q}$ via language templates

10:return ${\mathcal{E}}$

##  Appendix B User Study Details

Questionnaire on explanation goodness metrics. Participants were instructed to
rate on a 5-point Likert scale (1 - strongly disagree, 5 - strongly agree)
about the following statements, which were adapted from Hoffman et al. [2018].

  * •

The explanations help me _understand_ how the team of robots completes the
search and rescue mission.

  * •

The explanations are _satisfying_.

  * •

The explanations are sufficiently _detailed_.

  * •

The explanations are sufficiently _complete_ , that is, they provide me with
all the needed information to answer the questions.

  * •

The explanations are _actionable_ , that is, they help me know how to answer
the questions.

  * •

The explanations let me know how _reliable_ the robot team is for completing
the mission.

  * •

The explanations let me know how _trustworthy_ the robot team is for
completing the mission.

User interface and sample questions. Figures 4-8 show examples of user
interface and questions presented to participants during the user study.

Figure 4: Question based on explanations for a “when” query. Figure 5:
Question based on explanations for a “why not” query. Figure 6: Question based
on explanations for a “what” query. Figure 7: Question based on policy
summarization (sequence chart). Figure 8: Question based on a policy
summarization (GIF animation:
https://github.com/kjboggess/IJCAI2022/blob/main/MissionGifExample.gif


# Robust nonequilibrium surface currents in the 3D Hofstadter model

Mark T. Mitchison<EMAIL_ADDRESS>School of Physics, Trinity College
Dublin, College Green, Dublin 2, Ireland  Ángel Rivas<EMAIL_ADDRESS>Miguel
A. Martin-Delgado<EMAIL_ADDRESS>Departamento de Física Teórica, Facultad de
Ciencias Físicas, Universidad Complutense, 28040 Madrid, Spain.  CCS-Center
for Computational Simulation, Campus de Montegancedo UPM, 28660 Boadilla del
Monte, Madrid, Spain.

###### Abstract

Genuinely two-dimensional robust crosscurrents — which flow against the
natural direction of heat flux — have been missing since the discovery of
their one-dimensional counterpart. We provide a setup to realize them on a
cubic three-dimensional (3D) lattice hosting a Hofstadter model coupled to two
heat baths with different temperatures. We show that these currents exhibit
dissipative robustness: they are stable against the presence of impurities and
tilting of the gauge field in certain nonequilibrium configurations. Moreover,
we find protected boundary currents with genuinely 3D robustness, i.e. they
are only stable if tunnelling can occur in all three spatial directions. The
model also presents generic surface currents, which are robust for both
bosonic and fermionic systems. We identify the underlying qualitative
mechanism responsible for the robustness of the surface currents and the
crucial role played by certain discrete symmetries.

##  I Introduction

The fascinating world of exotic quantum phenomena in condensed matter physics
has traditionally been associated with low-dimensional systems, i.e. those
existing in one and two dimensions [1, 2, 3]. This has been true not only for
quantum phase transitions, but also for transport phenomena [4, 5, 6]. The
advent of topological quantum computing on the one hand [7, 8, 9, 10], and
topological insulators and superconductors on the other [11, 12, 13], provides
the opportunity to find captivating new quantum properties in systems in three
dimensions and even formally in higher dimensions [14, 15]. A particularly
interesting avenue of research in this direction concerns the emergence of
such phenomena in open quantum systems far from equilibrium [16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27], since these may exhibit novel behaviour that
cannot occur in closed systems.

In this paper, we address an open problem in the study of transport in three-
dimensional lattices that are out of thermodynamic equilibrium due to coupling
to thermal baths at different temperatures. The subject of our study will be
the current created by such a temperature gradient when the system has reached
the steady state, so that its entropy remains constant. For a two-dimensional
(2D) Hofstadter lattice model of bosons, a one-dimensional edge current was
recently found flowing in the opposite direction to the natural arrow of heat
flow, and for this reason is called a crosscurrent [20]. This exotic edge
current, which has also been observed in other 2D models [28] (see Refs. [29,
30, 31] for quasi-one-dimensional studies), is robustly protected by symmetry
properties with respect to the presence of point-like defects, and remains
stable for a wide range of reservoir coupling strengths [28]. Such
crosscurrents were first thought to be an intrinsically one-dimensional
phenomenon. Since then, the possibility of observing a truly two-dimensional
crosscurrent that is stable in the presence of dissipation has remained open.
Here, we fill this gap by constructing an explicit realization of a 2D
crosscurrent as a boundary current of a three-dimensional (3D) Hofstadter
lattice (see Fig. 1). Moreover, we demonstrate the existence of robust surface
currents that are stable with respect to defects in the presence of
dissipation, regardless whether they are crosscurrents or not.

Figure 1: Schematic of the 3D Hofstadter model coupled to hot and cold
reservoirs, with respective temperatures $T_{h}>T_{c}$, along the $X$-axis.
The coefficients $t_{X,Y,Z}$ indicate tunneling amplitudes along the three
spatial directions and each face of the unit cell is crossed by the magnetic
fluxes $b_{X}$, $b_{Y}$ and $b_{Z}$ depending on orientation of the gauge
field.

The boundary currents we consider here are descendants of the exotic surface
physics that appears in 3D topological insulators [11, 12, 13, 32, 33, 34, 35,
36, 37, 38], the main difference being that the latter are closed quantum
fermionic systems as opposed to the open bosonic and fermionic systems that we
consider. Typically, these 3D topological models can be divided into two
generic categories. One hosts a layered structure of connected 2D topological
nontrivial systems, in such a way that their exotic transport properties in 3D
are inherited from the nontrivial 2D properties of each layer. However, there
also exist systems with robust boundary states which are of genuinely 3D
origin. For time-reversal invariant (TRI) models, these two classes correspond
to the so-called “weak” and “strong” 3D topological insulators [12, 32],
respectively. For the case of 3D lattices with broken time-reversal (TRB)
symmetry, layered models are easy to construct, the 3D integer Hall effect
being the simplest instance [39, 40, 41, 42, 43, 44, 45]. However, observing
truly 3D robust boundary states for TRB systems is more challenging, and the
most studied model, the so-called Hopf insulator [46, 47, 48], presents
topologically protected boundary modes but its topological features are
unstable to the emergence of extra energy bands.

By extending the ideas previously developed in Refs. [20, 28] for 2D
topological systems coupled to thermal baths, in this work we show that it is
indeed possible to induce genuinely 3D robust boundary transport in TRB
systems. Namely, the robust surface currents we observe arise from a genuinely
three-dimensional effect and are not the result of stacking two-dimensional
models with corresponding one-dimensional robust edge currents. In this
regard, we obtain the following notable results:

  * (i)

We find robustness of boundary currents, allowing for a surface crosscurrent,
under perturbations created by the presence of a variety of defects in the 3D
lattice, including point-, surface- and volume-like defects.

  * (ii)

This robustness is induced by the fulfillment of certain symmetries in such
perturbations.

  * (iii)

For particular gauge field orientations, there are 2D layered symmetries which
can stabilize an unstable layered current (or vice versa), creating a 3D
stabilization effect.

  * (iv)

Moreover, there is a genuine 3D symmetry independent of any field orientation,
which protects the surface currents without any reference to their layered
structure.

In the next section, we introduce our 3D lattice model coupled to thermal
baths and detail the solution for its non-equilibrium steady state (NESS).
Following Ref. [28], we focus on the regime of weak system-bath coupling where
the boundary currents become dominant. In this regime, the NESS is well
approximated by the solution of a Lindblad master equation. In Sec. III, we
describe our results for the nonequilibrium current distribution obtained
within this framework. Sec. IV is devoted to our conclusions.

##  II Model

The Hofstadter model in three dimensions describes a cubic lattice of
$L_{X}\times L_{Y}\times L_{Z}$ sites in the presence of a gauge field,
governed by the Hamiltonian
$\hat{H}_{S}=\hat{H}_{0}+\hat{H}_{X}+\hat{H}_{Y}+\hat{H}_{Z}$, with
($\hbar=k_{B}=1$)

$\displaystyle\hat{H}_{0}$
$\displaystyle=\sum_{x,y,z}\omega_{0}\hat{a}^{\dagger}_{x,y,z}\hat{a}_{x,y,z},$
(1) $\displaystyle\hat{H}_{X}$
$\displaystyle=\frac{-t_{X}}{2}\sum_{x,y,z}\mathrm{e}^{\mathrm{i}2\pi
b_{Z}y}\hat{a}^{\dagger}_{x,y,z}\hat{a}_{x+1,y,z}+\rm h.c.,$ (2)
$\displaystyle\hat{H}_{Y}$
$\displaystyle=\frac{-t_{Y}}{2}\sum_{x,y,z}\mathrm{e}^{\mathrm{i}2\pi
b_{X}z}\hat{a}^{\dagger}_{x,y,z}\hat{a}_{x,y+1,z}+\rm h.c.,$ (3)
$\displaystyle\hat{H}_{Z}$
$\displaystyle=\frac{-t_{Z}}{2}\sum_{x,y,z}\mathrm{e}^{\mathrm{i}2\pi
b_{Y}x}\hat{a}^{\dagger}_{x,y,z}\hat{a}_{x,y,z+1}+\rm h.c.$ (4)

Here, $\hat{a}_{x,y,z}^{\dagger}$ creates a particle on lattice site
$(x,y,z)$, $t_{X,Y,Z}$ denote the hopping amplitudes in the three orthogonal
directions, and $\mathbf{b}=(b_{X},b_{Y},b_{Z})$ is an effective “magnetic”
flux vector on the three faces of the unit cell, as illustrated in Fig. 1, as
a result of the presence of the gauge field. Unless otherwise indicated, we
set $t_{X}=t_{Y}=t_{Z}=t$. The on-site energy shift $\omega_{0}$ is chosen to
be large enough to ensure that all single-particle eigenenergies are positive,
but its value is otherwise immaterial: we take $\omega_{0}=10t$ in the
following.

The system lattice is coupled to independent reservoirs, hot at one end,
$x=1$, and cold at the other end, $x=L_{X}$, with temperatures
$T_{h}=\beta_{h}^{-1}$ and $T_{c}=\beta_{c}^{-1}$ and chemical potentials
$\mu_{h}$ and $\mu_{c}$, respectively. We consider a system-reservoir
interaction of the bilinear form

$\displaystyle\hat{H}_{SB}=\sum_{j\in\partial
S}\sum_{q}g_{qj}(\hat{a}^{\dagger}_{j}\hat{b}_{qj}+\hat{b}_{qj}^{\dagger}\hat{a}_{j}),$
(5)

with $\partial S$ the boundary of the system, and $g_{qj}=0$ if
$j\notin\\{x=1\\}\cup\\{x=L_{X}\\}$. Here, $\hat{b}_{qj}^{\dagger}$ creates a
particle in mode $q$ of the reservoir coupled to the lattice site $j$. We
assume that system and reservoir operators, $\hat{a}$ and $\hat{b}$, satisfy
the same canonical commutation or anticommutation relations, depending on
whether we consider bosonic or fermionic systems, respectively. The total
Hamiltonian is given by

$\displaystyle\hat{H}=\hat{H}_{S}+\hat{H}_{SB}+\hat{H}_{B},$ (6)

with

$\displaystyle\hat{H}_{B}=\sum_{j\in\partial
S}\sum_{q}\Omega_{qj}\hat{b}^{\dagger}_{qj}\hat{b}_{qj}.$ (7)

For sufficiently long times the lattice system reaches a non-equilibrium
steady state (NESS). Assuming that the coupling between system and reservoir
is weak enough, and following the same steps as in [20] for the 2D case within
the Born-Markov and secular approximations [49, 50, 51], the NESS of the
Hofstadter system can be well approximated by the stationary solution of a
Gorini-Kossakowski-Lindblad-Sudarshan master equation

$\mathcal{L}\hat{\rho}=-\mathrm{i}[\hat{H},\hat{\rho}]+\mathcal{L}_{h}\hat{\rho}+\mathcal{L}_{c}\hat{\rho}=0.$
(8)

The dissipators for the hot and cold baths are given respectively by

$\displaystyle\mathcal{L}_{h}=\gamma\sum_{\alpha}s_{\alpha}\left(\bar{n}_{h}(\omega_{\alpha})\mathcal{D}[\hat{c}^{\dagger}_{\alpha}]+[1\mp\bar{n}_{h}(\omega_{\alpha})]\mathcal{D}[\hat{c}_{\alpha}]\right),$
(9)
$\displaystyle\mathcal{L}_{c}=\gamma\sum_{\alpha}r_{\alpha}\left(\bar{n}_{c}(\omega_{\alpha})\mathcal{D}[\hat{c}^{\dagger}_{\alpha}]+[1\mp\bar{n}_{c}(\omega_{\alpha})]\mathcal{D}[\hat{c}_{\alpha}]\right),$
(10)

with
$\mathcal{D}[\hat{L}]\bullet=\hat{L}\bullet\hat{L}^{\dagger}-\tfrac{1}{2}\\{\hat{L}^{\dagger}\hat{L},\bullet\\}$
and
$\bar{n}_{h,c}(\omega)=[\mathrm{e}^{\beta_{h,c}(\omega-\mu_{h,c})}\mp(-1)]^{-1}$,
where the minus sign is for fermions and the plus sign is for bosons.
Furthermore, $\hat{c}_{\alpha}=\sum_{j}U_{j\alpha}\hat{a}_{j}$ denotes a
canonical ladder operator of the diagonalised 3D Hofstadter Hamiltonian,
$\hat{H}_{S}=\sum_{\alpha}\omega_{\alpha}\hat{c}^{\dagger}_{\alpha}\hat{c}_{\alpha}$,
where $U_{j\alpha}$ is the unitary matrix comprising the one-particle
eigenvectors of $\hat{H}_{S}$. For simplicity, we have assumed that the
reservoir spectral density is constant within the frequency range of interest
(our results are independent of this choice). The coefficients $s_{\alpha}$
and $r_{\alpha}$ describe the dimensionless coupling strength of eigenmode
$\alpha$ to the hot and cold baths, respectively:

$\displaystyle
s_{\alpha}=\sum_{y=1}^{L_{Y}}\sum_{z=1}^{L_{Z}}|U_{(1,y,z),\alpha}|^{2},\quad
r_{\alpha}=\sum_{y=1}^{L_{Y}}\sum_{z=1}^{L_{Z}}|U_{(L_{X},y,z),\alpha}|^{2}.$
(11)

The current operators are defined via the continuity equation 111Here, the
continuity equation is derived by considering the Heisenberg equations of
motion generated by the microscopic Hamiltonian given in Eq. (6). for the
particle density $\hat{n}_{x,y,z}=\hat{a}^{\dagger}_{x,y,z}\hat{a}_{x,y,z}$
along the three spatial dimensions:

$\displaystyle\hat{J}^{X}_{x,y,z}$
$\displaystyle=\mathrm{i}\frac{t_{X}}{2}\mathrm{e}^{\mathrm{i}2\pi
b_{Z}y}\hat{a}_{x,y,z}^{\dagger}\hat{a}_{x+1,y,z}+\mathrm{h.c.}\,,$
$\displaystyle\hat{J}^{Y}_{x,y,z}$
$\displaystyle=\mathrm{i}\frac{t_{Y}}{2}\mathrm{e}^{\mathrm{i}2\pi
b_{X}z}\hat{a}_{x,y,z}^{\dagger}\hat{a}_{x,y+1,z}+\mathrm{h.c.}\,,$
$\displaystyle\hat{J}^{Z}_{x,y,z}$
$\displaystyle=\mathrm{i}\frac{t_{Z}}{2}\mathrm{e}^{\mathrm{i}2\pi
b_{Y}x}\hat{a}_{x,y,z}^{\dagger}\hat{a}_{x,y,z+1}+\mathrm{h.c.},$ (12)

so that

$\displaystyle\mathrm{i}[\hat{H}_{X},\hat{n}_{x,y,z}]$
$\displaystyle=\hat{J}^{X}_{x-1,y,z}-\hat{J}^{X}_{x,y,z}\,,$
$\displaystyle\mathrm{i}[\hat{H}_{Y},\hat{n}_{x,y,z}]$
$\displaystyle=\hat{J}^{Y}_{x,y-1,z}-\hat{J}^{Y}_{x,y,z}\,,$
$\displaystyle\mathrm{i}[\hat{H}_{Z},\hat{n}_{x,y,z}]$
$\displaystyle=\hat{J}^{Z}_{x,y,z-1}-\hat{J}^{Z}_{x,y,z}.$ (13)

The mean currents can be obtained by computing the correlation matrix
$\braket{\hat{c}^{\dagger}_{\alpha}\hat{c}_{\alpha^{\prime}}}$, which can be
done in a straightforward way by using the master equation (8). The NESS
correlation matrix takes the value [20, 28]

$\displaystyle\braket{\hat{c}^{\dagger}_{\alpha}\hat{c}_{\alpha^{\prime}}}=\delta_{\alpha\alpha^{\prime}}\frac{s_{\alpha}\bar{n}_{h}(\omega_{\alpha})+r_{\alpha}\bar{n}_{c}(\omega_{\alpha})}{s_{\alpha}+r_{\alpha}},$
(14)

so that only diagonal correlations in the eigenmode basis survive. In Ref.
[28] it was shown that, for a 2D model, the solution of the Lindblad master
equation gives an excellent approximation to the current distribution in the
weak-coupling regime. Similarly, therefore, here we expect Eq. (14) to yield
accurate values for the currents so long as $t_{X,Y,Z}\gg\gamma$. Our
approximations are also valid if one of them vanishes, say $t_{Z}=0$, so long
as the other two satisfy this condition, $t_{X,Y}\gg\gamma$.

Figure 2: Bosonic current distribution for a “tilted” field configuration with
$\mathbf{b}=(0,0.1,0.1)$. Arrows indicate the direction of the current, with
the length of the arrows proportional to the magnitude of the current. The
length of the arrows is scaled relative to the largest current value within
each plot. The right panel shows the topmost plane, demonstrating the
existence of 2D cross-currents flowing from the hot reservoir ($x=1$) to the
cold reservoir ($x=L_{X}$).

##  III Results

###  III.1 Symmetries and boundary currents

A perfect 3D Hofstadter lattice, without any defect, presents several
symmetries depending on the orientation of the magnetic flux vector
$\mathbf{b}$. We will be particularly interested in symmetries of $\hat{H}$
which interchange the hot and cold baths relative locations, as they have a
nontrivial action on the NESS. For instance, if $\mathbf{b}=(0,0,b_{z})$, the
Hamiltonian is invariant under the symmetries $\hat{\Theta}\hat{\Sigma}_{yz}$,
$\hat{\Theta}\hat{R}_{y}(\pi)$ and $\hat{\Theta}\hat{R}_{z}(\pi)$. Here,
$\hat{\Theta}$ is the time-reversal operation (the same as complex
conjugation, in this case), $\hat{\Sigma}_{yz}$ is a spatial reflection across
the $yz$-plane cutting the center of the lattice,
$\hat{\Sigma}_{yz}\hat{a}_{x,y,z}\hat{\Sigma}_{yz}^{\dagger}=\hat{a}_{L_{x}+1-x,y,z}$,
and $\hat{R}_{y,z}(\pi)$ are $\pi$-rotations about the $x$ and $y$ axes,
respectively, from the center of the lattice,
$\hat{R}_{y}(\pi)\hat{a}_{x,y,z}\hat{R}_{y}(\pi)^{\dagger}=\hat{a}_{L_{x}+1-x,y,L_{z}+1-z}$,
and
$\hat{R}_{z}(\pi)\hat{a}_{x,y,z}\hat{R}_{z}(\pi)^{\dagger}=\hat{a}_{L_{x}+1-x,L_{y}+1-y,z}$.
These are essentially 2D symmetries embedded in the 3D lattice, and they arise
for particular orientations of $\mathbf{b}$ or particular values of
$t_{X,Y,Z}$ but not for arbitrary directions of the magnetic field.
Nevertheless, for a general $\mathbf{b}$, the Hamiltonian is always invariant
under $\hat{\Theta}\hat{I}$, where $\hat{I}$ is a spatial inversion
transformation with respect to the center of the lattice
$\hat{I}\hat{a}_{\mathbf{r}}\hat{I}^{\dagger}=\hat{a}_{\mathbf{R}-\mathbf{r}}$,
with $\mathbf{r}=(x,y,z)$ and $\mathbf{R}=(L_{x}+1,L_{y}+1,L_{z}+1)$.

The master equation is also invariant under these symmetries, because the
Liouvillian (8) inherits them from the Hamiltonian within the Born-Markov-
secular (BMS) approximation. Of course, the total Liouvillian is not invariant
because it describes the evolution of an open system and these symmetries
involve the time-reversal operation, $\Theta$. Therefore, only a kind of weak
symmetry can be expected, in the sense of Ref. [53]. Each of the
aforementioned symmetries, $\hat{S}$, is a weak symmetry [53] for both the
Hamiltonian part $\mathcal{L}_{H}(\rho):=[H,\rho]$ and dissipative part
$\mathcal{L}_{D}(\rho):=\mathcal{L}(\rho)+\mathrm{i}\mathcal{L}_{H}(\rho)$ of
the Liouvillian: namely,
$\hat{S}\mathcal{L}_{H,D}(\rho)\hat{S}^{-1}=\mathcal{L}_{H,D}(\hat{S}\rho\hat{S}^{-1})$.
Note, however, that the total Liouvillian
$\mathcal{L}=-\mathrm{i}\mathcal{L}_{H}+\mathcal{L}_{D}$ is not invariant:
since $\hat{S}$ is anti-unitary, the Hamiltonian contribution flips sign.
Nevertheless, since $[\mathcal{L}_{H},\mathcal{L}_{D}]=0$ under the BMS
approximation, the NESS is annihilated by both $\mathcal{L}_{H}$ and
$\mathcal{L}_{D}$ individually. Therefore, a weak symmetry $\hat{S}$ of
$\mathcal{L}_{H}$ and $\mathcal{L}_{D}$ is sufficient to ensure that the NESS
is also invariant under $\hat{S}$.

Since the aforementioned weak symmetries interchange the reservoirs, the
invariance of $\hat{H}$ also ensures that the dimensionless couplings (11)
obey $s_{\alpha}=r_{\alpha}$. As a result, the mode occupations given by Eq.
(14) become independent of $s_{\alpha}$ and $r_{\alpha}$, being given simply
by the average of the reservoir distribution functions:

$n(\omega_{\alpha})=\frac{1}{2}\left[\bar{n}_{c}(\omega_{\alpha})+\bar{n}_{h}(\omega_{\alpha})\right].$
(15)

This is a function of the frequency $\omega_{\alpha}$ only, and leads to a
steady-state correlation matrix that is independent of the spatial orientation
of the reservoirs. In such a case, the induced current is highly suppressed in
the bulk and flows essentially on the boundary of the 3D lattice. To
understand this, let us resort to a semi-classical picture [54, 55, 56].

Under periodic boundary conditions, the 3D Hofstadter Hamiltonian can be
written in the form

$\hat{H}=\sum_{\mathbf{k}\in\mathrm{MBZ}}\mathbf{a}_{\mathbf{k}}^{\dagger}\cdot\mathbf{H}(\mathbf{k})\cdot\mathbf{a}_{\mathbf{k}},$
(16)

with $\mathbf{k}=(k_{X},k_{Y},k_{Z})$ the quasimomentum, ranging inside the
magnetic Bruillouin zone (MBZ),
$\mathbf{a}_{\mathbf{k}}^{\dagger}=(\hat{a}_{1\mathbf{k}}^{\dagger},\hat{a}_{2\mathbf{k}}^{\dagger},\ldots)$,
and a hermitian matrix $\mathbf{H}(\mathbf{k})$ whose dimension and structure
depend on $\mathbf{b}$. For instance, for rational fluxes
$\mathbf{b}=(p_{X}/q_{X},p_{Y}/q_{Y},p_{Z}/q_{Z})$, the number of eigenvalues
$\omega_{\alpha}(\mathbf{k})$ of $\mathbf{H}(\mathbf{k})$ (energy bands) is
given by the lowest common denominator of the three fractions [39, 40, 41,
42].

The action of a force $\mathbf{F}$ modifies the semiclassical equation for the
velocity $\mathbf{v}_{\alpha}$ of a Bloch wave-packet in the band
$\omega_{\alpha}(\mathbf{k})$ as

$\mathbf{v}_{\alpha}=\frac{\partial\omega_{\alpha}(\mathbf{k})}{\partial\mathbf{k}}-\mathbf{F}\times\bm{\mathcal{F}}\\!\\!_{\alpha}(\mathbf{k}).$
(17)

Here,
$\bm{\mathcal{F}}\\!\\!_{\alpha}(\mathbf{k})=\bm{\nabla}\times\bm{\mathcal{A}}_{\alpha}$
is the Berry curvature vector, defined in terms of the Berry connection,
$\bm{\mathcal{A}}_{\alpha}(\mathbf{k})=\mathrm{i}\langle
u_{\alpha}(\mathbf{k})|\bm{\nabla}_{\mathbf{k}}|u_{\alpha}(\mathbf{k})\rangle$,
with
$\mathbf{H}(\mathbf{k})|u_{\alpha}(\mathbf{k})\rangle=\omega_{\alpha}(\mathbf{k})|u_{\alpha}(\mathbf{k})\rangle$.

Under these conditions, we can generate a boundary by introducing a confining
potential $V(\mathbf{r})$ in the semiclassical equation (17) for
$\mathbf{v}_{\alpha}$. Assuming that the potential varies slowly on the scale
of the lattice, the energy of the wavepacket in the semi-classical
approximation is simply $\omega_{\alpha}(\mathbf{k})+V(\mathbf{r})$ [55].
Moreover, the potential induces a confining force
$\mathbf{F}=-\bm{\nabla}V(\mathbf{r})$, which approximately vanishes in the
bulk and becomes very large near the edge of the confined region, so that

$\displaystyle\mathbf{v}_{\alpha}^{\mathrm{(bulk)}}\simeq\frac{\partial\omega_{\alpha}(\mathbf{k})}{\partial\mathbf{k}},\quad\mathbf{v}_{\alpha}^{\rm(edge)}\simeq\bm{\nabla}V(\mathbf{r})\times\bm{\mathcal{F}}\\!\\!_{\alpha}(\mathbf{k}).$
(18)

Since $\mathbf{F}$ is normal and points inward from the confining boundary, a
circulating current is induced along this edge with direction given by the
orientation of the Berry curvature vector,
$\bm{\mathcal{F}}\\!\\!_{\alpha}(\mathbf{k})=\bm{\nabla}\times\bm{\mathcal{A}}_{\alpha}$,
which depends on the band $\alpha$, the direction of $\mathbf{b}$, and the
specific value of $\mathbf{k}$.

Thus, if there is an invariance under any of the symmetries which guarantee
that the NESS occupation number $n$ only depends on energy, we may estimate
the current density by [55]

$\displaystyle\mathbf{I}(\mathbf{r})=\sum_{\alpha}\int_{\mathrm{MBZ}}d^{3}\mathbf{k}\,n[\omega_{\alpha}(\mathbf{k})+V(\mathbf{r})]\mathbf{v}_{\alpha}(\mathbf{r},\mathbf{k}).$
(19)

Now, inside the confined region $V(\mathbf{r})\ll\omega_{\alpha}(\mathbf{k})$,
and

$\displaystyle\mathbf{I}(\mathbf{r})$
$\displaystyle\simeq\sum_{\alpha}\int_{\mathrm{MBZ}}d^{3}\mathbf{k}\,n[\omega_{\alpha}(\mathbf{k})]\left[\mathbf{v}_{\alpha}^{\mathrm{(bulk)}}+\mathbf{v}_{\alpha}^{\rm(edge)}\right].$
(20)

Since there is no privileged direction of
$\mathbf{v}_{\alpha}^{\mathrm{(bulk)}}$ on each energy shell (because of the
invariance under the aforementioned symmetries),
$\omega_{\alpha}(\mathbf{k})=\omega_{\alpha}(-\mathbf{k})$ and hence
$\frac{\partial\omega_{\alpha}(\mathbf{k})}{\partial\mathbf{k}}=-\frac{\partial\omega_{\alpha}(-\mathbf{k})}{\partial\mathbf{k}}$.
Therefore, the contribution of $\mathbf{v}_{\alpha}^{(\rm bulk)}$ to the
integrand is odd and integrates to zero. This cancellation of bulk currents is
referred to as the erasure effect [20, 28]. The remaining contribution to the
current close to the edge is thus given by

$\displaystyle\mathbf{I}(\mathbf{r})$
$\displaystyle\simeq\sum_{\alpha}\int_{\mathrm{M.B.Z.}}d^{3}\mathbf{k}\,n[\omega_{\alpha}(\mathbf{k})]\mathbf{v}_{\alpha}^{\rm(edge)}.$
(21)

The monotonically decreasing behaviour of $n(\omega)$ as a function of energy
[Eq. (15)] creates an imbalance of the contributions to
$\mathbf{I}(\mathbf{r})$ within the Brillouin zone in regions where
$\mathbf{v}_{\alpha}^{\rm(edge)}\propto\nabla V(\mathbf{r})\neq 0$. Thus, Eq.
(21) predicts a net circulating surface current. As a result, any magnetic
flux vector $\mathbf{b}$ non-parallel to the temperature gradient induces a
surface crosscurrent.

Fig. 2 illustrates the surface crosscurrents for a cuboid geometry in the
bosonic case. Here, and in all the following examples, we take the reservoir
temperatures to be $T_{h}=t$ and $T_{c}=0.01t$, and for bosons always set
$\mu_{\alpha}=0$. For visual clarity, we display relatively small systems but
large enough to avoid finite-size effects. Qualitatively similar results are
obtained for other temperature values and larger lattices. The left panel of
Fig. 2 demonstrates that all currents are restricted to the boundary of the
system and vanish in the bulk. The right panel of Fig. 2 zooms in on the
topmost boundary of the system, demonstrating that particles and energy flow
against the temperature gradient within a two-dimensional surface. The
occurrence of such two-dimensional crosscurrents represents the first main
result of this work.

###  III.2 Robustness under defects

Figure 3: Robustness of bosonic boundary currents to defects (gray circles),
(a–c) in the three-dimensional case with $t_{X}=t_{Y}=t_{Z}=t$ and (d–f) in
the two-dimensional case with $t_{Z}=0$ and $t_{X}=t_{Y}=t$. (a,d) Tilted
magnetic field, $\mathbf{b}=(0.1,0.1,0.1)$, with defects symmetric under the
reflection $\hat{\Sigma}_{xy}$. (b,e) Vertical magnetic field,
$\mathbf{b}=(0,0,0.1)$, with defects symmetric under the rotation
$\hat{R}_{y}(\pi)$. (c,f) Tilted magnetic field, $\mathbf{b}=(0.1,0.1,0.1)$,
with defects symmetric under the lattice inversion $\hat{I}$.

The erasure effect of the bulk currents sustains the appearance of the chiral
current in the symmetric situation, but it is spoilt in the nonsymmetric case
because the occupation number is no longer only a function of energy alone.
This situation arises when the lattice is not perfect, e.g. due to the
presence of defects. A defect is modeled here by a point-like impurity that
induces a large on-site energy shift, i.e., a term
$\Delta{\hat{a}}^{\dagger}_{\mathbf{r}}{\hat{a}}_{\mathbf{r}}$ added to the
Hamitonian, where $\Delta\gg\omega_{0},t$ and $\mathbf{r}$ is the position of
the impurity. Nevertheless, if the defect configuration complies with any of
the protecting symmetries, such that (15) remains valid, the surface current
will be stable.

For instance, in the case $t_{Z}=0$, the system reduces to a layered stack of
independent 2D Hofstadter lattices. The surface current is robust provided
that the defect configuration satisfies either of the two symmetries
$\hat{\Theta}\hat{\Sigma}_{yz}$ or $\hat{\Theta}\hat{R}_{z}(\pi)$ [20]. This
is so independently of the orientation of $\mathbf{b}$, as shown in Fig. 3(d).
However, if $\mathbf{b}$ is not orthogonal to the temperature gradient, any
nonzero $t_{Z}$ spoils the stability of currents 2D protected by
$\hat{\Theta}\hat{\Sigma}_{yz}$, as this operation ceases to be a symmetry of
the 3D Hofstadter Hamiltonian. This is illustrated in Fig. 3(a), where a
chaotic current pattern throughout the edges and bulk of the system is induced
by the defects. The 2D symmetry $\hat{\Theta}\hat{R}_{z}(\pi)$ is even more
fragile as any $\mathbf{b}$ not parallel to some lattice vector destabilizes
the boundary currents in 3D.

Remarkably, we also find the converse situation: starting from a stack of
independent 2D Hofstadter layers with unstable currents, we can make them
stable by switching on $t_{Z}$. This is shown in Figs. 3(b,e), where there are
two defects connected by the rotation $\hat{R}_{y}(\pi)$. Since
$\hat{\Theta}\hat{R}_{y}(\pi)$ is not a symmetry of the model for $t_{Z}=0$,
the defects destroy the edge currents [Fig. 3(e)]. However, for $\mathbf{b}$
aligned to the $Z$ direction, $\hat{\Theta}\hat{R}_{y}(\pi)$ becomes a
symmetry for $t_{Z}\neq 0$, and this produces a genuinely 3D stabilization
effect [Fig. 3(e)].

###  III.3 Robustness under defects and arbitrary magnetic tilting

Figure 4: Bosonic current distribution in the presence of bulk defects (gray
circles), with magnetic field $\mathbf{b}=(0,0.1,0)$.

(a)

(b)

(c)

Figure 5: Fermionic currents under a vertical magnetic field,
$\mathbf{b}=(0,0,0.4)$. (a) Edge currents emerge due to the erasure effect far
from half filling, with chemical potential $\mu=\omega_{0}-4t$. (b) The edge
currents are strongly disrupted near half-filling, $\mu=\omega_{0}$. (c) Same
as (b) but with weaker tunnelling in the $Z$ direction, $t_{Z}=0.1t$. The
currents are now concentrated near the boundary as the system becomes
quasi-2D.

In the previous examples, the protecting symmetries reduce in fact to 2D
symmetries in planes, which may lead to stable situations depending on the
orientation of $\mathbf{b}$. Therefore, those 2D symmetries are unstable under
tilting of $\mathbf{b}$ from directions orthogonal or parallel to those
planes. In order to obtain a surface current that is unconditionally symmetry-
protected in our 3D lattice, we must employ a genuine 3D symmetry. As
discussed in Sec. III.1, such a symmetry is $\Theta\hat{I}$. Perturbating
defects which comply with this symmetry do not spoil the surface current for
any orientation of $\mathbf{b}$ or values of $t_{X,Y,Z}\neq 0$. This is
illustrated in Fig. 3(c), where defects placed on opposite corners of the
cuboid do not destroy the boundary currents despite the genuinely 3D field
configuration, $\mathbf{b}=(0.1,0.1,0.1)$. Instead, the edge currents simply
detour around the impurity sites. When the tunnelling in the $Z$ direction is
switched off, however, the defects strongly disrupt the current pattern on the
planes to which they are confined [Fig. 3(f)].

We note that the point-like character of the defects employed so far does not
play any relevant role in the robustness of the currents. Surface and
volumentric defects can be introduced leading to stable currents provided that
the protection symmetry is satisfied; see Fig. 4, for example. Here, as in
Fig. 3, we observe that the introduction of defects in a symmetric
configuration leads to the emergence of counter-propagating currents that
“shield” the impurities. This can be understood as a manifestation of the
erasure effect [20, 28], since the defects effectively generate a new boundary
within the system.

###  III.4 Fermionic lattices

In the fermionic case, the situation is the same as for bosons at relatively
high temperatures $T_{c,h}\gtrsim t_{X,Y,Z}$, where particle exchange
statistics plays little role. In contrast, fermions and bosons behave very
differently when one of the reservoirs is at very low temperature, as in our
examples. In this case, the existence of surface currents depends strongly on
the value of the chemical potential, which, via the non-equilibrium
distribution function (15), selects the portions of the single-particle
spectrum that contribute significantly to the current pattern. In Fig. 5(a) we
show that, for $\mu\ll\omega_{0}$, the fermionic system behaves similarly to
the bosonic one: the current is fully localised on the surface of the system.
These fermionic surface currents enjoy the same robustness against symmetric
defects as in the bosonic case. This behaviour is explained by the erasure
effect described in Sec. III.1.

Close to half-filling, however, the situation is quite different as shown in
Figs. 5(b,c). Here, the terminology half-filling refers to the case
$\mu=\omega_{0}$, which would correspond to exactly half-filled bands at zero
temperature. If the tunnelling in one direction is significantly weaker, e.g.
$t_{Z}\ll t_{X,Y}$, the system behaves like a set of weakly coupled 2D layers,
each hosting edge modes near energy $\omega\approx\omega_{0}$. These modes are
the predominant carrier of current when the system is half-filled, leading to
currents localised near the boundaries as shown in Fig. 5(b) [20]. These
currents are robust against any distribution of impurities in the bulk of the
system, irrespective of their symmetries [28]. However, tunnelling between the
layers causes hybridisation of the 2D band structure [40, 41] — apart from at
very specific values of the field and tunnelling amplitudes [43]. This
hybridisation destroys the edge modes and thus fermionic surface currents
typically do not appear for isotropic tunnelling when one of the reservoirs is
at low temperature [Fig. 5(c)].

##  IV Conclusions

In summary, we have studied the distribution of currents within a 3D
Hofstadter lattice driven far from equilibrium by weakly coupled particle
reservoirs at very different temperatures. We have found dissipatively robust
surface currents, which are stable against the introduction of defects in
configurations that respect certain nonequilibrium symmetries. In particular,
we demonstrated the existence of surface crosscurrents flowing against the
temperature gradient within a 2D manifold on the boundary of the system. These
results represent the 3D generalization of the robust boundary currents in 2D
models reported in Refs. [20, 28].

Notably, however, we find new effects in 3D that cannot be understood in terms
of the quasi-2D physics of concatenated layers. Specifically, in the bosonic
case, we have found that the presence of tunnelling in all three directions
provides additional stability to certain defect configurations that would
otherwise destroy the boundary currents [Fig. 3]. We also found strikingly
different behaviour for fermions when one reservoir is at low temperature.
Here, the existence of three-dimensional tunnelling tends to destroy the edge
modes responsible for boundary currents [Fig. 5]. Remarkably, therefore, in
three dimensions these surface effects are more typical of the bosonic system
than the fermionic one, in contrast to the usual intuition from topological
physics in 2D.

Our results could be experimentally tested in photonic [57] or cold-atom [58,
59] systems, in which synthetic realisations of the Hofstadter Hamiltonian
have already been achieved, or even in solid-state systems where the 3D
quantum Hall effect was recently observed [45]. This would further augment the
panoply of exotic boundary phenomena that can be explored in higher-
dimensional quantum systems out of equilibrium.

###### Acknowledgements.

M. T. M. is supported by a Royal Society-Science Foundation Ireland University
Research Fellowship, and acknowledges funding from the ERC Starting Grant
ODYSSEY (Grant Agreement No. 758403) and the Engineering and Physical Sciences
Research Council-Science Foundation Ireland Joint Funding of Research project
QuamNESS. A. R. and M. A. M.-D. acknowledge financial support from the Spanish
MINECO grants MINECO/FEDER Projects FIS2017-91460-EXP, PGC2018-099169-B-I00
FIS-2018 and from CAM/FEDER Project No. S2018/TCS-4342 (QUITEMAD-CM). The
research of A. R. and M. A. M.-D. has been partially supported by the U.S.
Army Research Office through Grant No. W911NF-14-1-0103. Calculations were
performed on the GICC cluster at UCM and the Lonsdale cluster maintained by
the Trinity Centre for High Performance Computing, which was funded through
grants from Science Foundation Ireland.

## References

  * Auerbach [1994] A. Auerbach,  _Interacting Electrons and Quantum Magnetism_, Graduate texts in contemporary physics (Springer-Verlag, 1994).
  * Gonzalez  _et al._ [2008] J. Gonzalez, M. Martin-Delgado, G. Sierra, and A. Vozmediano,  _Quantum Electron Liquids and High-Tc Superconductivity_, Lecture Notes in Physics Monographs (Springer Berlin Heidelberg, 2008).
  * Ezawa [2008] Z. Ezawa,  _Quantum Hall Effects: Field Theoretical Approach And Related Topics (2nd Edition)_ (World Scientific Publishing Company, 2008).
  * Beenakker [1997] C. W. J. Beenakker, Rev. Mod. Phys. 69, 731 (1997).
  * Rammer [2018] J. Rammer,  _Quantum Transport Theory_ (CRC Press, 2018).
  * Nazarov and Blanter [2009] Y. Nazarov and Y. Blanter,  _Quantum Transport: Introduction to Nanoscience_ (Cambridge University Press, 2009).
  * Kitaev [2003] A. Kitaev, Annals of Physics 303, 2 (2003).
  * Bombin and Martin-Delgado [2006] H. Bombin and M. A. Martin-Delgado, Phys. Rev. Lett. 97, 180501 (2006).
  * Dennis  _et al._ [2002] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, Journal of Mathematical Physics 43, 4452 (2002).
  * Bombin and Martin-Delgado [2007a] H. Bombin and M. A. Martin-Delgado, Phys. Rev. Lett. 98, 160502 (2007a).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Ando [2013] Y. Ando, Journal of the Physical Society of Japan 82, 102001 (2013).
  * Bombin and Martin-Delgado [2007b] H. Bombin and M. A. Martin-Delgado, Phys. Rev. B 75, 075103 (2007b).
  * Qi  _et al._ [2008] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 78, 195424 (2008).
  * Diehl  _et al._ [2011] S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, Nature Phys. 7, 971 (2011).
  * Budich  _et al._ [2015] J. C. Budich, P. Zoller, and S. Diehl, Phys. Rev. A 91, 042117 (2015).
  * Iemini  _et al._ [2016] F. Iemini, D. Rossini, R. Fazio, S. Diehl, and L. Mazza, Phys. Rev. B 93, 115113 (2016).
  * Linzner  _et al._ [2016] D. Linzner, L. Wawer, F. Grusdt, and M. Fleischhauer, Phys. Rev. B 94, 201105 (2016).
  * Rivas and Martin-Delgado [2017] Á. Rivas and M. A. Martin-Delgado, Sci. Rep. 7, 6350 (2017).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Phys. Rev. X 9, 041015 (2019).
  * Song  _et al._ [2019] F. Song, S. Yao, and Z. Wang, Phys. Rev. Lett. 123, 170401 (2019).
  * Shavit and Goldstein [2020] G. Shavit and M. Goldstein, Phys. Rev. B 101, 125412 (2020).
  * Gau  _et al._ [2020] M. Gau, R. Egger, A. Zazunov, and Y. Gefen, Phys. Rev. Lett. 125, 147701 (2020).
  * Lieu  _et al._ [2020] S. Lieu, M. McGinley, and N. R. Cooper, Phys. Rev. Lett. 124, 040401 (2020).
  * McGinley and Cooper [2020] M. McGinley and N. R. Cooper, Nature Phys. 16, 1181 (2020).
  * Flynn  _et al._ [2021] V. P. Flynn, E. Cobanera, and L. Viola, Phys. Rev. Lett. 127, 245701 (2021).
  * Mitchison  _et al._ [2022] M. T. Mitchison, A. Rivas, and M. A. Martin-Delgado, Phys. Rev. Lett. 128, 120403 (2022).
  * Guo and Poletti [2016] C. Guo and D. Poletti, Phys. Rev. A 94, 033610 (2016).
  * Guo and Poletti [2017] C. Guo and D. Poletti, Phys. Rev. B 96, 165409 (2017).
  * Xing  _et al._ [2020] B. Xing, X. Xu, V. Balachandran, and D. Poletti, Phys. Rev. B 102, 245433 (2020).
  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007).
  * Fu and Kane [2008] L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Hasan and Moore [2011] M. Z. Hasan and J. E. Moore, Annual Review of Condensed Matter Physics 2, 55 (2011).
  * Moore and Balents [2007] J. E. Moore and L. Balents, Phys. Rev. B 75, 121306 (2007).
  * Roy [2009] R. Roy, Phys. Rev. B 79, 195322 (2009).
  * Hsieh  _et al._ [2008] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nature 452, 970 (2008).
  * Ortiz and Martin-Delgado [2016] L. Ortiz and M. Martin-Delgado, Annals of Physics 375, 193 (2016).
  * Montambaux and Kohmoto [1990] G. Montambaux and M. Kohmoto, Phys. Rev. B 41, 11417 (1990).
  * Hasegawa [1990] Y. Hasegawa, Journal of the Physical Society of Japan 59, 4384 (1990).
  * Kunszt and Zee [1991] Z. Kunszt and A. Zee, Phys. Rev. B 44, 6842 (1991).
  * Kohmoto  _et al._ [1992] M. Kohmoto, B. I. Halperin, and Y.-S. Wu, Phys. Rev. B 45, 13488 (1992).
  * Koshino  _et al._ [2001] M. Koshino, H. Aoki, K. Kuroki, S. Kagoshima, and T. Osada, Phys. Rev. Lett. 86, 1062 (2001).
  * Bernevig  _et al._ [2007] B. A. Bernevig, T. L. Hughes, S. Raghu, and D. P. Arovas, Phys. Rev. Lett. 99, 146804 (2007).
  * Tang  _et al._ [2019] F. Tang, Y. Ren, P. Wang, R. Zhong, J. Schneeloch, S. A. Yang, K. Yang, P. A. Lee, G. Gu, Z. Qiao, and L. Zhang, Nature 569, 537 (2019).
  * Moore  _et al._ [2008] J. E. Moore, Y. Ran, and X.-G. Wen, Phys. Rev. Lett. 101, 186805 (2008).
  * Deng  _et al._ [2013] D.-L. Deng, S.-T. Wang, C. Shen, and L.-M. Duan, Phys. Rev. B 88, 201105 (2013).
  * Kennedy [2016] R. Kennedy, Phys. Rev. B 94, 035137 (2016).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Oxford, 2002).
  * Rivas and Huelga [2012] A. Rivas and S. F. Huelga,  _Open quantum systems. An introduction_ (Springer, Heidelberg, 2012).
  * Gardiner and Zoller [2004] C. Gardiner and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer, Berlin, 2004).
  * Note [1] Here, the continuity equation is derived by considering the Heisenberg equations of motion generated by the microscopic Hamiltonian given in Eq. (6).
  * Buča and Prosen [2012] B. Buča and T. Prosen, New J. Phys. 14, 073007 (2012).
  * Chang and Niu [1996] M.-C. Chang and Q. Niu, Phys. Rev. B 53, 7010 (1996).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Price and Cooper [2012] H. M. Price and N. R. Cooper, Phys. Rev. A 85, 033620 (2012).
  * Hafezi  _et al._ [2013] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor, Nature Photonics 7, 1001 (2013).
  * Aidelsburger  _et al._ [2013] M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Phys. Rev. Lett. 111, 185301 (2013).
  * Miyake  _et al._ [2013] H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, Phys. Rev. Lett. 111, 185302 (2013).