# An Empirical Exploration in Quality Filtering of Text Data

Leo Gao
EleutherAI
<EMAIL_ADDRESS>
###### Abstract

While conventional wisdom suggests that more aggressively filtering data from
low-quality sources like Common Crawl always monotonically improves the
quality of training data, we find that aggressive filtering can in fact lead
to a decrease in model quality on a wide array of downstream tasks for a GPT-
like language model. We speculate that this is because optimizing sufficiently
strongly for a proxy metric harms performance on the true objective,
suggesting a need for more robust filtering objectives when attempting to
filter more aggressively. We hope this work leads to detailed analysis of the
effects of dataset filtering design choices on downstream model performance in
future work.

##  1 Introduction

As language models increase in size, the need for large, high-quality text
datasets has increased as well. Recent work in dataset construction for large
language models has centered largely on taking large internet corpora like
Common Crawl and employing some method of filtering using some proxy for
quality to extract a smaller, high quality training set (Wenzek et al., 2019;
Brown et al., 2020; Raffel et al., 2020; Yang et al., 2020). In particular, we
focus on shallow classifier-based quality filtering as in Brown et al. (2020)
because it provides a simple, continuous, and quantifiable way to adjust the
aggressiveness of filtering, and because this reflects the type of classifier
used in prior work.

While intuitively it may seem like the more data is discarded the higher
quality the remaining data will be, we find that this is not always the case
with shallow classifier-based filtering. Instead, we find that filtering
improves downstream task performance up to a point, but then decreases
performance again as the filtering becomes too aggressive.

We speculate that this decrease in performance is due to Goodhart’s law
(Goodhart, 1984), and specifically regressional Goodharting (Manheim and
Garrabrant, 2019):

Goodhart’s Law. Any observed statistical regularity will tend to collapse once
pressure is placed upon it for control purposes. (Goodhart, 1984)

In other words, optimizing a metric that is a proxy for a desired outcome
tends to invalidate the proxy. By optimizing too strongly for the classifier’s
score by discarding too many low-scoring documents, the documents that are
kept are consistently biased towards the ones with features superficially
resembling the high quality data in a way that satisfies the classifier,
rather than truly high quality data.

Figure 1: Average accuracy across all 13 tasks for various different filtering
ratios using a shallow quality classifier.222The average is taken across all
task accuracies, with each task weighted equally. The error bars in this plot
represent standard error and are computed by
$\mathrm{se_{mean}}=n^{-1}\sqrt{\sum\mathrm{se}_{i}^{2}}$, where
$\mathrm{se}_{i}$ represents the standard error for each individual task.The
amount of data post-filtering is held constant. Although filtering improves
performance at first, discarding more data can actually reduce accuracy, due
to misalignment between filtering classifier objective and text quality.

##  2 Related work

The recent proliferation of ever larger language models has led to increasing
demands on training data (Radford et al., 2018, 2019; Gokaslan and Cohen,
2019; Rosset, 2019; Shoeybi et al., 2019; Devlin et al., 2019; Liu et al.,
2019; Raffel et al., 2020; Brown et al., 2020; Zeng et al., 2021). This data
is increasingly derived from internet corpora like Common Crawl (Radford et
al., 2019; Ortiz Suárez et al., 2019; Wenzek et al., 2020; Conneau et al.,
2020; Brown et al., 2020; Gao et al., 2020; Raffel et al., 2020).

However, the quality of raw Common Crawl data is often insufficient to be
directly used. To combat this, many existing works use some kind of proxy for
quality, like a classifier between known high quality data and low quality
data (Brown et al., 2020; Gao et al., 2020; Zeng et al., 2021), hand-crafted
heuristics (Yang et al., 2020; Raffel et al., 2020), or keeping only documents
with perplexity scores that fall in some middle quantile of an existing
language model (Wenzek et al., 2020). Brown et al. (2020) in particular filter
extremely aggressively using their classifier, discarding about 98.7% of their
data.

Previous work has shown that models trained on heuristic-filtered datasets
perform better on downstream tasks (Raffel et al., 2020). However, Gao et al.
(2020) show that a perplexity-filtered CC-derived dataset actually performs
worse than unfiltered CC on certain tasks. Brown et al. (2020) do not provide
any detailed analysis, but claim better quality for filtered data as evaluated
through loss on held out sets of “generative text samples."

##  3 Downstream Evaluation Experiment

Figure 2: Plots of results for all downstream tasks explored in this paper.
Higher is better on all metrics except LAMBADA perplexity (first plot in the
third row), where lower is better.

To evaluate the effect of different degrees of filtering, we create a series
of training sets with a controlled filtering methodology but with different
hyperparameter settings to result in varied filtering ratios. We filter using
the same method used in Brown et al. (2020), with a Pareto-distribution
thresholded filtering method and a shallow CommonCrawl-WebText classfier. In
this method, rather than using a hard threshold, the threshold
$\tau\sim\mathrm{Pareto}(\alpha)$ is sampled from a Pareto distribution, such
that each document is kept if $\tau>1-\text{{score}}$, where $\alpha$ is a
hyperparameter that controls the permissivity of the filter (see Table 1). In
effect, this relaxes the filter when compared to a hard threshold and allows
some low-scoring data to be kept.

$\alpha$ | Fraction Discarded
---|---
1 | 0.4107
2 | 0.6351
3 | 0.7610
4 | 0.8329
5 | 0.8761
6 | 0.9026
7 | 0.9198
8 | 0.9315
Table 1: Percentage of discarded documents of various settings using our
classifier.

As none of the data or models used in Brown et al. (2020) has been made
public, we instead use the same type of fasttext (Joulin et al., 2017)
classifier between unfiltered Common Crawl and OpenWebText2 as used in Gao et
al. (2020).

We use GPT-Neo (Black et al., 2021) to train a series of models on each
training set and evaluate on downstream tasks using the EleutherAI LM
evaluation harness (Gao et al., 2021). Each model is 1.3 billion parameters,
has a GPT-2 architecture (Radford et al., 2019) with the same model
hyperparameters as the GPT-3-XL setting in Brown et al. (2020), and is trained
for 25k iterations with a batch size of 256.

To ensure that the effect is not confined to any specific task. we evaluate on
a series of many downstream tasks. We use zero-shot prompting with no task-
specific fine tuning and with prompting inspired by Brown et al. (2020) for
many tasks. In total, we evaluate on ANLI Round 3 (Nie et al., 2020), BoolQ
(Clark et al., 2019), CommitmentBank (de Marneffe et al., 2019), COPA (Gordon
et al., 2012), Hellaswag (Zellers et al., 2019), LAMBADA (Paperno et al.,
2016), MathQA (Amini et al., 2019), MultiRC (Khashabi et al., 2018),
OpenbookQA (Mihaylov et al., 2018), PiQA (Bisk et al., 2019), PubmedQA (Jin et
al., 2019), SciQ (Welbl et al., 2017), and Winogrande (Sakaguchi et al.,
2019). Error bars in all evaluation task plots indicate standard error with
respect to instances of the evaluation task.

For the training data, we create 40 GB filtered chunks of the Common Crawl
data for each value of $\alpha\in\\{1,2,3,4,5,8\\}$; in other words, different
amounts of raw Common Crawl data are consumed for different $\alpha$ to
produce the same fixed 40GB size result. For reference, Brown et al. (2020)
filter even more aggressively than we do, discarding about 98.7% of their
data. The 40GB size is chosen because it is approximately the size of
OpenWebText, which is representative of the amount of data usually used to
train models of this size.

###  3.1 Results

Of the tasks evaluated, several tasks remained near chance or had very high
variance, resulting in no clear trend. Of the remainder, an absolute majority
exhibited an initial increase in performance and then a decrease in
performance after the amount of documents discarded surpassed a threshold that
varied by task. Additionally, for almost all tasks the most filtered model was
not the best performing. Some tasks like BoolQ exhibit little clear trend. Not
all tasks have the same optimal $\alpha$—compare PiQA and LAMBADA—and some
tasks like PubmedQA show a much more sudden decrease in accuracy. For results
on all tasks, see Figure 2.

###  3.2 Analysis

We hypothesize that this decline in performance is because of misalignment
between the classifier objective, intended to be a proxy for quality, and
actual document quality. For instance, a classifier to distinguish WebText2
from Common Crawl, as in GPT-3, would also exclude domains of text data not
found as often in WebText2.

We also hypothesize that the difference in optimal $\alpha$ between different
tasks is because the characteristics of the different types of data that help
the most with each task are over/underdiscarded to a different extent due to
spurious correlations with the quality metric. As such, we do not expect the
exact thresholds to transfer to other tasks, classifiers, or datasets. This is
an expected consequence of Goodharting, because the degree to which different
types of text data correlate with the features learned by the classifier is
mostly spurious.

##  4 Domain Misalignment Experiment

To test the hypothesis that the misalignment of the objective leads to the
exclusion of non-OpenWebText2-like data, we train a fasttext classifier to
classify between BookCorpus2 (Gao et al., 2020) and OpenWebText (Gokaslan and
Cohen, 2019), and compute the mean BookCorpus2-probability of each training
set. If the classification model is favoring OpenWebText-like data over
generally high-quality data, then as filtering increases in intensity, the
proportion of BookCorpus2-like data should decrease as the data consists
increasingly of OpenWebText-like text. Conversely, if the classification model
is robustly favoring high quality text, then as filtering increases in
intensity, the proportion of BookCorpus2-like data should increase, as low-
quality text looks nothing like BookCorpus2. We also repeat this experiment
for Pubmed Abstracts.

