#  Developing cooperative policies for multi-stage reinforcement learning
tasks

Jordan Erskine1, Chris Lehnert2 Manuscript received: December, 23, 2021;
Revised March, 26, 2022; Accepted April, 19, 2022.This paper was recommended
for publication by Editor Jens Kober upon evaluation of the Associate Editor
and Reviewers’ comments. [Note that the Editor is the Senior Editor who
communicated the decision; this is not necessarily the same as the Editor-in-
Chief.]1J. Erskine and C. Lehnert are with the Queensland University of
Technology (QUT), Brisbane, Australia and affiliated with Queensland Centre of
Robotics (QCR<EMAIL_ADDRESS><EMAIL_ADDRESS>Object Identifier (DOI): see top of this page.

###### Abstract

Many hierarchical reinforcement learning algorithms utilise a series of
independent skills as a basis to solve tasks at a higher level of reasoning.
These algorithms don’t consider the value of using skills that are cooperative
instead of independent. This paper proposes the Cooperative Consecutive
Policies (CCP) method of enabling consecutive agents to cooperatively solve
long time horizon multi-stage tasks. This method is achieved by modifying the
policy of each agent to maximise both the current and next agent’s critic.
Cooperatively maximising critics allows each agent to take actions that are
beneficial for its task as well as subsequent tasks. Using this method in a
multi-room maze domain and a peg in hole manipulation domain, the cooperative
policies were able to outperform a set of naive policies, a single agent
trained across the entire domain, as well as another sequential HRL algorithm.

###### Index Terms:

Reinforcement Learning

##  I Introduction

Many of the struggles that Reinforcement Learning (RL) faces have been
addressed within the field of Heirarchical Reinforcement Learning (HRL) [1].
Complex domains such as ant [2], humanoid [3] and swimmer [4] have all had
high levels of success with hierarchical methods while non-hierarchical
methods have struggled to make progress. Hierarchical methods take advantage
of the ability of many tasks to be abstracted, solving the problem of high
level reasoning separately from low level control [5].

Many HRL approaches involve developing a series of skills (low level policies)
that are used by a high-level controller to solve tasks [6]. However, these
skills cannot always be utilised by a high level controller to optimally solve
a task [7]. For example, for the task of cutting a cake, a high level
controller may employ the skill of grasping a knife, followed by manoeuvring
the knife above the cake, and finally lowering the knife to cut the cake. A
skill trained to grasp a knife may choose to grasp the knife by the blade, a
behaviour appropriate for other tasks, such as passing a knife safely to a
human. This method of grasping however would not allow for effective cake
cutting. The solution to this problem that is proposed in this paper is to use
policies that are aware of the next policy’s goal.

We propose a method of training agents to accommodate for consecutive agents.
This method allows a collection of agents to work together more cohesively to
complete multi-stage tasks. In our proposed approach, each agent is
incentivised to cooperate with the next agent by training the agent’s policy
network to produce actions that maximise both the current agent’s critic and
the next agent’s critic, weighted by an introduced parameter, the cooperative
ratio. By incorporating the next agent’s critic, the current agent can
continue to achieve its own goal while also producing a solution that is
beneficial for the next agent.

The contributions produces by this paper are:

  * •

A novel method, Cooperative Consecutive Policies (CCP), which enables agents
to learn behaviours that maximise reward for their own task while
accommodating for subsequent tasks, improving performance for learning multi-
stage tasks.

  * •

Two case studies in a continuous state/action space maze domain and a robotic
manipulation domain. In these experiments, CCP outperformed a set of naive
policies (trained to greedily maximise subtask reward), a single agent trained
end-to-end on the same task, as well as sequential HRL baseline, using the
transition policies method.

  * •

An ablation study on the effect of varying the cooperative ratio. This
includes a study on the effect of the cooperative ratio on success rate, as
well as the ability to use different cooperative ratios for different
cooperative policies.

Figure 1: A high level view of the Cooperative Consecutive Policies method
works. Each agent acts within its own subtask, and is informed by the next
agent about how to act to assist in completing the subsequent task. In the
maze domain, the next agent can give information about where the dead ends are
and how to avoid them.  Figure 2: The structure of a cooperative policy
implementation in a N subtask environment. Each policy is trained using it’s
own critic and the subsequent agent’s critic, though the last agent is only
trained relative to its own critic. Each critic is trained using data from the
current and previous subtasks, utilising the policy from the agent
corresponding to that subtask.

##  II Related Work

Hierarchical Reinforcement Learning (HRL) has been a long standing field in RL
[8]. HRL methods capitalise on the inherent structure that is present in many
tasks. An important benefit of using HRL is that learning a series of smaller,
simpler skills is easier and faster to learn than learning to solve a single,
more complex task [9] [10] [11] .

Some hierarchical reinforcement learning approaches have been designed to deal
with sequential tasks [12] [13]. Typically they have the form of using a meta-
controller in conjunction with a series of lower-level policies for each
subtask [14] [15] [16] [17], which allows for more generalisation, as the
series of subtasks can be combined in more versatile ways. These methods
follow one of two approaches; manual definition of each subtask, or learning
of the subtasks.

Many HRL algorithms seek to learn skills autonomously. These skills are simple
behaviours that can be utilised by the high level meta-controller [18] [19]
[20] [21]. This approach ensures that the skills are useful for the overall
task, but this tends to be very slow to train.

Manually defining subtasks is an another approach to HRL. This approach
involves decomposing a task into subtasks manually, defining the reward signal
and termination signal for each subtask [18] [22]. Separate agents are then
trained to solve these subtasks independently. The fact that these agents are
not incentivised to assist in solving the overall task means that without
careful engineering of subtasks, suboptimal solutions to tasks are likely
[23]. Our method is designed to overcome the difficulties of engineering
subtasks by enabling the agents to cooperate towards completing the task,
loosening the requirements of careful design. There are other methods that
seek to do a solve a similar problem by learning transition policies between
previously learnt subtasks [24] [25], but these transition policies may still
struggle if the solution to one subtask is too far from an adequate starting
position for the subsequent subtask.

##  III Method

This paper introduces the Cooperative Consecutive Policies (CCP) method. In
the CCP method, a cooperative policy takes into account the subsequent agent’s
critic. The critic informs the policy of how good certain actions and states
are to complete its task effectively. Using both the subsequent and current
critic allows a policy to act to solve the current subtask in a way that
allows an effective solution in future subtasks.

The CCP method is not an algorithm in and of itself. It is an algorithmic
change that can be applied to any actor critic algorithm, assuming there are
consecutive agents used to solve a task. The experiments described in this
paper were implemented using the SAC algorithm [26] [27].

###  III-A Procedure

We consider a modified MDP formulation for solving the problem of picking
optimal actions to solve a task. We assume an environment that involves a task
that is decomposed into a series of $N$ subtasks. At each timestep $t$ an
agent can take an action $a_{t}$ from the current state $s_{t}$, which results
in the environment evolving to the next state $s_{t+1}$, producing a
transition signal $U(s_{t})\in[1,N]$ that determines what subtask is currently
active, and produces a series of $N$ reward signals $r_{n,t}|n\in[1,N]$ that
correspond to each subtask.

The CCP method requires creating $N$ agents, one for each subtask, that each
consist of a policy $\pi_{n}|n\in[1,N]$ with parameters $\theta_{n}$ and a
critic $Q_{n}|n\in[1,N]$ with parameters $\beta_{n}$. This method requires
that the subsequent agent is known to the current agent. Each critic maximises
the discounted sum of future rewards $r_{n}$ from its subtask. Each policy,
rather than maximising the critic that corresponds to their own subtask,
instead maximises a convex sum of the current and subsequent critics. The
standard critic loss signal (Equation 1) and the updated policy loss signal
(Equation 2) are,

$\nabla_{\beta_{n}}=\frac{1}{M}\sum_{0}^{M}(Q_{n}(s,a)-(r_{n}+\gamma
Q_{n}(s^{\prime},\pi_{n}(s^{\prime})))^{2}$ (1)

$\nabla_{\theta_{n}}=\frac{1}{M}\sum^{M}_{0}(-C(Q_{n},Q_{n+1}))$ (2)

Where $M$ is the number of samples in a batch $b$ sampled from replay buffer
$B_{n}$, $\alpha$ is the entropy maximisation term, and $C$ is a convex
combination of the current and subsequent critics, as defined by:

$C(Q_{n},Q_{n+1})=\eta\hat{Q}_{n}(s,\pi_{n}(s))+(1-\eta)\hat{Q}_{n+1}(s,\pi_{n}(s))$
(3)

Where $\eta$ is the cooperative ratio. This ratio affects how much the current
policy acts with respect to the subsequent critic, and is a number between 0
and 1. A cooperative ratio closer to 1 is incentivised to maximise the current
critic’s estimate, whereas a cooperative ratio closer to 0 is incentivised to
maximise the subsequent critic’s estimate. Each Q function is normalised
across the batch as denoted by $\hat{Q}$ using:

$\hat{Q}(s,a)=\frac{Q(s,a)-\displaystyle\min_{s^{\prime}\in
b}Q(s^{\prime},\pi(s^{\prime}))}{\displaystyle\max_{s^{\prime}\in
b}Q(s^{\prime},\pi(s^{\prime}))-\displaystyle\min_{s^{\prime}\in
b}Q(s^{\prime},\pi(s^{\prime}))}$ (4)

This normalisation is done for each critic separately and is recalculated for
each batch. It bounds the critics output to values between 0 and 1. This is
required to compare the current and subsequent critics, which can produce
diverse ranges of values during training. The full algorithmic approach
involves simultaneously gathering data (as shown in Algorithm 1) and training
on that data (as shown in Algorithm 2). Refer to appendix (Section VIII) for
mathematical analysis of this method.

This implementation of CCP is designed using SAC as the base algorithm. To
apply CCP to other algorithms is done by modifying the algorithm’s policy
update to using $C$ in place of the traditional critic evaluation, and
ensuring each agent is updated using the correct buffer $B$, as shown in
Figure 2.

Environment with N subtasks and associated reward signals $r_{(1,...,N)}$;

For each subtask initialise an agent $A_{n}$, including a policy $\pi_{n}$
with parameters $\theta_{n}$, a critic $Q_{n}$ with parameters $\beta_{n}$,
and a replay buffer $B_{n}$;

while  _timestep $<$ maxTimestep_ do

$s,n\xleftarrow[]{}$ reset environment;

while  _not done_ do

$a$ $\sim\pi_{n}(s)$;

$s^{\prime},r_{(1,...,N)},done\xleftarrow[]{}$ environment step with $a$;

record $(s,a,r_{1,...,N},s^{\prime},done)$ in $B_{n}$;

$s\xleftarrow[]{}s^{\prime}$;

$n\xleftarrow[]{}U(s^{\prime})$;

end while

end while

Algorithm 1 Gathering Data

Set of N agents $A_{n}$, each with policy $\pi_{n}$ with parameters
$\theta_{n}$, a critic $Q_{n}$ with parameters $\beta_{n}$, and a replay
buffer $B_{n}$;

cooperative ratio $\eta$;

discount factor $\gamma$;

entropy maximisation term $\alpha$;

for  _n in (1,…,N)_ do

sample minibatch $b$ from $B_{n}$ of M samples
$\xrightarrow{}(s,a,r_{(1,...,N)},s^{\prime},d)$;

for  _j in (n, n+1)_ do

$y(r_{j},s^{\prime},d)=r_{j}+(1-d)\gamma(Q_{j}(s^{\prime},a^{\prime})-\alpha$
log($\pi_{n}(a^{\prime}|s^{\prime})))$, $a^{\prime}\sim\pi_{n}(s^{\prime})$;

$\nabla_{\beta_{j}}=\frac{1}{M}\sum(Q_{j}(s,a)-y(r_{j},s^{\prime},d))^{2}$;

end for

$a^{\prime}\sim\pi_{n}(s)$;

for  _j in (n, n+1)_ do

Calculate $\hat{Q}_{n}$ across minibatch $b$ according to Equation 4;

end for

$C(Q_{n},Q_{n+1})=\eta\hat{Q}_{n}(s,a^{\prime})+(1-\eta)\hat{Q}_{n+1}(s,a^{\prime})$;

$\nabla_{\theta_{n}}$ =
$\frac{1}{M}\sum(\alpha\log\pi_{n}(a^{\prime}|s)-C(Q_{n},Q_{n+1}))$;

end for

Algorithm 2 Cooperative Training

##  IV Experimental studies

Figure 3: A 3 room example of the Maze domain. The first agent begins in the
starting area and produces actions to navigate the first room. As the agent
enters the next room, the next takes charge to navigate that room. The agent’s
scan range is shown in purple.

###  IV-A Environments

Two domains were used to test the efficacy of the CCP method: the maze domain,
and the peg in hole domain.

####  IV-A1 Maze Navigation

In the maze domain, a series of consecutive rooms were created. Each room had
two paths to exit into the next room, one of which leads to a dead end. Each
room was considered a subtask with the goal being to exit the room, with a
reward signal dependent on horizontal position, starting at 0 on the left of
the room and linearly increasing to 1 at the right of the room. The overall
task in this domain was to travel through all the rooms to get from one side
to the other. The rooms are designed such that an optimal subtask solution
would not lead to an optimal overall solution, as the reward signal for each
task increases as the position moves to the right, not towards the correct
path. This domain allows for easy extension in terms of the number of subtasks
that are required to be able to solve the overall task.

The agent in this domain uses two continuous actions: linear velocity and
angular velocity. The agent uses an observation of the environment that
includes a laser scan as well as a global position in the maze, all of which
are continuous measurements. This domain is shown in Figure 3.

####  IV-A2 Peg in hole manipulation

The second domain used was the peg in hole domain . In this domain, the task
is to use a Panda arm to insert a peg into a hole. This domain is decomposed
into two subtasks: grasping and inserting. The grasping subtask requires the
gripper to establish a grasp on the peg. This subtask uses an exponentially
decaying reward based on distance between the centre of the peg and the
fingers of the gripper. The second subtask is to move the peg towards and then
into the hole. The reward for this subtask is also an exponentially decaying
reward, this time based on the distance between the centre of the peg and the
centre of the hole. The task is considered a success if the peg is fully
inserted into the hole.

This domain was chosen as it represents a real world example of how completing
one subtask can affect the completion of a second subtask. If the peg is
grasped by the thin section, then the insertion subtask cannot be completed
optimally, as the peg cannot be fully inserted into the hole. This domain is
more complex than the maze domain, due to the higher dimensionality of the
state space and the complexity of grasping dynamics

The agent in this domain has two continuous actions: end effector velocity
along the axis of actuation and gripper force. The state that the agent
observes is the positions, orientations and velocities of the peg, hole and
end effector.

Both domains in this paper use relatively simple reward signals. Though more
complex, and potentially more informative, reward signals could be
constructed, these simple reward signals were used on purpose to show that
using a imperfect reward signal can be overcome using the CCP method.
Engineering a reward signal that can effectively avoid globally suboptimal
behaviours is an expensive and sometimes intractable problem, and being able
to solve a task without it is a valuable quality.

###  IV-B Algorithms

Four different methods were evaluated in these domains:

  * •

CSAC: The CSAC method is an implementation of Soft Actor Critic (SAC) [26][27]
that utilises the CCP methodology of incentivising cooperative solving of
sequential tasks, making a set of Cooperative SAC agents.

  * •

Naive SAC: The naive method uses a separate SAC agent for each subtask,
without any communication between agents. Each agent is attempting to maximise
solely its own reward signal. This method represents the naive approach of
treating each agent independently.

  * •

SAC: The SAC method involves using a single end-to-end SAC agent trained to
perform across the whole domain, utilising a reward signal that is the
combined reward signal from all the subtasks. This method represents the
standard RL approach to solving a task.

  * •

TP: The Transition Policies method [24] used as a baseline. This method trains
primitive policies to complete each subtask, and then trains transition
policies to move from subtask termination states to states that are good
initialisations for the subsequent subtask.

Hyperparameters and further implementation details are listed in Appendix 2
(Section VII).

###  IV-C Studies

####  IV-C1 Optimal cooperative ratio

The first study tested the performance of the four different methods in both
the maze domain and the peg in hole domain. The CCP method requires the
selection of a cooperative ratio to determine the trade-off between current
and subsequent critics. In the maze domain three different length mazes were
tested using a sweep across the cooperative ratio parameter to determine its
effects on learning performance, with a similar sweep conducted in the peg in
hole domain. The best results from these sweeps were then compared to the
results from using the other methods. This study also investigated the
sensitivity of the method’s performance with respect to the cooperative ratio.

####  IV-C2 Independent cooperative ratios

The previous experiment used the same cooperative ratio for each agent within
each experiment. An investigation was conducted to determine whether using
independent cooperative ratios for each agent would improve the performance,
conducted in the 3 room maze domain, in which there are two cooperative
agents. A parameter sweep was conducted across the cooperative ratio for each
cooperative agent to investigate it’s effect on learning performance.

##  V Results

Figure 4: Success rate in the 2, 3 and 4 room maze domains. Each plot shows the average success rate across several different random seeds, where success is defined as reaching the end of the maze. The domains use a cooperative ratio of 0.1, 0.1 and 0.9 respectively. Figure 5: Success rate in the peg in hole domain. The plot shows the average success rate across several different random seeds. This uses a cooperative ratio of 0.9 TABLE I: Success rate of each method in maze and peg in hole domains Method |  Cooperative Ratio | Maze - Subtasks | Peg in hole
---|---|---|---
2 | 3 | 4
CSAC | 0.1 | 1.00 | 0.62 | 0.00 | 0.11
0.2 | 1.00 | 0.23 | 0.00 | 0.11
0.3 | 1.00 | 0.07 | 0.00 | 0.11
0.4 | 0.31 | 0.00 | 0.00 | 0.17
0.5 | 0.05 | 0.06 | 0.00 | 0.21
0.6 | 0.17 | 0.04 | 0.00 | 0.34
0.7 | 0.42 | 0.40 | 0.17 | 0.31
0.8 | 0.51 | 0.16 | 0.50 | 0.37
0.9 | 0.06 | 0.21 | 0.79 | 0.40
SAC | 1.00 | 0.00 | 0.00 | 0.11
Naive | 0.54 | 0.20 | 0.14 | 0.13
TP | 0.43 | 0.50 | 0.40 | -

###  V-A Optimal cooperative ratio

####  V-A1 Maze navigation

The results for the 2 room experiment within the maze domain are presented in
Figure 4 and are summarised in Table I. This experiment shows that the
cooperative and naive agents both learn a successful policy in a similar time
period, whereas the single agent policy took more than 3 times as long to
reach a similar performance. The naive agents, though they reached a high
level of success quickly, had a decaying performance. This is due to the fact
that each agent learned a solution to their domain quickly, which included
travelling to the further door. As each agent refined its solution, the
shorter path, which is suboptimal overall, was used more frequently, reducing
the performance of the overall task. This experiment shows that decomposing a
task into subtasks is beneficial in terms of training speed, shown by the
relative training speed of the cooperative and naive policies compared to the
single agent. This experiment also shows that just decomposing a task into
subtasks and then treating them as entirely separate problems can lead to
suboptimal or decaying solutions.

Figures 4 and Table I show the success rates of the three different agent
types in the 3 and 4 room mazes, where success is defined as reaching the end
of the maze. In both of these domains, the single agent was not able to learn
a successful policy at all within the training window. This shows that the
local minima, the dead end paths, proved too difficult for the end-to-end SAC
agent to overcome within the time frame. The cooperative agents were able to
achieve a consistently higher performance than the naive agents. This
demonstrates the value of splitting a task into subtasks and then solving them
cooperatively.

Across all room configurations the TP baseline had a suboptimal level of
performance. This method learns a set of transition policies that attempts to
manipulate the agent from the termination state of a subtask to a good
starting state for the subsequent subtask. Due to the way the doorways are
arranged within the maze, a subtask that concludes at a dead end is too far
from a good starting state for the next subtask for a transition policy to be
able to rectify.

####  V-A2 Peg in hole manipulation

Figure 5 shows the results of using each of the 3 methods in the peg in hole
domain. CSAC is able reach a far higher level of performance in this domain
compared to both the single agent SAC method and the naive method. The low
success rate in this domain is due to the difficulty of using contact physics.
If the agent grasps the peg in the wrong way or applies too much force, the
peg can be pushed into a pose that is unreachable for the robot, making the
episode unsolvable. Despite this, CSAC is able to complete the task 40% of the
time, far greater than the other two methods which only achieved a success
rate of 11% and 13% for single SAC and naive SAC respectively.

The SAC method has a low performance in this domain due to the fact that it is
encouraged to move the peg to the hole as quickly as possible. This requires
that the peg be grasped by the narrow section and dragged to the hole. Due to
this poor grasp location, the task cannot be fully solved. In the case of the
naive SAC, the reward signal that the first agent is getting is to move close
to the peg, not necessarily to grasp it. This means that the first agent will
hover close to the peg to maximise it’s reward. The CCP method avoids both
these behavioural pitfalls and completes the task with a relatively high level
of success.

Table I shows the different success rates in each domain using a range of
cooperative ratios. These results show that the cooperative ratio is a
sensitive variable, and it’s optimal value is task dependent. This shows that
is important to use a sweep across the variable to ensure a good performance.

###  V-B Independent cooperative ratios

The results of the individual cooperative ratio experiments are shown in
Figure 6. This figure shows the average success rate of a set of cooperative
policies in a 3 room domain. In this experiment the first and second policies
used independent cooperative ratios, represented along the bottom two axes of
Figure 6. Using the individually tuned cooperative ratios, a higher level of
performance was found in the 3 room domain compared to using the same
cooperative ratio for both policies. Using cooperative ratios of 0.1 and 1.0
for the first two policies respectively result in the highest level of success
(78%), outperforming the policies trained using a shared cooperative ratio
(54%).

Figure 6: The measured success rate in the 3 room maze domain with different
cooperative ratios for each cooperative policy. The success rate is an average
success rate across the last 10 epochs of each of the 10 different randomly
initialised iterations of each configuration. Blue represents a low success
rate and red represents a high success rate.

Some insights can be gathered from analysing the effect on success rate when
changing the cooperative ratio by looking at Table I and Figure 6. In all of
the shared ratio experiments, using a cooperative ratio around 0.5 produces a
low level of performance. In the maze domain, the optimal behaviour for one
room compared to the next room involves utilising different paths. Trying to
achieve both objectives equally appears to lead to a poor performance overall,
as the policy can’t decide which objective to prioritise. This trend can be
seen in the individual cooperative ratio experiments when comparing the
cooperative ratios for the first policy. Using a cooperative ratio of 0.5 for
either policy leads to lower performance as seen in the shared ratio
experiments.

The independent ratio experiments also show that using a low cooperative ratio
for the first policy has a higher success rate than using a high cooperative
ratio, similar to results in the 2 room results in the shared ratio
experiments. This likely shows that in the first subtask it is harder to avoid
the dead-end without increased direction from the subsequent critic.

##  VI Conclusion

This paper introduces the CCP method for cooperatively solving multi-stage
tasks. This method was tested using the SAC algorithm (implementation called
CSAC) in two different domains, the maze domain and the peg in hole domain,
and was compared against three other methods; a SAC agent trained end-to-end
across the whole domain, a set of naive agents trained to solve each subtask
greedily, and a baseline HRL algorithm for sequential tasks, the Transition
Policies algorithm.

The CCP method outperformed each other method in the maze domain, as
summarised in Table I and shown in Figure 4. In the simplest domain (2 room
maze), CSAC converged on a solution 4 times faster than the single agent and
was able to maintain a high level of performance that the naive policies were
not able to maintain, while the TP baseline was unable to solve the domain. In
the more complex domains (3 and 4 room mazes), the cooperative policies had a
consistently higher level of performance than the naive policies and TP
baseline, whereas the single agent was not able to find any solution to the
task within 3 million training steps. Similar results were found in the peg in
hole domain (Figure 5), in which the algorithm using CSAC had a success rate
approximately 30% higher than the other methods.

Currently CCP requires that the task be decomposed into subtasks before
training. There are works that seek to learn the decomposition of a task [18],
and these methods could be combined with CCP for a more generalisable
algorithm. Additionally, the cooperative ratio variable is required to be
defined for each subtask. A proposed solution to this is to utilise meta-
learning [10] to tune this parameter. This method also only provides direction
based on the subsequent subtask’s critic. If a specific action needs to be
taken to solve the overall task and the subsequent agent isn’t aware of it,
then CCP will still fail. Future works will seek to address these issues and
make this method more robust.

## References

  * [1] O. Nachum, H. Tang, X. Lu, S. Gu, H. Lee, and S. Levine, “Why does hierarchy (sometimes) work so well in reinforcement learning?” _arXiv preprint arXiv:1909.10618_ , 2019.
  * [2] O. Nachum, S. Gu, H. Lee, and S. Levine, “Near-optimal representation learning for hierarchical reinforcement learning,” _arXiv preprint arXiv:1810.01257_ , 2018.
  * [3] X. B. Peng, M. Chang, G. Zhang, P. Abbeel, and S. Levine, “Mcp: Learning composable hierarchical control with multiplicative compositional policies,” _arXiv preprint arXiv:1905.09808_ , 2019.
  * [4] S. Li, R. Wang, M. Tang, and C. Zhang, “Hierarchical reinforcement learning with advantage-based auxiliary rewards,” _arXiv preprint arXiv:1910.04450_ , 2019.
  * [5] M. Wulfmeier, D. Rao, R. Hafner, T. Lampe, A. Abdolmaleki, T. Hertweck, M. Neunert, D. Tirumala, N. Siegel, N. Heess _et al._ , “Data-efficient hindsight off-policy option learning,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 11 340–11 350.
  * [6] A. S. Vezhnevets, S. Osindero, T. Schaul, N. Heess, M. Jaderberg, D. Silver, and K. Kavukcuoglu, “Feudal networks for hierarchical reinforcement learning,” in _Proceedings of the 34th International Conference on Machine Learning-Volume 70_. JMLR. org, 2017, pp. 3540–3549.
  * [7] O. Nachum, S. S. Gu, H. Lee, and S. Levine, “Data-efficient hierarchical reinforcement learning,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 3303–3313.
  * [8] M. Al-Emran, “Hierarchical reinforcement learning: a survey,” _International journal of computing and digital systems_ , vol. 4, no. 02, 2015.
  * [9] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [10] K. Frans, J. Ho, X. Chen, P. Abbeel, and J. Schulman, “Meta learning shared hierarchies,” _arXiv preprint arXiv:1710.09767_ , 2017.
  * [11] S. Iqbal and F. Sha, “Actor-attention-critic for multi-agent reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 2961–2970.
  * [12] J. Oh, S. Singh, H. Lee, and P. Kohli, “Zero-shot task generalization with multi-task deep reinforcement learning,” in _Proceedings of the 34th International Conference on Machine Learning-Volume 70_. JMLR. org, 2017, pp. 2661–2670.
  * [13] R. T. Icarte, T. Klassen, R. Valenzano, and S. McIlraith, “Using reward machines for high-level task specification and decomposition in reinforcement learning,” in _International Conference on Machine Learning_ , 2018, pp. 2107–2116.
  * [14] P.-L. Bacon, J. Harb, and D. Precup, “The option-critic architecture,” in _Thirty-First AAAI Conference on Artificial Intelligence_ , 2017.
  * [15] A. Azarafrooz and J. Brock, “Hierarchical soft actor-critic: Adversarial exploration via mutual information optimization,” _arXiv preprint arXiv:1906.07122_ , 2019.
  * [16] A. Gupta, V. Kumar, C. Lynch, S. Levine, and K. Hausman, “Relay policy learning: Solving long-horizon tasks via imitation and reinforcement learning,” _arXiv preprint arXiv:1910.11956_ , 2019.
  * [17] C. Tessler, S. Givony, T. Zahavy, D. Mankowitz, and S. Mannor, “A deep hierarchical approach to lifelong learning in minecraft,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 31, no. 1, 2017.
  * [18] J. Andreas, D. Klein, and S. Levine, “Modular multitask reinforcement learning with policy sketches,” in _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ser. ICML’17. JMLR.org, 2017, pp. 166–175.
  * [19] T. Haarnoja, K. Hartikainen, P. Abbeel, and S. Levine, “Latent space policies for hierarchical reinforcement learning,” _arXiv preprint arXiv:1804.02808_ , 2018.
  * [20] D. Esteban, L. Rozo, and D. G. Caldwell, “Hierarchical reinforcement learning for concurrent discovery of compound and composable policies,” _arXiv preprint arXiv:1905.09668_ , 2019.
  * [21] A. Sharma, S. Gu, S. Levine, V. Kumar, and K. Hausman, “Dynamics-aware unsupervised discovery of skills,” _arXiv preprint arXiv:1907.01657_ , 2019\.
  * [22] S. Forestier, Y. Mollard, and P.-Y. Oudeyer, “Intrinsically motivated goal exploration processes with automatic curriculum learning,” _arXiv preprint arXiv:1708.02190_ , 2017.
  * [23] C. Florensa, Y. Duan, and P. Abbeel, “Stochastic neural networks for hierarchical reinforcement learning,” _arXiv preprint arXiv:1704.03012_ , 2017.
  * [24] Y. Lee, S.-H. Sun, S. Somasundaram, E. S. Hu, and J. J. Lim, “Composing complex skills by learning transition policies,” in _International Conference on Learning Representations_ , 2018.
  * [25] Y. Lee, J. J. Lim, A. Anandkumar, and Y. Zhu, “Adversarial skill chaining for long-horizon robot manipulation via terminal state regularization,” _arXiv preprint arXiv:2111.07999_ , 2021.
  * [26] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” _arXiv preprint arXiv:1801.01290_ , 2018.
  * [27] T. Haarnoja, A. Zhou, K. Hartikainen, G. Tucker, S. Ha, J. Tan, V. Kumar, H. Zhu, A. Gupta, P. Abbeel _et al._ , “Soft actor-critic algorithms and applications,” _arXiv preprint arXiv:1812.05905_ , 2018.
  * [28] Vitchyr, “vitchyr/rlkit.” [Online]. Available: https://github.com/vitchyr/rlkit
  * [29] J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz, “Trust region policy optimization,” in _International conference on machine learning_. PMLR, 2015, pp. 1889–1897.

##  VII Appendix 1

The implementation used in this paper was based on the RLkit implementation of
SAC [28]. This implementation utilises an epoch based approach. All methods
use the hyperparameters shown in Table II. These hyperparameters were adapted
from the RLkit SAC implementation [28]. For each experiment, 5 random seeds
were used per method and domain, with 3 used for the transition policies
experiments.

It is recommended that if this algorithm is recreated, the batch size is at
least as large as in this paper, to ensure that the convex combination of
critics is effective across each batch. It is also recommended that the
discount factor $\gamma$ is not increased above 0.95, to ensure the method
doesn’t become unstable when estimating too far into the future.

The transition policies baseline was trained using the TRPO[29] algorithm
within the implementation provide with its publication [24]. The plots in this
paper that show the performance of the TP algorithm only represent the
training of the transition policies themselves. The primitive policies that
learned to solve each subtask were trained separately, and then each run that
the TP algorithm underwent used this same set of primitive policies.

TABLE II: Hyperparameters used in Maze Domain Hyperparameter | Value
---|---
Gamma ($\gamma$) | 0.95
Replay buffer length | 1e6
Batch size | 256
Maximum episode length | 1000
Soft target update factor | 0.005
Timesteps per epoch | 5000
Training loops per epoch | 1000

##  VIII Appendix 2

This section proves that using a convex combination of two critic functions
can be used for the purposes of training a policy to solve a task.

###  VIII-A Preliminaries

Consider a task that has an action space $A$, a state space $S$, and a reward
signal $R$. This task can be decomposed into a series of $N$ subtasks such
that each subtask has an action space $A_{n}\subseteq A$, a state space
$S_{n}\subseteq S$, and a reward signal $R_{n}\in[0,1]$. Each subtask reward
signal is defined by,

$R_{n}(s,a)=\begin{cases}\in[0,1],&\text{if }s\in S_{n}\\\ 0,&\text{else if
}s\in S_{m}|m\in[1,n)\\\ 1,&\text{else if }s\in S_{m}|m\in(n,N]\end{cases}$
(5)

The reward signal $R$ for the overall task is defined as,

$R(s,a)=\sum_{n=1}^{N}R_{n}(s,a)$ (6)

which is equivalent to,

$R(s,a)=n-1+R_{n}(s,a)\quad|\quad s\in S_{n},\forall a\in A$ (7)

The environment also contains a subtask transition function $U(s)\in[1,N)$
that denotes the current subtask. This function will only transition from
subtask $n$ to $n+1$ when $R_{n}\geq 1-e$, where $e$ is some sensitivity
metric.

A Q-function $Q$ approximates the expected discounted sum of future rewards,

$\displaystyle Q(s_{t},a_{t})$ $\displaystyle=R(s_{t},a_{t})+\gamma
Q(s_{t+1},\pi(s_{t+1}))$ (8)
$\displaystyle=\sum_{i=t}^{\infty}\gamma^{i-t}R(s_{i},a_{i}),\quad
a_{i}\sim\pi(s_{i})$ (9)

where $\pi$ is a policy chosen to maximise the Q function and the environment
transitions deterministically.

###  VIII-B Theorem

A convex combination of two normalised Q functions can be used in the place of
the sum of two Q functions for the purposes of picking an optimal action. The
original maximisation objective can be recovered using the correct weighting.

A standard policy $\pi$ selects actions that maximise a Q function such that,

$Q(s,\pi(s))=\max_{a\in A}Q(s,a),$ (10)

This can be decomposed into two Q functions according to Lemma 1.

####  VIII-B1 Lemma 1

A Q function for a task involving multiple sequential subtasks can be be
approximated as the sum of two Q functions, each corresponding to the current
and subsequent subtasks

Assume that for each subtask there exists a Q function $Q_{n}$. Assume that
the cooperative reward signal $r^{coop}_{n}$ for the current subtask $n$ is
given by

$r^{coop}_{n}=\eta r_{n}+(1-\eta)r_{n+1}$ (11)

where $r_{n}$ is the reward signal from the environment for subtask $n$, and
$\eta\in[0,1]$ is a weighting parameter.

Using this assumed reward signal, the cooperative Q function $C_{n}$ for a
subtask $n$ can be found.

$\displaystyle C_{n}^{\pi}(s,a)$
$\displaystyle=\mathbb{E}\left[\sum^{\infty}_{t=0}\gamma^{t}r^{coop}_{n}|\pi,s,a\right]$
(12)
$\displaystyle\begin{split}&=\eta\mathbb{E}\left[\sum^{\infty}_{t=0}\gamma^{t}r_{n}|\pi,s,a\right]+...\\\
&\qquad(1-\eta)\mathbb{E}\left[\sum^{\infty}_{t=0}\gamma^{t}r_{n+1}|\pi,s,a\right]\\\
\end{split}$ (13) $\displaystyle=\eta
Q_{n}^{\pi}(s,a)+(1-\eta)Q^{\pi}_{n+1}(s,a)\blacksquare$ (14)

This formulation of a cooperative Q function uses the weighted sum of two Q
functions for consecutive subtasks. The subsequent subtask will always exhibit
a smaller Q value, due to the rewards of the subtask being further away
through time. The above formulation for a cooperative Q function uses the
weighting variable $\eta$ to represent both this natural weighting of time
discounting of rewards, as well as the weighting that represents balancing
between the current and subsequent subtasks. To remove this natural time
discounting, a normalised cooperative Q function $\hat{C}_{n}$ is introduced.

$\hat{C}_{n}=\eta\hat{Q}_{n}(s,a)+(1-\eta)\hat{Q}_{n+1}(s,a),$ (15)

where $\hat{Q}$ is a normalisation such that,

$\hat{Q}(s,a)=\frac{Q(s,a)-\displaystyle\min_{a\in
A}Q(s,a)}{\displaystyle\max_{a\in A}Q(s,a)-\displaystyle\min_{a\in A}Q(s,a)}$
(16)

To simplify this operation, the following functions are defined,

$Ran(Q)=\max_{a\in A}Q(s,a)-\min_{a\in A}Q(s,a)$ (17) $Min(Q)=\min_{a\in
A}Q(s,a)$ (18)

Using these it can be shown that using the correct $\eta$ can lead to the same
action.

$\displaystyle\max_{a\in A}(\hat{C}_{n})=\max_{a\in
A}(\eta\hat{Q}_{n}+(1-\eta)\hat{Q}_{n+1})$ (19) $\displaystyle=\max_{a\in
A}\left(\eta\frac{Q_{n}-Min(Q_{n})}{Ran(Q_{n})}+(1-\eta)\frac{Q_{n+1}-Min(Q_{n+1})}{Ran(Q_{n+1})}\right)$
(20) $\displaystyle=\max_{a\in
A}\left(\frac{\eta}{Ran(Q_{n})}Q_{n}+\frac{(1-\eta)}{Ran(Q_{n+1})}Q_{n+1}\right)$
(21)

The evaluations $\max(\hat{C}_{n})$ and $\max(Q)$ produce the same actions
when the coefficients of the convex operation are equivalent. Therefore,

$\displaystyle\frac{\eta}{Ran(Q_{n})}$
$\displaystyle=\frac{(1-\eta)}{Ran(Q_{n+1})}$ (22) $\displaystyle\eta
Ran(Q_{n+1})$ $\displaystyle=(1-\eta)Ran(Q_{n})$ (23)
$\displaystyle\eta(Ran(Q_{n+1})+Ran(Q_{n}))$ $\displaystyle=Ran(Q_{n})$ (24)
$\displaystyle\eta$
$\displaystyle=\frac{Ran(Q_{n})}{(Ran(Q_{n+1})+Ran(Q_{n}))}$ (25)

This definition of $\eta$ is also bounded between 0 and 1. Using this $\eta$,
the original maximisation action can be recovered. $\blacksquare$


# Collisionless Magnetorotational Turbulence in Pair Plasmas:
Steady-state Dynamics, Particle Acceleration, and Radiative Cooling

Fabio Bacchini<EMAIL_ADDRESS>Centre for mathematical Plasma
Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B,
B-3001 Leuven, Belgium  Royal Belgian Institute for Space Aeronomy, Solar-
Terrestrial Centre of Excellence, Ringlaan 3, 1180 Uccle, Belgium  Vladimir
Zhdankin  Department of Physics, University of Wisconsin-Madison, Madison, WI
53706, USA  Center for Computational Astrophysics, Flatiron Institute, 162
Fifth Avenue, New York, NY 10010, USA  Evgeny A. Gorbunov  Centre for
mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven,
Celestijnenlaan 200B, B-3001 Leuven, Belgium  Gregory R. Werner  Center for
Integrated Plasma Studies, Department of Physics, University of Colorado, 390
UCB, Boulder, CO 80309-0390, USA  Lev Arzamasskiy  School of Natural Sciences,
Institute for Advanced Study, Princeton, NJ 08544, USA  Mitchell C. Begelman
JILA, University of Colorado and National Institute of Standards and
Technology, 440 UCB, Boulder, CO 80309-0440, USA  Department of Astrophysical
and Planetary Sciences, University of Colorado, 391 UCB, Boulder, CO
80309-0391, USA  Dmitri A. Uzdensky  Center for Integrated Plasma Studies,
Department of Physics, University of Colorado, 390 UCB, Boulder, CO
80309-0390, USA

###### Abstract

We present 3D fully kinetic shearing-box (SB) simulations of pair-plasma
magnetorotational turbulence with unprecedented macro-to-microscopic scale
separation. We retrieve the expected fluid-model behavior of turbulent
magnetic fields and angular-momentum transport, and observe fundamental
differences in turbulent fluctuation spectra linked with plasma heating. For
the first time, we provide a definitive demonstration of nonthermal particle
acceleration in kinetic magnetorotational turbulence agnostically of SB
initial conditions, by means of a novel strategy exploiting synchrotron
cooling.

_Introduction._ — Recent years have witnessed a surge of interest in the
plasma and astrophysics communities toward high-energy environments such as
the surroundings of supermassive black holes (SMBHs). Groundbreaking
observations with the Event Horizon Telescope [1, 2] and Gravity [3] hinted at
the existence of (quasi)stationary plasma structures (e.g. disks) around
SMBHs, underlining the importance of modeling plasma flows in strong-gravity
environments. Most efforts in this direction have relied on
magnetohydrodynamic (MHD) simulations [4, 5, 6, 7, 8, 9, 10], often resorting
to relatively simple prescriptions to include limited microscopic (i.e.
kinetic) effects [11, 12, 13, 14, 15]. However, it is commonly accepted that
plasmas around low-luminosity (LL) SMBHs are collisionless, and that kinetic
plasma physics can shape the global dynamics, thereby invalidating basic MHD
assumptions. For example, accretion onto LL SMBHs is subject to enhanced
angular momentum transport (AMT) by beyond-MHD collisionless turbulence.
Moreover, electron-ion plasmas in disks likely exist in two-temperature states
where relativistic electrons produce most of the radiation (e.g. via
synchrotron cooling, SC). Finally, relativistic electron energies could result
from nonthermal particle acceleration (NTPA) via turbulent dissipation of
magnetic fields. MHD models cannot completely capture these effects, rendering
fully kinetic simulations of accreting plasmas necessary.

Due to computational constraints, only local Particle-in-Cell (PIC) kinetic
simulations of turbulent accretion (e.g. considering individual disk-plasma
patches) are possible to date [16, 17, 18, 19, 20] (with recent, important
progress for large-scale/global modeling [21, 22, 23, 24, 25]). The shearing-
box (SB) paradigm is a reliable tool for such local simulations [26, 27, 28,
29, 30, 31, 32] focusing on the magnetorotational instability (MRI) [33, 34,
35], whose nonlinear development drives turbulence. However, the handful of
existing PIC-SB MRI simulations suffer from severe limitations [16, 17, 18,
19, 20]. Of particular importance is scale ordering: for realism, macroscopic
(MRI–MHD) scales should be well separated from kinetic plasma scales, but PIC
simulations often employ rescaled parameters to limit computing costs. This
can dramatically alter the MRI development and produce large deviations in
overall dynamics between MHD and PIC-SB runs. Recently, we demonstrated that
“large-enough” 3D PIC simulations retrieve the expected MHD-like dynamics
[20], attaining a saturated, nonlinear turbulent state. However, it is still
essential to explore larger scale separations to definitively determine the
MHD-to-kinetic transition of MRI turbulence in magnetic and bulk-velocity
power spectra. Limited scale separation in earlier 3D works also caused the
unrealistically slow growth of kinetic (especially mirror) instabilities
regulating pressure anisotropy, thereby affecting overall AMT.

Importantly, previous PIC-SB MRI simulations [16, 17, 18, 19, 20] did not
distinguish NTPA in the initial transient stages (including linear MRI growth)
that precede fully developed MRI turbulence. These “preturbulent” stages
develop macroscopic (system-scale) current sheets that undergo magnetic
reconnection, releasing large amounts of energy and producing considerable
NTPA in PIC simulations; however, these early stages are artificial (i.e. only
related to the development of the MRI in the SB paradigm) and hence irrelevant
for realistic accretion disks, and should therefore be excluded. Indeed, we
assert that the capability of MRI turbulence to accelerate particles _has not
been assessed at all thus far_ in a manner that clearly distinguishes
turbulent energization from earlier artificial injection mechanisms.

In this Letter, we advance the frontier of collisionless-MRI studies by
performing the largest (to date) 3D PIC-SB simulations. We employ
unprecedented scale separations to obtain PIC results retrieving expected MHD
dynamics, and allowing us to study detailed microphysical effects potentially
affecting the global accretion flow. Moreover, by imposing radiative cooling
to suppress NTPA during (only) the preturbulence stages, we distinctly
characterize NTPA solely due to MRI turbulence. We also show that cooling has
negligible effect on the MRI evolution aside from extinguishing the high-
energy nonthermal particles; this is the first exploration of radiation-
reaction (RR) effects in PIC-SB simulations.

_Numerical Model and Setup._ — We conduct 3D PIC-SB simulations with Zeltron
[36], employing the kinetic SB with orbital advection (KSB-OA) [20] to model a
local Cartesian ($x\equiv\,$radial direction from the central SMBH,
$y\equiv\,$toroidal, $z\equiv\,$vertical) sector of an accretion disk111Since
the publication of the original KSB-OA method [20], we have improved the SB
implementation in Zeltron with a more stable and less diffusive treatment of
the shearing-periodic boundary conditions.. The local plasma patch represented
by the simulation box has orbital frequency $\Omega_{0}$ and is subjected to a
background force corresponding to a (linearized) Keplerian shearing-velocity
profile
$\mbox{\boldmath{$v$}}_{\mathrm{s}}(x)=-(3/2)\Omega_{0}x\mbox{\boldmath{$\hat{y}$}}$.
We apply shearing-periodic boundary conditions at $x$-boundaries, while $y$\-
and $z$-boundaries are purely periodic [20]. Our initial setup consists of a
thermal pair plasma for simplicity ($m=m_{i}=m_{e}$), with uniform per-species
number density $n_{0}=n_{i,0}=n_{e,0}$, nonrelativistic uniform (normalized)
temperature $\theta_{0}\equiv k_{\mathrm{B}}T_{0}/(mc^{2})\ll 1$ (with
$T_{0}=T_{i,0}=T_{e,0}$), and zero (fluctuating) bulk speed. A weak, purely
vertical magnetic field $\mbox{\boldmath{$B$}}_{0}=(0,0,B_{0})$ initially
threads the plasma. The free parameters of our model are therefore the
physical box size $L_{x}\times L_{y}\times L_{z}$, the initial temperature
$\theta_{0}$, the initial Alfvén-to-light speed ratio (including both particle
species) $v_{\mathrm{A},0}/c\equiv B_{0}/\sqrt{8\pi n_{0}mc^{2}}$, and the
initial cyclotron-to-rotational frequency ratio
$\omega_{\mathrm{c},0}/\Omega_{0}$ where $\omega_{\mathrm{c},0}\equiv
eB_{0}/(mc)$. On the basis of previous work [20], we choose
$L_{z}=2\lambda_{\mathrm{MRI}}$ (with $\lambda_{\mathrm{MRI}}\equiv 2\pi
v_{\mathrm{A},0}/\Omega_{0}$ the most-unstable MRI wavelength),
$\theta_{0}=1/128$, $v_{\mathrm{A},0}/c\simeq 0.007$, and
$\omega_{\mathrm{c},0}/\Omega_{0}=30$ to obtain the nonrelativistic length-
scale ordering
$c/\omega_{\mathrm{p}}<\rho_{\mathrm{c},0}<\lambda_{\mathrm{MRI}}<L_{z}\leq
H_{0}$; here, $\omega_{\mathrm{p}}\equiv\sqrt{8\pi e^{2}n_{0}/m}$ is the total
plasma frequency, $\rho_{\mathrm{c},0}\equiv\sqrt{\theta_{0}}mc^{2}/(eB_{0})$
is the (initial) typical particle gyroradius, and
$H_{0}\equiv\sqrt{\theta_{0}}c/\Omega_{0}$ is the disk scale height. These
parameters thus determine a separation
$\rho_{\mathrm{c},0}/(c/\omega_{\mathrm{p}})\simeq 12.5$,
$\lambda_{\mathrm{MRI}}/\rho_{\mathrm{c},0}\simeq 15$, and
$H_{0}/\lambda_{\mathrm{MRI}}=L_{z}/\lambda_{\mathrm{MRI}}=2$ between the
macroscopic and microscopic length scales involved. Our choice of parameters
results in an initial plasma $\beta_{0}\equiv
2n_{0}k_{\mathrm{B}}T_{0}/[B_{0}^{2}/(8\pi)]=2\theta_{0}/(v_{\mathrm{A},0}/c)^{2}=2(\rho_{\mathrm{c},0}/(c/\omega_{\mathrm{p}}))^{2}\simeq
312$. Additionally, we choose $L_{x}=2L_{z}$ such that the shearing-velocity
offset between $x$-boundaries ($\pm L_{x}/2$) is $(3/2)\Omega_{0}L_{x}<c$, and
we set $L_{y}/L_{z}=4$. We employ 54 particles per cell with grid spacing
$\Delta x\simeq(c/\omega_{\mathrm{p}})/2$.

To model SC in MRI turbulence, we add an RR force to the right-hand side of
the KSB-OA particle equations of motion [20]. As in several other works
including SC in PIC [37, 36, 38, 39, 40, 41], we write this term as

$\left(\frac{\textrm{d}\mbox{\boldmath{$u$}}}{\textrm{d}t}\right)_{\mathrm{RR}}=-\kappa\gamma\left[\left(\mbox{\boldmath{$E$}}+\frac{\mbox{\boldmath{$u$}}}{\gamma}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$B$}}\right)^{2}+\left(\frac{\mbox{\boldmath{$u$}}}{\gamma}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$E$}}\right)^{2}\right]\mbox{\boldmath{$u$}},$
(1)

