# Human Reactions to Incorrect Answers from Robots
††thanks: This material is based upon work supported by the National Science
Foundation under Grants #DUE-2142360, #IUSE-2142428, and #IUSE-1730568. All
study activities were supervised by UNR’s IRB office.

Ponkoj Chandra Shill  Department of CSE
University of Nevada
Reno, USA
<EMAIL_ADDRESS>Md. Azizul Hakim  Department of CSE
University of Nevada
Reno, USA
<EMAIL_ADDRESS>Muhammad Jahanzeb Khan  Department of CSE
University of Nevada
Reno, USA
<EMAIL_ADDRESS>Bashira Akter Anima  Department of CSE
University of Nevada
Reno, USA
<EMAIL_ADDRESS>
###### Abstract

As robots grow more and more integrated into numerous industries, it is
critical to comprehend how humans respond to their failures. This paper
systematically studies how trust dynamics and system design are affected by
human responses to robot failures. The three-stage survey used in the study
provides a thorough understanding of human-robot interactions. While the
second stage concentrates on interaction details, such as robot precision and
error acknowledgment, the first stage collects demographic data and initial
levels of trust. In the last phase, participants’ perceptions are examined
after the encounter, and trust dynamics, forgiveness, and propensity to
suggest robotic technologies are evaluated. Results show that participants’
trust in robotic technologies increased significantly when robots acknowledged
their errors or limitations to participants and their willingness to suggest
robots for activities in the future points to a favorable change in
perception, emphasizing the role that direct engagement has in influencing
trust dynamics. By providing useful advice for creating more sympathetic,
responsive, and reliable robotic systems, the study advances the science of
human-robot interaction and promotes a wider adoption of robotic technologies.

###### Index Terms:

Human-Robot Interaction(HRI), Error Acknowledgment by Robots, Trust Dynamics
in HRI, ANOVA.

##  I Introduction

The field of Human-Robot Interaction (HRI) has witnessed remarkable
advancements in recent years, propelling robots beyond mere tools and into the
realm of intelligent agents capable of complex interactions and decision-
making processes [1]. Robots are no longer confined to controlled
environments; they are increasingly integrated into various sectors such as
healthcare, manufacturing, retail, and customer service, interacting directly
with humans. As robots take on roles that involve conveying information, the
question of how humans react to errors made by these machines becomes
paramount [2].

One of the significant challenges in this context is understanding how humans
perceive and respond to robot errors, particularly in situations where robots
are expected to provide accurate information or guidance [3]. Errors made by
robots can lead to a range of outcomes, from mistrust and frustration to
potentially hazardous consequences, depending on the context. Investigating
human reactions to incorrect or misleading information provided by robots is
not merely an academic exercise but holds crucial implications for the design
and acceptance of robotic systems [4].

Figure 1: Participants are interacting with the robot

Several recent studies have shed light on various facets of human-robot
interaction and responses to robot errors. In [5], the researchers highlight
the significance of multimodal cues in detecting and addressing conversational
failures in HRI. This research underscores the complexity of human responses
to errors in robot communication, emphasizing the role of non-verbal cues in
this process.

Despite these valuable contributions to the field of HRI, there remains a gap
in understanding how these insights can be applied to real-world scenarios and
the broader implications for trust and acceptance of robotic technologies.
This research project aims to address this gap by conducting a comprehensive
study on human responses to robot errors and their implications for human-
robot trust dynamics and system design.

Based on the objectives of the study, the following hypotheses are proposed:

  * H1:

The extent to which a robot acknowledges its limitations or errors can
significantly enhance human trust and emotional acceptance.

  * H2:

Humans are more forgiving of informational errors when they come from familiar
sources like other humans or internet search engines, compared to robots.

This study is conducted through a three-stage survey designed to capture a
comprehensive view of human-robot interactions. By systematically
investigating these aspects, this research aims to bridge a critical knowledge
gap in HRI. It contributes to a deeper understanding of human-robot trust
dynamics and provides insights for developing more effective, reliable, and
user-friendly robotic systems. This study is poised to enhance the field of
HRI by offering a nuanced perspective on how human reactions to robotic errors
can inform the design of more empathetic, responsive, and context-aware
robots.

Figure 2: Research Methodology

##  II Related Work

The interaction between humans and robots, particularly in scenarios involving
errors, has been the subject of extensive research. Studies have explored
various facets of human-robot interaction (HRI), ranging from the
psychological impacts of robot errors to the nuances of multimodal human
responses. This section discusses key findings from previous research and
contextualizes how our study contributes to the existing body of knowledge,
addressing gaps in the literature.

Previous research has extensively examined how humans perceive and react to
robot errors. For instance, studies have shown that human tolerance and
reactions to robot mistakes vary based on the context and nature of the
interaction. Yasuda and Matsumoto (2013) explored the psychological impact on
humans when a robot makes errors, suggesting that errors can significantly
influence human trust and perception of robots [6].

Traeger et al. delved into the dynamics of robots admitting mistakes and its
impact on human-robot conversations. Their findings indicated that robot error
admission could lead to more effective and constructive interactions, thereby
influencing the trust and reliability in human-robot relationships [7].

Roya Salek Shahrezaie et al. investigates how homophily, the tendency to
connect with similar others, influences human-robot interactions (HRI) and
trust. Using two NAO robots, it examines whether sharing interests affects
interaction quality. Results support the idea that similarity fosters trust in
HRI. The related work section discusses previous studies on homophily in HRI,
emphasizing its relevance for designing socially assistive robots. [8]

Kontogiorgos et al. conducted a systematic cross-corpus analysis of human
reactions to robot conversational failures. This study analyzed multimodal
behavioral responses, including facial expressions and acoustic features, to
understand how humans react to different types of robot failures [5].

Hayes et al. examined nonverbal behaviors in LfD settings, focusing on how
participants nonverbally communicated with robots during teaching tasks. This
work emphasized the importance of gestures and facial expressions in conveying
feedback to robots, particularly in scenarios where the robot committed errors
during the learning process [9].

Salem et al. investigated the effects of different types of errors on trust in
human-robot cooperation. Their findings underscored that the context of the
interaction, such as the task’s criticality, significantly affects how humans
perceive and react to robot errors [10].

Stiber et al. highlighted how individual differences in expressiveness affect
responses to robot behaviors, suggesting that cultural background and personal
traits can influence the perception and tolerance of robotic errors [11].

Our research extends these foundational studies by focusing on a detailed
examination of nonverbal human responses in a specific robotic interaction
context. We aim to fill the gap in understanding the nuanced ways humans react
to robots in task-oriented environments, particularly in scenarios involving
errors and miscommunications. Our study offers new insights into the
subtleties of human behavior in HRI, contributing to the development of more
adaptive and sensitive robotic systems that can effectively navigate complex
social interactions.

Figure 3: Nao

##  III Experimental Setup

System Configuration: In this experiment, we used the NAO robot, a humanoid
platform renowned for its agility and interactive capabilities depicts in Fig.
3. Nao was positioned directly across from the participants in a lab
environment. It started the session by greeting them and encouraging them to
ask questions from a predetermined set of ten questions designed to evaluate
its conversational abilities. Nao is equipped with advanced sensors and a 1.91
GHz Quad core Atom E3845 Processor and 4 GB DDR3 ram , which enables real-time
data collection and efficient processing of complex interactions. This
combination of features facilitates a comprehensive assessment of human-robot
dialogue. It also highlights Nao’s suitability for intricate interaction
studies, allowing it to navigate the nuances of human engagement effectively
through computational intelligence.

Hardware and Software Integration: The Nao robot was interconnected with a
personal computer serving as the speech recognition and response generation
processing unit. This integration was facilitated by the Robot Operating
System (ROS), which provided a flexible framework for controlling the robot.
The speech interaction mechanism was engineered such that the participant’s
spoken queries were captured by Nao’s built-in microphone, transmitted to the
computer, and processed using an advanced speech recognition application.

Speech Recognition and Natural Language Processing: Our speech recognition
application leveraged the ”speech_recognition” Python library to convert
spoken language into text form. Leveraging the Natural Language Toolkit
(NLTK), we performed various natural language processing (NLP) tasks such as
tokenization, lemmatization, and stop-word analysis. These NLP techniques
enabled the robot to comprehend human language more effectively by allowing it
to discern the context and significance of the queries posed by the
participants.

The NLP application was intricately designed to sift through the participant’s
speech, searching for specific keywords or phrases that aligned with any of
the ten pre-selected questions. This identification process was crucial for
the system’s ability to accurately understand and categorize the participant’s
inquiries. Upon successfully matching the participant’s input to one of these
predetermined questions, the application then proceeded to generate a relevant
answer. This answer was formulated as text and sent back to the Nao robot.
Utilizing its integrated speaker system, the Nao robot converted the text-
based response into audible speech, effectively communicating the answer back
to the participant. This seamless integration of speech recognition, natural
language processing, and robotic vocalization exemplified the advanced
capabilities of the system to facilitate interactive and meaningful human-
robot dialogues.