We chose BookCorpus2 and Pubmed Abstracts because of their similarity in
distribution to LAMBADA and PubmedQA respectively, in the hopes of observing a
similarity between the task evaluation curves and the data domain curves.

###  4.1 Results

As seen in Figure 3, the fraction of BookCorpus2-like data remains mostly
constant until around 0.6, after which it declines sharply. A similar pattern
is observed with Pubmed Abstracts, albeit with an earlier drop (Figure 4).

The BookCorpus2-like data curve’s drop precedes the LAMBADA performance drop
by about 0.2. Similarly, the Pubmed Abstracts drop also precedes the
PubmedQA’s main drop slightly.

###  4.2 Analysis

Figure 3: Fraction of documents in filtered Common Crawl classified as
BookCorpus2-like by a shallow classifier trained to distinguish OpenWebtext
and BookCorpus2. Note that this plot has a different x-axis scale from the
task evaluation plots.

The decrease in Pubmed Abstracts and BookCorpus2 like data as filtering
increases in aggressiveness supports the hypothesis that part of the problem
is that text domains not similar to OpenWebText2 are being discarded.

Our main hypothesis for why the domain data content starts decreasing before
the evaluation metric performance does is that these tasks are sufficiently
different in distribution to the respective datasets.

Figure 4: Fraction of documents in filtered Common Crawl classified as
PubmedAbstracts-like by a shallow classifier trained to distinguish
OpenWebtext and PubmedAbstracts. Note that this plot has a different x-axis
scale from the task evaluation plots.

##  5 Limitations

This work is intended to show that the common assumption that more aggressive
data filtering is better is not always true, and thus focuses on one
particular classifier used in the real world as an illustrative example.
Depending on the type of classifier, the training data used for the
classifier, and the downstream task, this effect may not be relevant in
certain settings. We leave an exhaustive exploration of the contribution of
these various factors to future work.

##  6 Conclusion

In this paper, we explored the effect of filtering the training data using a
shallow model trained on a proxy for quality on downstream language model
performance. We showed that increasing the aggressiveness of filtering against
this signal actually decreases model performance past a certain point, and
speculate that this is due to Goodhart’s law, as the misalignment between
proxy and true objective becomes more significant with increased optimization
pressure. We hope that this work leads to more careful analysis of the effects
of filtering in future language modeling work.

## Acknowledgements

The author would like to thank TPU Research Cloud for providing the
computational resources for the training, and CoreWeave for providing the
computational resources for data processing and evaluation.

The author would also like to thank Stella Biderman, Sid Black, Charles
Foster, Eric Hallahan, Kyle McDonell, Jason Phang, and Laria Reynolds for
providing feedback on the manuscript.

## References

  * Amini et al. (2019) Aida Amini, Saadia Gabriel, Shanchuan Lin, Rik Koncel-Kedziorski, Yejin Choi, and Hannaneh Hajishirzi. 2019.  MathQA: Towards interpretable math word problem solving with operation-based formalisms.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 2357–2367, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Bisk et al. (2019) Yonatan Bisk, Rowan Zellers, Ronan Le Bras, Jianfeng Gao, and Yejin Choi. 2019.  Piqa: Reasoning about physical commonsense in natural language.  _arXiv preprint arXiv:1911.11641_.
  * Black et al. (2021) Sid Black, Gao Leo, Phil Wang, Connor Leahy, and Stella Biderman. 2021.  GPT-Neo: Large Scale Autoregressive Language Modeling with Mesh-Tensorflow.  If you use this software, please cite it using these metadata.
  * Brown et al. (2020) Tom B. Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M. Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei. 2020.  Language models are few-shot learners.  _arXiv preprint arXiv:2005.14165_.
  * Clark et al. (2019) Christopher Clark, Kenton Lee, Ming-Wei Chang, Tom Kwiatkowski, Michael Collins, and Kristina Toutanova. 2019.  BoolQ: Exploring the surprising difficulty of natural yes/no questions.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 2924–2936, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Conneau et al. (2020) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2020.  Unsupervised cross-lingual representation learning at scale.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 8440–8451, Online. Association for Computational Linguistics.
  * de Marneffe et al. (2019) Marie-Catherine de Marneffe, Mandy Simons, and Judith Tonhauser. 2019.  The commitmentbank: Investigating projection in naturally occurring discourse.  _Proceedings of Sinn und Bedeutung_ , 23(2):107–124.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_. Association for Computational Linguistics.
  * Gao et al. (2020) Leo Gao, Stella Biderman, Sid Black, Laurence Golding, Travis Hoppe, Charles Foster, Jason Phang, Horace He, Anish Thite, Noa Nabeshima, Shawn Presser, and Connor Leahy. 2020.  The pile: An 800gb dataset of diverse text for language modeling.  _arXiv preprint arXiv:2101.00027_.
  * Gao et al. (2021) Leo Gao, Jonathan Tow, Stella Biderman, Sid Black, Anthony DiPofi, Charles Foster, Laurence Golding, Jeffrey Hsu, Kyle McDonell, Niklas Muennighoff, Jason Phang, Laria Reynolds, Eric Tang, Anish Thite, Ben Wang, Kevin Wang, and Andy Zou. 2021.  A framework for few-shot language model evaluation.
  * Gokaslan and Cohen (2019) Aaron Gokaslan and Vanya Cohen. 2019.  Openwebtext corpus.  http://Skylion007.github.io/OpenWebTextCorpus.
  * Goodhart (1984) C. A. E. Goodhart. 1984.  _Problems of Monetary Management: The UK Experience_, pages 91–121. Macmillan Education UK, London.
  * Gordon et al. (2012) Andrew Gordon, Zornitsa Kozareva, and Melissa Roemmele. 2012.  SemEval-2012 task 7: Choice of plausible alternatives: An evaluation of commonsense causal reasoning.  In _*SEM 2012: The First Joint Conference on Lexical and Computational Semantics – Volume 1: Proceedings of the main conference and the shared task, and Volume 2: Proceedings of the Sixth International Workshop on Semantic Evaluation (SemEval 2012)_ , pages 394–398, Montréal, Canada. Association for Computational Linguistics.
  * Jin et al. (2019) Qiao Jin, Bhuwan Dhingra, Zhengping Liu, William W. Cohen, and Xinghua Lu. 2019\.  Pubmedqa: A dataset for biomedical research question answering.  _arXiv preprint arXiv:1909.06146_.
  * Joulin et al. (2017) Armand Joulin, Edouard Grave, Piotr Bojanowski, and Tomas Mikolov. 2017.  Bag of tricks for efficient text classification.  In _Proceedings of the 15th Conference of the European Chapter of the Association for Computational Linguistics: Volume 2, Short Papers_ , pages 427–431. Association for Computational Linguistics.
  * Khashabi et al. (2018) Daniel Khashabi, Snigdha Chaturvedi, Michael Roth, Shyam Upadhyay, and Dan Roth. 2018.  Looking beyond the surface: A challenge set for reading comprehension over multiple sentences.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 252–262, New Orleans, Louisiana. Association for Computational Linguistics.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  RoBERTa: A robustly optimized BERT pretraining approach.  _arXiv preprint arXiv:1907.11692_.
  * Manheim and Garrabrant (2019) David Manheim and Scott Garrabrant. 2019.  Categorizing variants of goodhart’s law.  _arXiv preprint arXiv:1803.04585_.
  * Mihaylov et al. (2018) Todor Mihaylov, Peter Clark, Tushar Khot, and Ashish Sabharwal. 2018.  Can a suit of armor conduct electricity? a new dataset for open book question answering.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 2381–2391, Brussels, Belgium. Association for Computational Linguistics.
  * Nie et al. (2020) Yixin Nie, Adina Williams, Emily Dinan, Mohit Bansal, Jason Weston, and Douwe Kiela. 2020.  Adversarial nli: A new benchmark for natural language understanding.  _arXiv preprint arXiv:1910.14599_.
  * Ortiz Suárez et al. (2019) Pedro Javier Ortiz Suárez, Benoît Sagot, and Laurent Romary. 2019.  Asynchronous pipelines for processing huge corpora on medium to low resource infrastructures.  Proceedings of the Workshop on Challenges in the Management of Large Corpora (CMLC-7) 2019. Cardiff, 22nd July 2019, pages 9 – 16, Mannheim. Leibniz-Institut für Deutsche Sprache.
  * Paperno et al. (2016) Denis Paperno, Germán Kruszewski, Angeliki Lazaridou, Ngoc Quan Pham, Raffaella Bernardi, Sandro Pezzelle, Marco Baroni, Gemma Boleda, and Raquel Fernández. 2016.  The LAMBADA dataset: Word prediction requiring a broad discourse context.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1525–1534, Berlin, Germany. Association for Computational Linguistics.
  * Radford et al. (2018) Alec Radford, Karthik Narasimhan, Time Salimans, and Ilya Sutskever. 2018.  Improving language understanding with unsupervised learning.  _Technical report, OpenAI_.
  * Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019.  Language models are unsupervised multitask learners.  _OpenAI Blog_ , 1(8):9.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu. 2020.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _arXiv preprint arXiv:1910.10683_.
  * Rosset (2019) C Rosset. 2019.  Turing-NLG: A 17-billion-parameter language model by Microsoft.  _Microsoft Blog_.
  * Sakaguchi et al. (2019) Keisuke Sakaguchi, Ronan Le Bras, Chandra Bhagavatula, and Yejin Choi. 2019.  Winogrande: An adversarial winograd schema challenge at scale.  _arXiv preprint arXiv:1907.10641_.
  * Shoeybi et al. (2019) Mohammad Shoeybi, Mostofa Patwary, Raul Puri, Patrick LeGresley, Jared Casper, and Bryan Catanzaro. 2019.  Megatron-LM: Training multi-billion parameter language models using gpu model parallelism.  _arXiv preprint arXiv:1909.08053_.
  * Welbl et al. (2017) Johannes Welbl, Nelson F. Liu, and Matt Gardner. 2017.  Crowdsourcing multiple choice science questions.  In _Proceedings of the 3rd Workshop on Noisy User-generated Text_ , pages 94–106, Copenhagen, Denmark. Association for Computational Linguistics.
  * Wenzek et al. (2020) Guillaume Wenzek, Marie-Anne Lachaux, Alexis Conneau, Vishrav Chaudhary, Francisco Guzmán, Armand Joulin, and Edouard Grave. 2020.  CCNet: Extracting high quality monolingual datasets from web crawl data.  In _Proceedings of the 12th Language Resources and Evaluation Conference_ , pages 4003–4012, Marseille, France. European Language Resources Association.
  * Wenzek et al. (2019) Guillaume Wenzek, Marie-Anne Lachaux, Alexis Conneau, Vishrav Chaudhary, Francisco Guzmán, Armand Joulin, and Edouard Grave. 2019.  Ccnet: Extracting high quality monolingual datasets from web crawl data.  _arXiv preprint arXiv:1911.00359_.
  * Yang et al. (2020) Zhilin Yang, Zihang Dai, Yiming Yang, Jaime Carbonell, Ruslan Salakhutdinov, and Quoc V. Le. 2020.  Xlnet: Generalized autoregressive pretraining for language understanding.  _arXiv preprint arXiv:1906.08237_.
  * Zellers et al. (2019) Rowan Zellers, Ari Holtzman, Yonatan Bisk, Ali Farhadi, and Yejin Choi. 2019.  HellaSwag: Can a machine really finish your sentence? In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 4791–4800, Florence, Italy. Association for Computational Linguistics.
  * Zeng et al. (2021) Wei Zeng, Xiaozhe Ren, Teng Su, Hui Wang, Yi Liao, Zhiwei Wang, Xin Jiang, ZhenZhang Yang, Kaisheng Wang, Xiaoda Zhang, Chen Li, Ziyan Gong, Yifan Yao, Xinjing Huang, Jun Wang, Jianfeng Yu, Qi Guo, Yue Yu, Yan Zhang, Jin Wang, Hengtao Tao, Dasen Yan, Zexuan Yi, Fang Peng, Fangqing Jiang, Han Zhang, Lingfeng Deng, Yehong Zhang, Zhe Lin, Chao Zhang, Shaojie Zhang, Mingyue Guo, Shanzhi Gu, Gaojun Fan, Yaowei Wang, Xuefeng Jin, Qun Liu, and Yonghong Tian. 2021\.  Pangu-$\alpha$: Large-scale autoregressive pretrained chinese language models with auto-parallel computation.  _arXiv preprint arXiv:2104.12369_.


