\subsection{#1}}
\newcommand{\myparagraph}[1]{\paragraph{#1}}
\newcommand{\newmyparagraph}[1]{\paragraph{#1}}
\newcommand{\vspace{-5pt}\item}{\vspace{-5pt}\item}
\newcommand{\mybf}[1]{\vspace{4pt}\noindent\textbf{#1}}
\newcommand{\mycell}[1]{\begin{tabular}{@{}l@{}l}#1\end{tabular}}
\newcommand{\mycellc}[1]{\begin{tabular}{@{}c@{}l}#1\end{tabular}}
\newcommand\crule[3][black]{\textcolor{#1}{\rule{#2}{#3}}}

\newtheorem{example}{Example}
\newtheorem{theorem}{Theorem}

\pdfinfo{
/TemplateVersion (IJCAI.2021.0)
}

\title{Temporal and Object Quantification Networks}

\author{
Jiayuan Mao$^{1*}$\and
Zhezheng Luo$^{1*}$\and
Chuang Gan$^2$\and
Joshua B. Tenenbaum$^1$\and
Jiajun Wu$^3$\and\\
Leslie Pack Kaelbling$^1$\And
Tomer D. Ullman$^4$\\
\affiliations
$^1$Massachusetts Institute of Technology\quad\\
$^2$MIT-IBM Watson AI Lab\quad\\
$^3$Stanford University\quad\\
$^4$Harvard University
}

\begin{document}

\maketitle

\footnotetext{\hspace{-7pt}*: equal contribution; order determined by a coin toss.\par Email: {\{jiayuanm,ezzluo\}@mit.edu}\par Project page: \url{http://toqnet.csail.mit.edu}}

\begin{abstract}

We present {\it Temporal and Object Quantification Networks\xspace} (TOQ-Nets\xspace), a new class of neuro-symbolic networks with a structural bias that enables them to learn to recognize complex relational-temporal events. This is done by including reasoning layers that implement finite-domain quantification over objects and time.   The structure allows them to generalize directly to input instances with varying numbers of objects in temporal sequences of varying lengths.
We evaluate TOQ-Nets\xspace on input domains that require recognizing event-types in terms of complex temporal relational patterns.
We demonstrate that TOQ-Nets\xspace can generalize from small amounts of data to scenarios containing more objects than were present during training and to temporal warpings of input sequences.
\end{abstract}
 %
\section{Introduction}

\begin{figure*}[t]
    \includegraphics[width=0.9\textwidth]{fig/ltl_teaser_new2.pdf}
    \vspace{-10pt}
    \caption{(a) An input sequence composed of relational states: each column represents the state of an entity that changes over time. A logic formula describes a complex concept or feature that is true of this temporal sequence using object and temporal quantification. The sequence is segmented into three stages: throughout the first stage, \crule[red]{1em}{1em} holds for at least one entity, until the second stage, in which each entity is always either \crule[blue]{1em}{1em} or \crule[green]{1em}{1em}, until the third stage, in which  \crule[yellow]{1em}{1em} eventually becomes true for at least one of the entities. (b) Such events can be described using first-order linear temporal logic expressions.}
    \label{fig:teaser}
    \vspace{-5pt}
\end{figure*}
Every day, people interpret events and actions in terms of concepts, defined over temporally evolving relations among agents and objects~\cite{zacks2007event,stranger1996perception}. For example, in a soccer game, people can easily recognize when one player has control of the ball, when a player {\it passes} the ball to another player, or when a player is {\it offsides}.  Although it requires reasoning about intricate relationships among sets of objects over time, this cognitive act is effortless, intuitive, and fast. It also generalizes directly over different numbers and arrangements of players, and detailed timings and trajectories.
In contrast, most computational representations of sequential concepts are based on fixed windows of space and time, and lack the ability to perform relational generalization.

In this paper, we develop generalizable representations for learning complex activities in time sequences from realistic data.  As illustrated in \fig{fig:teaser}, we can describe complex events with a first-order linear temporal logic~\citeb{FO-LTL;}{pnueli1977temporal} formula,
which allows us to flexibly decompose an input sequence into stages that satisfy different criteria over time.  Object quantifiers ($\forall$ and $\exists$) are used to specify conditions on sets of objects that define each stage.
Such representations immediately support generalization to situations with a varying number of objects, and sequences with different time warpings.

More concretely, the variety of complicated spatio-temporal trajectories that {\it high pass} can refer to in a soccer game can be described in these terms:
in a high pass from player $A$ to teammate $B$, $A$ is {\it close} to the ball ($\mathit{distance}(A, \mathit{ball}) < \theta_1$) and moving ($\mathit{velocity}(A) > \theta_2$) {\it until} the ball moves over the ground ($\mathit{z_\textit{position}}(\mathit{ball}) > \theta_3$), which is in turn  {\it until} teammate $B$ gets control of the ball ($\mathit{teammate}(A, B) \wedge \mathit{distance}(B, \mathit{ball}) < \theta_1$).
Beyond modeling human actions in physical environments, these structures can be applied to events in any time sequence of relational states, e.g\onedot} \def\Eg{E.g\onedot, characterizing particular offensive or defensive maneuvers in board games such as chess or in actual conflicts, or detecting a process of money-laundering amidst financial transaction records.

In this paper, we propose a neuro-symbolic approach to learning to recognize temporal relational patterns, called {\it Temporal and Object Quantification Networks\xspace} (TOQ-Nets\xspace), in which we design structured neural networks with an explicit bias that represents finite-domain quantification over both entities and time.
A TOQ-Net\xspace is a multi-layer neural network whose inputs are the properties of agents and objects and their relationships, which may change over time. Each layer in the TOQ-Net\xspace performs either {\it object} or {\it temporal} quantification.

The key idea of TOQ-Nets\xspace is to use {\it tensors} to represent relations between objects, and to use tensor pooling operations over different dimensions to realize temporal and object quantifiers (\textbf{\textit{G}} and $\forall$).  For example, the colorful matrix in \fig{fig:teaser}(a) can be understood as representing a unary color property of a set of 8 entities over a sequence of 8 time steps. Representing a sequence of relations among objects over time would require a 3-dimensional tensor.
Crucially, the design of TOQ-Nets\xspace allows {\em the same network weights} to be applied to domains with different numbers of objects and time sequences of different lengths.
By stacking object and temporal quantification operations, TOQ-Nets\xspace can easily learn to represent higher-level sequential concepts based on the relations between entities over time, starting from low-level sensory input and supervised with only high-level class labels.

There are traditional symbolic learning or logic synthesis methods that construct first-order or linear temporal logic expressions from accurate symbolic data~\cite{neider2018learning,camacho2018finite,chou2020explaining}.  TOQ-Nets\xspace take a different approach and can learn from noisy data by backpropagating gradients, which allows them to start with a general perceptual processing layer that is directly fed into logical layers for further processing.

We evaluate TOQ-Nets\xspace on two perceptually and conceptually different benchmarks: trajectory-based sport event detection and human activity recognition, demonstrating several important contributions.
First, TOQ-Nets\xspace outperform both convolutional and recurrent baselines for modeling temporal-relational concepts across benchmarks. Second, by exploiting  temporal-relational features learned through supervised learning, TOQ-Nets\xspace achieve strong few-shot generalization to novel actions. Finally, TOQ-Nets\xspace exhibit strong generalization to scenarios with more entities than were present during training and are robust w.r.t\onedot time warped input trajectories. %
These results illustrate the power of combining symbolic representational structures with learned continuous-parameter representations to achieve robust, generalizable interpretation of complex relational-temporal events.

 %
\section{TOQ-Nets\xspace}

\begin{figure*}[!tp]
    \includegraphics[width=0.97\textwidth]{fig/model-overview.pdf}
    \vspace{-5pt}
    \caption{A TOQ-Net\xspace contains three modules: (i) an input feature extractor, (ii) relational reasoning layers, and (iii) temporal reasoning layers.
    %
    To illustrate the model's representational power, we show that logical forms of increasing complexity can be realized by stacking multiple layers.}
    \label{fig:model-overview}
\end{figure*}

 %

The input to a TOQ-Net\xspace is a tensor representation of the properties of all entities at each moment in time. For example, in a soccer game, the input encodes the position of each player and the ball, as well as their team membership at each step of an extended time interval. The output is a label of the category of the sequence, such as the type of soccer play it contains.

The first layer of a TOQ-Net\xspace (\fig{fig:model-overview} (i)) extracts temporal features for each entity with an {\it input feature extractor} that focuses on entity features within a fixed and local time window. These features may be computed, e.g., by a convolutional neural network
or a bank of parametric feature templates.
The output of this step is a collection of nullary, unary, and binary relational features over time for all entities. Throughout the paper we will assume that all output tensors of this layer are binary-valued, but it can be extended directly to real-valued functions.
This input feature extractor is task-specific and is not the focus of this paper.

Second, these temporal-relational features go through several {\em relational reasoning} layers ({\sc rrl}s), detailed in \sect{sect:relational}, each of which performs linear transformations, sigmoid activation, and {\it object quantification} operations. The linear and sigmoid functions allow the network to realize learned Boolean logical functions, and the object quantification operators can realize quantifiers.  Additional {\sc rrl}s enable deeper nesting of quantified expressions, as illustrated in \fig{fig:model-overview}. All operations in these layers are performed for all time steps in parallel.

Next, the {\sc rrl}s perform a final quantification, computing for each time step a set of a nullary features that are passed to the \textit{temporal reasoning layers} ({\sc trl}s), as detailed in \sect{sect:temporal}. Each {\sc trl} performs linear transformations, sigmoid activation, and {\it temporal quantification}, allowing the model to realize a subset of linear temporal logic~\cite{pnueli1977temporal}.
As with {\sc rrl}s, adding more {\sc trl}s enables the network to realize logical forms with more deeply nested temporal quantifiers.

In the last layer, all object and time information is projected into a set of features of the initial time step,  which summarize the temporal-relational properties of the entire trajectory (e.g., ``the kicker eventually scores''), and fed into a final softmax unit to obtain classification probabilities for the sequence.

It is important to understand the representational power of this model. The {\it input transformation layer} learns basic predicates and relations that will be useful for defining more complex concepts, but no specific predicates or relations are built into the network in advance.
The relational reasoning layers build quantified expressions over these basic properties and relations, and might construct expressions that could be interpretable as ``the player is close to the ball.'' Finally, the temporal reasoning layer applies temporal operations to these complex expressions, such as ``the player is close to the ball until the ball moves with high speed.'' Critically, {\it none} of the symbolic properties or predicates are hand defined---they are all constructed by the initial layer in order to enable the network to express the concept it is being trained on.

 TOQ-Nets\xspace are not fully first-order: all quantifiers operate only over the finite domain of the input instance, and can be seen as ``short-hand'' for finite disjunctions or conjunctions over objects or time points.  In addition, the depth of the logical forms it can learn is determined by the fixed depth of the network.
 However, our goal is not to fully replicate temporal logic, but to bring ideas of object and temporal quantification into neural networks, and to use them as structural inductive biases to build models that generalize better from small amounts of data to situations with varying numbers of objects and time courses.

\mysubsection{Temporal-Relational Feature Representation}
TOQ-Nets\xspace use tensors as internal representations between layers;  they represent, all at once, the values of all predicates and relations grounded on all objects at all time points.  The operations in a TOQ-Net\xspace are vectorized, operating in parallel on all objects and times, sometimes expanding the dimensionality via outer products, and then re-projecting into smaller dimensions via max-pooling.
This processing style is analogous to representing an entire graph using an adjacency matrix and using matrix operations to compute properties of the nodes or of the entire graph.
In TOQ-Nets\xspace, the input to the network, as well as the feature output of intermediate layers, is represented as a tuple of three tensors.

Specifically, we use a vector of dimension $D_0$ to represent
aspects of the state that are global and do not depend on any specific object at each time $t$.  We use a matrix of shape $N \times D_1$ to represent the unary properties of each entity at time $t$, where $N$ is the number of entities and $D_1$ is the hidden dimension size. Similarly, we use a tensor of shape $N \times N \times D_2$ to represent the relations between each pair of entities at time step $t$.
As a concrete example, illustrated in \fig{fig:model-overview}, the number of entities $N$ is the total number of players plus one (the ball). For each entity $x$ and each time step, the inputs are their 3D position, type (ball or player) and team membership. The TOQ-Net\xspace outputs the action performed by the target player. Since there are only entity features, the input trajectory is encoded with a ``unary'' tensor of shape $T \times N \times D_1$, where $T$ is the length of the trajectory. That is, there are no nullary or binary inputs in this case.

\begin{figure*}[!t]
    \centering
    \includegraphics[width=\textwidth]{fig/model.pdf}
    \caption{Illustration of (i) relational reasoning layers and (ii) temporal reasoning layers. We provide two illustrative running traces. (i) The first relational reasoning layer takes unary predicates $q_1$ and $q_2$ as input and its output $Q_1$ is able to represent $q_1 \wedge q_2$. The $\max(Q_1, dim=0)$ in layer 2 can represent $\exists x.~q_1(x, t) \wedge q_2(x, t)$.
    (ii) Assume $P_K$ encodes the occurance of events $e_1$ and $e_2$ at each time step. The first temporal reasoning layer can realize {\it always} $e_2$ with a temporal pooling from time step $3$ to time step $T$. In the second temporal reasoning layer, the temporal pooling summarizes that $e_1$ holds true from time step 1 to 2. Thus, the NN should be able to realize $e_1$ {\it until} ({\it always} $e_2$).
    }
    \vspace{-5pt}
    \label{fig:model}
\end{figure*}
\mysubsection{Relational Reasoning Layers}
\label{sect:relational}
Our Relational reasoning layers ({\sc rrl}s) follow prior work on Neural Logic Machines \cite{dong2018neural}, illustrated in \fig{fig:model}~(i).
Consider a specific time step $t$. At each layer $l$, the input to a neural logic layer is a 3-tuple $(P_{l-1}, Q_{l-1}, R_{l-1})$, which corresponds to nullary, unary, and binary features respectively. Their shapes are $D_0$, $N \times D_1$, and $N \times N \times D_2$. The output is another 3-tuple $(P_{l}, Q_{l}, R_{l})$, given by
\vspace{-5pt}
\begin{adjustwidth}{0em}{0em}
\small
\begin{eqnarray*}
    P_{l} & = & \mathrm{NN}_P \left( \mathrm{Concat}\left[ P_{l-1} ;\max(Q_{l-1}, \mathrm{dim}=0) \right] \right),\\
    Q_{l} & = & \mathrm{NN}_Q \left( \mathrm{Concat}\left[ Q_{l-1} ; \max(R_{l-1}, \mathrm{dim}=0) ;
    \right.\right.\\
    & & \left.\left.\max(R_{l-1}, \mathrm{dim}=1) ; \mathrm{expand}(P_{l-1}, \mathrm{dim}=1) \right] \right),\\
    R_{l} & = & \mathrm{NN}_R \left( \mathrm{Concat}\left[ R_{l-1} ; \mathrm{expand}(Q_{l-1}, \mathrm{dim}=0) ;
    \right.\right.\\
    & & \left.\left.\mathrm{expand}(Q_{l-1}, \mathrm{dim}=1) \right] \right).
\end{eqnarray*}
\end{adjustwidth}
where $\mathrm{NN}_{*}$ are single fully-connected layers with sigmoid activations. For unary and binary features, $\mathrm{NN}_Q$ and $\mathrm{NN}_R$ are applied along the feature dimension. That is, we apply the same linear transformation to the unary features of all entities. A different linear transformation is applied to the binary features of all entity pairs.
$\mathrm{Concat}[\cdot~;~\cdot]$ is the concatenation operation, applied to the last dimension of the tensors (the feature dimension).
$\max$, also called the ``reduced'' max operation, takes the maximum value along the given axis of a tensor. The $\mathrm{expand}$ operation, also called ``broadcast,'' will duplicate the input tensor $N$ times and stack them together along the given axis. {\sc rrl}s are applied identically to the input features at every time step $t$. That is, we use the same neural network weights in a {\sc rrl} for all time steps in parallel.

{\sc rrl}s are motivated by relational logic rules in a finite and fully grounded universe.
The $\max$ reduction operations implement a differentiable version of an existential quantifier over the finite universe of individuals, given that the truth values of the propositions are represented as values in $(0.0, 1.0)$. Because preceding and subsequent {\sc rrl}s
can negate propositions as needed, we omit explicit implementation of finite-domain universal quantification, although it could be added by including analogous $\min$ reductions. Thus, as illustrated in \fig{fig:model} (i), given input features $\mathit{q_1}(x, t)$ and $\mathit{q_2}(x, t)$, we can realize the formula $\exists x.~\mathit{q_1}(x, t) \wedge \mathit{q_2}(x, t)$ by stacking two such layers.

Throughout the paper we have been using only nullary, unary, and binary features, but the proposed framework itself can be easily extended to higher-order relational features. From a graph network~\cite{bruna2013spectral,kipf2016semi,battaglia2018relational} point of view, one can treat these features as the node and edge embeddings of a fully-connected graph and the relational reasoning layers as specialized graph neural network layers for realizing object quantifiers.

\mysubsection{Temporal Reasoning Layers}
\label{sect:temporal}
Temporal reasoning layers ({\sc trl}s) perform {\it quantification} operations similar to relational reasoning layers, but along the {\it temporal} rather than the object dimension. The first {\sc trl} takes as input the summarized event representation produced by the $K$-th relational reasoning layer, $P_{K}^{(t)}$ for at all time steps $t$, as a matrix of shape $T \times D$.
Each {\sc trl} is computed as
\vspace{-1em}
\begin{adjustwidth}{0em}{0em}
\begin{equation}
    P_{K+l}^{(t)}  = \max_{t' > t} \mathrm{NN}_l \left( \mathrm{Concat}\left[  P_{K+l-1}^{(t')}; \max_{t \le t'' < t'} P_{K+l-1}^{(t'')} \right] \right).
\label{eq:temporallayer}
\end{equation}
\end{adjustwidth}
We will walk through the subexpressions of this formula.
\begin{enumerate}
    \item $P_{K+l-1}$ is the output tensor of the previous temporal reasoning layer, of shape $T \times C$, where $T$ is the number of time steps, and $C$ is the number of feature channels. Each entry in this tensor $P_{K+l-1}^{(t)}\left[i\right]$ can be interpreted as: event $i$ happens at time $t$.
    \item $\left(\max_{t \le t'' < t'} P_{K+l-1}^{(t'')}\right)$, abbreviated as $Q^{(t, t')}_{k+l-1}$ in the following text, is a vector of shape $C$. Its entry $Q^{(t, t')}_{k+l-1}[i]$, where $t \le t'$, represents the concept that event $i$ happens some time between $t$ and $t'$. {\bf Importantly, } together with the preceding and subsequent neural network operations, which can realize negation operations, it also allows us to describe the event that $i$ holds true for all time steps between $t$ and $t'$.
    \item $\mathit{NN}_l$ is a fully connected neural network with sigmoidal activation, which gets uniformly applied to all time steps $t$ and a future time step $t' > t$. Its input is composed of two parts: the events that happen at $t'$, i.e\onedot} \def\Ie{I.e\onedot, $P_{K+l-1}^{(t)}$, and the events that happen between $t$ and $t'$, summarized with temporal quantification operations, i.e\onedot} \def\Ie{I.e\onedot, $\left(\max_{t \le t'' < t'} P_{K+l-1}^{(t'')}\right)$.
    \item The outer-most max pooling operation $\max_{t'>t}$ enumerates over all $t' > t$, and test whether the condition specified by $\mathit{NN}_l$ holds for at least one such $t'$.
\end{enumerate}

A special case is the first temporal reasoning layer. It takes $P_{K}$ as its input, which is the output of last relational reasoning layer. Thus, the first temporal reasoning layer implements:
\vspace{-1em}
\begin{adjustwidth}{0em}{0em}
\begin{eqnarray*}
    P_{K+1}^{(t)}  = \mathrm{NN}_1 \left(\max_{t \le t'' < T} P_{K}^{(t'')} \right),
\end{eqnarray*}
\end{adjustwidth}
where $T$ is the sequence length. Note that there is no enumeration for a future time step $t' > t$ involved. In addition, for all temporal layers, we add residual connections by concatenating their inputs with the outputs.

Next, let's consider a running example illustrating how a TOQ-Net\xspace can recognize the event that: event $p$ holds true until event $q$ becomes always true. In LTL, this can be written as $p \textit{\textbf{U}} (\textit{\textbf{G}} q)$. Using the plain first-order logic (FOL) language, we can describe it as: $\exists~t. \left[ \forall t'.~(0 \le t' < t) \implies p(t') \right) \wedge \left( \forall t'.~(t' \ge t) \implies q(t') \right]$.

For simplicity, we consider a tensor representation for two events $p(t)$ and $q(t)$, where $p(t) = 1$ if it happens at time step $t$ and $p(t) = 0$ otherwise. Given the input sequence of length 4 in \tbl{tab:temporalexample} ($p(t)$ and $q(t)$), the first layer is capable of computing the following four properties for each time step $t$: $\textit{\textbf{G}}p$ ($\textit{\textbf{always}}~p$), which is true if $p$ holds true for all future time steps starting from $t$, $\textit{\textbf{F}}p$ ($\textit{\textbf{eventually}}~p$), which is true if $p$ is true for at least one time step starting from $t$, and similarly, $\textit{\textbf{G}}q$ ($\textit{\textbf{always}}~q$) and $\textit{\textbf{F}}q$ ($\textit{\textbf{eventually}}~q$). Overall, together with residual connections, the first layer can recognize six useful events: $p(t)$, $q(t)$ (from residual connection), $\textit{\textbf{G}}p$, $\textit{\textbf{F}}p$, $\textit{\textbf{G}}q$, and $\textit{\textbf{F}}q$ (by temporal quantification, i.e\onedot} \def\Ie{I.e\onedot pooling operations along the temporal dimension).

\begin{table}[!tp]
    \centering
    \setlength{\tabcolsep}{7pt}
    \small
    \begin{tabular}{l|cc|cccc|c}
    \toprule
        & \multicolumn{2}{c}{Input} &  \multicolumn{4}{|c|}{1-st Layer} &  \multicolumn{1}{c}{2-nd Layer}\\
        $t$ & $p(t)$ & $q(t)$ & \textit{\textbf{G}}$p$ & \textit{\textbf{F}}$p$ & \textit{\textbf{G}}$q$ & \textit{\textbf{F}}$q$ & $p$~\textit{\textbf{U}}~(\textit{\textbf{G}}$q$) \\ \midrule
        1 & T & T & F & T & F & T & T \\
        2 & T & F & F & T & F & T & T \\
        3 & F & T & F & F & T & T & F \\
        4 & F & T & F & F & T & T & F \\
    \bottomrule
    \end{tabular}
    \caption{A running example of different temporal quantification formulas that TOQ-Nets\xspace can realize. For clarity, we use T and F for True/False. In the actual computation, they are ``soft'' Boolean values ranges in $[0, 1]$. \textbf{G} means \textit{\textbf{always}}; \textbf{F} means \textit{\textbf{Eventually}}; \textbf{U} means \textit{\textbf{until}}.}
    \label{tab:temporalexample}
\end{table}

\begin{algorithm}[!tp]
\caption{An example temporal structure that the second temporal reasoning layer can recognize.}
\label{alg:secondlayer}
\begin{algorithmic}[1]
\Require{$p(t)$, $q(t)$, $(\textit{\textbf{G}}p)(t)$, $(\textit{\textbf{F}}p)(t)$, $(\textit{\textbf{G}}q)(t)$, and $(\textit{\textbf{F}}q)(t)$}
\Ensure{$\left(p\textit{\textbf{U}} (\textit{\textbf{G}} q))\right)(t)$, which is {\it true} if $p$ holds {\it true} from time step $t$ until $q$ becomes always {\it true}.}

\For{$t \gets 1$ to $T$}
  \For{$t' \gets t + 1$ to $T$}
    \If{$\forall t'' \in [t, t').~p(t'')$ and $(\textit{\textbf{G}}q)(t')$}
      \State $\left(p\textit{\textbf{U}} (\textit{\textbf{G}} q))\right)(t) \gets \mathit{true}$
    \EndIf
  \EndFor
\EndFor
\end{algorithmic}
\end{algorithm}

The second layer can realize the computation depicted in Algorithm~\ref{alg:secondlayer}. For every time step $t$, it enumerates all $t' > t$ and computes the output based on
\begin{enumerate}
    \item the events at $t'$ (represented as $P^{(t')}_{K+l-1}$ in \eqn{eq:temporallayer}, concretely the $(\textit{\textbf{G}}q)(t')$ in the Algorithm~\ref{alg:secondlayer} example), and
    \item the state of another event between $t$ and $t'$ (represented as $\left(\max_{t \le t'' < t'} P_{K+l-1}^{(t'')}\right)$ in \eqn{eq:temporallayer}, concretely the $\forall t''\in[t,t').~p(t'')$ in the Algorithm~\ref{alg:secondlayer} example).
\end{enumerate}

From the perspective of First-Order Linear Temporal Logic (FO-LTL), stacking multiple temporal reasoning layers enables us to realize FO-LTL formulas such as:
\[ p_1~\textit{\textbf{U}}~p_2~\textit{\textbf{U}}~p_3~\textit{\textbf{U}}~\cdots~\textit{\textbf{U}}~p_k, \]
which is interpreted as $p_1$ holds true until $p_2$ becomes true and $p_2$ holds true from that until $p_3$ becomes true $\cdots$, and
\[ \textit{\textbf{F}} p_1~\textit{\textbf{XF}}~p_2~\textit{\textbf{XF}}~p_3~\textit{\textbf{XF}}~\cdots~\textit{\textbf{XF}}~p_k, \]
which is interpreted as $p_1$ eventually becomes true and after that $p_2$ eventually becomes true and after that $\cdots$, and in addition, formulas with interleaved until and eventually quantifiers. Here, $\textit{\textbf{XF}}$ is a composition of the ne\textit{\textbf{X}}t operator and the \textit{\textbf{F}}inally operator in LTL. Meanwhile, as described so far, TOQ-Nets\xspace can only nest object quantification inside temporal quantification, so it can represent
\textbf{\textit{always}} $\exists x.~q_1(x) \wedge q_2(x)$, but not $\exists x.~\textbf{\textit{always}}\; q_1(x) \wedge q_2(x)$. This can be solved by interleaving relational and temporal reasoning layers.

It is important to notice that, by using object and temporal pooling operations together with trainable neural networks to realize logic formulas with object and temporal quantifiers, the idea itself generalizes to a broader set of FO-LTL formulas. We design TOQ-Nets\xspace to model only a subset of FO-LTL formulas, because they can be computed efficiently (with only $O(T^2)$ space) and they are expressive enough for the type of data we are trying to model.

 %
\section{Experiments}

We compare our model with other approaches to object-centric temporal event detection in this section,
and include an application of TOQ-Nets\xspace to concept learning over robot object-manipulation trajectories in the supplementary material.
The setups and metrics focus on data efficiency and generalization.

\mysubsection{Baseline Approaches}
We compare TOQ-Nets\xspace against five baselines. The first two are spatial-temporal graph convolutional neural networks \citeb{STGAN;}{yan2018spatial} and its variant STGCN-MAX, which models entity relationships with graph neural networks and models temporal structure with temporal-domain convolutions. The third is STGCN-LSTM, which uses STGCN layers for entity relations but LSTM~\cite{hochreiter1997long} for temporal structures. The last two baselines are based on space-time graphs: Space-Time Graph~\cite{wang2018videos} and Non-Local networks~\cite{wang2018non}. We provide details about our implementation and how we choose the model configurations in the supplementary material.

\mysubsection{Trajectory-Based Soccer Event Detection}

We start our evaluation with an event-detection task in soccer games. The task is to recognize the action performed by a specific player at specific time step in a soccer game trajectory.

\begin{table}[!tp]
\vspace{0pt}
    \centering\small
    \setlength{\tabcolsep}{10.5pt}
    %
    \begin{tabular}{lccc}
    \toprule
        Model & Reg. & Few-Shot  & Full  \\ \midrule
        %
        %
        %
        %
        %
        %
        STGCN & 73.2$_{\pm 1.6}$ & 26.0$_{\pm 5.7}$ & 62.8$_{\pm 0.6}$ \\
        STGCN-MAX & 73.6$_{\pm 1.5}$ & 28.6$_{\pm 5.0}$ & 63.6$_{\pm 0.7}$ \\
        STGCN-LSTM & 72.7$_{\pm 1.4}$ & 23.8$_{\pm 5.9}$ & 61.9$_{\pm 0.6}$ \\
        Space-Time & 74.8$_{\pm 1.5}$ & 31.7$_{\pm 6.1}$ & 65.2$_{\pm 0.6}$ \\
        Non-Local & 76.5$_{\pm 2.4}$ & 45.0$_{\pm 6.3}$ & 69.5$_{\pm 2.4}$ \\
        TOQ-Net\xspace (ours) & \textbf{87.7$_{\pm 1.3}$} & \textbf{52.2$_{\pm 6.3}$} & \textbf{79.8$_{\pm 0.8}$} \\
        %
        %
    \bottomrule
    \end{tabular}
    \vspace{-5pt}
        \captionof{table}{Results on the soccer event dataset. Different columns correspond to different action sets (the regular, few-shot, and full action sets). The performance is measured by per-action (macro) accuracy, averaged over nine few-shot splits. The ${\pm}$ values indicate standard errors. TOQ-Net\xspace significantly outperforms all baseline methods on the few-shot action set.}
        \label{tab:gfootball}
\end{table}

    %

\newmyparagraph{Dataset and setup.}
We collect training and evaluation datasets based on the gfootball simulator\footnote{\url{https://research-football.dev/}}, which provides a physics-based 3D football simulation. It also provides AI agents that can be used to generate random plays.
The simulator provides the 3D coordinates of the ball and the players as well as the action each player is performing at each time step. There are in total 13 actions defined in the simulator, including \textit{movement, ball\_control, trap, short\_pass, long\_pass, high\_pass, header, shot, deflect, catch, interfere, trip} and \textit{sliding}.
We exclude \textit{header} and \textit{catch} actions, as they never appear in AI games. We also exclude \textit{ball\_control} and \textit{movement}, since they just mean the agent is moving (with or without the ball). Thus, in total, we have nine action categories. We run the simulator with AI-controlled players to generate plays, and formulate the task as classifying the action (9-way classification) of a specific player at a specific time step given a temporal context (25 frames). For each action, we have generated 5,000 videos, expect for \textit{sliding}, for which we generated 4,000 videos because it is rare in the AI games. Among the generated examples, 62\% (2,462 or 3,077) are used for training, 15\% are used for validation, and 23\% are used for testing. Each trajectory is an 8-fps replay clip that contains 17 frames (about two seconds). There is a single ``target'' player in each trajectory. The action label of the trajectory is the action performed by this target player at frame \#9. We randomly split all actions into two categories: seven ``regular'' actions, for which all game plays are available, and two ``few-shot'' actions, for which only 50 clips are available during training.

\begin{table*}
    \centering\small
    \begin{tabular}{lccccccc}
    \toprule
        Model   & 3v3 & 4v4 & 6v6 & 6v6 \scriptsize{(Time Warp)} & 8v8 & 11v11  \\ \midrule
STGCN &\cellcolor{yellow!25}{\mycellc{40.7$_{\pm 1.0}$~~{\scriptsize (-40.4\%)}}}  & \cellcolor{yellow!25}{\mycellc{63.2$_{\pm 4.9}$~~{\scriptsize (-7.4\%)}}}  & {\mycellc{68.2$_{\pm 2.8}$}} & \cellcolor{blue!5}{\mycellc{52.8$_{\pm 7.0}$~~{\scriptsize (-22.6\%)}}}  & \cellcolor{yellow!25}{\mycellc{55.4$_{\pm 3.3}$~~{\scriptsize (-18.8\%)}}}  & \cellcolor{yellow!25}{\mycellc{44.4$_{\pm 2.1}$~~{\scriptsize (-34.9\%)}}} \\
STGCN-MAX & \cellcolor{yellow!25}{\mycellc{47.4$_{\pm 3.2}$~~{\scriptsize (-33.7\%)}}}  & \cellcolor{yellow!25}{\mycellc{68.8$_{\pm 2.0}$~~{\scriptsize (-3.8\%)}}}  & {\mycellc{71.5$_{\pm 1.9}$}} & \cellcolor{blue!5}{\mycellc{56.5$_{\pm 4.5}$~~{\scriptsize (-21.0\%)}}}  & \cellcolor{yellow!25}{\mycellc{59.1$_{\pm 0.7}$~~{\scriptsize (-17.3\%)}}}  & \cellcolor{yellow!25}{\mycellc{45.6$_{\pm 2.5}$~~{\scriptsize (-36.2\%)}}} \\
STGCN-LSTM & \cellcolor{yellow!25}{\mycellc{39.7$_{\pm 1.1}$~~{\scriptsize (-43.1\%)}}}  & \cellcolor{yellow!25}{\mycellc{60.4$_{\pm 0.2}$~~{\scriptsize (-13.5\%)}}}  & {\mycellc{69.8$_{\pm 0.1}$}} & \cellcolor{blue!5}{\mycellc{30.6$_{\pm 0.6}$~~{\scriptsize (-56.1\%)}}}  & \cellcolor{yellow!25}{\mycellc{55.8$_{\pm 2.0}$~~{\scriptsize (-20.0\%)}}}  & \cellcolor{yellow!25}{\mycellc{44.1$_{\pm 0.7}$~~{\scriptsize (-36.8\%)}}} \\
Space-Time & \cellcolor{yellow!25}{\mycellc{29.0$_{\pm 1.6}$~~{\scriptsize (-60.4\%)}}}  & \cellcolor{yellow!25}{\mycellc{53.5$_{\pm 3.2}$~~{\scriptsize (-27.0\%)}}}  & {\mycellc{73.3$_{\pm 0.3}$}}  & \cellcolor{blue!5}{\mycellc{70.7$_{\pm 0.3}$~~{\scriptsize (-3.5\%)}}} & \cellcolor{yellow!25}{\mycellc{33.9$_{\pm 2.8}$~~{\scriptsize (-53.7\%)}}}  & \cellcolor{yellow!25}{\mycellc{15.2$_{\pm 1.8}$~~{\scriptsize (-79.3\%)}}} \\
Non-Local & \cellcolor{yellow!25}{\mycellc{45.9$_{\pm 5.1}$~~{\scriptsize (-41.2\%)}}}  & \cellcolor{yellow!25}{\mycellc{70.7$_{\pm 5.3}$~~{\scriptsize (-9.5\%)}}}  & {\mycellc{78.1$_{\pm 5.8}$}} & \cellcolor{blue!5}{\mycellc{77.7$_{\pm 5.0}$~~{\scriptsize (-0.5\%)}}}  & \cellcolor{yellow!25}{\mycellc{58.5$_{\pm 10.8}$~~{\scriptsize (-25.1\%)}}}  & \cellcolor{yellow!25}{\mycellc{41.8$_{\pm 13.6}$~~{\scriptsize (-46.5\%)}}} \\
TOQ-Net\xspace & \cellcolor{yellow!25}{\mycellc{\bf 77.4$_{\pm 3.5}$~~{\scriptsize (-12.4\%)}}}  & \cellcolor{yellow!25}{\mycellc{\bf 88.3$_{\pm 0.7}$~~{\scriptsize (-0.0\%)}}}  & {\mycellc{\bf 88.4$_{\pm 0.6}$}} & \cellcolor{blue!5}{\mycellc{\bf 86.9$_{\pm 0.4}$~~{\scriptsize (-1.7\%)}}}  & \cellcolor{yellow!25}{\mycellc{\bf 81.3$_{\pm 1.7}$~~{\scriptsize (-8.0\%)}}}  & \cellcolor{yellow!25}{\mycellc{\bf 77.1$_{\pm 1.7}$~~{\scriptsize (-12.8\%)}}} \\
    \bottomrule
    \end{tabular}
    \vspace{-5pt}
    \caption{Results on generalization to scenarios with more agents and temporally warped trajectories on the soccer event dataset.
    All models are trained only on 6v6 games. The standard errors indicated by the $\pm$ signs are computed with three random seeds.
    %
    %
    }
    \label{tab:gfootball-generalization}
\end{table*}

    %
    %
    %
    %
    %
    %
    %
    %
    %
    %
    %
    %

\newmyparagraph{Input features.}
Each trajectory is represented with 7 time-varying unary predicates, including the 3D coordinate of each player and the ball and four extra predicates defining the type of each entity $x$: $\mathrm{IsBall}(x)$, $\mathrm{IsTargetPlayer}(x)$, $\mathrm{SameTeam}(x)$, $\mathrm{OpponentTeam}(x)$, where $\mathrm{SameTeam}(x)$ and $\mathrm{OpponentTeam}(x)$ indicates whether $x$ is of the same team as the target player. We also add a temporal indicator function which is a Gaussian function centered at frame of interest with variance $\sigma^2=25$.

\newmyparagraph{Results.}
\tbl{tab:gfootball} shows the result. Our model significantly outperforms all the baselines in all three action settings,  suggesting that our model is able to discover a set of useful features at both input and intermediate levels and use them to compose new action classifiers from only a few examples.

\newmyparagraph{Generalization to more players.} Due to their object-centric design, TOQ-Nets\xspace can generalize to soccer games with a varying number of agents. After training on 6v6 soccer games (i.e\onedot} \def\Ie{I.e\onedot, 6 players on each team), we evaluate the performance of different models on games with different numbers of players: 3v3, 4v4, 8v8, and 11v11.
For each action we have generated, on average, 1,500 examples for testing. \tbl{tab:gfootball-generalization} summarizes the result and the full results are provided in the supplementary material.
Comparing the columns highlighted in yellow, we notice a significant performance drop for all baselines while TOQ-Net\xspace performs the best.
By visualizing data and predictions, we found that misclassifications of instances of {\it shot} as {\it short pass} contribute most to the performance degradation of our model when we have more players. Specifically, the recall of {\it shot} drops from 97\% to 60\%. In soccer plays with many agents, a shot is usually unsuccessful and a player from another team steals the ball in the end. In such scenarios, TOQ-Net\xspace tends to misclassify such trajectories as a {\it short pass}. Ideally, this issue should be addressed by understanding actions based on agents' goals instead of the actual outcome~\cite{intille2001recognizing}. We leave this extension as a future direction.

\newmyparagraph{Generalization to temporally warped trajectories.} Another crucial property of TOQ-Nets\xspace is to recognize actions based on their sequential order in the input trajectory, instead of binding features to specific time steps in the trajectory. To show this, we test the performance of different models on time warped trajectories. Each test trajectory has a length of 25, and each trajectory is labeled by the action performed by the target player at any time step between the 6-th and 19-th frame. We ensure that the target player performs only one action during the entire input trajectory. Thus, the label is unambiguous.
The results are shown in \tbl{tab:gfootball-generalization}. Specifically, our test set consists of 25-frame trajectories, and the action may occur at anytime between the 6th and the 19th frame. By comparing rows with and without time warping, we notice a $60\%$ performance drop for STGCN, STGCN-MAX, and STGCN-LSTM. In contrast, TOQ-Nets\xspace still have reasonable performance. Note that Space-Time and Non-Local model have almost no performance drop against time warping because they are completely agnostic to temporal ordering.

\mysubsection{Extension to Real-World Datasets}
The proposed TOQ-Net\xspace can also be extended to other real-world datasets. %
These examples further illustrate the robustness of TOQ-Net to temporal variations in activities.

\subsubsection{Toyota Smarthome}
Toyota Smarthome~\cite{das2019toyota} is a dataset that contains videos of humans performing everyday activities such as ``walk'', ``take pills'', and ``use laptop''. It also comes with 3D-skeleton detections. There are around 16.1k videos in the dataset, and 19 activity classes. The videos’ length varies between a few seconds to 3 minutes. We subsample 30 frames for each video. We split frames into training (9.9k), validation (2.5k), and testing (3.6k). We treat human joints as entities. The input is then the position of joints, the velocity of joints, limb lengths, and joint angles. We evaluated our model and STGCN on a 19-way classification task. We also test model performance on time-warped sequences by accelerating the trajectories by two times.

Our model achieves a comparable accuracy to STGCN on the standard classification task: (42.0\% vs. 43.0\%). Importantly, on the generalization test to time-warped sequences, our model has only a 0.8\% performance drop (41.2\%), while STGCN drops 10.7\% (32.3\%). This indicates that the temporal structures learned by TOQ-Net\xspace improve model generalization to varying time courses.

\subsubsection{Volleyball Activity}
The volleyball dataset~\cite{msibrahi2016deepactivity} contains 4830 video clips collected from 55 youtube volleyball videos. They are labeled with 8 group activities (e.g. “left spike” and “right pass”). Each video contains 20 frames with the labeled group activity performed at the 10-th frame. The dataset also includes annotations for players, including the bounding box, the indicator of whether the player is involved in the group activity, and the individual action such as ``setting'', ``digging'', and ``spiking''. We use the manual annotations (processed by an MLP) as the input features. We train models to classify the video into one of the eight group activities, following the original split, i.e., 24, 15, and 16 of 55 videos are used for training, validation, and testing.

On the standard classification task, TOQ-Net achieves a comparable performance with STGCN (73.3\% vs. 73.6\%). When we perform time warping on the input sequences, STGCN’s performance drops by more than 25.0\% (39.5\% on temporally shifted trajectories and 48.6\% on 2$\times$ quick motion trajectories), while our model drops only 3\% (70.3\% on temporally shifted trajectories and 70.7\% on quick motion trajectories). This again shows the generalization ability of TOQ-Net w.r.t\onedot varying time courses, and the robustness of learned temporal structures.

\iffalse
\mysubsection{Skeleton-based Activity Recognition}
\jiayuan{Given that we may drop this dataset in the end, I will stop editing it}
The second benchmark that we use for evaluation is Toyota Smarthome, a skeleton-based action recognition dataset. It contains real-world videos recorded in an apartment where human subjects will perform different household actions.

\begin{table}[]
    \caption{Results on the Toyota Smarthome dataset for skeleton-based action recognition. Different columns correspond to different action set (the regular actions set, the few-shot action set, and the full action set). In each cell we show both of the micro average accuracy and the macro average accuracy. Our model TOQ-Net\xspace achieves a comparable results with the baselines on the regular action set and surpass them by a significant margin in the few-shot learning setting. The results with STGCN as the feature extractor shows that our model can be integrated with convolutional neural networks as local feature extractors to perform temporal-relational reasoning.}
    \vspace{5pt}
    \centering\small
    \begin{tabular}{llccc}
    \toprule
        Model & Feature & \mycellc{Reg.\\Micro/Macro} & \mycellc{Few-Shot\\Micro/Macro}  & \mycellc{Full \\Micro/Macro}  \\ \midrule
        STGCN~\cite{yan2018spatial} & * & 40.6 / 38.9 & 23.9 / 21.6 & 35.8 / 36.1 \\
        \midrule
        %
        %
        %
        TOQ-Net\xspace & Human-defined feature & ??.? / ??.? & ??.? / ??.? & ??.? / ??.? \\
        TOQ-Net\xspace & MLP feature           & ??.? / ??.? & ??.? / ??.? & ??.? / ??.? \\
        TOQ-Net\xspace & STGCN feature         & ??.? / ??.? & ??.? / ??.? & ??.? / ??.? \\
    \bottomrule
    \end{tabular}
    \label{tab:smarthome}
\end{table}
\myparagraph{Dataset and setup.}
Toyota Smarthome dataset contains 19 actions in total. In each video, there is a single person performing different household actions such as drinking water or using a laptop. We only use the 3D skeleton of the person as the model input.

Following the same evaluation metric as in the soccer play dataset, we split the action set into two subsets: regular actions and few-shot actions. For each regular action, on the average 500 video clips will be available for the model while for few-shot actions, there will be only 25 videos. We report the classifications on both the regular action set and the full action set.

\myparagraph{Baselines.}
We again compare our model with the STGCN model, which was originally designed for skeleton-based action recognition. We follow the original STGCN design where the edges in the graph are pre-defined human joint connections (e.g\onedot} \def\Eg{E.g\onedot, the limbs of the human). For TOQ-Net\xspace, we input the binary joint connectivity as a binary feature in the model input.

\myparagraph{Results.}
The results are summarized in \tbl{tab:smarthome}. TOQ-Net\xspace outperforms the baseline with a significant margin, especially on the full action classification benchmark.
\fi %
\section{Related Work}

\newmyparagraph{Action concept representations and learning.}
First-order and linear temporal logics~\citeb{LTL;}{pnueli1977temporal}
have been used for analyzing sporting events~\cite{intille1999framework,intille2001recognizing} and activities of daily living~\cite{tran2008event,brendel2011probabilistic} in logic-based reasoning frameworks. However, these frameworks require extra knowledge to annotate relationships between low-level, primitive actions and complex ones, or performing search in a large combinatorial space for candidate temporal logic rules~\cite{penning2011nsca,lamb2007connectionist}. By contrast, TOQ-Nets\xspace enable end-to-end learning of complex action descriptions from sensory input with only high-level action-class labels.

\newmyparagraph{Temporal and relational reasoning.}
This paper is also related to work on using data-driven models for modeling relational and temporal structure, such as LTL ~\cite{neider2018learning,camacho2018finite,chou2020explaining}, Logical Neural Networks~\cite{riegel2020logical}, ADL description languages~\cite{intille1999framework}, and hidden Markov models~\cite{tang2012learning}. These models need human-annotated action descriptions and symbolic state variables (e.g\onedot} \def\Eg{E.g\onedot, {\it pick up} $x$ means a state transition from not {\it holding} $x$ to {\it holding} $x$), and dedicated inference algorithms such as graph structure learning. In contrast, TOQ-Nets\xspace have an end-to-end design, and can be integrated with other neural networks.
People have also used structural representations to model object-centric temporal concepts with graph neural networks~\citeb{GNNs;}{yan2018spatial}, recurrent neural networks~\citeb{RNNs;}{msibrahi2016deepactivity}, and integrated GNN-RNN architectures \cite{deng2016sim,qi2018stag}.
TOQ-Nets\xspace use a similar relational representation, but different models for temporal structures.

 %
\section{Conclusion and Discussion}
The design of TOQ-Nets\xspace suggests multiple research directions. For example, the generalization of the acquired action concepts to novel object kinds, such as from {\it opening fridges} to {\it opening bottles}, needs further exploration. Meanwhile, TOQ-Nets\xspace are based on physical properties, e.g\onedot} \def\Eg{E.g\onedot 6D poses. Incorporating representations of mental variables such as goals, intentions, and beliefs can aid in action and event recognition~\cite{baker2017rational,zacks2001perceiving,vallacher1987people}.

In summary, we have presented TOQ-Nets\xspace, a neuro-symbolic architecture for learning to describe complex state sequences with quantification over both entities and time.  TOQ-Nets\xspace use tensors to represent the time-varying properties and relations of different entities, and use tensor pooling operations over different dimensions to realize temporal and object quantifiers.
TOQ-Nets\xspace generalize well to scenarios with varying numbers of entities and time courses.

\section*{Acknowledgements} We thank Tom\'as Lozano-P\'erez for helpful discussions and suggestions. We gratefully acknowledge support from the Center for Brains, Minds and Machines (NSF STC award CCF1231216), ONR MURI N00014-16-1-2007, NSF grant 1723381, AFOSR grant FA9550-17-1-0165, ONR grant N00014-18-1-2847, the Honda Research Institute, MIT-IBM Watson Lab, SUTD Temasek Laboratories, the Samsung Global Research Outreach (GRO) Program, and the Stanford Institute for Human-Centered AI (HAI). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of our sponsors.

\clearpage

{
\bibliographystyle{named}

\subsection{#1}}
\newcommand{\myparagraph}[1]{\paragraph{#1}}
\newcommand{\newmyparagraph}[1]{\paragraph{#1}}
\newcommand{\vspace{-5pt}\item}{\vspace{-5pt}\item}
\newcommand{\mybf}[1]{\vspace{4pt}\noindent\textbf{#1}}
\newcommand{\mycell}[1]{\begin{tabular}{@{}l@{}l}#1\end{tabular}}
\newcommand{\mycellc}[1]{\begin{tabular}{@{}c@{}l}#1\end{tabular}}
\newcommand\crule[3][black]{\textcolor{#1}{\rule{#2}{#3}}}

\title{Supplementary Material for\\Temporal and Object Quantification Networks}

\author{
Jiayuan Mao$^{1*}$\and
Zhezheng Luo$^{1*}$\and
Chuang Gan$^2$\and
Joshua B. Tenenbaum$^1$\and
Jiajun Wu$^3$\and\\
Leslie Pack Kaelbling$^1$\And
Tomer D. Ullman$^4$\\
\affiliations
$^1$Massachusetts Institute of Technology\quad\\
$^2$MIT-IBM Watson AI Lab\quad\\
$^3$Stanford University\quad\\
$^4$Harvard University
}

\begin{document}

\maketitle

We start from inrtoducing a new dataset  we collected for model evaluation, namely the RLBench dataset (\sectapp{app:rlbench}). We then provide implementation details of different models in \sectapp{app:implementation}. Finally, in \sectapp{app:additional}, we present additional ablation studies and visualizations.


\section*{\sc \bf Supplemental Material}
\subsection{Steady-state solution}
\label{app:solution}

In this section we briefly detail the solution for the non-equilibrium steady state using the quantum Langevin formalism. We consider a general quadratic Hamiltonian of the form $\hat{H}_{\rm tot} = \hat{H}_S + \hat{H}_{SB} + \hat{H}_B$, with
\begin{align}\label{H_S}
	\hat{H}_S  & = \sum_{j,k} H_{jk} \ensuremath{\hat{a}^\dagger}_j\a_k, \\
	\label{H_SB}
	\hat{H}_{SB} & = \sum_{j\in \partial S} \sum_q g_{qj} \left (\ensuremath{\hat{a}^\dagger}_j \b_{qj} + \ensuremath{\hat{b}^\dagger}_{qj} \a_j \right ), \\
	\label{H_B}
	\hat{H}_B & =  \sum_{j\in \partial S} \sum_q \Omega_{qj} \ensuremath{\hat{b}^\dagger}_{qj} \b_{qj}.
\end{align}
In Eq.~\eqref{H_S}, $\a^\dagger_j$ creates a boson or fermion on site $j$ of the system and $H_{jk}$ are the elements of a hermitian single-particle Hamiltonian matrix, $\mathbf{H}$. Similarly, the canonical operators $\ensuremath{\hat{b}^\dagger}_{qj}$ create a particle in mode $q$ of the reservoir connected to site $j$, with corresponding frequency $\Omega_{qj}$ and coupling $g_{qj}$ to a site $j$ lying on the boundary of the system, $\partial S$.

Following the standard procedure, we formally solve the equations of motion for $\b_{qj}$ in the Heisenberg picture and substitute the result into the Heisenberg equation for $\a_j$, obtaining
\begin{equation}\label{QLE_general}
	\ensuremath{\mathrm{i}} \partial_t\hat{a}_j(t) =  \sum_k H_{jk} \a_k(t) + \int_{t_0}^\infty\mathrm{d} t'\, \chi_j(t-t') \a_j(t') + \hat{\xi}_j(t).
\end{equation}
Above, we have introduced the retarded memory kernel $\chi_j(t)$ and the noise operator $\hat{\xi}_j(t)$, which for $j\in \partial S$ are given by
\begin{align}
	\label{noise_op}
	\hat{\xi}_j(t) & = \sum_q g_{qj} \ensuremath{\mathrm{e}}^{-\ensuremath{\mathrm{i}} \Omega_{qj}(t-t_0)} \b_{qj}(t_0), \\
	\label{memory_kernel}
	\chi_j(t) & = -\ensuremath{\mathrm{i}} \theta(t)\sum_q g^2_{qj} \ensuremath{\mathrm{e}}^{-\ensuremath{\mathrm{i}} \Omega_{qj}t},
\end{align}
with $\theta(t)$ the unit step function, while for $j\notin \partial S$ we define $	\hat{\xi}_j(t) = 0 = \chi_j(t) $. Note that the memory kernel is nothing but the retarded Green function of the noise operator, since $\theta(t-t')\langle [\hat{\xi}_j(t),\hat{\xi}^\dagger_k(t')]_\pm\rangle = \ensuremath{\mathrm{i}}  \delta_{jk}\chi_j(t-t')$, where the minus (plus) sign indicates the (anti-)commutator for bosons (fermions).

Since we are interested in the steady state, we take the limit $t_0\to -\infty$ in Eq.~\eqref{QLE_general} and solve in the Fourier domain. The solution can be expressed in the compact form
\begin{equation}\label{a_solution_vector}
	\tilde{\mathbf{a}}(\omega) = \mathbf{G}(\omega) \cdot \tilde{\boldsymbol{\xi}}(\omega),
\end{equation}
by defining vectors $\tilde{\mathbf{a}}= (\tilde{a}_1, \tilde{a}_2,\ldots  )^T$ and $\tilde{\bm{\xi}}(\omega) = (\tilde{\xi}_1, \tilde{\xi}_2,\ldots  )^T$ of Fourier-transformed operators and introducing the (matrix-valued) retarded Green function
\begin{equation}\label{G_ret}
	\mathbf{G}(\omega) = \left [ \omega\mathbf{1} - \mathbf{H} -\mathbf{\Sigma}(\omega) \right ]^{-1}.
\end{equation}
Here, $\mathbf{1}$ is the identity matrix and the retarded self-energy matrix $\mathbf{\Sigma}(\omega)$ is diagonal with elements $\Sigma_{jk}(\omega) = \delta_{jk} \tilde{\chi}_j(\omega)$ given by the Fourier transform of the memory kernel. The latter decomposes as
\begin{equation}\label{memory_kernel_decomp}
	\tilde{\chi}_j(\omega) = \fint \frac{\mathrm{d}\omega'}{2\pi} \frac{\gamma_j(\omega')}{\omega-\omega'} - \ensuremath{\mathrm{i}} \frac{\gamma_j(\omega)}{2},
\end{equation}
where $ \fint$ denotes a principal-value integral and we introduced the spectral density
\begin{equation}\label{spectral_density}
	\gamma_j(\omega) = 2\pi \sum_q g_{qj}^2\delta(\omega-\Omega_{qj}).
\end{equation}

The non-equilibrium steady state of the system depends on the statistical properties of the noise operator via Eq.~\eqref{a_solution_vector}. Assuming that the initial state of the reservoirs (at time $t=t_0$) is thermal and uncorrelated, we find the noise spectrum
\begin{equation}\label{noise_spectrum}
	\langle \tilde{\xi}^\dagger_j(\omega) \tilde{\xi}_k(\omega')\rangle = 2\pi\delta(\omega-\omega')\delta_{jk} \gamma_j(\omega) f_j(\omega),
\end{equation}
where $f_j(\omega) = (\ensuremath{\mathrm{e}}^{\beta_j(\omega-\mu_j)} \pm 1 )^{-1}$ is the reservoir distribution function, characterized by an inverse temperature $\beta_j$ and chemical potential $\mu_j$, and expectation values are taken with respect to the state at time $t=t_0$. It is now straightforward to write down the two-point correlation function of the system in its steady state, $C_{jk}(t) = \langle \ensuremath{\hat{a}^\dagger}_k(t) \a_j(0)\rangle$, as
\begin{equation}\label{two_point_function}
	\mathbf{C}(t) = \int\frac{\mathrm{d}\omega	}{2\pi} \ensuremath{\mathrm{e}}^{\ensuremath{\mathrm{i}}\omega t} \mathbf{G} \cdot (\mathbf{\Gamma} \circ \mathbf{F}) \cdot \mathbf{G}^\dagger.
\end{equation}
Here we have suppressed frequency arguments and defined diagonal matrices $\mathbf{\Gamma} = \ensuremath{\mathrm{i}} (\mathbf{\Sigma} - \mathbf{\Sigma}^\dagger)$ and $F_{jk} = \delta_{jk} f_j$, while $\circ$ denotes the elementwise (Hadamard) product. The steady-state correlation matrix $\mathbf{C} = \mathbf{C}(0)$ is found simply by setting $t=0$ in Eq.~\eqref{two_point_function}.

The setup considered in the main text features hot baths coupled to sites on the boundary at $x=1$ and cold baths coupled to the other boundary at $x=L_X$. Moreover, the wide-band limit approximation is made, which amounts to replacing the spectral densities $\gamma_j(\omega)$ by frequency-independent constants. This approximation is valid so long as $\gamma_j(\omega)$ does not vary significantly over the frequency range within which the integrand of Eq.~\eqref{two_point_function} differs appreciably from zero.  Assuming symmetric coupling, $\gamma_j = \gamma$, one then recovers Eq.~(4). In this case, $\mathbf{\Gamma} = \mathbf{\Gamma}_h + \mathbf{\Gamma}_c$, where the self-energy $\mathbf{\Gamma}_h$ ($\mathbf{\Gamma}_c$) of the hot (cold) bath is given by a diagonal matrix, with $[\mathbf{\Gamma}_{h,c}]_{jk} = \gamma \delta_{jk}$ whenever $j$ corresponds to a boundary site with $x=1$ ($x=L_X$) and $[\mathbf{\Gamma}_{h,c}]_{jk} = 0$  otherwise. Likewise, we have $f_j = \bar{n}_h$ ($f_j = \bar{n}_c$) for $j$ corresponding to $x=1$ ($x=L_X$).

\subsection{Definition of currents}
\label{app:currents}

The current densities are defined via the continuity equation for the particle density, $\hat{n}_{x,y} = {\bf\a}^\dagger_{x,y}
\cdot {\bf \a}_{x,y}$, as
\begin{align}\label{continuity_equations}
	\ensuremath{\mathrm{i}} [\hat{H}_X, \hat{n}_{x,y}] & = \hat{J}^X_{x-1,y} - \hat{J}^X_{x,y} \notag \\
	\ensuremath{\mathrm{i}} [\hat{H}_Y, \hat{n}_{x,y}] 	& =  \hat{J}^Y_{x,y-1} - \hat{J}^Y_{x,y}.
\end{align}
Here $\hat{J}^X_{x,y}$ denotes the current density operator describing particle flow from site $(x,y)$ to its nearest neighbour in the $x$ direction, i.e.~$(x,y) \to (x+1,y)$. Similarly, $\hat{J}^Y_{x,y}$ is the corresponding current in the $y$ direction, $(x,y) \to (x,y+1)$. Using the canonical (anti-)commutation relations, it is straightforward to derive
\begin{align}
	\label{current_op_X}
	\hat{J}^X_{x,y} & = \frac{t_X}{2\ensuremath{\mathrm{i}} } {\bf \a}^\dagger_{x+1,y}\cdot (\sigma_z+\ensuremath{\mathrm{i}}\sigma_y)\cdot {\bf \a}_{x,y} + {\rm h.c.},\\
	\label{current_op_Y}
	\hat{J}^Y_{x,y} & = \frac{t_Y}{2\ensuremath{\mathrm{i}} } {\bf \a}^\dagger_{x,y+1}\cdot (\sigma_z+\ensuremath{\mathrm{i}}\sigma_x)\cdot {\bf \a}_{x,y} + {\rm h.c.}
\end{align}
The expectation values of these one-body observables can then be read off directly from the elements of the steady-state correlation matrix.

The total particle current flowing, say, into the cold bath is represented by the observable $\hat{J}_{\rm tot} = \ensuremath{\mathrm{i}}[\hat{H}_{\rm tot}, \hat{N}_c]$, where $\hat{N}_c = \sum_{j\in c}\sum_q \hat{b}_{qj}^\dag \hat{b}_{qj}$ is the particle-number operator of the cold bath. Following a similar calculation to Sec.~\ref{app:solution}, we find the standard Landauer formula for the mean steady-state current~\cite{Ryndyk2016}
\begin{equation}\label{J_tot}
J_{\rm tot} = \int\frac{\mathrm{d}	\omega}{2\pi} \,\mathcal{T}(\omega) \left [\bar{n}_h(\omega) - \bar{n}_c(\omega)\right ],
\end{equation}
where $\mathcal{T}(\omega) = \mathrm{Tr}\left [ \mathbf{\Gamma}_c\cdot \mathbf{G} \cdot\mathbf{\Gamma}_h\cdot \mathbf{G}^\dagger\right ]$ is the transmission function. Particle conservation implies that $\sum_y J^X_{x,y} = J_{\rm tot}$ in the NESS.

\subsection{Particle-hole symmetry}
\label{app:PH_symmetry}

In the fermionic setting, the system obeys a generalized particle-hole (GPH) symmetry \cite{Chiu2016} when $\mu=\omega_0$. It is simplest to first analyse this situation assuming that $\omega_0=0$.  In this case,  the fermionic Hamiltonian [Eqs.~(1)--(3)] is invariant under the combined flavor swap and particle hole transformation,  $\hat{\Upsilon}\mathbf{\hat{a}}_{x,y}^\dagger\hat{\Upsilon}^{\dagger} = (\sigma_x \mathbf{\hat{a}}_{x,y})^{\rm T}$ or, explicitly,  $\hat{\Upsilon}(\hat{a}^\dagger_{x,y,\uparrow}, \hat{a}^\dagger_{x,y,\downarrow})\hat{\Upsilon}^{\dagger} = (\hat{a}_{x,y,\downarrow}, \hat{a}_{x,y,\uparrow})$. Meanwhile, the current operators defined in Eqs.~\eqref{current_op_X} and~\eqref{current_op_Y} are odd under this transformation, i.e.~$\hat{\Upsilon}\hat{J}^{X,Y}_{x,y} \hat{\Upsilon}^{\dagger} = -\hat{J}^{X,Y}_{x,y}$. It follows that any state $\hat{\rho}$ that is GPH-invariant, satisfying $[\hat{\rho},\hat{\Upsilon}] = 0$, has vanishing currents. Indeed, schematically we have $\braket{\hat{J}} = \mathrm{Tr}[\hat{J}\hat{\rho}] = -\mathrm{Tr}[\hat{\Upsilon}\hat{J}\hat{\Upsilon}^{\dagger}\hat{\rho}] = -\braket{\hat{J}}$, where $\hat{J}$ stands for any current operator that is GPH-odd.

In general, the correlation matrix can be written in block form as
\begin{equation}\label{C_block}
	\mathbf{C} = \left ( \begin{matrix}
		\mathbf{C}_{\uparrow \uparrow} & \mathbf{C}_{\uparrow \downarrow}\\		\mathbf{C}_{\downarrow \uparrow} & \mathbf{C}_{\downarrow \downarrow}
	\end{matrix}\right ),
\end{equation}
where $\mathbf{C}_{\uparrow \uparrow}$ is the correlation matrix for the $\uparrow$ flavor, $\mathbf{C}_{\uparrow \downarrow} = \mathbf{C}_{\downarrow \uparrow}^\dagger $ describes coherences between the different flavors, etc. The condition for GPH-invariance, $[\hat{\rho},\hat{\Upsilon}] = 0$, implies that the fermionic correlation matrix satisfies
\begin{equation}\label{PH_symmetry_C}
	\mathbf{C} = \mathbf{1} - \mathbf{S}\cdot\mathbf{C}^*\cdot \mathbf{S},
\end{equation}
where $\mathbf{S}$ is a matrix that swaps the flavors, and $\mathbf{C}^*$ is the complex conjugate of $\mathbf{C}$. With the basis order implicit in Eq.~\eqref{C_block}, we have
\begin{equation}\label{S_matrix}
	\mathbf{S} = \left ( \begin{matrix}
		 \mathbf{0} & \mathbf{1}\\	\mathbf{1} & \mathbf{0}
	\end{matrix}\right ),
\end{equation}
where $\mathbf{0}$ is the zero matrix. The GPH symmetry of the many-body Hamiltonian translates to the condition
\begin{equation}\label{PH_single_particle_H}
	\mathbf{H} = -\mathbf{S}\cdot\mathbf{H}^*\cdot\mathbf{S}.
\end{equation}
In turn, this implies the following constraint on the retarded Green function:
\begin{equation}\label{PH_Green_function}
	\mathbf{G}(-\omega) = - \mathbf{S}\cdot\left [\mathbf{G}(\omega)\right ]^*\cdot\mathbf{S}.
\end{equation}
Note that here we assume the self-energy is also GPH-symmetric, i.e.~$\mathbf{\Sigma}(\omega) = -\mathbf{S}\cdot[\mathbf{\Sigma}(-\omega)]^*\cdot\mathbf{S}$. In particular, this applies when the system-bath coupling is identical for both flavor states, so that $[\mathbf{\Sigma},\mathbf{S}] =0$, and when the spectral densities are symmetric, $\gamma_j(\omega) = \gamma_j(-\omega)$.

At the GPH-symmetric point, $\mu=0$, the distribution functions satisfy $\bar{n}_{h,c}(\omega) =  1 - \bar{n}_{h,c}(-\omega)$. Thus, reversing the dummy integration variable as $\omega \to -\omega$ in Eq.~(4), we obtain
\begin{align}
	\label{C_PH_derivation}
\mathbf{C} & = \int\frac{\mathrm{d}\omega}{2\pi} \, \mathbf{G}(\omega) \cdot \Big [\mathbf{ \Gamma}_h \bar{n}_h(\omega) + \mathbf{\Gamma}_c \bar{n}_c(\omega)\Big]\cdot \mathbf{G}^\dagger(\omega)\notag \\
 & = \int\frac{\mathrm{d}\omega}{2\pi} \, \mathbf{G}(\omega) \cdot \Big [\mathbf{ \Gamma} - \mathbf{ \Gamma}_h \bar{n}_h(-\omega) - \mathbf{\Gamma}_c \bar{n}_c(-\omega)\Big]\cdot\mathbf{G}^\dagger(\omega)\notag \\
& =  \mathbf{1} - \int\frac{\mathrm{d}\omega}{2\pi} \, \mathbf{G}(-\omega) \cdot \Big[ \mathbf{ \Gamma}_h \bar{n}_h(\omega) + \mathbf{\Gamma}_c \bar{n}_c(\omega)\Big]\cdot\mathbf{G}^\dagger(-\omega)\notag \\
& =  \mathbf{1} - \mathbf{S}\cdot\mathbf{C}^*\cdot \mathbf{S}.
\end{align}
Above, we have made use of Eq.~\eqref{PH_Green_function} as well as the identity $\mathbf{G}\cdot \mathbf{\Gamma}\cdot \mathbf{G}^\dagger = \ensuremath{\mathrm{i}} \left (\textbf{G} - \mathbf{G}^\dagger\right )$ and the property
\begin{equation}\label{commutation_Gret}
\int\frac{\mathrm{d}\omega}{2\pi }\ensuremath{\mathrm{i}}  \left (\textbf{G} - \mathbf{G}^\dagger \right ) = \mathbf{1},
\end{equation}
which follows from the (anti-)commutation relations~\cite{Kamenev2009}. Since the NESS is the unique Gaussian state with correlation matrix given by Eq.~\eqref{C_PH_derivation}, it follows that the state is GPH-symmetric when $\mu= \omega_0 = 0$. This implies that all currents vanish. Consistent with this observation, GPH symmetry implies that the transmission function is even, $\mathcal{T}(\omega) = \mathcal{T}(-\omega)$, and therefore the integral in Eq.~\eqref{J_tot} is odd and equates to zero.

The same GPH symmetries and associated vanishing currents are obtained when $\mu = \omega_0 \neq 0$. This can be shown simply by shifting the frequency variable $\omega\to \omega +\mu$ in Eq.~\eqref{C_PH_derivation}. The final result then follows from the GPH symmetry of the shifted Hamiltonian $\mathbf{H}' = \mathbf{H}-\omega_0\mathbf{1}$. We have numerically confirmed the GPH symmetry of the NESS in the wide-band limit when $\mu=\omega_0$, both in weak- and strong-coupling regimes.

\subsection{Weak-coupling limit}
\label{app:weak_coupling}

Following Ref.~\cite{Dhar2012}, we show that in the weak-coupling limit the NESS reduces to an energy-diagonal density matrix that coincides with the solution of the secular-Born-Markov master equation~\cite{Breuer2002,RivasHuelga,Rivas2017}. The system Hamiltonian can be written in its spectral representation
\begin{equation}\label{single_particle_H_diag}
	\mathbf{H} = \mathbf{U}\cdot\mathbf{\Omega}\cdot\mathbf{U}^\dagger,
\end{equation}
where $\mathbf{\Omega} = {\rm diag}\{\omega_\alpha\}$ and $\mathbf{U}$ comprise the eigenvalues and eigenvectors of $\mathbf{H}$. The many-body Hamiltonian in this basis reads as
\begin{equation}\label{many_body_H_diag}
\hat{H} = \sum_{\alpha} \omega_\alpha \ensuremath{\hat{c}^\dagger}_\alpha \c_\alpha, \qquad \c_\alpha = \sum_{j} U^*_{j\alpha}\a_j	.
\end{equation}

Now consider the retarded Green function~\eqref{G_ret} in the wide-band limit, $\mathbf{\Sigma} = -\ensuremath{\mathrm{i}} \mathbf{\Gamma}/2 \propto \gamma$. As $\gamma\to 0$, $\mathbf{G}(\omega)$ becomes approximately diagonal in the energy representation. Following identical arguments to Ref.~\cite{Dhar2012}, we obtain the steady-state correlation matrix in the energy eigenbasis as
\begin{equation}\label{diagonal_NESS}
\braket{\ensuremath{\hat{c}^\dagger}_\alpha \c_{\alpha'}} = \delta_{\alpha\alpha'} \frac{s_\alpha \bar{n}_h(\omega_\alpha) + r_\alpha \bar{n}_c(\omega_\alpha)}{s_\alpha + r_\alpha},
\end{equation}
where $s_\alpha$ and $r_\alpha$ describe the dimensionless coupling strength of eigenmode $\alpha$ to the hot and cold baths, respectively:
\begin{align}
\label{s_alpha}
& s_\alpha = \gamma^{-1}\sum_{j,k}U^*_{j\alpha} [\Gamma_h]_{jk}	U_{k\alpha} = \sum_{y=1}^{L_Y}\sum_{\sigma=\uparrow,\downarrow} |U_{(1,y,\sigma),\alpha}|^2,\\
\label{r_alpha}
& r_\alpha = \gamma^{-1}\sum_{j,k}U^*_{j\alpha} [\Gamma_c]_{jk}	U_{k\alpha} = \sum_{y=1}^{L_Y}\sum_{\sigma=\uparrow,\downarrow} |U_{(L_X,y,\sigma),\alpha}|^2.
\end{align}

We now consider the Lindblad master equation obtained within the Born-Markov and secular approximations. Following the standard derivation~\cite{Breuer2002,RivasHuelga} yields the master equation
\begin{equation}\label{Lindblad}
	\dot{\hat{\rho}} = -\ensuremath{\mathrm{i}} [\hat{H}, \hat{\rho}] + \mathcal{L}_h\hat{\rho} + \mathcal{L}_c\hat{\rho},
\end{equation}
with dissipators for the hot and cold baths given respectively by
\begin{align}
	\label{hot_dissipator}
	\mathcal{L}_h = \gamma \sum_\alpha s_\alpha \left ( \bar{n}_h(\omega_\alpha)\mathcal{D}[\ensuremath{\hat{c}^\dagger}_\alpha] + [1\mp \bar{n}_h(\omega_\alpha)]\mathcal{D}[\c_\alpha]  \right ),\\
		\label{cold_dissipator}
	\mathcal{L}_c = \gamma \sum_\alpha r_\alpha \left ( \bar{n}_c(\omega_\alpha)\mathcal{D}[\ensuremath{\hat{c}^\dagger}_\alpha] + [1\mp \bar{n}_c(\omega_\alpha)]\mathcal{D}[\c_\alpha]  \right ),
\end{align}
where $\mathcal{D}[\hat{L}]\bullet  = \hat{L}\bullet \hat{L}^\dagger - \tfrac{1}{2}\{\hat{L}^\dagger \hat{L},\bullet\}$ is a Lindblad superoperator and the plus (minus) sign is for bosons (fermions). It is easy to check that the steady-state solution of Eq.~\eqref{Lindblad}, $\dot{\hat{\rho}}=0$, is a Gaussian state with correlation matrix given exactly by Eq.~\eqref{diagonal_NESS}.

As is well known, the Lindblad equation~\eqref{Lindblad} can be written as
\begin{equation}\label{Lindblad_nonHermitian}
\dot{\hat{\rho}} = -\ensuremath{\mathrm{i}} \hat{H}_{\rm eff} \hat{\rho} + \ensuremath{\mathrm{i}} \hat{\rho} \hat{H}_{\rm eff}^\dagger + \mathcal{J}_h\hat{\rho} + \mathcal{J}_c\hat{\rho},
\end{equation}
where the non-Hermitian Hamiltonian is
\begin{align}\label{H_nonHermitian}
\hat{H}_{\rm eff} & = \hat{H} - \frac{\ensuremath{\mathrm{i}}}{2} \sum_{\alpha} \gamma s_\alpha \left[ \bar{n}_h(\omega_\alpha) + (1\mp 2\bar{n}_{h}(\omega_\alpha)) \hat{c}^\dagger_\alpha \hat{c}_\alpha\right ] \notag \\
& \qquad - \frac{\ensuremath{\mathrm{i}}}{2} \sum_{\alpha} \gamma r_\alpha \left[ \bar{n}_c(\omega_\alpha) + (1\mp 2\bar{n}_{c}(\omega_\alpha)) \hat{c}^\dagger_\alpha \hat{c}_\alpha\right ],
\end{align}
and the superoperators $\mathcal{J}_{c,h}$ represent ``quantum jumps'' due to exchange of particles with the baths~\cite{Plenio1998},
\begin{align}
\label{jump_superoperators}
\mathcal{J}_h\hat{\rho} & = \sum_\alpha \gamma s_\alpha \left [\bar{n}_h(\omega_\alpha) \hat{c}_\alpha^\dagger \hat{\rho}\hat{c}_\alpha + \left (1 + \bar{n}_h(\omega_\alpha)\right ) \hat{c}_\alpha \hat{\rho}\hat{c}_\alpha^\dagger  \right ] ,\\
\mathcal{J}_c\hat{\rho} & = \sum_\alpha \gamma r_\alpha \left [\bar{n}_c(\omega_\alpha) \hat{c}_\alpha^\dagger \hat{\rho}\hat{c}_\alpha + \left (1 + \bar{n}_c(\omega_\alpha)\right ) \hat{c}_\alpha \hat{\rho}\hat{c}_\alpha^\dagger  \right ] .
\end{align}
By including both the non-Hermitian Hamiltonian and the jump terms, we ensure a completely positive and physically consistent description of the full many-body state of the non-equilibrium lattice system, valid at arbitrary times. This is in contrast to non-Hermitian models that neglect the jumps (valid for rare post-selected trajectories of continuously measured systems) or the effective non-Hermitian Hamiltonians appearing in single-particle Green functions, the topological properties of which have recently been of intense interest; e.g., see Ref.~\cite{MartinezAlvarez2018} and references therein.

\begin{figure*}
	\begin{minipage}{0.31\linewidth}
		\includegraphics[width=\linewidth]{impurity_fermion_counter}
	\end{minipage}\quad
	\begin{minipage}{0.31\linewidth}
		\includegraphics[width=\linewidth]{impurity_fermion_highT}
	\end{minipage}\quad
	\begin{minipage}{0.31\linewidth}
		\includegraphics[width=\linewidth]{impurity_fermion_trivial}
	\end{minipage}
	\caption{Current patterns for the same asymmetric impurity distribution shown in Fig.~4(c) but with different parameters. (a,b) Fermionic currents in the topologically nontrivial regime, $m=t$, with on-site energy $\omega_0=10t$, chemical potential $\mu = \omega-0.5t$, and (a) low temperatures, $T_h=t$ and $T_c=0.01t$, or (b) higher temperatures, $T_h=100t$ and $T_c=t$. (c) Same as (a) but showing bosons (red) and fermions (blue) in a topologically trivial regime, $m=3t$. \label{fig:fermion_impurity_quiver}}
\end{figure*}

\subsection{Weak-coupling distribution functions}
\label{app:distribution}

\begin{figure}[b]
	\begin{minipage}{0.495\linewidth}
		\includegraphics[width=\linewidth]{dist_symm}
	\end{minipage}
	\begin{minipage}{0.495\linewidth}
		\includegraphics[width=\linewidth]{dist_asymm}
	\end{minipage}
	\caption{Distribution functions in the NESS for bosons (red) and fermions (blue). (a)~Symmetric impurity configuration shown in Fig.~4(b), with $T_h=t$, $T_c=0.01t$, $\omega_0=10t$ and $\mu=\omega_0-0.1t$. (b)~Asymmetric impurity configuration shown in Fig.~4(c), but with $\mu=\omega_0-0.5t$ as in Fig.~\ref{fig:fermion_impurity_quiver}. Circles show the data for $T_h=t$ and $T_c=0.01t$ [Fig.~\ref{fig:fermion_impurity_quiver}(a)], while light blue squares show the result for fermions at high temperature, $T_h=100t$ and $T_c = t$ [Fig.~\ref{fig:fermion_impurity_quiver}(b)]. The values of $n(\omega_\alpha)$ for bosons are multiplied by $10^3$ to be visible on the same axes. \label{fig:distributions}}
\end{figure}

We now investigate in more detail the form of the weak-coupling distribution function, $n(\omega_\alpha) = \braket{\hat{c}_\alpha^\dagger\hat{c}_\alpha}$, defined by Eq.~\eqref{diagonal_NESS}.  When the reservoirs are at the same temperature, $\bar{n}_h=\bar{n}_c = \bar{n}$, we obtain simply $ n(\omega_\alpha) = \bar{n}(\omega_\alpha)$. This is a smooth, monotonically decreasing function of frequency that is independent of the mode-specific couplings $s_\alpha$ and $r_\alpha$.  As a consequence, the equilibrium currents---as seen, for example, in Fig.~2(c)---are unaffected by the addition of defects in any configuration~\cite{Rivas2017}.

Out of equilibrium, an analogous result is obtained in the clean QWZ model described by Eqs.~(1)--(3) due to the symmetries $\hat{\Pi}\hat{\Theta}\hat{\Sigma}_y$ and $\hat{\Pi}\hat{R}_\pi$. In particular, the symmetry $\hat{\Pi}\hat{\Theta}\hat{\Sigma}_y$ implies that $|U_{(1,y,\sigma),\alpha}|^2 = |U_{(L_X,y,\sigma),\alpha}|^2 $ and therefore $s_\alpha = r_\alpha$ identically. Similarly, the symmetry under $\hat{\Pi}\hat{R}_\pi$ means that $|U_{(1,y,\sigma),\alpha}|^2 = |U_{(L_X,L_Y+1-y,\sigma),\alpha}|^2$, and again $s_\alpha = r_\alpha$ upon summing over $y$ in Eqs.~\eqref{s_alpha} and~\eqref{r_alpha}. As a consequence, the mode occupations~\eqref{diagonal_NESS} become independent of $s_\alpha$ and $r_\alpha$, being given simply by the average of the reservoir distribution functions:
\begin{equation}\label{mode_occ_symmetry}
	n(\omega_\alpha) = \frac{1}{2}\left [\bar{n}_c(\omega_\alpha) + \bar{n}_h(\omega_\alpha) \right ].
\end{equation}
This is again a function of the frequency $\omega_\alpha$ only, and therefore leads to stable current patterns as in the equilibrium case~\cite{Rivas2017}. Any distribution of defects that respects one of the above discrete symmetries will also yield $s_\alpha = r_\alpha$, leading to similar behaviour.

Note that the relevant symmetries are determined both by the Hamiltonian and by the configuration of the baths driving the system out of equilibrium. For example, the transformation $\hat{\Theta}\hat{\Sigma}_x$ is a symmetry of the QWZ Hamiltonian, where $\hat{\Sigma}_x = \hat{R}_{\pi} \hat{\Sigma}_y$ is a reflection about the $x$ axis. However, because this reflection does not interchange the hot and cold baths it is irrelevant for the NESS. Therefore, defect distributions that are only symmetric under $\hat{\Sigma}_x$ destroy the bosonic edge currents whenever $T_c\neq T_h$.

Examples of the nonequilibrium distribution functions obtained for a symmetric defect configuration are shown in Fig.~\ref{fig:distributions}(a), corresponding to the system depicted in Fig.~4(b). As expected, the distributions are well behaved functions of frequency. Conversely, Fig.~\ref{fig:distributions}(b) shows that when the defect configuration does not respect either of the discrete symmetries, as in Figs.~4(c) or~\ref{fig:fermion_impurity_quiver}, the mode occupations are erratically varying functions of frequency. This disrupts the erasure effect and leads to the destruction of the boundary current patterns, in general.

Fig.~\ref{fig:distributions} also shows that the fermionic distributions exhibit a sudden change at the chemical potential whenever $T_c$ is sufficiently small. This sharp feature explains the absence of counter-currents around the impurities in Figs.~4(b,c). In this case, the chemical potential is just below the particle-hole symmetric point $\mu=\omega_0$. This leads to a single ``hole'' that occupies an edge state, visible in Fig.~\ref{fig:distributions} as the isolated point just below $\omega_\alpha=\omega_0$. The edge current pattern observed in Fig.~4(b) is primarily due to this hole state which, occupying an edge mode, is robust against any bulk perturbation. Conversely, when the chemical potential is shifted down to $\mu = \omega_0- 0.5t$, several holes are created in other edge states [Fig.~\ref{fig:distributions}(b), blue circles], including edge states localized around the impurities. Consequently, we see the appearance of counter-currents shielding the impurities in Fig.~\ref{fig:fermion_impurity_quiver}(a).

At higher temperature, the sharp feature in $n(\omega_\alpha)$ disappears [Fig.~\ref{fig:distributions}(b), light blue squares]. The current distribution is then determined mainly by bulk states and thus depends intimately on the symmetries of the impurity configuration. If the configuration is asymmetric, as in Fig.~\ref{fig:fermion_impurity_quiver}(b), the boundary currents are strongly disrupted.

In the topologically trivial regime, there are no edge states and both bosons and fermions behave very similarly with respect to the introduction of defects. In particular, asymmetric impurity distributions tend to inhibit the formation of clean boundary current patterns. An example is shown in Fig.~\ref{fig:fermion_impurity_quiver}(c).

\begin{figure*}
	\includegraphics[width=0.99\linewidth]{semiclassical}
	\caption{The emergence of edge currents in the weak coupling limit can be understood by a semiclassical argument. (a) The factor $n[\omega_-(\mathbf{k})]-n[\omega_+(\mathbf{k})]$ in Eq. \eqref{semi-classical_edge_current} for the bosonic case with $\omega_0=10t$, $T_h=t$, $T_c=0.01t$ and $m=2.1t$. (b) The Berry curvature ${\mathcal{F}\!\!}_-(\mathbf{k})$, which has a balanced positive and negative contributions over B.Z. (the zero value plane is depicted in light blue color) as for  $m>2t$ the model is in a trivial phase and ${\mathcal{F}\!\!}_-(\mathbf{k})$ integrates out to zero. (c) The product $\{n[\omega_-(\mathbf{k})]-n[\omega_+(\mathbf{k})]\}{\mathcal{F}\!\!}_-(\mathbf{k})$, which breaks such a balance, leads to a nonvanishing integral in Eq. \eqref{semi-classical_edge_current}. (d) The semiclassical current pattern generated according \eqref{semi-classical_edge_current} for the potential $V(x,y)=(x/8)^{10}+(y/8)^{10}$.  \label{fig:semiclassical}}
\end{figure*}

\subsection{Semiclassical picture}
\label{app:semiclassical}

The behaviour in the weak coupling can be understood by using a semiclassical picture \cite{Xiao2010,Price2012}.  Under periodic boundary conditions,  the QWZ Hamiltonian can be written in the form
\begin{equation}
\hat{H}=\sum_{\mathbf{k}\in \mathrm{B.Z.}} \mathbf{a}_{\mathbf{k}}^\dagger\cdot \mathbf{H}(\mathbf{k})\cdot \mathbf{a}_{\mathbf{k}},
\end{equation}
with $\mathbf{k}=(k_X,k_Y)$ the quasimomentum and
\begin{equation}
	\label{single-particle-Bloch-H}
\mathbf{H}(\mathbf{k})=\omega_0\mathbf{1}+\mathbf{h}(\mathbf{k})\cdot\bm{\sigma},
\end{equation}
where
\begin{equation}
\mathbf{h}(\mathbf{k})=[t_Y\sin(k_Y),t_X\sin(k_X),m+t_X\cos(k_X)+t_Y\cos(k_Y)].
\end{equation}
Note that we have taken a slightly different Hamiltonian with respect to the original reference \cite{Qi2006}.  The two energy bands of the model are then given by
\begin{equation}
\omega_{\pm}(\mathbf{k})=\omega_0\pm |\mathbf{h}(\mathbf{k})|,
\end{equation}
and for the isotropic case considered throughout this work  ($t_X=t_Y=t>0$),
\begin{equation}
|\mathbf{h}(\mathbf{k})|=t\sqrt{\sin^2(k_X)+\sin^2(k_Y)+[m/t+\cos(k_X)+\cos(k_Y)]^2}.
\end{equation}
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})$ (assuming the single-particle Hamiltonian~\eqref{single-particle-Bloch-H} is gapped,  i.e.  $m\neq0,\pm 2t$, and $m$ is small enough to justify the semiclassical approximation) as
\begin{equation}\label{Vsem}
\mathbf{v}_\alpha=\frac{\partial \omega_\alpha(\mathbf{k})}{\partial\mathbf{k}}-\mathbf{F}\times \bm{\mathcal{F}}\!\!_\alpha(\mathbf{k}).
\end{equation}
Here, $\bm{\mathcal{F}}\!\!_\alpha(\mathbf{k})=[0,0,{\mathcal{F}\!\!}_\alpha(\mathbf{k})]$, with ${\mathcal{F}\!\!}_\alpha(\mathbf{k})$ the Berry curvature of the band, which, for a two-band model, can be written in the form
\begin{align}
{\mathcal{F}\!\!}_\pm(\mathbf{k})&=\mp\frac{1}{2}\mathbf{n}(\mathbf{k})\cdot\left[\frac{\partial\mathbf{n}(\mathbf{k})}{\partial k_X}\times \frac{\partial\mathbf{n}(\mathbf{k})}{\partial k_Y}\right],
\end{align}
with $\mathbf{n}(\mathbf{k}):=\mathbf{h}(\mathbf{k})/|\mathbf{h}(\mathbf{k})|$. For the QWZ Hamiltonian we find
\begin{equation}
{\mathcal{F}\!\!}_\pm(\mathbf{k})=\frac{\pm\big[\cos(k_X)+\cos(k_Y)+\tfrac{m}{t}\cos(k_X)\cos(k_Y)\big]}{2\big\{\!\sin^2(k_X)+\sin^2(k_Y)+\big[\tfrac{m}{t}+\cos(k_X)+\cos(k_Y)\big]^2\big\}^\frac32}.
\end{equation}
We can now argue for the presence of an edge circulation under open boundary conditions by introducing a confining potential $V(x,y)$ in the semiclassical equation for $\mathbf{v}_\pm$ \eqref{Vsem}.  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(x,y)$~\cite{Xiao2010}. Moreoever, the potential induces a confining force $\mathbf{F}=-\bm{\nabla}V(x,y)$,  which approximately vanishes in the bulk and becomes very large in the boundary of the confined region,
\begin{align}
	\label{edge-bulk-velocity}
\mathbf{v}_\alpha^\mathrm{(bulk)}\simeq  \frac{\partial \omega_\alpha(\mathbf{k})}{\partial\mathbf{k}}, \quad \mathbf{v}_\alpha^\mathrm{(edge)}\simeq \bm{\nabla}V(x,y)\times \bm{\mathcal{F}}\!\!_\alpha(\mathbf{k}).
\end{align}
Since $\mathbf{F}$ is normal and points inward from the confined boundary,  a circulating current is induced along this boundary with direction given the sign of the Berry curvature, which depends on the band $\alpha$, the sign of $m$, and the specific value of $\mathbf{k}$: clockwise for ${\mathcal{F}\!\!}_\alpha(\mathbf{k})>0$ and counterclockwise for ${\mathcal{F}\!\!}_\alpha(\mathbf{k})<0$. Thus, for the NESS at the weak coupling regime,  we may estimate the current density (current per unit width) by \cite{Xiao2010}
\begin{align}
	\label{semiclassical_current}
\mathbf{I}(x,y)=\sum_{\alpha=\pm}\int_{\mathrm{B.Z.}}d^2\mathbf{k}\, n[\omega_\alpha(\mathbf{k})+V(x,y)] \mathbf{v}_\alpha(x,y,\mathbf{k}),
\end{align}
where the integral extends over the first Brillouin zone. Here we have assumed that the NESS occupation number $n$ only depends on energy. As discussed in Sec.~\ref{app:distribution}, this is the case under a $\hat{\Pi}\hat{\Theta}\hat{\Sigma}_y$ or a $\hat{\Pi}\hat{R}_\pi$ symmetric configuration.

Now, if we consider regions far from the boundary, so that $V(x,y)\ll \omega_\alpha(\mathbf{k})$, we may approximate Eq.~\eqref{semiclassical_current} by
\begin{align}
	\mathbf{I}(x,y)&\simeq\int_{\mathrm{B.Z.}}d^2\mathbf{k}\, \{n[\omega_-(\mathbf{k})]-n[\omega_+(\mathbf{k})]\} \mathbf{v}_-,
\end{align}
using the fact that $\mathbf{v}_+= -\mathbf{v}_-$. Since there is no privileged direction of $\mathbf{v}_-^{\mathrm{(bulk)}}$ on each energy shell: $\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. We conclude that
\begin{align}
	\label{semi-classical_edge_current}
	\mathbf{I}(x,y)&\simeq\int_{\mathrm{B.Z.}}d^2\mathbf{k}\, \left \{n[\omega_-(\mathbf{k})]-n[\omega_+(\mathbf{k})]\right \} \mathbf{v}_-^{(\rm edge)},
\end{align}
which clearly vanishes in the bulk of the system where $V(x,y)=0$ (see Eq.~\eqref{edge-bulk-velocity}).  Nevertheless, since $n(\omega)$ is a monotonically decreasing function of energy [Eq.~\eqref{mode_occ_symmetry}], there is an imbalance of the contributions to $\mathbf{I}(x,y)$ within the Brillouin zone in regions where $\mathbf{v}_\alpha^{\rm edge} \propto \nabla V(x,y) \neq 0$. Thus, Eq.~\eqref{semi-classical_edge_current} predicts a net chiral edge current,  see Fig. \ref{fig:semiclassical}. This erasure effect sustains the appearance of the chiral current in the symmetric situation, but it is spoilt in the nonsymmetric case,  as the occupation number is no longer only a function of the energy alone.

\section{Results}

	\begin{figure*}
		\begin{minipage}{0.31\linewidth}
			\includegraphics[width=\linewidth]{quiver_strong_combi}
		\end{minipage}\quad
		\begin{minipage}{0.31\linewidth}
			\includegraphics[width=\linewidth]{quiver_weak_combi}
		\end{minipage}\quad
		\begin{minipage}{0.31\linewidth}
			\includegraphics[width=\linewidth]{quiver_eq_combi}
		\end{minipage}
		\caption{Current distributions on a lattice of size $L = 8$. Red arrows show boson currents in a topologically trivial phase ($m=3t$), blue arrows show fermion currents in a topologically nontrivial regime ($m=t$). Arrow sizes indicate the magnitude of the currents scaled relative to the largest value within each plot. (a,b) Nonequilibrium current profile with different temperatures $T_h =  t$ and $T_c= 0.01t$, on-site energy $\omega_0=10t$, chemical potential $\mu = \omega_0-0.01t$ for fermions (we always set $\mu=0$ for bosons), and coupling strength (a)~$\gamma = 0.5t$ and (b)~$ \gamma = 0.005t$. (c)~Equilibrium case with temperatures $T_c=T_h=t$, weak coupling, $\gamma = 0.005t$, and other parameters identical to (a,b). \label{fig:quiver_coupling}}
	\end{figure*}

	\noindent {\em Results.---}In the following examples, we focus on the symmetric case with $L_X = L_Y = L$ and $t_X = t_Y=t$, so that $|m|<2t$ defines the topologically non-trivial phase. We also fix $\omega_0=10t$ and consider relatively low temperatures, $T_{c,h}\lesssim t$, to accentuate the role of band topology and particle exchange statistics.

	\begin{figure}[b]
		\begin{minipage}{0.495\linewidth}
			\includegraphics[width=\linewidth]{edge_bulk_currents_bosons}
		\end{minipage}
		\begin{minipage}{0.495\linewidth}
			\includegraphics[width=\linewidth]{edge_bulk_currents_fermions}
		\end{minipage}
		\caption{Average edge currents, $J_{\rm edge}$ (solid lines), and bulk currents, $J_{\rm bulk}$ (dashed lines), as a function of the system-bath coupling strength, $\gamma$, for bosons (left panel) and fermions (right panel) in topologically non-trivial ($m/t=1$) and trivial ($m/t=2.1,3$) regimes, and with $T_h= t$, $T_c=0.01t$, $\omega_0=10t$, and $\mu = \omega_0+0.1t$. \label{fig:edge_bulk_currents}}
	\end{figure}

	Figs.~\ref{fig:quiver_coupling}(a,b) plot the nonequilibrium current distributions for two different values of the system-reservoir coupling, $\gamma$. Red arrows show the currents for a bosonic system with $m=3t$ (similar results are obtained for $|m|<2t$). Remarkably, the currents become progressively localized on the boundary of the system as $\gamma$ is reduced, even though the band topology is trivial. These edge currents also arise in the equilibrium case, $T_h=T_c$, as shown in Fig.~\ref{fig:quiver_coupling}(c). Qualitatively similar results are obtained for fermions, as shown for a nontrivial phase ($m=t$) by the blue arrows in Fig.~\ref{fig:quiver_coupling}. The direction of fermionic particle flow depends on the chemical potential; for the parameters in Figs.~\ref{fig:quiver_coupling}(a,b), the thermoelectric induced current flows in the opposite direction to the temperature gradient. Conversely, when $\mu=\omega_0$, all currents vanish due to particle-hole symmetry~\cite{SM}.

	In order to quantify the emergence of boundary currents more precisely, we define the average edge and bulk currents
	\begin{align}\label{edge_and_bulk_current}
		J_{\rm edge} &= \frac{1}{2L_X}\sum_{x} \left( J^{X}_{x,L_Y} - J^X_{x,1} \right ), 	\qquad J_{\rm bulk} = \frac{1}{L_Y}\sum_{y} J^{X}_{x,y}.
	\end{align}
	The choice of $x$-coordinate in the definition of $J_{\rm bulk}$ is arbitrary due to particle-number conservation; we take $x=\lfloor L_X/2\rfloor$. The total current flowing between the two reservoirs is given by $J_{\rm tot} = L_Y J_{\rm bulk}$.

	We plot the edge and bulk currents in Fig.~\ref{fig:edge_bulk_currents} for bosonic and fermionic systems as a function of the system-reservoir coupling strength, $\gamma$. In the bosonic case, both bulk and edge currents increase with $m$ and are thus larger in the topologically trivial phase. The situation is reversed for fermions, with the largest edge currents obtained in the topologically nontrivial phase. In all cases, the bulk currents are proportional to $\gamma$ and thus vanish as $\gamma\to 0$, while the edge currents remains invariant for a wide range of values of $\gamma$ and persist in the weak-coupling limit. These quantitative results therefore confirm the qualitative picture of Fig.~\ref{fig:quiver_coupling}, i.e., edge currents arise for small coupling irrespectively of band topology or particle statistics. In topologically trivial phases, edge currents begin to dominate once $\gamma$ becomes comparable to the level spacing of the single-particle Hamiltonian. However, in the topologically nontrivial phase for fermions near half filling,  currents remain localized on the edges even for strong coupling, $\gamma\sim  t$.

	\begin{figure*}
		\begin{minipage}{0.31\linewidth}
			\includegraphics[width=\linewidth]{impurity_edge}
		\end{minipage}\quad
		\begin{minipage}{0.31\linewidth}
			\includegraphics[width=\linewidth]{impurity_sym}
		\end{minipage} \quad
		\begin{minipage}{0.31\linewidth}
			\includegraphics[width=\linewidth]{impurity_asym}
		\end{minipage}
		\caption{Effect of impurities on the current distribution in the weak-coupling limit, with impurity positions marked by filled black circles. Red arrows show boson currents and blue arrows show fermion currents. In all plots, $m=t$, $T_h= t$, $T_c=0.01t$, $\omega_0=10t$ and $\mu = \omega_0 - 0.1t$.\label{fig:quiver_impurities}}
	\end{figure*}

	Finally, we investigate whether these effects withstand the presence of imperfections, a property which defines the notion of dissipatively robust currents. Specifically, we examine how the currents change when static impurities are added to the lattice. Impurities are modelled by a large on-site energy shift, i.e.,~a term $\Delta {\bf\a}^\dagger_{x,y}
	\cdot {\bf \a}_{x,y}$ added to the Hamitonian, where $\Delta \gg \omega_0, t$ and $(x,y)$ are the coordinates of the impurity site~\cite{Rivas2017}. We focus hereafter on the weak-coupling limit where the edge currents are most prominent in the absence of impurities. In this limit, an analytical expression for the NESS~\cite{Dhar2012} can be obtained, which is identical to the solution of the Lindblad equation derived under the Born-Markov and secular approximations~\cite{SM}.

	In Fig.~\ref{fig:quiver_impurities} we show three examples of how impurities affect the nonequilibrium current distribution in the topologically nontrivial phase. We observe substantial differences between bosonic and fermionic systems when one of the reservoirs is at low temperature. The bosonic edge currents are preserved only when the impurity distribution is invariant under one of the symmetries~$\hat{\Pi}\hat{\Theta}\hat{\Sigma}_y$ or $\hat{\Pi}\hat{R}_\pi$, which correspond to the purely spatial symmetries $\hat{\Sigma}_y$ and $\hat{R}_\pi$ since the defects satisfy $\hat{\Pi}$ and $\hat{\Theta}$ automatically. In the case of impurities placed on the edge (and assuming the symmetries are respected in the bosonic case), the currents simply detour around the impurity sites [Fig.~\ref{fig:quiver_impurities}(a)]. If impurities are placed in the bulk of a bosonic system in an appropriately symmetric way, counter-currents shield the impurity sites by circulating in the opposite direction to the edge currents~[Fig.~\ref{fig:quiver_impurities}(b)]. In the absence of these symmetries, the bosonic edge currents are strongly disrupted [Fig.~\ref{fig:quiver_impurities}(c)]. For fermions in a topological phase, the situation is starkly different when at least one reservoir is at low temperature: impurities placed in the bulk have no effect whatsoever on the current distribution~[Figs.~\ref{fig:quiver_impurities}(b,c)]. In the topologically trivial phase or at high temperature, both bosonic and fermionic edge currents enjoy the same symmetry-protected robustness as  bosons in the nontrivial phase.

\noindent {\em Discussion.---}The question posed in the introduction on the nontopological origin of the edge currents can now be answered. In the weak-coupling limit, nonequilibrium coherences between energy eigenstates become small~\cite{Dhar2012} and $J_{\rm bulk}$, which is directly proportional to these coherences~\cite{Mitchison2018}, is negligible in comparison to the contributions from individual energy eigenstates. Moreover, the single-particle eigenmodes in the QWZ model present a nonzero Berry curvature $\bm{\mathcal{F}}_\alpha$ for $m\neq 0$~\cite{Qi2006}, which breaks time-reversal symmetry and endows the eigenmodes with a particular chirality. This can be understood in a semi-classical picture~\cite{Haldane2004}, where a wavepacket propagates with velocity $\mathbf{v}_\alpha = \partial_{\mathbf{k}} \omega_\alpha +  \nabla V \times \bm{\mathcal{F}}_\alpha$, with $\mathbf{k}$ the wavevector and $V(x,y)$ a confining potential describing the edge of the system~\cite{Xiao2010}. Assuming that the non-equilibrium distribution function only depends on energy, the net current due to the derivative term $\partial_{\mathbf{k}} \omega_\alpha$ can be shown to vanish in the bulk due to cancellations between positive and negative wavevectors~\cite{SM}. However, the contribution of the axial force $\nabla V \times \bm{\mathcal{F}}_\alpha$ is non-zero near the boundary, thus generating a net chiral edge current~\cite{SM}. This erasure effect arises whenever there is broken time-reversal symmetry and a smooth, energetically monotonic distribution of chiral eigenmodes~\cite{Rivas2017}; it is thus independent of exchange statistics or band topology. Intuitively, the effect can be understood in terms of circulating currents that cancel in the bulk but add constructively on the boundary [Fig.~\ref{fig:QWZ}(b)].

A chiral current between two thermal reservoirs unavoidably entails a local cross-current phenomenon:  on one edge there are particles flowing from the cold bath to the hot bath. In bosonic systems, this implies a ``violation'' of the second law of thermodynamics within a one-dimensional subsystem, i.e.,~energy flows against the temperature gradient along one edge even though the total rate of entropy production is positive. This effect was reported in Refs.~\cite{Rivas2017,Xing2020} for bosonic lattices governed by the Hofstadter Hamiltonian. Now, we see that the key requisite for this behavior is not the nontrivial band topology of the lattice, but rather the nonzero Berry curvature in a regime of weak system-bath coupling [Fig.~\ref{fig:quiver_coupling}(b,c)].

The edge currents are robust to impurities if the defect distribution satisfies either of the symmetries $\hat{\Pi}\hat{\Theta}\hat{\Sigma}_y$ or $\hat{\Pi}\hat{R}_\pi$. Note that the relevant nonequilibrium symmetries are determined both by the Hamiltonian and by the configuration of the baths. These symmetries leave the NESS invariant and lead to a steady-state correlation matrix that is independent of the spatial orientation of the reservoirs. Under these conditions, the nonequilibrium distribution function is simply the average of the reservoir distributions~\cite{SM}, i.e.,~$n(\omega_\alpha) = \tfrac{1}{2}\left [ \bar{n}_h(\omega_\alpha) + \bar{n}_c(\omega_\alpha)\right ]$ is the population of the eigenmode with frequency $\omega_\alpha$. This smooth distribution function yields boundary currents due to the erasure effect in the bulk. In the absence of symmetry, however, the contribution of each mode depends not only on energy but also on the spatial profile of the corresponding wavefunction. This creates an erratically varying distribution function that destroys the erasure effect~\cite{SM}.

For fermions at low temperature $T_c$, $n(\omega_\alpha$) undergoes a sharp change near $\omega_\alpha \approx \mu$. In a topologically nontrivial phase near half filling, this feature acts as a filter that populates only one or two edge states~\cite{SM}, which are topologically protected against perturbations in the bulk. The corresponding boundary current is thus completely unaffected by impurities [Figs.~\ref{fig:quiver_impurities}(b,c)]. Remarkably, this holds for arbitrarily large $T_h$, providing an instance of a topologically protected property that is also stable under dissipation. The number and character of edge modes that contribute to the current pattern is determined by the chemical potential~\cite{SM}.

In summary, our nonperturbative analysis has shown that robust boundary currents emerge at weak coupling in the QWZ lattice driven out of equilibrium by a thermal gradient. Our results reveal a novel dissipative mechanism for symmetry-protected edge transport, which arises in bosonic and fermionic systems with and without nontrivial band topology. This makes quantum simulators based on either fermionic or bosonic degrees of freedom attractive candidates to observe this exotic boundary physics.

\textit{Acknowledgements.---}M.~T.~M. acknowledges funding from the ERC Starting  Grant ODYSSEY (Grant  Agreement No.~758403) and the EPSRC-SFI 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 Lonsdale cluster maintained by the Trinity Centre for High Performance Computing. This cluster was funded through grants from SFI.

	\bibliographystyle{apsrev4-2}