where $u$ is the spatial part of the particle four-velocity,
$\gamma\equiv\sqrt{1+u^{2}/c^{2}}$ is the particle Lorentz factor, $E$ is the
electric field, and $\kappa$ is a numerical coefficient chosen such that RR
maintains a quasithermal state at $\gamma-1\sim 1$. The MRI in our setup grows
out of numerical noise; we evolve the initial state for a total time
$t_{\mathrm{end}}=15P_{0}$ ($P_{0}\equiv 2\pi/\Omega_{0}$).

Figure 1: Properties of mesoscale MRI turbulence. $(a1\mbox{--}3)$: Component-
wise magnetic-field distribution at $t=8.1P_{0}$; $(b)\mbox{--}(c)$: Average
magnetic-field amplification, plasma-$\beta$, and gyroradius throughout the
run; $(d)\mbox{--}(e)$: Average pressure anisotropy over time and
$(\beta_{\|},p_{\perp}/p_{\|})$-distribution with respect to mirror/firehose
approximate thresholds during the nonlinear stage; $(f)$: $\alpha$-model
stresses (contributing to AMT) over time. Figure 2: $(a)$: Toroidal and
poloidal magnetic-field power spectra versus $k_{\perp}$ averaged over
$t\in[7,8]P_{0}$ during the nonlinear stage; $(b)$: Likewise for compressive
and solenoidal velocity-fluctuation spectra. Figure 3: Particle energy
distributions throughout the MRI nonlinear stage, $(a)$ without SC and $(b)$
with SC activated for the whole run.