# Gap Opening in Protoplanetary Disks: Gas Dynamics from Global Non-ideal MHD
Simulations with Consistent Thermochemistry

Xiao Hu (胡晓)1,2, Zhi-Yun Li1, Lile Wang (王力乐)3, Zhaohuan Zhu (朱照寰)4,5 and
Jaehan Bae2
1Department of Astronomy, University of Virginia, Charlottesville, VA 22904,
USA
2Department of Astronomy, University of Florida, Gainesville, FL 32608, USA
3The Kavli Institute for Astronomy and Astrophysics, Peking University,
Beijing 100084, China
4Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 S.
Maryland Parkway, Las Vegas, NV 89154, USA
5Nevada Center for Astrophysics, University of Nevada, Las Vegas, 4505 South
Maryland Parkway, Las Vegas, NV 89154, USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>0000-0003-3201-4549
0000-0002-6540-7042 0000-0003-3616-6822 0000-0001-7258-770X

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Recent high angular resolution ALMA observations have revealed numerous gaps
in protoplanetary disks. A popular interpretation has been that planets open
them. Most previous investigations of planet gap-opening have concentrated on
viscous disks. Here, we carry out 2D (axisymmetric) global simulations of gap
opening by a planet in a wind-launching non-ideal MHD disk with consistent
thermochemistry. We find a strong concentration of poloidal magnetic flux in
the planet-opened gap, where the gas dynamics are magnetically dominated. The
magnetic field also drives a fast (nearly sonic) meridional gas circulation in
the denser disk regions near the inner and outer edges of the gap, which may
be observable through high-resolution molecular line observations. The gap is
more ionized than its denser surrounding regions, with a better magnetic
field-matter coupling. In particular, it has a much higher abundance of
molecular ion HCO+, consistent with ALMA observations of the well-studied AS
209 protoplanetary disk that has prominent gaps and fast meridional motions
reaching the local sound speed. Finally, we provide fitting formulae for the
ambipolar and Ohmic diffusivities as a function of the disk local density,
which can be used for future 3D simulations of planet gap-opening in non-ideal
MHD disks where thermochemistry is too computationally expensive to evolve
self-consistently with the magneto-hydrodynamics.

###### keywords:

accretion, accretion disks — magnetohydrodynamics (MHD) — planets and
satellites: formation — circumstellar matter — method: numerical

††pubyear: 2023††pagerange: Gap Opening in Protoplanetary Disks: Gas Dynamics
from Global Non-ideal MHD Simulations with Consistent Thermochemistry–Gap
Opening in Protoplanetary Disks: Gas Dynamics from Global Non-ideal MHD
Simulations with Consistent Thermochemistry

##  1 Introduction

Planet gap opening in protoplanetary disks (PPDs hereafter) is a classic
research topic in planet formation (Lin & Papaloizou, 1993; Armitage, 2010).
The gravity from the planet launches spiral density waves, which deposit
significant torque in the disk (Goldreich & Tremaine, 1979; Kley & Nelson,
2012). This torque is large enough to alter the disk structure, resulting in a
local low-density annulus (i.e., a planetary gap) along the planetary orbit
(Crida et al., 2006; Duffell & MacFadyen, 2013). This process affects both the
accretion of gas and solids and the orbital migration of the planet (Baruteau
et al., 2014; Kanagawa et al., 2015). These factors largely shape the outcome
of an exoplanet system, determining planet mass and orbital radius (Ida & Lin,
2004; Mordasini et al., 2012).

The critical mass of gap opening can be derived from the shear stability of
steep pressure gradient (Lin & Papaloizou, 1993) and the time scale for the
disk viscosity to refill it (e.g., Duffell & MacFadyen, 2013; Bitsch et al.,
2018; Hu et al., 2018). The process of heat dissipation is important for the
disk torque (Lega et al., 2014; Masset, 2017). The gap depth scaling relations
were obtained in semi-analytical work and simplified 2D (infinitely thin) disk
simulations, and they were shown to hold well when moving to higher dimensions
(Fung & Chiang, 2016). Three-dimensional (3D) simulations provide, in
addition, a characterization of the meridional flows, which turn out to be
amenable to ALMA observations and thus provide a powerful tool to probe the
planet-disk interaction (e.g., Teague et al., 2019; Pinte et al., 2020).

It is well known that disk dynamics and evolution are strongly affected by
magnetic fields. For example, magnetic fields have been shown to play
important roles in carving gaps by direct surface accretion stream (Suriano et
al., 2017) or radial magnetic flux redistribution (e.g., Suriano et al., 2018;
Hu et al., 2019; Hu et al., 2021; Cui & Bai, 2021) in protoplanetary disks,
and maintain meridional flows (Hu et al., 2022). In addition, they likely
dominate the angular momentum transport in PPDs through magnetized disk winds
(Blandford & Payne, 1982; Bai & Stone, 2013). Their dynamical importance
motivated numerical studies that combine planet-disk interaction with magnetic
fields and wind-driven accretion. Early simulations along this line usually
assumed ideal MHD and incorporated only a toroidal magnetic field without
launching a disk wind (Nelson & Papaloizou, 2003; Winters et al., 2003;
Baruteau et al., 2011). More recent local shearing-box MHD simulations have
included a net poloidal magnetic flux and found that the magnetic flux gets
concentrated into the planet-induced gap, making the gap deeper and wider due
to enhanced MRI turbulence within the gap (Zhu et al., 2013; Carballido et
al., 2017). However, besides being local, these studies are typically
unstratified in the vertical direction and thus incapable of wind-driven
angular momentum transport. As importantly, disks are weakly magnetized, so
non-ideal MHD effects are important. Gressel et al. (2013) included Ohmic
dissipation in a global simulation, but its extent is limited to 4.5 disk
scale heights, limiting the disk wind’s treatment. In this paper, we will go
beyond the previous work by extending the simulation domain to close to the
polar axis and including ambipolar diffusion to have comprehensive coverage of
the wind region, which is important for the scale of tens to hundreds of au
most accessible to ALMA observations.

Because of the complexity and especially the computational demand of self-
consistent thermochemistry evolution, we will carry out our global non-ideal
MHD simulations assuming axisymmetry as a first step, adopting a torque
profile to approximate the planet-disk interaction. This approximation will be
relaxed in future 3D simulations. The focus is on the gas kinematics in the
vicinity of the planet-opened gap, which is starting to be probed in
increasing detail by ALMA. Our approach complements the recent work of Aoyama
& Bai (2023), who carried out 3D simulations of gap-opening by planets in non-
ideal MHD disks but with a prescribed spatially-uniform ambipolar Elsasser
number inside the disk. However, a significant spatial variation of ambipolar
Elsasser number is expected given the large dynamic range of the density
involved in the gap opening process. We seek to capture this variation through
consistent thermochemistry that includes ionization.