Implementation of Error Handling and Apology Mechanisms: To investigate the
impact of robot errors on human trust and perception, our experimental
protocol deliberately included a mechanism for generating incorrect responses.
A randomized algorithm intermittently selected responses that were
intentionally erroneous, simulating real-world scenarios where robots may
misinterpret or incorrectly process human input. Upon issuing a mistaken
response, our system was designed to automatically detect the anomaly and
initiate an apology sequence. This sequence involved the Nao robot
acknowledging its error and expressing regret to the participant, either
through a direct apology or, in some cases, by attempting to correct the
misinformation with the accurate answer. This approach not only allowed us to
observe participant reactions to robot failures but also provided insights
into the dynamics of trust recovery and forgiveness in human-robot
interactions. By integrating these error-handling and apology mechanisms, our
study offers a comprehensive view of how such factors influence the overall
efficacy and acceptance of robotic technologies in interactive settings.

Following the completion of the ten-question interaction sequence, the Nao
robot expressed gratitude towards the participant, marking the end of the
session. This experimental setup provided a rich framework for studying human-
robot interaction dynamics, particularly focusing on the human response to
robot errors and the impact of such interactions on trust and perception of
robotic technologies.

##  IV Methodology

The research project’s methodology is meticulously planned to gauge
participants’ initial emotional responses and reflective thoughts following
interactions with a robot. A mixed-methods strategy is applied, integrating
qualitative and quantitative data via a three stage survey approach.

###  IV-A Three Stage Study

####  IV-A1 Preliminary Survey

This stage collects demographic data and baseline attitudes toward robots. We
conducted this survey before the robot interaction which involves a detailed
survey with specific questions designed to gather background information and
preconceptions about robots and automation from the participants. The survey
includes Age Group, Gender, Country, Prior Interaction with Robots, Perception
of Robots, and Trust in Automation. This survey sets a baseline for
understanding each participant’s initial stance towards robots and automated
systems, which is essential for analyzing any shifts in perception following
the robot interaction.

####  IV-A2 Interaction with Robot

This stage occurs immediately after the participants’ interaction with Nao. In
this stage, participants answer a series of survey questions about their
experience. They report the number of questions they asked Nao, evaluate the
correctness of Nao’s responses, indicate their awareness of any errors made by
the robot, and observe whether Nao apologized for its mistakes. This post-
interaction survey is crucial for capturing the participants’ immediate
reactions and perceptions of the robot’s performance, accuracy, and response
to errors. This stage is designed to quantitatively and qualitatively evaluate
the interaction, focusing on the robot’s accuracy and error-handling
capabilities.

####  IV-A3 Reflective Survey

Assesses the change in perceptions and trust levels post-interaction. This
stage is crucial in the methodology for gathering vital survey data post-
interaction with Nao. This stage delves into the participants’ nuanced
reactions and perceptions, focusing on several key areas such as Error Impact
on Trust, Trust Impact of Unacknowledged Error, Comparative Forgiveness (Human
vs. Robot), Effectiveness of Robot’s Apology, Importance of Acknowledging
Errors, Willingness to Rely on Robot Post-Error, Future Interaction Intent,
Perceived Reliability of Robots, and Recommendation Likelihood. This stage is
integral to the study as it seeks to understand the depth and breadth of human
responses to robotic performance, particularly in the context of trust and
reliability. The collected data will provide rich insights into the dynamics
of human-robot interaction and the psychological impacts of robotic errors and
their acknowledgments.

The methodology depicted in Fig. 2 maps out the process of a human-robot
interaction study. Initially, participants complete a survey gauging their
demographics, education, prior experience with robots, and trust in
automation. They then interact with NAO[12], in a series of interactions
during which they pose predetermined questions as shown in the Fig. 1.
Following the interaction, participants complete two further surveys that
capture their immediate and reflective responses to the robot’s performance.
The collected data from these surveys is methodically stored for in-depth
analysis to conclude the impact of robotic interactions on human trust and
perception. The design allows for a comprehensive analysis of the human-robot
interface, capturing the depth and evolution of participant responses.

###  IV-B Participant Selection

The participant group for this study comprises 31 individuals from different
countries between 21 and 64 years of age where 9 of them are female and 22 of
them are male. This group is largely inexperienced in interacting with robots
previously, as most participants(17 out of 31) are engaging in such
interactions for the first time. Their general perception of robots leans
towards viewing them as ”somewhat reliable and useful.” The trust in
automation technology among these participants varies, with a moderate level
of trust being the most common trend. These insights suggest that the group,
while diverse, shares an optimism about the potential of robotic systems
despite limited prior exposure. This demographic appears open to advancements
in human-robot interactions (HRI), emphasizing the significance of building
trust in automation for broader acceptance of future technologies. The
findings point towards the necessity of developing user-friendly interfaces
and ensuring transparency in robot functionalities, which could further
improve perceptions and trust in such technologies.

###  IV-C Interaction with Robot

In the robot interaction stage of the study, participants individually engage
with the Nao, a programmatically error-prone robot. During this phase, they
are presented with a predefined list of 10 general knowledge questions to ask
the Nao. These questions cover topics such as historical dates, geographical
facts, and common knowledge queries such as the independence day of the US,
capital of France, the first US president, shape of a stop sign, and so on.
Nao, equipped with a programmed response set, deliberately answers some
questions incorrectly. Notably, after each incorrect response, Nao is designed
to acknowledge its mistake and issue an apology to the participant. This
interaction phase is pivotal for observing and analyzing the participants’
real-time reactions to robotic errors and apologies, which is crucial for
understanding human-robot interaction dynamics. During this interaction, Nao
is programmed to answer some questions incorrectly and then acknowledge and
apologize for these errors, allowing for observation of participant reactions
to the robot’s mistakes and subsequent apologies.

##  V Results and Discussion

The survey data provides valuable insights into human reactions to robot
interactions, shedding light on several key aspects:

  * •

Perceived Reliability and Usefulness: Participants with limited prior
interaction with robots perceived them as somewhat reliable and useful. The
variation in trust levels among participants indicates a nuanced perspective,
about 58% of participants expressed moderate trust in robots before the
interactions with the robot.

  * •

Reaction to Errors: All of the participants observed errors in robot
responses. Interestingly, about 84% of participants reported feeling more
trusting when the errors were acknowledged and apologized for, highlighting
the potential impact of transparent communication on user trust.

  * •

Willingness to Rely on Robots: Despite noticing errors, many participants
remained willing to rely on robots for future tasks. This suggests a
recognition among users that robots, like humans, possess both strengths and
weaknesses.

  * •

Trust Post-Mistakes: Participants’ trust in robots either remained unchanged
or increased after witnessing mistakes. This finding suggests a growing
acceptance of robots, with many participants open to considering their use in
various domains.

Figure 4: Perception on Robots Before and After the Study

To support the first hypothesis (H1), asserting that ”The extent to which a
robot acknowledges its limitations or errors can significantly enhance human
trust and emotional acceptance” participants were asked to compare the
forgiveness between human and robot. The result shows that humans are more
forgiving of informational errors when they come from familiar sources like
other humans or internet search engines, compared to robots. The data reveals
that more than 25% of participants indicated they would forgive the robot more
than the humans, while approximately 55% were less forgiving to the robot.
Fig. 4 shows the perception on robots, before and after the study. So that we
can say our first hypothesis is supported.

To support the second hypothesis (H2), we compared the trust level on the
robot before and after the study. Experimental result shows that about 32% of
participants had low trust in automation before the experiment but after the
experiment, the trust level increased to around 84% and around 84% of
participants were likely to recommend the use of robots in high-stakes
domains. Fig. 5 shows the visual representation of these statistics. So that
we can say our second hypothesis is supported.

Figure 5: People are more forgiving to human than robot

To provide robust experimental evidence, a one-way ANOVA test [13] was
conducted to analyze the change in participants’ trust levels before and after
the experiment. The results, presented in Fig 6, demonstrate a highly
significant difference in mean trust levels. In this visualization, we observe
that prior to our experiment, the peak density of human trust in the robot
hovered around 0.5, but post-experiment, it noticeably increased to around
1.5. The one-way ANOVA test reveals a significant overall change in trust
levels across experimental groups, supported by a high F-statistic of 150.69
and a corresponding p-value of $3.16\times e^{-13}$.

Figure 6: Density Plot of Trust Levels Before and After the Experiment TABLE I: Contingency Table for Trust Levels Before and After Experiment | Trust/Reliability Level After Experiment
---|---
Trust Level Before Experiment |  Increases reliability |  No change |  Decreases reliability
Completely trusted (a) |  7 |  0 |  0
Moderately trusted (b) |  0 |  11 |  0
Slightly trusted (c) |  0 |  0 |  1

To compare the trust levels before and after the experiment, we created a
contingency table. Table I allows us to visualize the relationship between the
trust level before the experiment (rows) and the trust/reliability level after
the experiment (columns). Based on the table, we can confidently say that the
trust level increases after the experiment, indicating that the participants
perceived the robot as more reliable following the experimental intervention.

The chi-square test has been conducted on the contingency table, yielding the
following results:

  * •

Chi-Square Statistic: 38.0

  * •

p-value: $1.12\times 10^{-7}$

  * •

Degrees of Freedom: 4

  * •

Expected Frequencies:

$\begin{bmatrix}2.579&4.053&0.368\\\ 4.053&6.368&0.579\\\ 0.368&0.579&0.053\\\
\end{bmatrix}$

According to the chi-square statistic it is clear that there is a significant
association between the trust levels before and after the experiment.

This research provide the following additional insights from our participants
into specific aspects of human-robot interactions:

####  V-1 Trust Impact of Unacknowledged Errors