_Results._ — Via background shearing and small-scale perturbations, initially
vertical magnetic-field lines bend and $B$ acquires $x,y$-components. The
system’s evolution encompasses an initial stage, where small-scale
perturbations are excited; an exponential magnetic-field growth due to linear
MRI modes; and a nonlinear, saturated stage where primary-MRI channel flows
decay into a turbulent state. The channel decay is driven by parasitic
instabilities (e.g. tearing and kinking of current sheets at channel
interfaces) feeding off the primary MRI [34, 35]. These instabilities require
sufficiently large domains to develop, and can be suppressed if specific
conditions are not met [42]; our system size ensures that current-sheet
kinking and tearing efficiently destroy the primary MRI channels and initiate
the nonlinearly turbulent stage [20]. A typical spatial distribution of
magnetic-field components during the nonlinear stage is shown in Fig. 1
(panels $(a1\mbox{--}3)$), where we observe the chaotic configuration of the
MRI-driven turbulence. As outlined above, we inhibit NTPA until the turbulent
state is established by keeping SC active until right after the initial
channel disruption222Smaller simulations show that deactivating SC at later
times produces qualitatively similar results.. At $t\simeq 4P_{0}$, we remove
SC, such that particles can be energized by the MRI turbulence only, starting
from an approximately Maxwellian state (see below).