The paper is organized as follows. In §2, we describe the numerical methods
and simulation setup, especially the planet torque implementation and
validation through hydro-only setups. §3 analyzes the results from our
simulation, focusing on the magnetic effects on the gas kinematics in and
around the planet-opened gap and the relation between the ambipolar and ohmic
diffusivities and the disk density. We discuss the observation implications of
our results and conclude in §4.

##  2 Methods

###  2.1 Disk

We simulate the disk evolution with the combination of non-ideal MHD effects
using the higher-order Godunov MHD code `Athena++` (Stone et al., 2020), ray-
tracing radiative transfer for high-energy photons, and consistent
thermochemistry. For each MHD timestep, the non-equilibrium thermochemistry is
co-evolved in each zone throughout the simulation domain with a semi-implicit
method. In general, this numerical system is almost the same as Wang et al.
(2019) (WBG19 hereafter). We refer the reader to WBG19 and references therein
for the details of the initial conditions and the overall setup of
thermochemical reactions, which included 28 species: $\mathrm{e}^{-}$ (free
electrons), $\mathrm{H}^{+}$, H, $\mathrm{H_{2}}$, He, $\mathrm{He}^{+}$, O,
$\mathrm{O}^{+}$, OH, $\mathrm{OH}^{+}$, $\mathrm{H_{2}O}$, C,
$\mathrm{C}^{+}$, CO, CH, $\mathrm{CH}^{+}$, $\mathrm{HCO}^{+}$, Si,
$\mathrm{Si}^{+}$, SiO, $\mathrm{SiO}^{+}$, $\mathrm{SiOH}^{+}$, S,
$\mathrm{S}^{+}$, $\mathrm{HS}^{+}$, Gr, $\mathrm{Gr}^{+}$, $\mathrm{Gr}^{-}$.
Here Gr and Gr± denote neutral and singly-charged dust grains, respectively.
For the boundary conditions, we inherited the setups in Hu et al. (2021),
i.e., similar to WBG19 except for the toroidal field above the disk region
(viz. inside the wind region) at the inner radial boundary: we set
$B_{\phi}=-B_{r}|_{t=0}$ (the initial value of $r$ component) to suppress
magnetic instabilities there. Other hydrodynamic and field components are
identical to WBG19.

The dust grains are treated as single-sized carbonaceous grains co-moving with
the gas. The size is $a_{\mathrm{Gr}}=5$ Å. The basic properties of our model
are summarized in Table 1. The main parameter that determines the dust
thermochemical properties, dust grain cross section per hydrogen nuclei
$\sigma_{\mathrm{Gr}}/\mathrm{H}$, is set to $8\times
10^{-23}~{}\mathrm{cm}^{2}$ which corresponds to a dust-to-gas mass ratio of
$7\times 10^{-6}$. This is also the maximum $\sigma_{\mathrm{Gr}}/\mathrm{H}$
in Hu et al. (2021), which is 10 times lower than WBG19.

Table 1: Properties of Disk Model (§3)

Item | Value
---|---
Radial domain | $2~{}\mathrm{AU}\leq r\leq\ 100~{}\mathrm{AU}$
Latitudinal domain | $0.035~{}{\rm rad}\leq\theta\leq\pi/2~{}{\rm rad}$
Resolution |  $N_{\log r}=480$, $N_{\theta}=128$
Stellar mass | $1.0~{}M_{\odot}$
Initial mid-plane density | $8\times 10^{14}(R/\mathrm{AU})^{-2.2218}~{}m_{p}~{}\mathrm{cm}^{-3}$
Initial mid-plane plasma $\beta$ | $10^{4}$
Initial mid-plane temperature | $305(R/\mathrm{AU})^{-0.57}~{}\mathrm{K}$
Artificial heating profile$\dagger$ | $305(R/\mathrm{AU})^{-0.57}~{}\mathrm{K}$
Luminosities [photon $\mathrm{s}^{-1}$] |
$7~{}\mathrm{eV}$ (“soft” FUV) | $4.5\times 10^{42}$
$12~{}\mathrm{eV}$ (LW) | $1.6\times 10^{40}$
$3~{}\mathrm{keV}$ (X-ray) | $1.1\times 10^{38}$
Initial abundances [$n_{\mathrm{X}}/n_{\mathrm{H}}$] |
$\mathrm{H_{2}}$ | 0.5
He | 0.1
$\mathrm{H_{2}O}$ | $1.8\times 10^{-4}$
CO | $1.4\times 10^{-4}$
S | $2.8\times 10^{-5}$
SiO | $1.7\times 10^{-6}$
Dust/PAH properties |
$a_{\mathrm{Gr}}$ |  $5$ Å
$\sigma_{\mathrm{Gr}}/\mathrm{H}$ | $8\times 10^{-23}~{}{\rm cm^{2}}$

$\dagger$: The artificial heating profile indicates the temperature of dusts
in the mid-plane. Similar to WBG19, because the radiative transfer of diffuse
infrared radiation field is not calculated in this paper, we adopt this
profile of dust temperature as the floor of dust temperature.
---

###  2.2 Magnetic Diffusion

In a weakly ionized protoplanetary disk, the equation of motion for charged
species ($\mathrm{e}^{-}$, ions and charged grains $\mathrm{Gr}^{\pm}$) is set
by the balance between the Lorentz force and the neutral-ion drag force:

$Z_{j}e\left(\bm{E}^{{}^{\prime}}+\frac{\bm{v}_{j}}{c}\times\bm{B}\right)=\gamma_{j}\rho
m_{j}\bm{v}_{j}$ (1)

where, for a given charged particle type (denoted by $j$), $Z_{j}$ represents
its charge number (in units of charge $e$), $m_{j}$ is its mass, and
$\bm{v}_{j}$ is its drift velocity relative to the neutral background.
$\gamma_{j}$ is defined as $\langle\sigma v\rangle/(m+m_{j})$, where $m$ is
the average particle mass of the neutrals, and $\langle\sigma v_{j}\rangle$ is
the rate coefficient for momentum transfer between the charged particle and
neutrals. The electric field in the frame moving with the neutrals is
represented by $\bm{E}^{{}^{\prime}}$.

For very small grains, the interaction between charged grains and neutral
molecules is influenced by the electric field from induced electric dipoles in
the neutrals. This $r^{-4}$ electric potential results in a temperature-
independent collision rate coefficient $\langle\sigma v\rangle_{j}$(e.g.,
Draine, 2011). As grains grow larger, their geometrical cross-section becomes
more significant in interactions with neutrals. This grain size-driven
transition is reflected in the Hall parameter, which is the ratio of charged
particles’ gyrofrequency under the Lorentz force to their collision frequency
with neutrals:

$\beta_{j}=\frac{\absolutevalue{Z_{j}}eB}{m_{j}c}\frac{1}{\gamma_{j}\rho}$ (2)

In our simulation, we adopt the recipes in Bai (2011, 2014) to calculate the
collision coefficients between the charged grains and neutrals:

$\begin{split}\langle\sigma v\rangle_{\mathrm{Gr}}=\rm{max}&\bigg{[}1.3\times
10^{-9}\absolutevalue{Z_{\mathrm{Gr}}},\\\ 4\times
10^{-3}&\left(\frac{a_{\mathrm{Gr}}}{1\rm{\mu
m}}\right)^{2}\left(\frac{T}{100\rm{K}}\right)^{1/2}\bigg{]}\rm{~{}cm^{3}~{}s^{-1}}\end{split}$
(3)

So the transition from the electric potential-dominated cross-section to the
geometric cross-section is at $\sim nm$ scale. Given $T$=100K, any single
charged grain with $a_{\mathrm{Gr}}>5.7\times 10^{-8}~{}\mathrm{cm}$ needs to
consider the geometric effect when calculating the collision coefficient. This
affects the Ohmic, Hall and Pederson conductivities:

$\begin{split}&\sigma_{\mathrm{O}}=\frac{ec}{B}\sum_{j}n_{j}|Z_{j}|\beta_{j}\
,\\\
&\sigma_{\mathrm{H}}=\frac{ec}{B}\sum_{j}\frac{n_{j}Z_{j}}{1+\beta_{j}^{2}}\
,\\\
&\sigma_{\mathrm{P}}=\frac{ec}{B}\sum_{j}\frac{n_{j}|Z_{j}|\beta_{j}}{1+\beta_{j}^{2}}\
.\end{split}$ (4)

Here the summation index $j$ runs through all charged species, with $n_{j}$
and $Z_{j}$ donating, respectively, the number density and charge of
individual charged species. Using these conductivities, the general
expressions for Ohmic diffusivity $\eta_{\mathrm{O}}$, Hall diffusivity
$\eta_{\mathrm{H}}$ and Ambipolar diffusivity $\eta_{\mathrm{A}}$ – are (Bai,
2011; Wang et al., 2019).

$\begin{split}&\eta_{\mathrm{O}}=\frac{c^{2}}{4\pi}\frac{1}{\sigma_{\mathrm{O}}}\
,\quad\eta_{\mathrm{H}}=\frac{c^{2}}{4\pi}\frac{\sigma_{H}}{\sigma_{\mathrm{H}}^{2}+\sigma_{\mathrm{P}}^{2}}\
,\\\
&\eta_{\mathrm{A}}=\frac{c^{2}}{4\pi}\frac{\sigma_{\mathrm{P}}}{\sigma_{\mathrm{H}}^{2}+\sigma_{\mathrm{P}}^{2}}-\eta_{\mathrm{O}}\
,\end{split}$ (5)