Data indicates that 77% of participants experienced a decrease in trust when
errors were not acknowledged. This finding underscores the importance of
transparent communication in maintaining user trust in robotic systems.

####  V-2 Importance of Acknowledging Errors

All of the participants emphasized the significance of robots admitting
mistakes which means acknowledging errors appears to play a crucial role in
maintaining user trust.

####  V-3 Effeteness of Robot’s Apology

The 87% of participants expressed that after making a mistake, robot’s apology
is very effective to gain the trust.

####  V-4 Error Impact on Trust

Our analysis suggests that, for 62% of the participants, the errors have a
negative impact on trust. However, the extent of this impact varies among
participants, with some individuals displaying resilience to errors.

####  V-5 Perceived Reliability of Robots

Experimental result suggests that, after the interaction with the robot 84% of
participants perceive robots as reliable, aligning with the initial
expectation of participants finding robots somewhat reliable and useful.

####  V-6 Future Interaction Intent

We investigated the participants’ intent for future interactions with robots.
93% of participants expressed willingness to engage with robots in future
tasks, indicating a positive outlook on the integration of robots into various
domains.

##  VI Conclusion and Future Work

In this study, we found that participants’ degrees of confidence in robots
varied, with many exhibiting some degree of trust even in the case of those
who had not before interacted with robots. The experimental results support
both of our hypothesises (H1 and H2). The majority of the participants were
aware of the robot’s shortcomings, and their level of trust in it was
influenced by whether or not the robot admitted its mistakes. Even with fault
knowledge, there was a noticeable willingness to interact with robots in
subsequent assignments. Even when they made mistakes, participants tended to
believe that robots were either equally or even more trustworthy, indicating
that public acceptance of robots is growing. The study demonstrates that
people’s reactions to interacting with robots are nuanced, weighing doubt
against an appreciation of the potential benefits of robotics. In our upcoming
research, we will investigate how people’s long-term trust and dependence on
robots are impacted by the kind and frequency of errors made by them. We’ll
also investigate how people from various demographic groups and cultures
respond to engaging with robots. Lastly, to increase robot dependability and
user trust, we will strive to make them even more adept at identifying and
correcting mistakes. In future, we are planning to do the same experiment on
different types of robots to see how people react on different types of robot
and how different types of robots effect on human trust on automation.

## ACKNOWLEDGMENTS

We appreciate all of the help that was provided by Dr. David Feil-Seifer and
Dr. Monica Nicolescu for sharing their opinions with us. We also appreciate
our participants those who help us by participating in our survey.

## References

  * [1] H. Su, W. Qi, J. Chen, C. Yang, J. Sandoval, and M. A. Laribi, “Recent advancements in multimodal human–robot interaction,” Frontiers in Neurorobotics, vol. 17, 2023.
  * [2] A. M. Aroyo, J. de Bruyne, O. Dheu, E. Fosch-Villaronga, A. Gudkov, H. Hoch, S. Jones, C. Lutz, H. Sætra, M. Solberg, and A. Tamò-Larrieux, “Overtrusting robots: Setting a research agenda to mitigate overtrust in automation,” Paladyn, Journal of Behavioral Robotics, vol. 12, no. 1, pp. 423–436, 2021.
  * [3] N. Mirnig, G. Stollnberger, M. Miksch, S. Stadler, M. Giuliani, and M. Tscheligi, “To err is robot: How humans assess and act toward an erroneous social robot,” Frontiers in Robotics and AI, vol. 4, 2017.
  * [4] A. Weidemann and N. Rußwinkel, “The role of frustration in human–robot interaction – what is needed for a successful collaboration?,” Frontiers in Psychology, vol. 12, 2021.
  * [5] D. Kontogiorgos, M. Tran, J. Gustafson, and M. Soleymani, “A systematic cross-corpus analysis of human reactions to robot conversational failures,” ICMI ’21: Proceedings of the 2021 International Conference on Multimodal Interaction, 2021.
  * [6] H. Yasuda and M. Matsumoto, “Psychological impact on human when a robot makes mistakes,” in Proceedings of the 2013 IEEE/SICE International Symposium on System Integration, pp. 335–339, 2013.
  * [7] M. L. Traeger, “Robots that admit mistakes foster better conversation in humans.” Yale News, March 2020.
  * [8] R. Salek Shahrezaie, B. Akter Anima, and D. Feil-Seifer, “Human and robot bonding: A study on implying homophily in hri,” UNR, 2020.
  * [9] C. J. Hayes, M. Moosaei, and L. D. Riek, “Exploring implicit human responses to robot mistakes in a learning from demonstration task,” in 2016 25th IEEE International Symposium on Robot and Human Interactive Communication (RO-MAN), IEEE, Aug. 2016.
  * [10] M. Salem, G. Lakatos, F. Amirabdollahian, and K. Dautenhahn, “Would you trust a (faulty) robot? effects of error, task type and personality on human-robot cooperation and trust,” in Proceedings of the Tenth Annual ACM/IEEE International Conference on Human-Robot Interaction, HRI ’15, (New York, NY, USA), p. 141–148, Association for Computing Machinery, 2015.
  * [11] M. Stiber and C.-M. Huang, “Not all errors are created equal: Exploring human responses to robot errors with varying severity,” in Companion Publication of the 2020 International Conference on Multimodal Interaction, ICMI ’20 Companion, (New York, NY, USA), p. 97–101, Association for Computing Machinery, 2021.
  * [12] “Nao configuaration,” 2023.
  * [13] R. Bevans, “One-way anova — when and how to use it (with examples),” March 2020.