The right panels in Fig. 1 show the evolution of the physical properties of
MRI turbulence. Panel $(b)$ shows the evolution of each component of magnetic-
field energy, averaged over space. The exponential amplification of magnetic
fields, due to the primary MRI, can be observed until $t\sim 4P_{0}$. At that
time, the MRI channel flows decay into a sustained-turbulence state with a
statistically constant total magnetic energy. This quasi-steady state is due
to a balance between constant energy injection by the MRI and turbulent energy
cascade. When the turbulence reaches kinetic scales, plasma heating and high-
energy particle acceleration can occur. Panel $(c)$ shows that the volume-
averaged $\beta$ and average particle gyroradius slowly increase during the
nonlinear stage, as a consequence of turbulent heating333Note that $\beta\sim
10$ during the whole nonlinear stage, implying that turbulent models of disk
plasmas should consider a high-$\beta$ collisionless regime [43].. In this
turbulent state, we can evaluate AMT throughout the disk, e.g., via the
$\alpha$-parametrization [44], involving the Maxwell ($M_{xy}\equiv-
B_{x}B_{y}/(4\pi)$), Reynolds ($R_{xy}$), and anisotropic ($A_{xy}$) stresses.
Here, $R_{xy}\equiv mnU_{x}U_{y}/\Gamma^{2}$ is subdominant due to the lack of
macroscopically coherent flows in the turbulent accretion disk (with density
$n$ and bulk velocity
$\mbox{\boldmath{$U$}}\equiv(1/n)\int\mbox{\boldmath{$u$}}f(\mbox{\boldmath{$x$}},\mbox{\boldmath{$u$}},t)\textrm{d}^{3}\mbox{\boldmath{$u$}}$
calculated from the distribution function, $f$, and
$\Gamma\equiv\sqrt{1+U^{2}/c^{2}}$). We evaluate the anisotropic stress
$A_{xy}\equiv-(p_{\perp}-p_{\|})B_{x}B_{y}/B^{2}$ [45] by computing parallel
and perpendicular plasma pressure,
$p_{\|}\equiv\textbf{P}\mbox{\boldmath{$\,\mbox{\boldmath{$:$}}\,$}}\mbox{\boldmath{$B$}}\mbox{\boldmath{$B$}}/B^{2}$
and
$p_{\perp}\equiv\textbf{P}\,\mbox{\boldmath{$:$}}\,(\mathbb{I}-\mbox{\boldmath{$B$}}\mbox{\boldmath{$B$}}/B^{2})/2$,
where $\textbf{P}\equiv
m\int(\mbox{\boldmath{$u$}}\mbox{\boldmath{$u$}}/\gamma)f(\mbox{\boldmath{$x$}},\mbox{\boldmath{$u$}},t)\textrm{d}^{3}\mbox{\boldmath{$u$}}$
is the full pressure tensor; $A_{xy}$ does not exist in standard MHD
calculations with isotropic pressure. In panel $(d)$, we observe a positive
volume-averaged pressure anisotropy $\langle p_{\perp}-p_{\|}\rangle>0$ during
the entire nonlinear stage, indicating the presence of perpendicular-pressure
enhancement mechanisms (e.g., approximate conservation of the first adiabatic
invariant $p_{\perp}/(mnB)$ in growing fields). Perpendicular-pressure growth
is limited by kinetic effects, as shown in panel $(e)$ by the
$(\beta_{\|},p_{\perp}/p_{\|})$-distribution at a representative time
$t=8.1P_{0}$ during the turbulent stage. Here, we observe the presence of
mirror modes: as $p_{\perp}$ increases, approximate thresholds for mirror
instabilities (top dashed line) [46] are crossed, scattering particles
efficiently and limiting further increase in anisotropy. Similarly, the
firehose instability (bottom dashed line) [46] prevents the growth of $p_{\|}$
(due to turbulent reconnection and/or conservation of the adiabatic invariant
$p_{\|}B^{2}/(mn)^{3}$ [18] when $B$ is dissipated). Panel $(f)$ shows our
first-principles measurement of AMT: we observe that pressure anisotropy does
indeed substantially contribute to $\alpha$, via $\langle
A_{xy}\rangle\sim\langle M_{xy}\rangle$ at all times. The resulting
$\langle\alpha\rangle\equiv\langle R_{xy}+A_{xy}+M_{xy}\rangle/\langle
p\rangle\simeq 0.1$ throughout the nonlinear stage.