The non-ideal induction equation is:

$\displaystyle\frac{\partial{\bm{B}}}{\partial
t}=\nabla\times\left({\bm{v}}\times{\bm{B}}\right)-\frac{4\pi}{c}\nabla\times(\eta_{\mathrm{O}}{\bm{J}}$
(6)
$\displaystyle+\eta_{\mathrm{H}}{\bm{J}}\times{\bm{b}}+\eta_{\mathrm{A}}{\bm{J}}_{\bot}),$

where $\bm{v}$ is gas velocity, $\bm{B}$ is magnetic field,
${\bm{b}}={\bm{B}}/|B|$ is the unit vector representing field line direction.
$\bm{J}$ is current density vector, and $\bm{J}_{\bot}$ is the current
component perpendicular to the magnetic field. Note that the ambipolar
diffusion and Ohmic dissipation terms are included in our non-ideal MHD
simulations but not the Hall term.

###  2.3 Planet Torque

We use the same approach as Hallam & Paardekooper (2017) to replace the gap-
forming planet with a gap-forming one-dimensional torque density distribution
across the two-dimensional disk along the cylindrical radial direction. The
torque density profile is from D’Angelo & Lubow (2010):

$\Lambda=-\mathcal{F}\left(x,\beta,\zeta\right)\Omega^{2}R_{0}^{2}q^{2}\left(\frac{R_{0}}{H}\right)^{4},$
(7)

where $R_{0}$ is the planet’s orbital radius, $q$ is the planet-star mass
ratio, $\Omega$ is the orbital angular velocity, $H$ is the disk scale height
at $R_{0}$, $x=(R-R_{0})/H$, $\mathcal{F}$ is a dimensionless function and
$\beta$ and $\zeta$ are the surface density and temperature radial gradients
(with a factor of -1) respectively. $\mathcal{F}$ is found by fitting the
results of three-dimensional simulations (D’Angelo & Lubow, 2010):

$\mathcal{F}(x,\beta,\zeta)=\left(p_{1}e^{\left(-\frac{(x+p_{2})^{2}}{p_{3}^{2}}\right)}+p_{4}e^{\left(-\frac{(x-p_{5})^{2}}{p_{6}^{2}}\right)}\right)\\\
\times\textrm{tanh}(p_{7}-p_{8}x),$ (8)

where the fitting parameters $p_{1}$ through $p_{8}$ are constant for a given
$\beta$, $\zeta$. We chose $\beta=0.5$ and $\zeta=1$ from the table 1 of
D’Angelo & Lubow (2010), which are the closest fit to our disk parameter
($\beta=0.57$ and $\zeta=1$). We list the fitting parameters in Table 2. The
torque profile in the arbitrary units is shown in Figure 1, together with the
impulse approximation from Lin & Papaloizou (1986). Note the torque profile
here is for the disk midplane; for anywhere above the midplane, the strength
of the torque is multiplied by a height factor
$(R-R_{0})^{2}/((R-R_{0})^{2}+z^{2})$, with $z$ being the height above the
midplane.

Table 2: Values of the parameter p in Eq. 8. $p_{n}$ | Value
---|---
$p_{1}$ | $0.0297597$
$p_{2}$ | $1.09770$
$p_{3}$ | $0.938567$
$p_{4}$ | $0.0421186$
$p_{5}$ | $0.902328$
$p_{6}$ | $1.03579$
$p_{7}$ | $0.0981183$
$p_{8}$ | $4.68108$

This profile is initially tested in a 2.5D axisymmetric inviscid hydro-only
simulation that employs the same grid structure as our non-ideal MHD model,
which we will describe later. In Figure 2, we show that the adopted torque
profile can produce a similar gas velocity pattern as 3D hydrodynamics gap
opening simulations in the meridional plane. In particular, the planet drives
a dense midplane flow away from it near the midplane. At higher altitudes,
steep density gradients drive a flow into the gap towards the planet, forming
the characteristic “collapsing flow” evident in 3D hydro simulations (e.g.,
Fung & Chiang, 2016).

Figure 1: Torque profile (solid) described by Eq. 8, with the impulse
approximation (Lin & Papaloizou, 1986) plotted as dashed line for comparison.
Figure 2: Test of planetary torque in a hydro-only simulation. Plotted are the
density distribution (color map) and velocity field (vectors) near the planet-
opened gap, showing a meridional circulation pattern with flows moving away
from the planet (located at 10 au) near the midplane and collapsing flows
towards the planet at higher altitudes.

##  3 Model Results

We run our non-ideal MHD simulation to a time of $t=1000$ years, corresponding
to about 30 planet’s orbital periods at $10$ au. As expected, a gap is
gradually opened by the planetary torque111A caveat of our planet’s torque
profile is a potential pile-up of material near the gap edges since it lacks
spiral shocks from a real gap-opening planet that may transport angular
momentum farther from the planet. Since this study mainly focuses on gas
kinematics instead of density profile, we consider the order of unity
difference of gas surface density near the gap edge acceptable., as
illustrated in Fig. 3, where the column densities are plotted at different
times (panel a), together with the effective $\alpha$ parameter corresponding
to the magnetic stresses (panel b). The effective Shakura-Sunyaev $\alpha$ is
defined as the $R-\phi$ stress to pressure ratio:

$\alpha_{int}=\left|\left.\int_{0}^{z+}{B_{R}B_{\phi}}dz\right|\middle/\int_{0}^{z+}P_{\rm
gas}dz\right.$ (9)

where the integral goes from the disk midplane to the upper disk surface
($z^{+}\sim~{}4h$). The large values of $\alpha$ inside the gap point to the
dynamical importance of the magnetic field, as discussed in detail later. Most
of our analysis will be based on the snapshot at $t=800$ years; some analysis
will be conducted at an earlier time $t=483$ years.

Figure 3: Gap opening as illustrated by the evolution of the surface density
distribution over the first 1000 years (panel a). Panel (b) plots the
“effective” $\alpha$ parameter corresponding to the magnetic stresses, which
are particularly important inside the low-density gap.  Figure 4: Properties
of the fiducial model at a representative time of 800 years on the meridian
plane. Panel (a) shows the mass density $\rho$ scaled by a reference value
$\rho_{0}$ (initial midplane density at R=10 au) and poloidal magnetic field
lines (black); Panel (b) the radial mass flux per unit area in code units,
with the velocity unit vectors (white arrows) superposed; Panel (c) the
plasma-$\beta$; and Panel (d) the Elsasser number $Am$. See the supplementary
online material for an animated version of this figure.

In Fig. 4, we plot several quantities on the meridian plane at a
representative time of 800 years. Several broad features are worth noting.
First, as expected, a clear gap is created near the planet at 10 au, as seen
in the density map shown in panel (a). Second, the wind is driven off the
disk’s surface along the initially imposed large-scale poloidal magnetic field
lines, with a generally positive radial mass flux, as shown in panel (b). The
flow pattern inside the disk is more complex, however, with patches of inflow
(with a blue color in panel [b]) next to patches of outward expansion (colored
brown). Third, the poloidal field lines have a clear concentration in the low-
density gap, as shown in panel (a). The concentration continues into the wind.
The combination of flux concentration and low density in the gap means a large
drop in the plasma-$\beta$ compared to the surrounding disk regions, as
illustrated in panel (c). The low density enables a higher ionization level,
leading to a larger Elsasser number $Am$, as shown in panel (d). Indeed, the
$Am$ value approaches or even exceeds unity, much higher than the surrounding
disk regions, where $Am$ is well below unity. Therefore, the planet created a
gap where the magnetic field is expected to be more dynamically significant
and better coupled to the gas than its surroundings.

Figure 5: Face-on view of the fiducial model at a representative time of 800
years. Panel (a) shows a color map of the column density $\Sigma$; (b) the
vertical component of the magnetic field $B_{z,mid}$ (in Gauss) on the disk
midplane; (c) the plasma-$\beta$ on the disk midplane; and (d) the Elsasser
number $Am$ on the midplane. See the supplementary material online for an
animation of the figure.

To illustrate the magnetically significant gap more pictorially, we plot in
Fig. 5 the face-on view of the distributions of the column density $\Sigma$,
the vertical component of the magnetic field $B_{z,mid}$ on the midplane, and
the plasma-$\beta$ and Elsasser number $Am$ at the same location. The figure
clearly shows the strong correlation between the column density (panel a) and
the plasma-$\beta$ (c) and their anti-correlation with the field strength (b)
and the Elsasser number (d) that were already apparent in the earlier
meridional plot (Fig. 4) It also reveals a more subtle difference between the
dense ring immediately outside the gap (around $r\sim 12.5$ au) and that
inside it ($r\sim 8$ au). The former has a much weaker midplane field strength
(panel b) and a much lower Elsasser number (panel d) than the latter. The
lower Elsasser number may appear surprising since the outer ring is less dense
compared to the inner one and thus should be better ionized by high energy
radiation and cosmic rays. However, the Elsasser number depends on not only
the charge densities but also the magnetic field strength. In particular, a
weaker field reduces the Hall parameter (see equation 2), making it harder for
the charged particles (especially the more massive charged grains) to couple
to the field lines. We suspect that a feedback loop may be operating in the
outer ring where a lower Elsasser number leads to faster diffusion of the
magnetic field from the region, lowering the field strength, which, in turn,
leads to a weaker field-matter coupling and, thus, an even lower Elsasser
number.