11institutetext:  Institut für Astronomie und Astrophysik, Universität
Tübingen, Auf der Morgenstelle 10, D-72076, Germany 22institutetext: Astronomy
Unit, Department of Physics and Astronomy, Queen Mary University of London,
London E1 4NS,UK 33institutetext: Leibniz-Institut für Astrophysik Potsdam
(AIP), An der Sternwarte 16, 14482 Potsdam, Germany.
33email: jjoseph@aip.de,lucas.jordan@uni-tuebingen.de,{a.ziampras, g.a.turpin,
<EMAIL_ADDRESS>
# Measuring the numerical viscosity in simulations of protoplanetary disks in
Cartesian grids

The viscously spreading ring revisited

Jibin Joseph Measuring the numerical viscosity in simulations of
protoplanetary disks in Cartesian gridsMeasuring the numerical viscosity in
simulations of protoplanetary disks in Cartesian gridsMeasuring the numerical
viscosity in simulations of protoplanetary disks in Cartesian gridsMeasuring
the numerical viscosity in simulations of protoplanetary disks in Cartesian
grids Alexandros Ziampras Measuring the numerical viscosity in simulations of
protoplanetary disks in Cartesian gridsMeasuring the numerical viscosity in
simulations of protoplanetary disks in Cartesian gridsMeasuring the numerical
viscosity in simulations of protoplanetary disks in Cartesian gridsMeasuring
the numerical viscosity in simulations of protoplanetary disks in Cartesian
grids Lucas Jordan Measuring the numerical viscosity in simulations of
protoplanetary disks in Cartesian gridsMeasuring the numerical viscosity in
simulations of protoplanetary disks in Cartesian grids George A. Turpin
Measuring the numerical viscosity in simulations of protoplanetary disks in
Cartesian gridsMeasuring the numerical viscosity in simulations of
protoplanetary disks in Cartesian grids Richard P. Nelson Measuring the
numerical viscosity in simulations of protoplanetary disks in Cartesian
gridsMeasuring the numerical viscosity in simulations of protoplanetary disks
in Cartesian grids

(Accepted date: 24 July 2023)

###### Abstract

Context.  Hydrodynamical simulations solve the governing equations on a
discrete grid of space and time. This discretization causes numerical
diffusion similar to a physical viscous diffusion, whose magnitude is often
unknown or poorly constrained. With the current trend of simulating accretion
disks with no or very low prescribed physical viscosity, it becomes essential
to understand and quantify this inherent numerical diffusion, in the form of a
numerical viscosity.

Aims.  We study the behavior of the viscous spreading ring and the spiral
instability that develops in it. We then use this setup to quantify the
numerical viscosity in Cartesian grids and study its properties.

Methods.  We simulate the viscous spreading ring and the related instability
on a two-dimensional polar grid using PLUTO as well as FARGO, and ensure the
convergence of our results with a resolution study. We then repeat our models
on a Cartesian grid and measure the numerical viscosity by comparing results
to the known analytical solution, using PLUTO and Athena++.

Results.  We find that the numerical viscosity in a Cartesian grid scales with
resolution as approximately $\nu_{\mathrm{num}}\propto\Delta x^{2}$ and is
equivalent to an effective $\alpha\sim 10^{-4}$ for a common numerical setup.
We also show that the spiral instability manifests as a single leading spiral
throughout the whole domain on polar grids. This is contrary to previous
results and indicates that sufficient resolution is necessary in order to
correctly resolve the instability.

Conclusions.  Our results are relevant in the context of models where the
origin should be included in the computational domain, or when polar grids
cannot be used. Examples of such cases include models of disk accretion onto a
central binary and inherently Cartesian codes.

###### Key Words.:

Hydrodynamics – accretion, accretion disks – protoplanetary disks – Methods:
numerical

##  1 Introduction

Hydrodynamical simulations are a useful tool in studying a variety of
astrophysical processes. An example of such a process is the accretion flow of
a protoplanetary disk around a star. Accretion is thought to be achieved via
turbulence operating in parts of the disk (Lyra & Umurhan, 2019), that
transports angular momentum radially outwards similar to a viscosity (Lynden-
Bell & Pringle, 1974; Balbus & Papaloizou, 1999), or via a magnetic torque
that is applied on the material close to the surface of the disk as a stellar
wind expels gas from the disk surface layers (Bai & Stone, 2013) with bulk of
the disk remaining inviscid and laminar. The result, particularly in
protoplanetary disks, is a steady radial infall of gas that ultimately
depletes the disk over typical timescales of 1–10 Myr Haisch et al. (2001).

While turbulent angular momentum transport has been the traditional way of
modeling accretion, numerous recent observations of disks around T Tauri stars
(e.g., the DSHARP survey, Andrews et al., 2018) suggest very low
$\alpha\lesssim 10^{-4}$ in order for turbulence to be compatible with the
radial width of observed rings (Dullemond et al., 2018), the vertical
structure of mm grains (Dullemond et al., 2022) and the formation of rings and
gaps by embedded planets (Zhang et al., 2018). As a result, magnetic winds are
becoming a favored means of interpreting accretion in protoplanetary disks,
and numerical models of such disks tend towards the inviscid limit (e.g., Lega
et al., 2022).

Here, however, we run into a different problem: the numerical schemes used to
model protoplanetary disks introduce a certain amount of numerical diffusion.
The magnitude of this diffusion, is often unknown or poorly constrained. Thus,
calculating an upper limit to this non-physical diffusion, or equivalently a
”numerical viscosity”, is essential in ensuring that results are not affected
by the effects it can induce. This is particularly important for models using
Cartesian grids, which are primarily used when the central object should be
included in the simulation domain such as when modeling accretion patterns
around binary stars (e.g. Tiede et al., 2021), or for certain MHD codes (e.g.,
Fromang et al., 2006). Analogue issues also appear when an object is not
centered on a polar grid, which is often the case for circumplanetary disks
(e.g., Crida et al., 2009). Grid-noise due to the asymmetric nature of
Voronoi-mesh cells has also been identified as a source of numerical diffusion
in moving mesh codes (Zier & Springel, 2022). In these cases, the geometry of
the grid introduces a very high numerical viscosity and as a result models
require the execution of very computationally expensive, high-resolution
simulations.

In this work, we revisited the viscously spreading ring problem introduced
first by Lüst (1952) and again by Lynden-Bell & Pringle (1974), and analyzed
in 2D by Speith & Kley (2003) (henceforth SK03). We approached this with
numerical hydrodynamics simulations that first aimed to reproduce or improve
upon the results of SK03 by reanalyzing the azimuthal instability they
reported, and then performed calculations in Cartesian coordinates to quantify
the numerical viscosity of the codes PLUTO and Athena++.

The viscous ring problem and our physical and numerical setup are described in
Sect. 2. We compare our results with those of SK03 and point out the origin of
a numerical artifact in their study through a resolution analysis in Sect. 3.
We present our results for Cartesian models and our estimates of the numerical
viscosity in Sect. 4. We then discuss our findings in Sect. 5, and conclude in
Sect. 6.

##  2 Model setup

In this section, we revisit the setup for the viscously spreading ring problem
and describe our physics and numerical methods. We also provide a list of
models that were executed using different codes and in different coordinate
systems.

###  2.1 The viscous ring problem

Starting from an axisymmetric surface density distribution of a pressureless
($P=0$) ring of gas with mass $M_{\text{ring}}$ that is centered at distance
$R_{0}$ and is infinitely thin in the radial direction

$\Sigma_{\mathrm{ring}}(t=0,R)=\frac{M_{\mathrm{ring}}}{2\pi
R}\delta(R-R_{0}),$ (1)

the ring will spread along the radial direction $R$ following the analytical
solution (e.g., Lynden-Bell & Pringle, 1974; Pringle, 1981)

$\Sigma_{\mathrm{ring}}(\tau,x)=\frac{M_{\mathrm{ring}}}{\pi
R_{0}^{2}}\frac{1}{\tau
x^{1/4}}\,I_{\frac{1}{4}}\left(\frac{2x}{\tau}\right)\mathrm{exp}\left(-\frac{1+x^{2}}{\tau}\right),$
(2)

where $x=R/R_{0}$ and $\tau=12\nu tR_{0}^{-2}$ are normalized distance and
time quantities, $\nu$ is the kinematic viscosity of the gas, and
$I_{\frac{1}{4}}$ is the modified Bessel function of the order $1/4$. In
addition, the radial velocity of the ring is given by the following relation,

$u_{\mathrm{R,ring}}=-\frac{3}{\Sigma_{\mathrm{ring}}\,\sqrt{R}}\frac{\partial}{\partial
R}\left[\nu\Sigma_{\mathrm{ring}}\sqrt{R}\right]$ (3)

We note that, Eq. (2) is only an approximate solution to the hydrodynamic
equations (see Sect. 2.2), and assumes that the disk is massless with highly
supersonic and Keplerian azimuthal velocity. Furthermore, one assumes that the
kinematic viscosity, $\nu$, is much smaller than the specific angular momentum
$R^{2}\Omega$, an assumption that is violated close to the central object. A
time evolution of the viscous ring problem using our fiducial axisymmetric
numerical setup (described in Sect. 2.3) is shown in Fig. 1.

The problem was revisited by SK03, who showed that the ring is subject to a
viscosity-driven spiral instability in the azimuthal direction that results in
the development of a leading spiral arm that covers the full extent of the
computational domain (see Fig. 3). They performed a detailed stability
analysis and showed that unstable modes with wave-numbers $k$ develop where

$k^{2}>\frac{3}{R\Sigma_{0}}\frac{\partial\Sigma_{0}}{\partial R}$ (4)

with $\Sigma_{0}$ being the axisymmetric density profile. They further showed
that the spiral arm changes direction from leading to trailing at around the
peak of the ring (see second and third panels of Fig. 14), and verified their
results with two fundamentally different codes: the finite-difference upwind
grid code RH2D (Kley, 1989, 1999), and the smooth particle hydrodynamics (SPH)
code used by Flebbe et al. (1994).

###  2.2 Hydrodynamics

We consider the vertically-integrated hydrodynamics equations for a gas with
surface density $\Sigma$, velocity vector $\bm{u}$, and nearly zero pressure
$P$ at distance $R$,

$\begin{split}&\frac{\partial{\Sigma}}{\partial{t}}+\nabla\cdot(\Sigma\bm{u})=0,\\\
&\frac{\partial{(\Sigma\bm{u})}}{\partial{t}}+\nabla\cdot\left(\Sigma\bm{u}\otimes\bm{u}\right)=-\nabla
P-\Sigma\nabla\Phi_{\star}+\nabla\cdot\overline{\bm{{\upsigma}}}.\end{split}$
(5)

The gas orbits around a star with mass $M_{\star}$ such that the gravitational
potential is given by $\Phi_{\star}=-\text{G}M_{\star}/R$, with G being the
gravitational constant. The viscous stress tensor is represented with
$\overline{\bm{{\upsigma}}}$. We use a locally isothermal equation of state,
such that the gas has a (very small) sound speed
$c_{\text{s}}=\sqrt{P/\Sigma}=h\sqrt{\text{G}M_{\star}/R}=hR\Omega_{\text{K}}$,
where $h$ is the aspect ratio.

###  2.3 Numerical setup

For our standard setup, we used PLUTO version 4.3 (Mignone et al., 2007), a
finite-volume code with a second-order accurate scheme (HLLC, Toro et al.,
1994) and the van Leer flux limiter (van Leer, 1977). We followed the setup
laid out in SK03 to reproduce the original results. In this section, we
provide a detailed description of our model and its differences with the setup
in SK03.

To study the viscous ring and the related instability we used a polar grid
with computational domain $R\in[0.2,2]\,R_{0}$ and $\phi\in[0,2\pi)$. We used
both logarithmic ($\Delta R\propto R$) and arithmetic ($\Delta
R=\mathrm{const.}$) scaling in the radial direction for comparison reasons,
but focus on logarithmic setups first. We also enabled the FARGO transport
scheme (Masset, 2000) implemented in PLUTO by Mignone et al. (2012), while the
parabolic viscosity term is handled using the explicit time stepping method in
PLUTO. We used a Courant number of 0.4 in all simulations, unless specified
otherwise. The system was then evolved for various grid resolutions, the
results of which are discussed in Sect. 3.2. The Cartesian setup used to
measure the numerical viscosity is described in Sect. 4.2.

To evolve the viscous ring problem, we initialized the surface density using
Eq. (2) at $\tau_{0}=0.018$ and added a small, constant surface density floor
$\Sigma_{\mathrm{floor}}=10^{-7}\,\Sigma_{\mathrm{ref}}$, where
$\Sigma_{\mathrm{ref}}=\Sigma_{\mathrm{ring}}(\tau_{0},1)$. We chose a set of
code units where $\bar{\mathrm{G}}=\bar{M}_{\star}=\bar{R}_{0}=1$. This leaves
the kinematic viscosity as the only relevant physical constant in this setup,
which we chose to be $\bar{\nu}=10^{-5}$. For this viscosity, we get a viscous
ring spreading time (i.e., $\tau=1$) of $1326$ orbital periods ($P_{0}$) at
$R=R_{0}$. To facilitate development of the spiral instability, the initial
surface density distribution was seeded with a small amount of noise and
evolved up to $t=400\leavevmode\nobreak\ P_{0}$, or $\tau=0.3$.

We chose a strict outflow boundary condition in both radial directions such
that $\partial_{R}\Sigma=0$, $u_{R,\text{in}}=-|u_{R}(R_{\text{in}})|$,
$u_{R,\text{out}}=|u_{R}(R_{\text{out}})|$. The azimuthal velocity is fixed to
the Keplerian speed corrected for pressure support
$u_{\phi}=R\Omega_{\text{K}}\sqrt{1-h^{2}}$ at the boundaries. We explored
various values of the aspect ratio, shown in Appendix A, and chose $h=0.005$
for our simulations, as we found it to be low enough for the gas to act like a
pressureless fluid and high enough to ensure numerical stability in our
simulations. The RH2D setup in SK03 used a viscosity value of
$\bar{\nu}=4.77\times 10^{-5}$ and an initial viscous time $\tau_{0}=0.016$,
but this has no effect on the subsequent evolution.

To corroborate our results, we used the FARGO code (legacy version) (Masset,
2000). Similar to the RH2D code used in SK03, FARGO uses a finite difference
scheme as described in Stone & Norman (1992) with the second-order upwind
algorithm by van Leer (1977). Both codes therefore have similar behavior and
we used FARGO to compare our results to SK03.

##  3 Analysis of the viscous ring problem

In this section, we explore the spiral instability using logarithmically
spaced polar grids. In such grids, the radial cell size $\Delta R$ increases
with distance such that the ratio $\Delta R/H$ remains comparable throughout
the domain (in the case of $h=\mathrm{const.}$, this ratio is constant as
$H=hR$). This makes logarithmically spaced grids well-suited for general-
purpose disk models, where gas dynamics should be resolved over a scale height
$H$ using a reasonable number of cells. We briefly describe the general
behavior of the viscous spreading ring and show that our fiducial model is
numerically converged.

###  3.1 Time evolution of the viscous ring

We begin our analysis with a 1D model of the viscously spreading ring, using
465 cells logarithmic scaled in the radial direction. This resolution is
chosen such that the gas is resolved with at least one cell per scale height
$H=hR$. Even though the problem is formally defined for a pressureless fluid
($h=0$), Fig. 1 and Fig. 2 show that our model adequately reproduces the
analytical equations of Eq. (2) and Eq. (3), respectively.

Figure 1: Evolution of the one-dimensional viscously spreading ring. The
colors indicate different timestamps, solid lines represent simulation results
and dots follow the analytical solution at the corresponding time according to
Eq. (2). Figure 2: Similar to Fig. 1 showing the evolution of radial velocity
of the ring compared to the analytical solution Eq. (3).

To assess the robustness of the numerical methods used, we quantified
numerical errors in our results by measuring the average density-weighted
relative deviation from the analytical solution of both the surface density
$\Sigma$ and radial velocity $u_{R}$ on the grid, at various times. We find a
maximum deviation of $0.74\%$ for $\Sigma$ and $3.31\%$ for $u_{R}$. While
these values are largely influenced by regions close to the boundaries, where
deviations are the largest, for most of the domain (0.4–1.8 $R_{0}$) we find
maximum deviations of $0.15\%$ and $0.75\%$, for $\Sigma$ and $u_{R}$
respectively. In addition, we find that further reducing the aspect ratio $h$
does not substantially affect the measured errors.

We then expanded to a 2D $\\{R,\phi\\}$ domain, with a fiducial resolution of
$N_{R}\times N_{\phi}=465\times 256$. Here, we expect the same radial
spreading of the ring but also development of the viscosity-driven spiral
instability discussed in SK03. Figure 3 shows a time evolution of the ring
spreading and of the instability which manifests in the form of an outward-
propagating spiral arm. Weak radial wave-like perturbations are launched
immediately after the start of the simulation (see panel (a)).

The perturbations are not affected by changes to pressure or initial
conditions, but they do become weaker and move to smaller wavelengths with
lower viscosity. This viscosity dependency is in agreement with the viscous
overstability (e.g., Latter & Ogilvie, 2006) but a more thorough study of the
perturbations is not within the scope of this work.

These initial waves leave the domain through the outer boundary while the
spiral instability begins to develop near the inner boundary, with the spiral
becoming apparent at $\tau\approx 0.03$ (panel (b)). The spiral continues to
grow in amplitude and propagate outwards until it eventually spans the entire
radial extent of the disk (panels (c), (d)).

In contrast to the results in SK03 we do not see a flip of the direction of
the spiral from leading to trailing around the peak of the ring. Instead, we
find that this flip is merely a numerical artifact and does not appear in
simulations with a high enough radial resolution. We discuss this further in
Appendix B.

Figure 3: Snapshots of the normalized deviation from the analytical surface
density distribution for our fiducial model at different times. The dashed
green line marks the position of maximal surface density that indicates the
peak of the spreading ring. Initial waves (a) are seen as outward-moving rings
while the spiral instability develops as a spiral arm at the inner boundary
(b) and spreads through the whole domain (c,d). Figure 4: (a): Time evolution
of the normalized maximal deviation from the azimuthally averaged surface
density distribution for various numerical configurations. Inadequate
resolution in the $R$-direction can damp the viscous instability, with
convergence being achieved for our fiducial model with $N_{R}=465$. (b):
Comparison of the same quantity for different azimuthal resolutions and
$N_{R}=465$. The higher resolution of $N_{\phi}=1300$ azimuthal cells results
in square grid cells but does not influence the growth of the instability.

###  3.2 Resolution study

To ensure that we have properly resolved and developed the viscous
instability, we performed a resolution study by analyzing its growth phase.
For this, we measured the largest density deviation from the azimuthally
averaged profile at different radii and tracked this quantity over time. In
that sense, this quantity acts as a proxy for the time evolution of the
amplitude of the strongest spiral at a given radius, which hints at the growth
rate of the spiral instability. Figure 4 shows this evolution for various
resolutions. Alternatively, we also evaluated the amplitude of the azimuthal
Fourier modes of density rings at specific radii over time. This method yields
similar results to those obtained from measuring the density deviations and
are therefore not shown here.

After comparing our results against a model with square cells ($N_{R}\times
N_{\phi}=465\times 1300$), we found that increasing the azimuthal resolution
beyond 256 cells has no effect on the growth rate of the instability (see
panel (b) in Fig. 4) and therefore fixed the number of azimuthal cells to 256
in our resolution study. In the radial direction we found a resolution of 465
cells to be sufficient.

We confirmed our results using the FARGO code in Appendix B and achieved
numerical convergence with both codes at a resolution of $465\times 256$
cells. We conclude that this resolution is sufficient to fully resolve the
viscously spreading ring and the development and saturation of spiral arms due
to the spiral instability.

##  4 The viscous ring problem in Cartesian coordinates

Having profiled the viscous ring problem in a polar grid, we now switch to
Cartesian coordinates with our ultimate goal being the measurement of the
numerical viscosity of PLUTO for such a grid. For this set of models our
computational domain extends between $x,y\in[-2,2]\,R_{0}$, with a fiducial
resolution of $N_{x}\times N_{y}=1024\times 1024$ cells which corresponds to
approximately one cell per scale height at $R=\sqrt{x^{2}+y^{2}}=R_{0}$. We
note that, $\Delta x=\Delta y=\mathrm{const.}$ in these models, and therefore
our effective resolution is significantly lower at smaller radii.

For a fair comparison with our models on a polar grid, we retained the second-
order accurate HLLC scheme with a linear spatial reconstruction and van Leer
flux limiter for all simulations unless stated otherwise. In addition, we
emulated the inner radial boundary of a polar grid using a damping region,
where we damped $\Sigma$ to $\Sigma_{\mathrm{ring}}$ through Eq. (2) and
$u_{\mathrm{R}}$ to zero for $R<0.2\,R_{0}$. Finally, we set an outflow
boundary condition at all domain boundaries. To measure numerical viscosity in
PLUTO we used a similar setup but damped $\Sigma$ to $\Sigma_{\mathrm{floor}}$
following (de Val-Borro et al., 2006), and provide a comparison using Athena++
version 21.0 (Stone et al., 2008).

###  4.1 Viscous evolution

We first analyze the evolution of the ring with standard viscosity similar to
Sec. 3. The radial ring spreading and the radial velocity evolution can be
seen in Fig. 5 and Fig. 6 respectively. They match the analytical solution
well except for noticeable deviations in $u_{\mathrm{R}}$ close to the damping
region and boundaries. We again quantified the error by measuring the density-
weighted relative deviations in the domain excluding the damping region and
found a maximum deviation of $0.71\%$ for $\Sigma$ and of $15.74\%$ for
$u_{\mathrm{R,ring}}$. We note that, unlike the polar case, we find no
development of spirals in any Cartesian model. We suspect that the spiral
instability is suppressed as a result of the highly diffusive grid, further
discussed in Sec. 5. We see similar results at a lower and higher resolution
of $512^{2}$ and $2048^{2}$ cells, respectively.

Figure 5: Snapshots of slices through the $y=0$ plane showing the evolution of
surface density for our fiducial viscous Cartesian simulation. The colors
indicate different timestamps, solid lines represent simulation results and
dots follow the analytical solution at the corresponding time according to Eq.
(2). Unlike the 2D polar case, no spirals develop. Figure 6: Similar to Fig. 5
showing the evolution of radial velocity of the ring compared to the
analytical solution Eq. (3).

###  4.2 Inviscid models: estimates of numerical viscosity

The discretization in space and time inherent in the numerical schemes
employed by our codes results in numerical errors when solving the
hydrodynamics equations. In our analysis in Appendix C, we show that this
error has the form of a numerical viscosity and can lead to ring spreading,
even in the absence of any physical viscosity. In this section, we repeat the
models shown in Sect. 4.1 with an inviscid prescription ($\nu=0$), and analyze
the subsequent ring spreading in an attempt to quantify the numerical
viscosity of the methods in PLUTO.

Blue lines in Fig. 7 show the time evolution for our fiducial inviscid model
at a resolution of $1024\times 1024$. We found that the ring spreading is
indeed quite similar to viscous models, with the exception that the peak of
the ring is flatter compared to the analytical solution. We then extracted a
global estimate of the numerical viscosity $\nu_{\mathrm{num}}$ by fitting the
analytical solution Eq. (2) to our data. The fitting method is described in
Appendix D.

We repeated the simulation with a prescribed kinematic viscosity equal to the
measured numerical viscosity ($\nu=\nu_{\mathrm{num}}$) and once again with
$\nu=2\,\nu_{\mathrm{num}}$ , the time evolution of which is shown as orange
and green lines respectively in Fig. 7. We found that the
$\nu=\nu_{\mathrm{num}}$ model evolves twice as fast as our inviscid model,
and the $\nu=2\,\nu_{\mathrm{num}}$ model evolves thrice as fast. This
suggests that all models evolve as if the total viscosity were
$\nu_{\text{tot}}=\nu+\nu_{\mathrm{num}}$. We motivate our approach and
rationalize this result in Appendix C.

Figure 7: Slices of surface density profiles for our fiducial inviscid
Cartesian model ($\nu=0$, blue lines) at a given time $t_{0}$ against the
corresponding viscous models with $\nu=\nu_{\mathrm{num}}$ at $t_{0}/2$
(orange) and $\nu=2\,\nu_{\mathrm{num}}$ at $t_{0}/3$ (green). Solid and
dashed lines correspond to $t_{0}=224$ and 748 orbits at $R_{0}$,
respectively. The model with $\nu=\nu_{\mathrm{num}}$ evolves twice as fast as
the inviscid simulation and the $\nu=2\,\nu_{\mathrm{num}}$ model evolves
thrice as fast, indicating a viscous contribution by the non-negligible
numerical viscosity.

####  4.2.1 Resolution scaling of the numerical viscosity

In our analysis of the numerical viscosity for a first-order solver in
Appendix C, we show that the numerical viscosity scales as
$\nu_{\mathrm{num}}\propto\Delta x$ for a first-order method. As we used a
second-order solver for our tests, we expect the numerical viscosity to scale
as $\nu_{\mathrm{num}}\propto\Delta x^{2}$ in our models. To test this, we
conducted simulations for grid resolutions of $256^{2}$, $512^{2}$,
$1024^{2}$, $2048^{2}$ and $4096^{2}$ cells. The resulting values of
$\nu_{\mathrm{num}}$ as a function of cell count are listed in Table
LABEL:tab:numerical_viscosities and shown in Fig. 8. We then fit a power law
to the relation $\nu_{\mathrm{num}}(\Delta x)$ and found that the numerical
viscosity scales as $\nu_{\mathrm{num}}\propto\Delta x^{2.09}$. We also found
that reducing the maximum Courant number $C_{\mathrm{max}}\approx\Delta
t/\min\left(\nicefrac{{\Delta x}}{{u}}\right)$ from our nominal value of 0.4,
had no effect on the numerical viscosity, indicating that the dependence of
$\nu_{\mathrm{num}}$ on $\Delta x$, $\Delta t$ and the Courant number is more
complicated than it would be for a first-order scheme. While we find our
$\nu_{\mathrm{num}}\propto\Delta x^{2}$ estimate to roughly describe the
convergence, we also find that the convergence rate is not constant and
increases with resolution. We found similar results using the code Athena++
(see Appendix E).

Table 1: Numerical viscosity $\nu_{\mathrm{num}}$ with one standard deviation errors in code units for different grid types and resolutions. An equivalent $\alpha_{\text{num}}$ at $R_{0}$ and a realistic $h=0.05$ is also listed. Code | Grid | Resolution | $\bar{\nu}_{\mathrm{num}}$ |  $\alpha_{\text{num}}$ at $R_{0}$
---|---|---|---|---
PLUTO | Cartesian | $256\times 256$ | $(4.68\pm 0.19)\times 10^{-6}$ | $(1.87\pm 0.08)\times 10^{-3}$
PLUTO | Cartesian | $512\times 512$ | $(1.48\pm 0.01)\times 10^{-6}$ | $(5.92\pm 0.04)\times 10^{-4}$
PLUTO | Cartesian | $1024\times 1024$ | $(4.16\pm 0.03)\times 10^{-7}$ | $(1.67\pm 0.01)\times 10^{-4}$
PLUTO | Cartesian | $2048\times 2048$ | $(9.45\pm 0.11)\times 10^{-8}$ | $(3.78\pm 0.04)\times 10^{-5}$
PLUTO | Cartesian | $4096\times 4096$ | $(1.32\pm 0.08)\times 10^{-8}$ | $(5.29\pm 0.31)\times 10^{-6}$
PLUTO | polar | $465\times 256$ | $(3.89\pm 0.94)\times 10^{-15}$ | $(1.56\pm 0.38)\times 10^{-12}$
Athena++ | Cartesian | $256\times 256$ | $(7.31\pm 0.22)\times 10^{-6}$ | $(2.92\pm 0.08)\times 10^{-3}$
Athena++ | Cartesian | $512\times 512$ | $(2.04\pm 0.05)\times 10^{-6}$ | $(8.16\pm 0.19)\times 10^{-4}$
Athena++ | Cartesian | $1024\times 1024$ | $(4.26\pm 0.11)\times 10^{-7}$ | $(1.70\pm 0.04)\times 10^{-4}$
Athena++ | Cartesian | $2048\times 2048$ | $(1.11\pm 0.02)\times 10^{-7}$ | $(4.44\pm 0.08)\times 10^{-5}$
Figure 8: Numerical viscosity $\bar{\nu}_{\mathrm{num}}$ in code units and the
corresponding $\alpha$ at $R_{0}$ for $h=0.05$ at different resolutions,
extracted by fitting Eq. (2) to our Cartesian grid results for both Athena++
and PLUTO. The two codes agree very well for $N_{\textrm{cells}}\geq
1024^{2}$. See Appendix E for details on the Athena++ models.

For a realistic protoplanetary disk with $h\sim 0.05$, the numerical viscosity
$\bar{\nu}_{\mathrm{num}}=4.16\times 10^{-7}\sqrt{\mathrm{G}M_{0}R_{0}}$
extracted from our fiducial Cartesian model translates to an $\alpha$
viscosity parameter (Shakura & Sunyaev, 1973) of $\alpha\sim 1.67\times
10^{-4}/\sqrt{\bar{R}}$. More generally, we can write

$\alpha=\frac{\nu}{c_{\text{s}}H}\approx
400\times\bar{\nu}\,\left(\frac{0.05}{h}\right)^{2}\,\sqrt{\frac{R_{0}}{R}},$
(6)

##  5 Discussion

In this section, we compare our results on the spiral instability to those
found in SK03, and comment on the nature of numerical viscosity.

###  5.1 Comparison to previous results

Analyzing the spreading ring on a polar grid in Sect. 3, we could confirm the
main findings from SK03 that the viscous pressureless spreading ring is
subject to a viscous instability that manifests as a leading spiral arm. On
the other hand, we found their described flip from a leading to a trailing
spiral to be purely of numerical origin. As shown in Appendix B, we attribute
the spiral flip to an under-resolved inner disk; the quantity $\Delta R/R$
increases for small radii in an arithmetic grid like that of SK03, while a
constant smoothing length in their SPH simulations results in effectively the
same issue. In our adequately resolved models ($N_{R}\geq 465$), we find many
azimuthal spiral modes growing at the beginning of our simulation with their
amplitude decaying exponentially with mode number. In our arithmetic, low-
resolution analog of SK03, however, the growth of azimuthal Fourier mode
amplitudes is delayed, weaker and their amplitude decays faster with
increasing mode number (Fig. 13). We provide more details in Appendix B, but
conclude here that any effects beside the leading spiral found in SK03 was due
to their poor numerical resolution.

###  5.2 Differences between physical and numerical viscosity

In Sect. 3.1 we showed the development of the spiral instability for our
fiducial polar model (Fig. 3) and its growth (Fig. 4). Unlike the polar grid,
our Cartesian models did not develop this instability, even at comparatively
higher resolutions. This can be attributed to the lack of angular momentum
conservation, which results in a dramatically high numerical diffusivity. We
confirm this by analyzing the growth rates of the instability just like for
the polar case in Fig. 9, and show that azimuthal wavenumbers in Cartesian
runs saturate at levels that are six orders of magnitude weaker compared to
polar runs.

Furthermore, we interpolated our fiducial polar run with a fully developed
spiral onto a Cartesian grid and noticed that the spiral features of the
viscous instability completely vanished within a span of 2–4 orbits after
continuation. This suggests that even though numerical diffusion causes an
effect akin to physical viscosity, the properties of numerical and physical
viscosity are not the same. Both physical and numerical viscosity lead to a
viscous spreading of the initial ring distribution, whereas a physical
viscosity is necessary for the development of the spiral instability. The same
effect cannot be replicated by numerical diffusion alone.

Figure 9: Fourier amplitude normalized to azimuthally averaged density at
$R=0.5\,R_{0}$ for $m=1,2,3$ for the fiducial polar simulation (dashed lines)
and interpolated fiducial Cartesian simulation (solid lines). The spiral
instability does not develop, with amplitudes up to six orders of magnitude
weaker compared to the polar case.

###  5.3 Effect of characteristic limiting on numerical viscosity

A commonly used method to reduce diffusive numerical effects is characteristic
limiting (henceforth “CL”, implemented by Mignone et al., 2012), which
performs spatial reconstruction on characteristic variables instead of the
primitive variables in the system. We assessed the effect of this method by
comparing to our inviscid $512\times 512$ Cartesian simulation with otherwise
identical parameters. The low resolution was chosen to highlight the
substantial difference while using CL. As shown in Fig. 10, the simulation
with CL is less diffusive, evident by the lack of a flattened peak and an
overall less-diffused ring. On average, the model with CL had roughly half the
numerical viscosity as the standard inviscid $512\times 512$ model.

Figure 10: Inviscid Cartesian $512^{2}$ simulation utilizing characteristic
limiting (dashed line) compared to out standard inviscid models at resolution
$512^{2}$ and $1024^{2}$ (solid lines), after 150 orbits at $R_{0}$.
Characteristic limiting reduces numerical viscosity roughly by a factor of 2
at $512^{2}$ resolution. It also better maintains local features of the ring,
especially the density peak which is rather flattened by numerical viscosity
in the compared models. A dotted black line shows the fit using the analytical
ring profile for the $512^{2}$ simulation.

##  6 Conclusions

We revisited the viscous ring problem in SK03 using the finite-volume codes
PLUTO and Athena++ as well as the finite-difference code FARGO, in an attempt
to measure the numerical viscosity of hydrodynamical codes.

We first reproduced the viscous ring spreading in one dimension and the spiral
instability in two dimensions on a polar grid. We analyzed the growth of the
instability and showed that SK03 had insufficient spatial resolution to
properly resolve the instability growth.

We then evolved the viscous ring spreading in two dimensions on a Cartesian
grid, applying the same viscosity $\nu$ as used in the polar grid run. The
evolution matched the analytical solution to good accuracy but failed to
develop the spiral instability even at higher resolutions. We attributed the
absence of spiral development to the lack of angular momentum conservation and
subsequently high numerical diffusion inherent to Cartesian grids. Our
findings suggest that a high physical viscosity is not the only ingredient to
developing the spiral instability, but also a numerical setup with very good
angular momentum conservation and low numerical diffusivity.

The viscous spreading ring was then used to measure the numerical viscosity in
Cartesian grids. This was done by evolving the viscous ring for an inviscid
setup and fitting the simulated ring evolution with the analytical solution in
Eq. (2) as a function of time. As the ring evolution depends only on the
viscosity $\nu$, this method allowed us to extract the viscosity over which
the system evolved, which for an inviscid system is the numerical viscosity.
We then showed a scaling relation between numerical viscosity and resolution
that, for second-order methods, corresponds to
$\bar{\nu}_{\mathrm{num}}\approx 0.63\,\Delta x^{2}$. Translating our results
to a Shakura–Sunyaev $\alpha$ parameter, we found a relation $\alpha\approx
2\times 10^{-4}/\sqrt{\bar{R}}$ for our fiducial model with
$\bar{\nu}_{\mathrm{num}}\approx\ 4\times 10^{-7}$ and a realistic aspect
ratio $h=0.05$.

We highlight that our models all utilize second-order accurate schemes in both
space and time. Even though the effects of numerical viscosity can be
mitigated to some degree by using higher order spatial reconstruction and time
marching algorithms, further study needs to be done on that matter.

Our results show the existence of moderate diffusive effects in Cartesian
grids and quantify the resulting numerical viscosity for standard numerical
parameters and different grid resolutions. We also lay out a method that can
be used to quantify the numerical viscosity to good accuracy. This information
can be utilized to make informed decisions on how to measure and minimize
numerical diffusion in hydrodynamics simulations of accretion disks, and is
especially useful in the context of low-viscosity or even inviscid
configurations of circumbinary or protoplanetary disks on Cartesian grids, as
well as for inherently Cartesian codes.

###### Acknowledgements.

We dedicate this paper to the memory of our dear friend and mentor Willy Kley,
with whom we developed the idea for this paper. We thank him for his
fundamental work, counseling, and kindness. AZ and RPN are supported by STFC
grant ST/P000592/1, and RPN is supported by the Leverhulme Trust through grant
RPG-2018-418. GAT is supported by an STFC PhD studentship. JJ is co-funded by
the European Union (ERC, EPOCH-OF-TAURUS, 101043302). Views and opinions
expressed are however those of the author(s) only and do not necessarily
reflect those of the European Union or the European Research Council. Neither
the European Union nor the granting authority can be held responsible for
them. The authors acknowledge support by the High Performance and Cloud
Computing Group at the Zentrum für Datenverarbeitung of the University of
Tübingen, the state of Baden-Württemberg through bwHPC and the German Research
Foundation (DFG) through grant no INST 37/935-1 FUGG. This research utilised
Queen Mary’s Apocrita HPC facility, supported by QMUL Research-IT
(http://doi.org/10.5281/zenodo.438045). This work was performed using the
DiRAC Data Intensive service at Leicester, operated by the University of
Leicester IT Services, which forms part of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC
capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant
ST/R001014/1. DiRAC is part of the National e-Infrastructure. The plots in
this paper were prepared using the Python library matplotlib (Hunter, 2007).

## References

  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41
  * Bai & Stone (2013) Bai, X.-N. & Stone, J. M. 2013, ApJ, 769, 76
  * Balbus & Papaloizou (1999) Balbus, S. A. & Papaloizou, J. C. B. 1999, ApJ, 521, 650
  * Courant et al. (1952) Courant, R., Isaacson, E., & Rees, M. 1952, Communications on pure and applied mathematics, 5, 243
  * Crida et al. (2009) Crida, A., Baruteau, C., Kley, W., & Masset, F. 2009, A&A, 502, 679
  * de Val-Borro et al. (2006) de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529
  * Dullemond et al. (2018) Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46
  * Dullemond et al. (2022) Dullemond, C. P., Ziampras, A., Ostertag, D., & Dominik, C. 2022, A&A, 668, A105
  * Flebbe et al. (1994) Flebbe, O., Muenzel, S., Herold, H., Riffert, H., & Ruder, H. 1994, ApJ, 431, 754
  * Fromang et al. (2006) Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371
  * Haisch et al. (2001) Haisch, Karl E., J., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Kley (1989) Kley, W. 1989, A&A, 208, 98
  * Kley (1999) Kley, W. 1999, MNRAS, 303, 696
  * Latter & Ogilvie (2006) Latter, H. N. & Ogilvie, G. I. 2006, MNRAS, 372, 1829
  * Lega et al. (2022) Lega, E., Morbidelli, A., Nelson, R. P., et al. 2022, A&A, 658, A32
  * Lüst (1952) Lüst, R. 1952, Zeitschrift Naturforschung Teil A, 7, 87
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603
  * Lyra & Umurhan (2019) Lyra, W. & Umurhan, O. M. 2019, PASP, 131, 072001
  * Masset (2000) Masset, F. 2000, A&AS, 141, 165
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Mignone et al. (2012) Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7
  * Pringle (1981) Pringle, J. E. 1981, ARA&A, 19, 137
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  * Speith & Kley (2003) Speith, R. & Kley, W. 2003, A&A, 399, 395
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  * Stone & Norman (1992) Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753
  * Tiede et al. (2021) Tiede, C., Zrake, J., MacFadyen, A., & Haiman, Z. 2021, arXiv e-prints, arXiv:2111.04721
  * Toro et al. (1994) Toro, E. F., Spruce, M., & Speares, W. 1994, Shock waves, 4, 25
  * van Leer (1977) van Leer, B. 1977, Journal of Computational Physics, 23, 276
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47
  * Zier & Springel (2022) Zier, O. & Springel, V. 2022, arXiv e-prints, arXiv:2205.07916

##  Appendix A Effect of pressure scale height on the viscous ring spreading
and spiral instability

The analytical solution for the viscous ring spreading and the the stability
analysis of the spiral instability studied in SK03 assumes a disk without
pressure. For a locally isothermal system, this would imply an aspect ratio of
$h=0$ but the PLUTO solver becomes numerically unstable under this condition.
To ensure that the pressure in our simulations is small enough to properly
simulate the spiral instability we compared the maximal spiral amplitude time
averaged over $\tau=0.2$-0.3 for different values of aspect ratio $h$ in Fig.
11.

For aspect ratios $h>0.006$, the spiral instability is strongly damped and the
ring structure deviates substantially from the analytical solution. For aspect
ratios $h<0.006$, the instability becomes stronger, and is fully saturated at
our fiducial aspect ratio of $h=0.005$. We also tested and confirmed (not
shown) that the numerical viscosity measured in our Cartesian runs is fully
converged for $h=0.005$. Nevertheless, our fiducial aspect ratio is small
enough for the gas to be considered pressureless, and we found that the noise
we observed in the amplitude of the instability (see e.g., Fig. 4) fades away
for even lower aspect ratios and the time evolution becomes smoother like it
is in our FARGO simulations (see Fig. 12).

Figure 11: Time average of the maximal spiral amplitude at $R=0.5\,R_{0}$
between $\tau=0.2$–0.3 for different values of aspect ratio $h$. The spiral
instability is strongly damped for $h>0.006$ while spiral amplitudes converge
for $h<0.005$.

##  Appendix B Comparison of codes and the spiral flip in SK03

To test the robustness of our results on the spiral instability we performed a
similar suite of simulations using the FARGO code. Figure 12 shows the growth
of the spiral instability for different resolutions. We see that a resolution
of $465\times 256$ is sufficiently converged, with full convergence at
$1024\times 256$ and above. Even though PLUTO achieved convergence at a lower
resolution of $465\times 256$, FARGO shows slight improvements in the growth
onset and amplitudes for higher resolutions. Regardless, we can conclude that
the results are consistent between both codes. We note that just like PLUTO,
the spiral flip vanishes before it reaches the ring in FARGO simulations for
resolutions with more that 465 radial cells.

Figure 12: Time evolution of the normalized maximal deviation from the initial
surface density distribution for different FARGO simulations compared to the
fiducial PLUTO simulation.

To further understand the spiral flip results shown in SK03 we replicated
their setup in FARGO, with a $256\times 256$ arithmetically spaced polar grid.
Our simulation reproduces their Fig. 6 very well, but is run for longer. The
spiral flip pattern emerges after $\tau=0.04$ and gradually morphs into a
single leading spiral for $\tau>0.43$ (see Fig. 14). We conducted additional
simulations with higher resolutions and also on logarithmic grids. Each time
the resolution was increased on an arithmetic grid, the flip occurred further
inward and turns into a leading spiral earlier into the simulation. On a
logarithmic grid, the inner region is better resolved and the flip always
vanishes early into the simulation and never reaches the spreading ring, see
Fig. 14.

A comparison of the growth rates between the SK03 setup and a high resolution
FARGO run is shown in Fig. 13. As the grid is the only difference between the
simulations, the flip seems to be a numerical effect that is exacerbated when
using an arithmetic grid, which has a lower resolution at the inner boundary.
The constant smoothing length used for the SPH code in SK03 has the same
resolution effect as an arithmetic grid, which explains why their simulations
are in agreement. We conclude that the flipped spiral structure with a
trailing spiral observed at the inner domain in SK03 is merely a numerical
effect.

Figure 13: Time evolution of the first three Fourier azimuthal modes at
$R=R_{0}$. Solid lines represent a high-resolution ($2048\times 256$)
simulation using a logarithmic grid, and the dashed lines represent results
from an identical setup as in SK03, with a $256\times 256$ arithmetic grid.
Figure 14: Surface density deviation heatmaps similar to Fig. 3, for FARGO
simulations at a resolution of $256\times 256$ on an arithmetic grid (top) and
at $456\times 256$ on a logarithmic grid (bottom). Top: this setup reproduces
the results of SK03 (see Fig. 6 therein). The flip from a leading to trailing
spiral is visible in the middle two panels, before it disappears into a single
leading spiral. Bottom: this setup is equivalent to our fiducial PLUTO model,
and behaves very similar to our results in Fig. 3. Here, spirals develop
earlier and the initial flip is short-lived and closer to the inner boundary.

##  Appendix C Numerical diffusion as a consequence of advection

We consider the 1D advection equation of a quantity $q$ for a fluid moving at
constant velocity $u>0$

$\frac{\partial{q}}{\partial{t}}+u\frac{\partial{q}}{\partial{x}}=0.$ (7)

We can discretize this equation to a first-order upwind scheme (e.g., Courant
et al. 1952) where $q^{\prime}=q(t+\Delta t)$ on a grid where $i$ is the cell
index such that $q_{i-1}=q(x_{i}-\Delta x)$

$\frac{q_{i}^{\prime}-q_{i}}{\Delta t}+u\frac{q_{i}-q_{i-1}}{\Delta x}=0.$ (8)

By Taylor-expanding $q_{i}^{\prime}$ and $q_{i-1}$ to second order in time and
space respectively, and substituting a wave-like solution $q_{tt}=uq_{xx}$, we
arrive at a modified upwind equation:

$\frac{\partial{q}}{\partial{t}}+u\frac{\partial{q}}{\partial{x}}=D\frac{\partial^{2}q}{\partial
x^{2}},\qquad D=\frac{u\Delta x}{2}(1-C),\quad C\equiv u\frac{\Delta
t}{\Delta{x}}$ (9)

which corresponds to an advection–diffusion equation with diffusion
coefficient $D$. This term, while not physically motivated, allows the upwind
scheme to remain stable for Courant numbers $0<C<1$.

This approach does not exactly relate to our results since the gas velocity is
not constant and we further use second-order schemes with PLUTO. Nevertheless,
we can expect that our numerical scheme will give rise to a similarly-
motivated numerical diffusion term such that Eq. (5) effectively becomes

$\frac{\partial{(\Sigma\bm{u})}}{\partial{t}}+\nabla\cdot\left(\Sigma\bm{u}\otimes\bm{u}\right)=-\nabla
P-\Sigma\nabla\Phi_{\star}+\nabla\cdot(\overline{\bm{{\upsigma}}}+\overline{\bm{{\upsigma}}}_{\mathrm{num}}),$
(10)

with $\overline{\bm{{\upsigma}}}\propto\nu$ and
$\overline{\bm{{\upsigma}}}_{\mathrm{num}}\propto\nu_{\text{num}}$. Given that
our inviscid Cartesian models behave similarly to viscous models with
viscosity $\nu_{\text{num}}$ as far as ring spreading is concerned, it is
therefore unsurprising that a ring in a viscous model with
$\nu\sim\nu_{\text{num}}$ should spread as if the total viscosity were
$\nu_{\text{tot}}\approx\nu+\nu_{\text{num}}$.

We also note that, for a given Courant number $C$, the diffusion coefficient
in Eq. (9) is proportional to $\Delta x$ for the first-order method we
considered, and expect a $\Delta x^{2}$ scaling for a second-order method in
1D. Given that the viscous ring problem evolves on both the $x$ and $y$
directions in 2D, doubling the resolution would require increasing the number
of cells in both directions by a factor of 2 (thus increasing grid cell count
by a factor of 4). In doing so, we would expect a scaling
$\nu_{\mathrm{num}}\propto\Delta x^{2}$, which is verified in Fig. 8.

Finally, we found that doubling the resolution in only one direction (such as
$N_{x}\times N_{y}=512\times 1024$) results in a numerical viscosity estimate
that is much closer to the value expected for the low-resolution direction (in
this case, $N_{x}$) than to an estimate dictated by $\Delta x\times\Delta y$.
This rules out that the numerical diffusion we observe depends linearly on
$\Delta x$ and $\Delta y$ and suggests that it more accurately depends on
$\min(\Delta x,\Delta y)^{2}$, highlighting the second-order accuracy of our
solver.

##  Appendix D Measuring numerical viscosity

To measure the numerical viscosity in our inviscid simulations we utilize the
relation between the viscous time scale $\tau$ for the ring spreading and the
orbital time t at $R_{0}$ given by $\tau=12\nu t/R^{2}_{0}$. We note that for
a given viscosity $\nu$, this equation represents the slope–intercept form of
an equation for a line with slope $m=12\nu/R_{0}^{2}$ and intercept zero.

To do so, first we fit the analytical solution to a slice of the two-
dimensional data at time $t$ to compute $\tau$. Repeating this procedure for
all snapshots generates a curve for $\tau(t)$ shown in Fig. 15. We then fit a
straight line to the $\tau(t)$ curve for the last $\sim$ 180 orbits. The slope
of this straight line fit gives the numerical viscosity using the above
relation for the slope $m$.

Figure 15: The $\tau(t)$ curve for our inviscid 1024 model. The dashed line
shows the fit over the last 180 time units. The slope of this line is used to
extract the numerical viscosity.

##  Appendix E Inviscid Cartesian simulations with Athena++

For the suite of Cartesian models using Athena++ we implement an equivalent
setup to the PLUTO runs, as described in Sect. 2. A locally isothermal disk is
achieved by using an adiabatic equation of state combined with a thermal
relaxation timescale $\tau_{\textrm{cool}}=0.01\,\Omega_{\text{K}}^{-1}$. The
viscous ring sits on top of a constant density background with
$\Sigma_{\textrm{bgr}}=10^{-7}\,\Sigma_{\textrm{ref}}$, however the density
floor of the domain was set to
$\Sigma_{\textrm{floor}}=10^{-15}\,\Sigma_{\textrm{ref}}$ as interactions with
the density floor created numerical instabilities.

Our results are listed in Table LABEL:tab:numerical_viscosities and shown in
Fig. 8. Both codes agree very well for $N_{\textrm{cells}}\geq 1024^{2}$,
which corresponds to our fiducial resolution. The expected behavior of
$\nu_{\mathrm{num}}\propto\Delta x^{2}$ is also found with Athena++.