Fig. 2 (panel $(a)$) shows the isotropic power spectra of the poloidal
($B_{\mathrm{pol}}\equiv\sqrt{B_{x}^{2}+B_{z}^{2}}$) and toroidal ($B_{y}$)
magnetic field with respect to perpendicular wavenumbers $k_{\perp}$, averaged
over $t\in[7,8]P_{0}$ during the nonlinear stage. Our simulation allows us to
attain a decade-wide separation between MHD and kinetic scales, observable in
this figure: energy is injected around the most-unstable MRI wavenumber
$k_{\mathrm{MRI}}$; a distinct inertial range is identified by characteristic
power laws $\propto k_{\perp}^{-2}$ (for $B_{y}$) and $\propto
k_{\perp}^{-5/3}$ or $\propto k_{\perp}^{-3/2}$ (for $B_{\mathrm{pol}}$) at
$k_{\perp}\langle\rho_{\mathrm{c}}\rangle<1$. The spectral slopes in the
inertial range are compatible with earlier MHD and hybrid studies [47, 48]. We
then observe a spectral break around
$k_{\perp}\langle\rho_{\mathrm{c}}\rangle=1$, followed by a modification of
the spectra to power laws steeper than $\propto k_{\perp}^{-3}$.

Panel $(b)$ of Fig. 2 shows the isotropic power spectra of compressive
($\mbox{\boldmath{$V$}}_{\mathrm{compr}}$) and solenoidal
($\mbox{\boldmath{$V$}}_{\mathrm{sol}}$) bulk-velocity fluctuations,
calculated from
$\mbox{\boldmath{$V$}}\equiv\mbox{\boldmath{$U$}}/\Gamma=\mbox{\boldmath{$V$}}_{\mathrm{compr}}+\mbox{\boldmath{$V$}}_{\mathrm{sol}}$
such that
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$V$}}_{\mathrm{compr}}=\mbox{\boldmath{$0$}},\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$V$}}_{\mathrm{sol}}=0$.
We observe that the $\mbox{\boldmath{$V$}}_{\mathrm{sol}}$ spectrum possesses
a spectral slope $k_{\perp}^{-2}$ in the MHD range;
$\mbox{\boldmath{$V$}}_{\mathrm{compr}}$ follows a roughly similar scaling,
although numerical noise makes it harder to identify a clear spectral slope.
These spectra are fundamental for identifying the type of turbulent
fluctuations that contribute to plasma heating, including differential-heating
mechanisms of electron-ion plasmas in low-luminosity accretion flows [49, 50,
51]. In contrast with reduced-MHD expectations [51], we observe here that
solenoidal (predominantly Alfvénic) fluctuations strongly dominate, in terms
of spectral power, at all scales (including sub-Larmor scales, where, however,
numerical noise impacts the results).