###  3.1 Origins of Gap Magnetic Flux Concentration

To help understand the magnetic concentration in the gap, we examine the time
evolution of the vertical field component $B_{z}$ on the midplane. Under the
assumption $B_{R}=B_{\phi}=0$ (where $R$ is the cylindrical radius), which is
strictly true for an axisymmetric disk with a mirror symmetry for the upper
and lower hemispheres, as assumed here), the induction equation (Eq.2.2) can
be reduced to:

$\displaystyle\frac{\partial{B_{z}}}{\partial
t}=-\frac{1}{R}\frac{\partial\left(Rv_{R}B_{z}\right)}{\partial
R}+\frac{1}{R}\frac{\partial\left(R\left(\eta_{O}+\eta_{A}\right)\partial_{R}B_{z}\right)}{\partial
R}.$ (10)

The right side of the equation consists of two terms: the first term
represents advection, and the second term accounts for magnetic diffusion. We
have plotted these terms separately and their sum in Figure 6. To illustrate
the trend of magnetic field evolution, the data is the average of 50 snapshots
at relatively early times, from t=100 to t=150 years. The tendency for the
vertical field strength to increase in the gap around the planet is already
clear at these early times, with the magnetic diffusion causing $B_{z}$ to
increase with time faster than the advection causing $B_{z}$ to decrease with
time. This analysis indicates that magnetic diffusion is the root cause of the
magnetic flux concentration in the gap. Therefore, its proper treatment is
needed for quantifying the degree of flux concentration.

Figure 6: The blue lines correspond to the advection effect, while the orange
lines represent magnetic diffusion. The black lines depict the net effect
($\partial B_{z}/\partial t$) of these two phenomena. Solid lines signify
positive values, and dotted lines indicate negative values. The diffusion and
advection terms are directly from Eq. 10.

###  3.2 Gas Dynamics in and around the Gap

Figure 7: Panel (a): Gas poloidal velocity divided by the local sound speed,
with arrows indicating the flow direction; (b): Effective radial mass flux (in
code units), with positive indicating outward motion and the solid black lines
showing magnetic field lines; (c) Line integration contours (LIC) of the
poloidal streamlines showing the magnetically-mediated meridional gas
circulation near the planet, the area with a velocity above 600 m/s is not
saturated by color due to the properties of LIC; (d) Toroidal magnetic field
strength (color map, in Gauss) with the field lines (black lines), velocity
field (white vectors), and an isodensity contour (thick gray line) to help
mark the low-density gap.

To illustrate the gas dynamics in and around the planet-opened gap, we plot in
Fig. 7 the poloidal flow speed (color map) and flow directions (arrows) in a
meridional plane in a zoomed-in region around the gap. Part of the flow
pattern is similar to the non-magnetized model; the gas is repelled away from
the planet near the midplane and circulates back towards the planet at higher
altitudes. However, there are several significant differences. First, near the
inner edge of the gap, there is a fast (transonic) surface accretion stream.
It originates from about one scale height above the midplane, with an initial
speed of around a quarter of the sound speed. The stream accelerates as it
rises to a higher altitude and then slides inwards along the disk surface,
reaching a supersonic speed (Figure 7a). Interestingly, the accreting stream
(with a negative radial mass flux in panel b) is located close to a sharp kink
in the poloidal magnetic field lines (panel b). The spatial correlation is not
a coincidence; it is expected because the material near the tip of the kink
that opens outward rotates faster than the material along the same field line
but at larger radii, causing it to lose angular momentum and fall inward. Part
of the angular momentum is transported out of the disk along field lines to
the low-density wind. The remaining part is transported along field lines to
the denser disk region below the accretion stream, causing it to expand. This
magnetically mediated redistribution of angular momentum 222The phenomenon is
essentially the magnetically driven surface “avalanche flow” described in
Matsumoto et al. (1996) (see also Suriano et al., 2017, and references
therein). It is closely related to the well-known channel flows in early
simulations of the magneto-rotational instability (e.g., Hawley & Balbus,
1992). can be seen most clearly in Panel (d) of the Figure, where we plot the
toroidal magnetic field in the color map. Note the negative $B_{\phi}$ in the
accretion stream, indicating that it is magnetically braked. In contrast, the
region below the stream has a positive $B_{\phi}$, indicating that the
magnetic field is trying to spin up the material, causing it to expand
outward. The thermal pressure gradient can aid this outward expansion. The
positive magnetic torque and pressure gradient combine to overcome the
persistent negative planetary torque in the region that tends to remove its
angular momentum and drive it to accrete. Similarly, the surface accretion
stream may be helped by the negative planetary torque, which aids the negative
magnetic torque in removing its angular momentum and thus drives the stream to
accrete faster than it would otherwise. The net effect of the interplay
between the magnetic and planetary torques is the production of fast
meridional motions above the circulation cell driven directly by the planet
close to the midplane (see panel [c]).

The surface accretion stream and associated meridional circulations are
persistent features of our simulations (see the animated version of Figure 7).
It is not limited to the inner edge of the gap; it often occurs near the outer
edge as well, especially at earlier times. For example, we show in Fig. 8 the
same plots as in Fig. 7 but for an earlier time of $t=483$ years. Note the
double-kinked, “S-shaped” poloidal field line near the radius of 11 au outside
the planet, with a fast, nearly sonic infall around the upper field line kink
that opens outward and a fast expansion around the lower kink that opens
inward. As discussed earlier, this is the expected behavior of magnetically
mediated angular momentum redistribution, which can drive meridional gas
circulation up to the sound speed near both the inner and outer edges of the
planet-opened gap (see panel [a]).

To illustrate the role of magnetic fields on the gas meridional circulation
further, we plot in Fig. 9 the distribution of the magnetic torque (panel a)
at the same time shown in Fig. 8. As expected, regions near the gap where the
magnetic field is braking the disk rotation (with a negative torque) tend to
infall, while those with a positive magnetic torque tend to expand (compare
Fig. 9a with Fig. 8b). The correlation is not perfect, though, because of the
presence of a thermal pressure gradient and a persistent planetary torque
(Fig. 9b), which dominates the magnetic torque near the gap (inside the thick
black contours of Fig. 9c), particularly in the relatively high-density
regions close to the midplane inside the gap and near the inner and outer
edges of the gap and the surrounding rings. However, even in the regions with
higher planetary torque, the magnetic torque remains significant, with a value
typically smaller by a factor of only a few compared to the planetary torque
(see panel d). Interestingly, even in such regions, there are fast (sonic)
flows toward the planet despite the tendency for the planetary torque to push
them away from the planet. As mentioned earlier, these flows are likely helped
by strong pressure gradients near the gap edges, as in the hydro case.

It is worth noting that, at the (relatively early) time shown in Fig. 9, the
magnetic torque dominates the planetary torque in a roughly “V-shaped” region
directly above the planet (see panel c) because of a relatively low density
and strong magnetic field in this part of the gap (see Fig. 9). It is a region
completely dominated by the magnetic field, with a plasma-$\beta$ below unity
and an Elsasser number above unity. As the density in the gap continues to
drop at later times, this magnetically-dominated region extends all the way to
the disk midplane, forming a roughly “U-shaped” region that centers not on the
planet but to a location interior to it, where most of the poloidal magnetic
flux in the gap is concentrated (see, e.g., the time $t=1000$ years in the
animated version of Fig. 9 online).

The poloidal magnetic flux concentration in the gap results from the non-ideal
MHD effects (ambipolar diffusion and Ohmic dissipation) included in our
simulation. The reason is that the flux-to-mass ratio is conserved in the
ideal MHD limit, with a lower flux expected in the strongly mass-depleted gap.
This argument is consistent with our analysis in §3.1, although a complete
understanding of exactly how the needed redistribution of the magnetic flux
relative to the mass occurs is still lacking. From the evolution of the
poloidal magnetic field lines, we observed that some field lines are dragged
into the gap from larger radii by accretion flows and get stuck. Why the
dragged-in field lines stay in the gap is unclear, but it is consistent with
other non-ideal MHD disk simulations that include ambipolar diffusion where
the poloidal field strength strongly anti-correlates with the surface density
(e.g., Suriano et al., 2018; Riols et al., 2020; Cui & Bai, 2021; Hu et al.,
2022; Aoyama & Bai, 2023).

Figure 8: Same as Fig. 7 but at an earlier time of 483 years, showing
prominent magnetically-driven fast (sonic) meridional gas circulations near
both inner and outer edges of the planet-opened gap.  Figure 9: Torque (in
code units) distribution comparison at t=483. Panel (a) and (b) are magnetic
torque $\Lambda_{m}$ and planet torque $\Lambda_{p}$ plotted in the same
symmetric log scale; (c) shows the ratio of planet torque over magnetic
torque, and (d) illustrated the absolute value of the same ratio. The field
lines are black solid lines in all panels but (b), where they are color-coded
to indicate the magnetic flux within (to the inner boundary) each line and
have an offset to the black lines in the other panels. See the supplementary
material online for an animation of the figure.

###  3.3 Magnetic Diffusivity-Gas Density Relation

A major advantage of our 2D (axisymmetric) non-ideal MHD simulation enabled by
the prescribed planetary torque is that it allows us to evolve the disk
thermochemistry (including ionization) self-consistently. In contrast, most
existing disk simulations that include nonideal MHD processes either use
spatially dependent only diffusion profiles (e.g., Bai & Stone, 2017; Cui &
Bai, 2021) or simple power-laws that depend only on the local density (e.g.,
Suriano et al., 2018; Hu et al., 2022) or column density (e.g., Riols et al.,
2020). In this section, we approximate the numerically obtained magnetic
diffusivities as a function of the local density using analytic fitting
formulae. Such formulae will be useful for future 3D studies where it is too
computationally expensive to compute the thermochemistry self-consistently in
each cell and at each time step.

To illustrate the dependence of the magnetic diffusivities on the local
density, we plot in Figure 10 the time evolution (up to 1000 years) of the
disk mass density $\rho$, the Elsasser number $Am$, and the Ohmic diffusivity
$\eta_{O}$ as a function of the polar angle $\theta_{mid}$ away from the mid-
plane (in units of the ratio of the scale height $h$ over the radius $r$) at
the radius of $r=10$ au (where the planet is located). The sharp change of
$Am$ and $\eta_{O}$ at $\sim 4$ disk scale heights corresponds to the height
of UV penetration. The Ohmic diffusivity shown in Panel (c) appears to
correlate with the gas density shown in Panel (a), with the diffusivity
increasing nearly monotonically with the density. The change of $Am$ with the
density appears more complex. While $Am$ generally anti-correlates with
$\rho$, there is an $Am$ “valley” just below the disk surface, starting from a
height of $\sim$3h and reaching below $\sim$2h after 800 years. This trend
resembles the power-law prescription proposed by Suriano et al. (2018) where
the $Am$ inside the disk decreases with a decreasing density as a result of
balancing the recombination of the dominant ions and electrons with the cosmic
ray ionization.

To better understand the magnetic diffusion-density relation, we plot the
$Am-\rho$ and $\eta_{O}-\rho$ diagrams in Figure 11. Each track in the
diagrams corresponds to the time evolution of $Am$ and $\rho$ or $\eta_{O}$
and $\rho$ at a given height at $r=10$ au, with the height denoted by the
track’s color (the bluest corresponds to the midplane and the reddest marks
$\sim$ 5.5h, with an increment of 0.25 scale height). The Ohmic diffusivity in
Panel (b) has a simple dependence on the gas density: almost a single
monotonic track. It means $\eta_{O}$ barely depends on its vertical location.
We used two polynomials, $\mathcal{P}_{\rm O,1}$ and $\mathcal{P}_{\rm O,2}$,
to fit the Ohmic-density relation in the log-log space, with the coefficients
$k_{i}$ of the polynomials ($\mathcal{P}=\sum_{i=0}^{n}k_{i}x^{i}$) given in
Table 4).

The $Am-\rho$ diagram is a different story. In Panel (a) of Figure 11, the
tracks representing different heights from the midplane only start to converge
into a single track when the density drops below the initial midplane density
by about two orders of magnitude. The multiple ion species and charged grains
acting like “heavy ions” give the ambipolar diffusion much more complicated
behaviors than the Ohmic diffusivity. To have a better-performing look-up
table of $Am-\rho$, we ignored the $Am$ evolution in the first 100 yrs. The
fitting result is given in Table 3 and shown in Figure 11 as dashed lines.
This table333The table is available in both Python and C scripts at
https://github.com/astroxhu/diffusion-table could be used for future 3D
protoplanetary disk simulations covering a high dynamic gas density range,
especially near the midplane.

Figure 10: Time evolution of the gas density (panel a), $Am$ (panel b), and
$\eta_{O}$ (panel c) at different heights at the radius $r=10$ au. The white
contours in panel (b) indicate $Am=1$. Figure 11: Panel (a) shows the relation
between density and $Am$, while Panel (b) shows the relation between density
and $\eta_{O}$. The measurement is taken from the midplane (blue lines) to
$\sim 5$ scale heights (red lines) at an increment of 0.25 scale height. Table
3: Am-$\rho$ look-up table (§3.3)

$\log_{10}{\rho^{*}}$ | $\theta^{*}$ | $\log_{10}{Am}$
---|---|---
-0.5 to 0. | 0 to 4. | $(-2.04\theta^{*}_{1}-0.656)\times\log_{10}{\rho^{*}}-1.768$
-1.05 to -0.5 | 0 to 0.8 | $(0.945\theta^{*}-1.709)\times\log_{10}{\rho^{*}}+1.493\theta^{*}-2.295$
-1.05 to -0.5 | 0.8 to 4 | $(1.127\theta^{*}_{1}-0.953)\times\log_{10}{\rho^{*}}+\ 1.584\theta^{*}_{1}-1.1$
$\mathcal{F}(\log_{10}{\rho^{*}},\log_{10}{Am})$ to -1.05 | 0 to 0.8 | $(0.121-0.227\theta^{*}_{1})/(0.338-0.11\theta^{*}_{1})\times\log_{10}{\rho^{*}}+\ 0.5\theta^{*}_{1}+1.05\times(0.121-0.227\theta^{*}_{1})/(0.338-0.11\theta^{*}_{1})-0.5$
$\mathcal{F}(\log_{10}{\rho^{*}},\log_{10}{Am})$ to -1.05 | 0 to 0.8 | $(0.0756\theta^{*}_{2}-0.0606)/(0.25-0.0566\theta^{*}_{2})\times\log_{10}{\rho^{*}}+\ 0.4\theta^{*}_{2}-0.1+1.05\times(0.0756\theta^{*}_{2}-0.0606)/(0.25-0.0566\theta^{*}_{2})$
-2.5 to $\mathcal{F}(\log_{10}{\rho^{*}},\log_{10}{Am})$ | 0 to 0.8 | $(0.537\theta^{*}-0.702)/(0.587-0.421\theta^{*})\times\log_{10}{\rho^{*}}+\ 0.727\theta^{*}-0.621$
|  | $-(0.11\theta^{*}-1.388)*(0.537\theta^{*}-0.702)/(0.587-0.421\theta^{*})$
-2.5 to $\mathcal{F}(\log_{10}{\rho^{*}},\log_{10}{Am})$ | 0.8 to 4 | $(0.2864\theta^{*}_{2}-0.2724)/(0.2502-0.0496\theta^{*}_{2})\times\log_{10}{\rho^{*}}+\ 0.3244\theta^{*}_{2}-0.0394$
|  | $-(0.0566\theta^{*}_{2}-1.3)\times(0.2864\theta^{*}_{2}-0.2724)/(0.2502-0.0496\theta^{*}_{2})$
-2.5 to -0.5 | 0 to 4 |  Upper Limit : $\mathcal{P}_{\rm AD,1}$, see Table 4
$\log_{10}{\rho^{*}_{atm}}$ to -2.5 | 0 to 4 |  $\mathcal{P}_{\rm AD,2}$, see Table 4

Here $\rho^{*}=\rho/\rho_{\rm mid,init}$, $\theta^{*}=\theta/(h/r)$,
$\theta^{*}_{1}=\mathrm{min}{(\theta^{*},1)}$,
$\theta^{*}_{2}=\mathrm{min}{(\theta^{*},1.95)}.$
---

Table 4: Coefficients for polynomial fits of AD and Ohmic look-up table (§3.3)  $k_{i}$ | $\mathcal{P}_{\rm AD,1}$ | $\mathcal{P}_{\rm AD,2}$ | $\mathcal{P}_{\rm O,1}$ | $\mathcal{P}_{\rm O,2}$
---|---|---|---|---
$k_{6}$ | 0.26429362 | 0.134528857 |  |
$k_{5}$ | 2.308469 | 3.15915033 |  |
$k_{4}$ | 7.81413499 | 30.6754303 |  | -0.82772958
$k_{3}$ | 12.39250353 | 157.055386 |  | -11.48813613
$k_{2}$ | 7.90442493 | 445.96139 | -0.10996372 | -58.24798054
$k_{1}$ | -0.52179895 | 664.2113 | 0.49186494 | -124.86627111
$k_{0}$ | -1.52831614 | 404.541785 | 16.82894602 | -80.52351483

##  4 Discussion and Conclusions

Our consistent thermochemical modeling of the planet’s gap opening allows us
to determine various chemical species’ spatial distributions and evolution in
the disk and its surrounding environment. Of the modeled species, the most
directly relevant to observations are CO and HCO+. Their number densities are
shown in Fig. 12. As expected, the CO number density is lower in the gap than
in the denser surrounding regions because of gas depletion. However, the HCO+
number density is higher in the gap than in the surrounding regions despite a
lower gas density, indicating a much higher fractional abundance. The enhanced
HCO+ abundance is consistent with the ALMA observations of AS 209, one of the
protoplanetary disks with prominent gaps potentially opened by planets, which
show that HCO+ has an unusually high abundance in the radial range between
$\sim 50$ and $\sim 150$ au where the CO column density shows a depression
(Aikawa et al., 2021). It is also in agreement with the greatly enhanced DCO+
in the same region (Favre et al., 2019, see their Fig. 4).

Figure 12: Distributions of the number densities of CO (Panel a) and HCO+
(Panel b) at a representative time $t=800~{}$years, showing enhanced HCO+ in
the planet-opened gap where the CO density is reduced along with the gas.
Figure 13: Gas velocity structure presented the same general style as in Fig.
2 of Teague et al. (2019) to facilitate comparison with observations. Gas
velocity vectors in the meridional plane at the surface, with a vector for 0.1
$c_{s}$ shown for reference. The color map in the background is the gas
density.