Particle energy spectra $f(\gamma)\equiv\textrm{d}N/\textrm{d}\gamma$ are
shown in Fig. 3 for $t\in[4,15]P_{0}$, i.e., during the nonlinear stage. When
SC is turned off at $t=4P_{0}$, the distribution lacks strong nonthermal
features444In contrast with previous studies, where the initial distribution
of the turbulence stage was already strongly nonthermal.. We then observe the
effect of MRI turbulence on this initial distribution: Fig. 3 (panel $(a)$)
showcases continuous plasma heating, indicated by a shift of the distribution
peak toward higher energies over time. Concurrently, a steep nonthermal tail
$f(\gamma)\sim\gamma^{-4}$ develops and evolves toward
$f(\gamma)\sim\gamma^{-2.5}$ at the end of the run. At that point, particles
have accumulated around the highest possible energy determined by
$\rho_{\mathrm{c}}\sim L_{z}$ as also observed in previous MRI simulations
[20]. The power-law tail is less steep than in externally driven-turbulence
simulations of similar plasma beta [52]. For comparison, in panel $(b)$ we
show spectra for an equivalent simulation (same initial parameters except for
$\Delta x\simeq c/\omega_{\mathrm{p}}$), where SC remains active throughout
the nonlinear stage. Here, $f(\gamma)$ settles onto a low-energy Maxwellian
state that is maintained indefinitely, resulting from a balance between energy
gain from turbulence and energy loss from SC. This is consistent with previous
externally driven turbulence simulations of strongly radiative pair plasmas
[53, 54]. Despite the lack of NTPA, the radiative and nonradiative simulations
proceed almost indistinguishably, in terms of magnetic-field amplification,
average pressure anisotropy, and turbulent spectra (not shown), at least over
the simulated time scale $t\leq 15P_{0}$.

_Discussion._ — We have presented the largest (to date) 3D PIC-SB simulations
of collisionless pair-plasma MRI turbulence, with and without radiation-
reaction forces. Our unprecedentedly large runs are critical for determining
the mesoscale steady-state behavior of the 3D MRI from first principles. We
found that magnetic-field amplification and overall MRI dynamics, from the
initial linear instability to the quasi-stationary nonlinear turbulence,
broadly align with MHD and hybrid expectations: a saturated $\beta\sim 10$ is
maintained throughout the nonlinear stage, with plasma slowly heating due to
cascading turbulence from MHD down to kinetic scales. However, our fully
kinetic runs also add novel, fundamental information to the picture.