To better connect with observations, we follow Teague et al. (2019) and plot
the gas velocity structure in $R-z$ plane in Figure 13, where $R$ is the
cylindrical radius and $z$ the height above the midplane. For illustration
purposes, we picked a surface that has a constant column density above it
(from the top of the simulation domain). This figure is to be compared with
the Fig. 4 of Galloway-Sprietsma et al. (2023) for the AS 209 disk. Since our
simulation was not explicitly designed for the AS 209 system, a detailed match
is not to be expected. Nevertheless, we find velocity vectors pointing upward
and inward at the gap’s inner edge, which resembles the nearly sonic
meridional flow observed in the ${}^{12}{\rm CO}$-emitting region of the AS
209 disk. However, the current model does not capture the outflow in the low-
density gap region probed by 13CO, possibly because the adopted planetary
torque profile makes the gap region too evacuated (and thus too strongly
magnetically dominated) to launch an outflow efficiently.

Our work is most directly comparable to Aoyama & Bai (2023), who carried out
3D global simulations of gap opening by planets in non-ideal MHD disks with a
spatially constant ambipolar Elsasser number (either 1 or 3) inside the disk.
They find a concentration of the poloidal magnetic flux in the planet-opened
gap, which they believe is likely associated with the spiral density shocks at
the surface of their simulated disks. We also find a concentration of poloidal
magnetic flux in the planet-opened gap, despite a different (and more self-
consistent) treatment of the ambipolar diffusivity, indicating that it is a
rather general phenomenon. In particular, it does not depend on spiral density
shocks, which are absent in our 2D (axisymmetric) simulation. We believe it is
caused by the magnetic flux diffusing into the gap and being trapped because
of enhanced local magnetic coupling to the lower-density gas. However, the
details of the process remain to be determined.

In summary, we carried out 2D (axisymmetric) simulations of gap opening by a
planet embedded in a wind-launching non-ideal MHD disk with a prescribed
planetary torque but self-consistent thermochemistry. Our main conclusions are
as follows:

1\. There is a strong concentration of poloidal magnetic flux in the planet-
opened gap, where the magnetic field is much better coupled to the gas than in
the denser regions surrounding the gap and where the magnetic pressure becomes
comparable to, or even larger than, the gas pressure. Combined with similar
results recently obtained by Aoyama & Bai (2023) under very different
conditions, our finding indicates that the flux concentration and the
resulting magnetic domination of the gas dynamics are robust features of the
planet-opened gaps in non-ideal MHD disks.

2\. Magnetic fields are also dynamically important in the denser regions
surrounding the planet-opened gap. They drive fast infall and expansion
motions up to the local sound speed through magnetic torque-induced angular
momentum redistribution, thereby controlling the meridional gas circulation
near the inner and outer edges of the gap.

3\. The plant-opened gap has a much higher abundance of molecular ion HCO+
than its denser surrounding regions, consistent with high-resolution ALMA
observations of the protoplanetary disk AS 209 with prominent rings and gaps.
The magnetically-induced fast accretion stream near the disk surface at the
inner edge of the gap is also consistent with sonic meridional flows probed by
the AS 209 12CO observations.

4\. We numerically obtained fitting formulae for the ambipolar and Ohmic
diffusivities as a function of the local density that can be used for future
3D simulations of planet gap-opening in non-ideal MHD disks where
thermochemistry is too computationally expensive to evolve self-consistently
with the magneto-hydrodynamics.

## Acknowledgements

X.H. acknowledges support from the University of Virginia through VICO
(Virginia Initiative on Cosmic Origins) and NSF AST-1815784. ZYL is supported
in part by NASA 80NSSC20K0533 and NSF AST-1910106. Our simulations are made
possible by an XSEDE allocation (AST200032). L.W. acknowledges the computation
resources provided by the KIAA. Z.Z. acknowledges support from the National
Science Foundation under CAREER Grant Number AST-1753168. Figures in this
paper were made with the help of `Matplotlib` (Hunter, 2007) and `NumPy`
(Harris et al., 2020).

## Data Availability

The data from the simulations will be shared on reasonable request to the
corresponding authors.

## References

  * Aikawa et al. (2021) Aikawa Y., et al., 2021, ApJS, 257, 13
  * Aoyama & Bai (2023) Aoyama Y., Bai X.-N., 2023, ApJ, 946, 5
  * Armitage (2010) Armitage P. J., 2010, Astrophysics of Planet Formation. Cambridge, UK: Cambridge University Press
  * Bai (2011) Bai X.-N., 2011, ApJ, 739, 50
  * Bai (2014) Bai X.-N., 2014, ApJ, 791, 72
  * Bai & Stone (2013) Bai X.-N., Stone J. M., 2013, ApJ, 769, 76
  * Bai & Stone (2017) Bai X.-N., Stone J. M., 2017, ApJ, 836, 46
  * Baruteau et al. (2011) Baruteau C., Fromang S., Nelson R. P., Masset F., 2011, A&A, 533, A84
  * Baruteau et al. (2014) Baruteau C., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. pp 667–689 (arXiv:1312.4293), doi:10.2458/azu˙uapress˙9780816531240-ch029
  * Bitsch et al. (2018) Bitsch B., Morbidelli A., Johansen A., Lega E., Lambrechts M., Crida A., 2018, A&A, 612, A30
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  * Carballido et al. (2017) Carballido A., Matthews L. S., Hyde T. W., 2017, MNRAS, 472, 3277
  * Crida et al. (2006) Crida A., Morbidelli A., Masset F., 2006, Icarus, 181, 587
  * Cui & Bai (2021) Cui C., Bai X.-N., 2021, MNRAS, 507, 1106
  * D’Angelo & Lubow (2010) D’Angelo G., Lubow S. H., 2010, ApJ, 724, 730
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. Princeton University Press
  * Duffell & MacFadyen (2013) Duffell P. C., MacFadyen A. I., 2013, ApJ, 769, 41
  * Favre et al. (2019) Favre C., et al., 2019, ApJ, 871, 107
  * Fung & Chiang (2016) Fung J., Chiang E., 2016, ApJ, 832, 105
  * Galloway-Sprietsma et al. (2023) Galloway-Sprietsma M., et al., 2023, arXiv e-prints, p. arXiv:2304.03665
  * Goldreich & Tremaine (1979) Goldreich P., Tremaine S., 1979, ApJ, 233, 857
  * Gressel et al. (2013) Gressel O., Nelson R. P., Turner N. J., Ziegler U., 2013, ApJ, 779, 59
  * Hallam & Paardekooper (2017) Hallam P. D., Paardekooper S. J., 2017, MNRAS, 469, 3813
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hawley & Balbus (1992) Hawley J. F., Balbus S. A., 1992, ApJ, 400, 595
  * Hu et al. (2018) Hu X., Tan J. C., Zhu Z., Chatterjee S., Birnstiel T., Youdin A. N., Mohanty S., 2018, ApJ, 857, 20
  * Hu et al. (2019) Hu X., Zhu Z., Okuzumi S., Bai X.-N., Wang L., Tomida K., Stone J. M., 2019, ApJ, 885, 36
  * Hu et al. (2021) Hu X., Wang L., Okuzumi S., Zhu Z., 2021, ApJ, 913, 133
  * Hu et al. (2022) Hu X., Li Z.-Y., Zhu Z., Yang C.-C., 2022, MNRAS, 516, 2006
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ida & Lin (2004) Ida S., Lin D. N. C., 2004, ApJ, 616, 567
  * Kanagawa et al. (2015) Kanagawa K. D., Muto T., Tanaka H., Tanigawa T., Takeuchi T., Tsukagoshi T., Momose M., 2015, ApJ, 806, L15
  * Kley & Nelson (2012) Kley W., Nelson R. P., 2012, ARA&A, 50, 211
  * Lega et al. (2014) Lega E., Crida A., Bitsch B., Morbidelli A., 2014, MNRAS, 440, 683
  * Lin & Papaloizou (1986) Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846
  * Lin & Papaloizou (1993) Lin D. N. C., Papaloizou J. C. B., 1993, in Levy E. H., Lunine J. I., eds, Protostars and Planets III. p. 749
  * Masset (2017) Masset F. S., 2017, MNRAS, 472, 4204
  * Matsumoto et al. (1996) Matsumoto R., Uchida Y., Hirose S., Shibata K., Hayashi M. R., Ferrari A., Bodo G., Norman C., 1996, ApJ, 461, 115
  * Mordasini et al. (2012) Mordasini C., Alibert Y., Georgy C., Dittkrist K. M., Klahr H., Henning T., 2012, A&A, 547, A112
  * Nelson & Papaloizou (2003) Nelson R. P., Papaloizou J. C. B., 2003, MNRAS, 339, 993
  * Pinte et al. (2020) Pinte C., et al., 2020, ApJ, 890, L9
  * Riols et al. (2020) Riols A., Lesur G., Menard F., 2020, A&A, 639, A95
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4
  * Suriano et al. (2017) Suriano S. S., Li Z.-Y., Krasnopolsky R., Shang H., 2017, MNRAS, 468, 3850
  * Suriano et al. (2018) Suriano S. S., Li Z.-Y., Krasnopolsky R., Shang H., 2018, MNRAS, 477, 1239
  * Teague et al. (2019) Teague R., Bae J., Bergin E. A., 2019, Nature, 574, 378
  * Wang et al. (2019) Wang L., Bai X.-N., Goodman J., 2019, ApJ, 874, 90
  * Winters et al. (2003) Winters W. F., Balbus S. A., Hawley J. F., 2003, ApJ, 589, 543
  * Zhu et al. (2013) Zhu Z., Stone J. M., Rafikov R. R., 2013, ApJ, 768, 143