For the first time, we obtain an MHD-like, fully kinetic ordering of stresses
contributing to AMT in collisionless MRI turbulence. Braginskii-MHD
simulations [55, 56] have shown that substantial AMT can occur in anisotropic-
pressure fluid models. Those simulations necessarily employ ad-hoc closures
from kinetic plasma theory to limit pressure-anisotropy growth. Enhanced AMT
was also observed in hybrid (fluid electrons/kinetic ions) SB simulations
[47]. These earlier models provided an expectation for the hierarchy of the
different stresses contributing to AMT, showing that $A_{xy}\leq M_{xy}$. All
previous 2D/3D PIC-MRI studies [16, 17, 18, 19, 20]) showed AMT enhancement,
but often with a reversed ordering $M_{xy}\leq A_{xy}$ especially during the
late nonlinear stage, where $M_{xy}\ll A_{xy}$ was typically observed. This
was likely due to limited separation between macro- and microscopic scales.
The present study employs a scale separation
$\omega_{\mathrm{c},0}/\Omega_{0}$ twice the largest employed among all
previous 3D works [20]; this proved sufficient to recover the expected fluid
trend $A_{xy}\lesssim M_{xy}$. We also showed that pressure anisotropy,
regulating $A_{xy}$, remains well-controlled by fast mirror and firehose
instabilities, which are naturally included in our first-principles model. Our
results could be used as an input, e.g., for (global/local) Braginskii-MHD
simulations of accretion, replacing ad-hoc mirror/firehose thresholds with
PIC-derived functions of local plasma conditions.

Exploiting our large scale separation, we analyzed the MRI-turbulent magnetic
and velocity-fluctuation spectra from macroscopic to microscopic scales, over
two decades in $k_{\perp}$. Magnetic spectra show a clear steepening of
spectral slopes, from $\sim-2$ (for $B_{y}$) and $\sim-5/3$ or $\sim-3/2$ (for
$B_{\mathrm{pol}}$) at large scales, aligning with MHD expectations [48], to
slopes $<-3$ at kinetic scales. Presently, our results cannot be compared with
analytic expectations for spectral slopes, since a theory of MRI turbulence at
sub-Larmor scales has not been fully developed. Our findings are therefore
important for comparison with such analytic endeavors in the future. Standard
pair-plasma-turbulence simulations, for which firmer analytic expectations
exist, have shown a steepening of magnetic spectra at kinetic scales, broadly
aligned with our results [52]. Intriguingly, our velocity-fluctuation spectra
are instead in direct contrast with earlier reduced-MHD MRI-turbulence
estimates, which showed a dominance of compressive over solenoidal
fluctuations [51], while we observe the opposite. Importantly, those previous
MHD studies considered a substantially different setup with a strong,
unidirectional toroidal field throughout the disk. This is not the case for
our runs, where a strong mean $B_{y}$ is absent. A thorough setup comparison
deserves a dedicated study to be pursued elsewhere.

Finally, the main result of this study is the first definitive demonstration
that MRI turbulence can efficiently accelerate particles to nonthermal
energies. Our quantification is, for the first time, agnostic of initial
nonthermal particle injection, i.e. we completely separated NTPA due to
steady-state MRI turbulence from “spurious” NTPA due to initial MRI stages,
where nonthermal particles are preaccelerated by magnetic reconnection
associated with the collapse of transient macroscopic channel flows. We
suppressed this effect by including SC at preturbulence stages, and turning SC
off upon reaching the nonlinear stage. We then followed the evolution of
particle energy spectra from a “clean”, near-Maxwellian state, measuring
significant NTPA creating high-energy power-law tails. This also marks the
first time that radiation reaction is included self-consistently in kinetic
SB-MRI simulations. For demonstration, we also conducted a full simulation
where SC is kept active at all times; there, a stationary thermal distribution
arises from the balance between MRI energy injection and RR backreaction, as
expected. Our addition of RR forces is fundamental for future, more realistic
electron-ion semirelativistic MRI scenarios [E. Gorbunov et al., in prep.],
where electron SC must be included in simulations to correctly interpret SMBH
observations [1, 2].

###### Acknowledgements.

F.B. would like to thank Benjamin Crinquand for discussions and guidance on
the implementation of radiation-reaction forces. F.B. also thanks Lorenzo
Sironi, Giuseppe Arrò, Mario Riquelme, and Astor Sandoval for useful
discussions throughout the development of this work. F.B. acknowledges support
from the FED-tWIN programme (profile Prf-2020-004, project “ENERGY”) issued by
BELSPO. The authors acknowledge support from a NASA ATP grant 80NSSC20K0545.
Additional support was provided by National Science Foundation grants AST
1806084 and AST 1903335 to the University of Colorado. Support for L.A. was
provided by the Institute for Advanced Study. Support for V.Z. at the Flatiron
Institute is provided by the Simons Foundation. This work was performed in
part at Aspen Center for Physics, which is supported by National Science
Foundation grant PHY-1607611. This work made use of substantial computational
resources provided by several entities. We acknowledge support by the VSC
(Flemish Supercomputer Center), funded by the Research Foundation – Flanders
(FWO) and the Flemish Government – department EWI. An award of computer time
was provided by the Innovative and Novel Computational Impact on Theory and
Experiment (INCITE) program, using resources (namely Theta) of the Argonne
Leadership Computing Facility, which is a DOE Office of Science User Facility
supported under contract DE-AC02-06CH11357. Additional computer time was
provided by the Extreme Science and Engineering Discovery Environment (XSEDE),
which is supported by National Science Foundation grant number ACI-1548562
[57], and by the Frontera computing project [58] at the Texas Advanced
Computing Center (TACC, www.tacc.utexas.edu). We acknowledge TACC for
providing HPC resources on both Stampede2 and Frontera. Additional resources
(namely Pleiades) supporting this work were provided by the NASA High-End
Computing (HEC) Program through the NASA Advanced Supercomputing (NAS)
Division at Ames Research Center.

## References

  * EHT Collaboration [2019] EHT Collaboration, ApJL 875, L5 (2019).
  * EHT Collaboration [2022] EHT Collaboration, ApJL 930, L16 (2022).
  * GRAVITY Collaboration  _et al._ [2019] GRAVITY Collaboration, R. Abuter, A. Amorim, M. Bauböck, J. P. Berger, H. Bonnet, W. Brandner, Y. Clénet, V. Coudé Du Foresto, P. T. de Zeeuw, J. Dexter, G. Duvert, A. Eckart, F. Eisenhauer, N. M. Förster Schreiber, P. Garcia, F. Gao, E. Gendron, R. Genzel, O. Gerhard, S. Gillessen, M. Habibi, X. Haubois, T. Henning, S. Hippler, M. Horrobin, A. Jiménez-Rosales, L. Jocou, P. Kervella, S. Lacour, V. Lapeyrère, J. B. Le Bouquin, P. Léna, T. Ott, T. Paumard, K. Perraut, G. Perrin, O. Pfuhl, S. Rabien, G. Rodriguez Coira, G. Rousset, S. Scheithauer, A. Sternberg, O. Straub, C. Straubmeier, E. Sturm, L. J. Tacconi, F. Vincent, S. von Fellenberg, I. Waisberg, F. Widmann, E. Wieprecht, E. Wiezorrek, J. Woillez, and S. Yazici, A & A 625, L10 (2019).
  * Gammie  _et al._ [2003] C. F. Gammie, J. McKinney, and G. Tóth, ApJ 589, 444 (2003).
  * Narayan  _et al._ [2012] R. Narayan, A. Sa̧dowski, R. F. Penna, and A. K. Kulkarni, Monthly Notices of the Royal Astronomical Society 426, 3241 (2012).
  * Porth  _et al._ [2019] O. Porth, K. Chatterjee, R. Narayan, C. F. Gammie, Y. Mizuno, P. Anninos, J. G. Baker, M. Bugli, C. kwan Chan, J. Davelaar, L. D. Zanna, Z. B. Etienne, P. C. Fragile, B. J. Kelly, M. Liska, S. Markoff, J. C. McKinney, B. Mishra, S. C. Noble, H. Olivares, B. Prather, L. Rezzolla, B. R. Ryan, J. M. Stone, N. Tomei, C. J. White, Z. Younsi, K. Akiyama, A. Alberdi, W. Alef, K. Asada, R. Azulay, A.-K. Baczko, D. Ball, M. Baloković, J. Barrett, D. Bintley, L. Blackburn, W. Boland, K. L. Bouman, G. C. Bower, M. Bremer, C. D. Brinkerink, R. Brissenden, S. Britzen, A. E. Broderick, D. Broguiere, T. Bronzwaer, D.-Y. Byun, J. E. Carlstrom, A. Chael, S. Chatterjee, M.-T. Chen, Y. Chen, I. Cho, P. Christian, J. E. Conway, J. M. Cordes, B. C. Geoffrey, Y. Cui, M. D. Laurentis, R. Deane, J. Dempsey, G. Desvignes, S. S. Doeleman, R. P. Eatough, H. Falcke, V. L. Fish, E. Fomalont, R. Fraga-Encinas, B. Freeman, P. Friberg, C. M. Fromm, J. L. Gómez, P. Galison, R. García, O. Gentaz, B. Georgiev, C. Goddi, R. Gold, M. Gu, M. Gurwell, K. Hada, M. H. Hecht, R. Hesper, L. C. Ho, P. Ho, M. Honma, C.-W. L. Huang, L. Huang, D. H. Hughes, S. Ikeda, M. Inoue, S. Issaoun, D. J. James, B. T. Jannuzi, M. Janssen, B. Jeter, W. Jiang, M. D. Johnson, S. Jorstad, T. Jung, M. Karami, R. Karuppusamy, T. Kawashima, G. K. Keating, M. Kettenis, J.-Y. Kim, J. Kim, J. Kim, M. Kino, J. Y. Koay, M. K. Patrick, S. Koyama, M. Kramer, C. Kramer, T. P. Krichbaum, C.-Y. Kuo, T. R. Lauer, S.-S. Lee, Y.-R. Li, Z. Li, M. Lindqvist, K. Liu, E. Liuzzo, W.-P. Lo, A. P. Lobanov, L. Loinard, C. Lonsdale, R.-S. Lu, N. R. MacDonald, J. Mao, D. P. Marrone, A. P. Marscher, I. Martí-Vidal, S. Matsushita, L. D. Matthews, L. Medeiros, K. M. Menten, I. Mizuno, J. M. Moran, K. Moriyama, M. Moscibrodzka, C. Müller, H. Nagai, N. M. Nagar, M. Nakamura, G. Narayanan, I. Natarajan, R. Neri, C. Ni, A. Noutsos, H. Okino, T. Oyama, F. Özel, D. C. M. Palumbo, N. Patel, U.-L. Pen, D. W. Pesce, V. Piétu, R. Plambeck, A. PopStefanija, J. A. Preciado-López, D. Psaltis, H.-Y. Pu, V. Ramakrishnan, R. Rao, M. G. Rawlings, A. W. Raymond, B. Ripperda, F. Roelofs, A. Rogers, E. Ros, M. Rose, A. Roshanineshat, H. Rottmann, A. L. Roy, C. Ruszczyk, K. L. J. Rygl, S. Sánchez, D. Sánchez-Arguelles, M. Sasada, T. Savolainen, F. P. Schloerb, K.-F. Schuster, L. Shao, Z. Shen, D. Small, B. W. Sohn, J. SooHoo, F. Tazaki, P. Tiede, R. P. J. Tilanus, M. Titus, K. Toma, P. Torne, T. Trent, S. Trippe, S. Tsuda, I. van Bemmel, H. J. van Langevelde, D. R. van Rossum, J. Wagner, J. Wardle, J. Weintroub, N. Wex, R. Wharton, M. Wielgus, G. N. Wong, Q. Wu, K. Young, A. Young, F. Yuan, Y.-F. Yuan, J. A. Zensus, G. Zhao, S.-S. Zhao, Z. Zhu, and T. E. H. T. Collaboration), The Astrophysical Journal Supplement Series 243, 26 (2019).
  * Tomei  _et al._ [2019] N. Tomei, L. Del Zanna, M. Bugli, and N. Bucciantini, Monthly Notices of the Royal Astronomical Society 491, 2346 (2019).
  * Dexter  _et al._ [2020] J. Dexter, A. Tchekhovskoy, A. Jiménez-Rosales, S. M. Ressler, M. Bauböck, Y. Dallilar, P. T. de Zeeuw, F. Eisenhauer, S. von Fellenberg, F. Gao, R. Genzel, S. Gillessen, M. Habibi, T. Ott, J. Stadler, O. Straub, and F. Widmann, Monthly Notices of the Royal Astronomical Society 497, 4999 (2020).
  * Ripperda  _et al._ [2020] B. Ripperda, F. Bacchini, and A. A. Philippov, The Astrophysical Journal 900, 100 (2020).
  * Ripperda  _et al._ [2022] B. Ripperda, M. Liska, K. Chatterjee, G. Musoke, A. A. Philippov, S. B. Markoff, A. Tchekhovskoy, and Z. Younsi, The Astrophysical Journal Letters 924, L32 (2022).
  * Ressler  _et al._ [2015] S. M. Ressler, A. Tchekhovskoy, E. Quataert, M. Chandra, and C. F. Gammie, MNRAS 454, 1848 (2015).
  * Chael  _et al._ [2018] A. Chael, M. Rowan, R. Narayan, M. Johnson, and L. Sironi, MNRAS 478, 4 (2018).
  * Dexter  _et al._ [2021] J. Dexter, N. Scepi, and M. C. Begelman, ApJL 919, L20 (2021).
  * Hankla  _et al._ [2022] A. M. Hankla, N. Scepi, and J. Dexter, Monthly Notices of the Royal Astronomical Society 515, 775 (2022).
  * Scepi  _et al._ [2022] N. Scepi, J. Dexter, and M. C. Begelman, MNRAS 511, 3536 (2022).
  * Riquelme  _et al._ [2012] M. A. Riquelme, E. Quataert, P. Sharma, and A. Spitkovsky, ApJ 755, 50 (2012).
  * Hoshino [2013] M. Hoshino, ApJ 773, 118 (2013).
  * Hoshino [2015] M. Hoshino, PRL 114, 061101 (2015).
  * Inchingolo  _et al._ [2018] G. Inchingolo, T. Grismayer, N. F. Loureiro, R. A. Fonseca, and L. O. Silva, ApJ 859, 149 (2018).
  * Bacchini  _et al._ [2022] F. Bacchini, L. Arzamasskiy, V. Zhdankin, G. Werner, M. Begelman, and D. Uzdensky, ApJ 938, 86 (2022).
  * Parfrey  _et al._ [2019] K. Parfrey, A. Philippov, and B. Cerutti, PRL 122, 035101 (2019).
  * Crinquand  _et al._ [2020] B. Crinquand, B. Cerutti, A. Philippov, K. Parfrey, and G. Dubus, PRL 124, 145101 (2020).
  * El Mellah  _et al._ [2022] I. El Mellah, B. Cerutti, B. Crinquand, and K. Parfrey, A&A 663, A169 (2022).
  * Galishnikova  _et al._ [2023] A. Galishnikova, A. Philippov, E. Quataert, F. Bacchini, K. Parfrey, and B. Ripperda, PRL 130, 115201 (2023).
  * Sandoval  _et al._ [2023] A. Sandoval, M. Riquelme, A. Spitkovsky, and F. Bacchini, arXiv:2308.12348 (2023).
  * Hill [1878] G. W. Hill, American Journal of Mathematics 1, 5 (1878).
  * Hawley  _et al._ [1995] J. F. Hawley, C. F. Gammie, and S. A. Balbus, ApJ 440, 742 (1995).
  * Matsumoto and Tajima [1995] R. Matsumoto and T. Tajima, ApJ 445, 767 (1995).
  * Stone  _et al._ [1996] J. M. Stone, J. F. Hawley, C. F. Gammie, and S. A. Balbus, ApJ 463, 656 (1996).
  * Sano and Inutsuka [2001] T. Sano and S. Inutsuka, ApJ 561, L179 (2001).
  * Bai and Stone [2013] X. Bai and J. M. Stone, ApJ 767, 1 (2013).
  * Hirai  _et al._ [2018] K. Hirai, Y. Katoh, N. Terada, and S. Kawai, ApJ 853, 174 (2018).
  * Balbus and Hawley [1991] S. A. Balbus and J. F. Hawley, ApJ 376, 214 (1991).
  * Goodman and Xu [1994] J. Goodman and G. Xu, ApJ 432, 213 (1994).
  * Pessah and Goodman [2009] M. E. Pessah and J. Goodman, ApJ 698, L72 (2009).
  * Cerutti  _et al._ [2013] B. Cerutti, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, ApJ 770, 147 (2013).
  * Tamburini  _et al._ [2011] M. Tamburini, F. Pegoraro, A. Di Piazza, C. Keitel, T. Liseykina, and A. Macchi, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 653, 1 (2011).
  * Cerutti  _et al._ [2016] B. Cerutti, A. Philippov, and A. Spitkovsky, MNRAS 457, 3 (2016).
  * Comisso and Sironi [2021] L. Comisso and L. Sironi, Phys. Rev. Lett. 127, 255102 (2021).
  * Hakobyan  _et al._ [2023] H. Hakobyan, B. Ripperda, and A. A. Philippov, The Astrophysical Journal Letters 943, L29 (2023).
  * Schoeffler  _et al._ [2023] K. M. Schoeffler, T. Grismayer, D. Uzdensky, and L. O. Silva, MNRAS 523, 3812 (2023).
  * Latter  _et al._ [2009] H. N. Latter, P. Lesaffre, and S. A. Balbus, MNRAS 394, 715 (2009).
  * Arzamasskiy  _et al._ [2023] L. Arzamasskiy, M. W. Kunz, J. Squire, E. Quataert, and A. A. Schekochihin, Phys. Rev. X 13, 021014 (2023).
  * Shakura and Sunyaev [1973] N. I. Shakura and R. A. Sunyaev, A&A 500, 33 (1973).
  * Kulsrud [1983] R. M. Kulsrud, Handbook of plasma physics 1, 115 (1983).
  * Kunz  _et al._ [2014] M. W. Kunz, A. A. Schekochihin, and J. M. Stone, PRL 112, 205003 (2014).
  * Kunz  _et al._ [2016] M. W. Kunz, J. M. Stone, and E. Quataert, PRL 117, 235101 (2016).
  * Walker  _et al._ [2016] J. Walker, G. Lesur, and S. Boldyrev, MNRAS 457, L39 (2016).
  * Kawazura  _et al._ [2020] Y. Kawazura, A. A. Schekochihin, M. Barnes, J. M. TenBarge, Y. Tong, K. G. Klein, and W. Dorland, PRX 10, 041050 (2020).
  * Zhdankin [2021] V. Zhdankin, The Astrophysical Journal 922, 2 (2021).
  * Kawazura  _et al._ [2022] Y. Kawazura, A. A. Schekochihin, M. Barnes, W. Dorland, and S. A. Balbus, JPP 88, 905880311 (2022).
  * Zhdankin  _et al._ [2017] V. Zhdankin, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, PRL 118, 055103 (2017).
  * Zhdankin  _et al._ [2020] V. Zhdankin, D. A. Uzdensky, G. R. Werner, and M. C. Begelman, Monthly Notices of the Royal Astronomical Society 493, 1 (2020).
  * Grošelj  _et al._ [2023] D. Grošelj, H. Hakobyan, A. M. Beloborodov, L. Sironi, and A. A. Philippov, arXiv:2301.11327 (2023).
  * Sharma  _et al._ [2006] P. Sharma, G. W. Hammett, E. Quataert, and J. M. Stone, ApJ 637, 952 (2006).
  * Kempski  _et al._ [2019] P. Kempski, E. Quataert, J. Squire, and M. W. Kunz, Monthly Notices of the Royal Astronomical Society 486, 3 (2019).
  * Towns  _et al._ [2014] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, R. Roskies, J. R. Scott, and N. Wilkens-Diehr, Comput. Sci. Eng. 16, 62 (2014).
  * Stanzione  _et al._ [2020] D. Stanzione, J. West, R. T. Evans, T. Minyard, O. Ghattas, and D. K. Panda, “Frontera: The evolution of leadership computing at the national science foundation,” (Association for Computing Machinery, New York, NY, USA, 2020) p. 106–111.