# DCMN+: Dual Co-Matching Network for Multi-choice Reading Comprehension

Shuailiang Zhang,1,2,3 Hai Zhao,1,2,3 Yuwei Wu,1,2,3 Zhuosheng Zhang,1,2,3
Xi Zhou,4 Xiang Zhou4
1Department of Computer Science and Engineering, Shanghai Jiao Tong University
2Key Laboratory of Shanghai Education Commission for Intelligent Interaction
and Cognitive Engineering, Shanghai Jiao Tong University, Shanghai, China
3MoE Key Lab of Artificial Intelligence, AI Institute, Shanghai Jiao Tong
University, Shanghai, China
4CloudWalk Technology, Shanghai, China
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Corresponding author. This paper was partially supported by National Key
Research and Development Program of China (No. 2017YFB0304100), Key Projects
of National Natural Science Foundation of China (U1836222 and 61733011).

###### Abstract

Multi-choice reading comprehension is a challenging task to select an answer
from a set of candidate options when given passage and question. Previous
approaches usually only calculate question-aware passage representation and
ignore passage-aware question representation when modeling the relationship
between passage and question, which cannot effectively capture the
relationship between passage and question. In this work, we propose dual co-
matching network (DCMN) which models the relationship among passage, question
and answer options bidirectionally. Besides, inspired by how humans solve
multi-choice questions, we integrate two reading strategies into our model:
(i) passage sentence selection that finds the most salient supporting
sentences to answer the question, (ii) answer option interaction that encodes
the comparison information between answer options. DCMN equipped with the two
strategies (DCMN+) obtains state-of-the-art results on five multi-choice
reading comprehension datasets from different domains: RACE, SemEval-2018 Task
11, ROCStories, COIN, MCTest.

## Introduction

Machine reading comprehension (MRC) is a fundamental and long-standing goal of
natural language understanding which aims to teach machine to answer question
automatically according to given passage (?; ?; ?; ?). In this paper, we focus
on multi-choice MRC tasks such as RACE (?) which requests to choose the right
option from a set of candidate answers according to given passage and
question. Different from MRC datasets such as SQuAD (?) and NewsQA (?) where
the expected answer is usually in the form of a short span from the given
passage, answer in multi-choice MRC is non-extractive and may not appear in
the original passage, which allows rich types of questions such as commonsense
reasoning and passage summarization, as illustrated by the example in Table 1.

Passage: Runners in a relay race pass a stick in one direction. However,
merchants passed silk, gold, fruit, and glass along the Silk Road in more than
one direction. They earned their living by traveling the famous Silk Road. …
The Silk Road was made up of many routes, not one smooth path. They passed
through what are now 18 countries. The routes crossed mountains and deserts
and had many dangers of hot sun, deep snow and even battles…
---
Question: The Silk Road became less important because _ .
A. it was made up of different routes
B. silk trading became less popular
C. sea travel provided easier routes
D. people needed fewer foreign goods
Table 1:  An example passage with related question and options from RACE
dataset. The ground-truth answer and the evidence sentences in the passage are
in bold. Figure 1: The framework of our model. P-Passage, Q-Question,
O-Option.

Pre-trained language models such as BERT (?) and XLNet (?) have achieved
significant improvement on various MRC tasks. Recent works on MRC may be put
into two categories, training more powerful language models or exploring
effective applying pattern of the language models to solve specific task.
There is no doubt that training a better language model is essential and
indeed extremely helpful (?; ?) but at the same time it is time-consuming and
resource-demanding to impart massive amounts of general knowledge from
external corpora into a deep language model via pre-training (?; ?). For
example, training a 24-layer transformer (?) requires 64 TPUs for 4 days. So
from the practical viewpoint, given limited computing resources and a well
pre-trained model, can we improve the machine reading comprehension during
fine-tuning instead of via expensive full pre-training? This work starts from
this viewpoint and focuses on exploring effective applying pattern of language
models instead of presenting better language models to furthermore enhance
state-of-the-art multi-choice MRC. We will show the way to use a strong pre-
trained language model may still have a heavy impact on MRC performance no
matter how strong the language model itself is.

To well handle multi-choice MRC problem, an effective solution has to
carefully model the relationship among the triplet of three sequences, passage
(P), question (Q) and answer candidate options (A) with a matching module to
determine the answer. However, previous unidirectional matching strategies
usually calculate question-aware passage representation and ignore passage-
aware question representation when modeling the relationship between passage
and question (?; ?; ?).

To alleviate such an obvious defect in modeling the {P, Q, A} triplet from
existing work, we propose dual co-matching network (DCMN) which
bidirectionally incorporates all the pairwise relationships among the {P, Q,
A} triplet. In detail, we model the passage-question, passage-option and
question-option pairwise relationship simultaneously and bidirectionally for
each triplet, exploiting the gated mechanism to fuse the representations from
two directions. Besides, we integrate two reading strategies which humans
usually use into the model. One is passage sentence selection that helps
extract salient evidence sentences from the given passage, and then matches
evidence sentences with answer options. The other is answer option interaction
that encodes comparison information into each option. The overall framework is
shown in Figure 1. The output of pre-trained language model (i.e. BERT (?) and
XLNet (?)) is used as the contextualized encoding. After passage sentence
selection and answer option interaction, bidirectional matching
representations are built for every pairwise relationship among the {P, Q, A}
triplet.

Our model achieves new state-of-the-art results on the multi-choice MRC
benchmark challenge RACE (?). We further conduct experiments on four
representative multi-choice MRC datasets from different domains (i.e.,
ROCStories (?), SemEval-2018 Task 11 (?), MCTest (?), COIN Shared Task 1 (?))
and achieve the absolute improvement of 4.9% and 2.8% in average accuracy from
directly fine-tuned BERT and XLNet, respectively, which indicates our method
has a heavy impact on the MRC performance no matter how strong the pre-trained
language model itself is.

## Our Proposed Model

The illustration of our model is shown in Figure 1. The major components of
the model are Contextualized Encoding, Passage Sentence Selection, Answer
Option Interaction and Bidirectional Matching. We will discuss each component
in detail.

### Task Definition

For the task of multi-choice reading comprehension, the machine is given a
passage (P), a question (Q), and a set of answer candidate options (A) to
select the correct answer from the candidates, where
$\textbf{P}=\\{\bf{p_{1},p_{2},...,p_{n}}\\}$ is the passage composed of $n$
sentences, $\textbf{A}=\\{\bf{A_{1},A_{2},...,A_{m}}\\}$ is the option set
with $m$ answer candidates.

N sent |  % on RACE |  % on COIN |  Passage |  Question
---|---|---|---|---
1 |  65 |  76 |  Soon after, the snack came out. I then opened the chips and started to enjoy them, before enjoying the soda. I had a great little snack… |  What else did the person enjoy?
2 |  22 |  10 |  She lived in the house across the street that had 9 people and 3 dogs and another cat living in it. She didn’t seem very happy there, especially with a 2 year old that chased her and grabbed her. The other people in the house agreed… |  What did the 2 year old’s mom own?
$>=$3 |  13 |  14 |  When I was hungry last week for a little snack and a soda, I went to the closest vending machine. I felt that it was a little overpriced, but being as though I needed something… |  What’s the main idea of this passage?

Table 2:  Analysis of the sentences in passage required to answer questions on
RACE and COIN. 50 examples from each dataset are sampled randomly. N sent
indicates the number of sentences required to answer the question. The
evidence sentences in the passage are in emphasis and the correct answer is
with bold.

### Contextualized Encoding

In this work, pre-trained language models are used as the encoder of our model
which encodes each token in passage and question into a fixed-length vector.
Given an encoder, the passage, the question, and the answer options are
encoded as follows:

$\begin{split}\textbf{H}^{p}=Encode(\textbf{P})&,\textbf{H}^{q}=Encode(\textbf{Q})\\\
\textbf{H}^{a}=&Encode(\textbf{A})\end{split}$ (1)

where $Encode(\cdot)$ returns the last layer output by the encoder, which can
be well pre-trained language models such as BERT (?) and XLNet (?), as using
transformer as the contextualized encoder has shown to be very powerful in
language representation (?; ?; ?; ?). $\textbf{H}^{p}\in R^{|P|\times l}$,
$\textbf{H}^{q}\in R^{|Q|\times l}$, and $\textbf{H}^{a}\in R^{|A|\times l}$
are sequence representation of passage, question and answer option,
respectively. $|P|$, $|Q|$, $|A|$ are the sequence length, respectively. $l$
is the dimension of the hidden state.

### Passage Sentence Selection

Existing multi-choice MRC models learn the passage representation with all the
sentences in one-shot, which is inefficient and counter-intuitive. To explore
how many sentences are necessarily required to answer the question, we
randomly extract 50 examples from the development set of RACE and COIN, as
shown in Table 2. Among all examples, 87% questions on RACE and 86% on COIN
can be answered within two evidence sentences. From this observation, the
model should be extremely beneficial if focusing on a few key evidence
sentences.

To select the evidence sentences from the passage
$\textbf{P}=\\{\bf{p_{1}},p_{2},..,p_{i},..,p_{n}\\}$, this module scores each
sentence $p_{i}$ with respect to the question Q and answer option A in
parallel. The top $K$ scored sentences will be selected. This module shares
the encoder with the whole model. For each
{$\bf{p_{i}},\textbf{Q},\textbf{A}$} triplet, $\textbf{H}^{p_{i}}\in
R^{|p_{i}|\times l}$, $\textbf{H}^{q}$, and $\textbf{H}^{a}$ are all
representations offered by the encoder. Here we introduce two methods to
compute the score of the triplet based on the representations.

  * •

Cosine score: The model computes word-by-word cosine similarity between the
sentence and question-option sequence pair.

$\begin{split}\textbf{D}^{pa}&=Cosine(\textbf{H}^{a},\textbf{H}^{p_{i}})\in
R^{|A|\times|p_{i}|}\\\
\textbf{D}^{pq}&=Cosine(\textbf{H}^{q},\textbf{H}^{p_{i}})\in
R^{|Q|\times|p_{i}|}\\\ \bar{\textbf{D}}^{pa}&=MaxPooling(\textbf{D}^{pa})\in
R^{|A|}\\\ \bar{\textbf{D}}^{pq}&=MaxPooling(\textbf{D}^{pq})\in R^{|Q|}\\\
score&=\frac{\sum_{k=1}^{|A|}\bar{\textbf{D}}^{pa}_{k}}{|A|}+\frac{\sum_{k=1}^{|Q|}\bar{\textbf{D}}^{pq}_{k}}{|Q|}\end{split}$
(2)

where $\textbf{D}^{pa}$, $\textbf{D}^{pa}$ are the distance matrices and
$\textbf{D}^{pa}_{ij}$ is the cosine similarity between the $i$-th word in the
candidate option and the $j$-th word in the passage sentence.

  * •

Bilinear score: Inspired by (?), we compute the bilinear weighted distance
between two sequences, which can be calculated as follows:

$\begin{split}\alpha&=SoftMax(\textbf{H}^{q}W_{1})\in R^{|Q|\times l}\\\
\textbf{q}&=\alpha^{T}\textbf{H}^{q}\in R^{l}\\\
\bar{\textbf{P}}_{j}&=\textbf{H}^{p_{i}}_{j}W_{2}\textbf{q}\in
R^{l},j\in[1,|p_{i}|]\\\
\hat{\textbf{P}}^{pq}&=Max(\bar{\textbf{P}}_{1}\bar{\textbf{P}}_{2},...,\bar{\textbf{P}}_{|p_{i}|})\in
R^{l}\end{split}$ (3)

where $W_{1}$, $W_{2}\in R^{l\times l}$ are learnable parameters,
$\hat{\textbf{P}}^{pq}$ is the bilinear similarity vector between the passage
sentence and question. Similarly, the vector $\hat{\textbf{P}}^{pa}$ between
the passage sentence and answer can be calculated with the same procedure. The
final score can be computed as follows:

$score=W_{3}^{T}\hat{\textbf{P}}^{pq}+W_{4}^{T}\hat{\textbf{P}}^{pa}$ (4)

where $W_{3}$, $W_{4}\in R^{l}$ are learnable parameters.

After scoring each sentence, top $K$ scored sentences are selected and
concatenated together as an updated passage $\textbf{P}_{s}$ to replace
original full passage. So the new sequence triplet is {$\textbf{P}_{s}$, Q, A}
and the new passage is represented as $\textbf{H}^{ps}$.

### Answer Option Interaction

Human solving multi-choice problem may seek help from comparing all answer
options. For example, one option has to be picked up not because it is the
most likely correct, but all the others are impossibly correct. Inspired by
such human experience, we introduce the comparison information among answer
options so that each option is not independent of the other. Here we build
bilinear representations between any two options. Gated mechanism (?) is used
to fuse interaction representation into the original answer option
representations.

The encoder encodes each answer option $\textbf{A}_{i}$ as
$\textbf{H}^{a_{i}}$. Then the comparison vector between option
$\textbf{A}_{i}$ and $\textbf{A}_{j}$ can be computed as follows:

$\begin{split}\textbf{G}&=SoftMax(\textbf{H}^{a_{i}}W_{5}{\textbf{H}^{a_{j}}}^{T})\in
R^{|A_{i}|\times|A_{j}|}\\\
\textbf{H}^{a_{i,j}}&=ReLU(\textbf{G}\textbf{H}^{a_{j}})\in R^{|A_{i}|\times
l}\end{split}$ (5)

where $W_{5}\in R^{l\times l}$ is one learnable parameter, G is the bilinear
interaction matrix between $A_{i}$ and $A_{j}$, $\textbf{H}^{a_{i,j}}$ is the
interaction representation. Then gated mechanism is used to fuse interaction
representation into the original answer option representations as follows:

$\begin{split}\hat{\textbf{H}}^{a_{i}}&=[\\{\textbf{H}^{a_{i,j}}\\}_{j\neq
i}]\in R^{|A_{i}|\times(m-1)l}\\\
\bar{\textbf{H}}^{a_{i}}&=\hat{\textbf{H}}^{a_{i}}W_{6}\in R^{|A_{i}|\times
l}\\\ g&=\sigma(\bar{\textbf{H}}^{a_{i}}W_{7}+\textbf{H}^{a_{i}}W_{8}+b)\\\
\textbf{H}^{o_{i}}&=g*\textbf{H}^{a_{i}}+(1-g)*\bar{\textbf{H}}^{a_{i}}\end{split}$
(6)

where $W_{7}$, $W_{8}\in R^{l\times l}$ and $W_{6}\in R^{(m-1)l\times l}$ are
learnable parameters, $\hat{\textbf{H}}^{a_{i}}$ is the concatenation of all
the interaction representations. $g\in R^{|A_{i}|\times l}$ is a reset gate
which balances the influence of $\bar{\textbf{H}}^{a_{i}}$ and
$\textbf{H}^{a_{i}}$, and $\textbf{H}^{o_{i}}$ is the final option
representation of $\textbf{A}_{i}$ encoded with the interaction information.
At last, we denote
$\textbf{O}=\\{\textbf{H}^{o_{1}},\textbf{H}^{o_{2}},...,\textbf{H}^{o_{m}}\\}$
as the final answer option representation set fused with comparison
information across answer options.

### Bidirectional Matching

The triplet changes from {P, Q, A} to {$\textbf{P}_{s}$, Q, O} with passage
sentence selection and answer option interaction. To fully model the
relationship in the {$\textbf{P}_{s}$, Q, O} triplet, bidirectional matching
is built to get all pairwise representations among the triplet, including
passage-answer, passage-question and question-answer representation. Here
shows how to model the relationship between question-answer sequence pair as
an example and it is the same for the other two pairs.

Bidirectional matching representation between the question $\textbf{H}^{q}$
and answer option $\textbf{H}^{o}$ can be calculated as follows:

$\begin{split}\textbf{G}^{qo}&=SoftMax(\textbf{H}^{q}W_{9}{\textbf{H}^{o}}^{T})\\\
\textbf{G}^{oq}&=SoftMax(\textbf{H}^{o}W_{10}{\textbf{H}^{q}}^{T})\\\
\textbf{E}^{q}&=\textbf{G}^{qo}\textbf{H}^{o},\textbf{E}^{o}={\textbf{G}^{oq}}\textbf{H}^{q}\\\
\textbf{S}^{q}&=ReLU(\textbf{E}^{q}W_{11})\\\
\textbf{S}^{o}&=ReLU(\textbf{E}^{o}W_{12})\end{split}$ (7)

where $W_{9}$, $W_{10}$, $W_{11}$, $W_{12}\in R^{l\times l}$ are learnable
parameters. $\textbf{G}^{qo}\in R^{|Q|\times|O|}$ and $\textbf{G}^{oq}\in
R^{|O|\times|Q|}$ are the weight matrices between question and answer option.
$\textbf{E}^{q}\in R^{|Q|\times l},\textbf{E}^{o}\in R^{|A|\times l}$
represent option-aware question representation and question-aware option
representation, respectively. The final representation of question-answer pair
is calculated as follows:

$\begin{split}\textbf{S}^{q\\_o}&=MaxPooling(\textbf{S}^{q})\\\
\textbf{S}^{o\\_q}&=MaxPooling(\textbf{S}^{o})\\\
g&=\sigma(\textbf{S}^{q\\_o}W_{13}+\textbf{S}^{o\\_q}W_{14}+b)\\\
\textbf{M}^{q\\_o}&=g*\textbf{S}^{o\\_q}+(1-g)*\textbf{S}^{o\\_q}\end{split}$
(8)

where $W_{13},W_{14}\in R^{l\times l}$ and $b\in R^{l}$ are three learnable
parameters. After a row-wise max pooling operation, we get the aggregation
representation $\textbf{M}^{q}\in R^{l}$ and $\textbf{M}^{o}\in R^{l}$. $g\in
R^{l}$ is a reset gate. $\textbf{M}^{q\\_o}\in R^{l}$ is the final
bidirectional matching representation of the question-answer sequence pair.

Passage-question and passage-option sequence matching representation
$\textbf{M}^{p\\_q},\textbf{M}^{p\\_o}\in R^{l}$ can be calculated in the same
procedure from Eq.(7) to Eq.(8). The framework of this module is shown in
Figure 1.

### Objective Function

With the built matching representations
$\textbf{M}^{p\\_q},\textbf{M}^{p\\_o},\textbf{M}^{q\\_o}$ for three sequence
pairs, we concatenate them as the final representation $\textbf{C}\in R^{3l}$
for each passage-question-option triplet. We denote the representation
$\textbf{C}_{i}$ for each {$P_{s},Q,O_{i}$} triplet. If $A_{k}$ is the correct
option, then the objective function can be computed as follows:

$\begin{split}\textbf{C}&=[\textbf{M}^{p\\_q};\textbf{M}^{p\\_o};\textbf{M}^{q\\_o}]\\\
L(A_{k}|P,Q)&=-log{\frac{\text{exp}(V^{T}\textbf{C}_{k})}{\sum_{j=1}^{m}{\text{exp}(V^{T}\textbf{C}_{j})}}}\end{split}$
(9)

where $V\in R^{3l}$ is a learnable parameter and $m$ is the number of answer
options.

## Experiments

### Dataset

We evaluate our model on five multi-choice MRC datasets from different
domains. Statistics of these datasets are detailed in Table 3. Accuracy is
calculated as $acc=N^{+}/N$, where $N^{+}$ and $N$ are the number of correct
predictions and the total number of questions. Some details about these
datasets are shown as follows:

Task |  Domain |  #o |  #p |  #q
---|---|---|---|---
RACE |  general |  4 |  27,933 |  97,687
SemEval |  narrative text |  2 |  2,119 |  13,939
ROCStories |  stories |  2 |  3472 |  3472
MCTest |  stories |  4 |  660 |  2,640
COIN |  everyday scenarios |  2 |  _ |  5,102
Table 3:  Statistics of multi-choice machine reading comprehension datasets.
#o is the average number of candidate options for each question. #p and #q are
the number of documents and questions in the dataset.

  * •

RACE (?): RACE consists of two subsets: RACE-M and RACE-H respectively
corresponding to middle school and high school difficulty levels, which is
recognized as one of the largest and most difficult datasets in multi-choice
reading comprehension.

  * •

SemEval-2018 Task11 (?): Multi-choice questions should be answered based on
narrative texts about everyday activities.

  * •

ROCStories (?): This dataset contains 98,162 five-sentence coherent stories in
the training dataset, 1,871 four-sentence story contexts along with a right
ending and a wrong ending in the development and test datasets, respectively.

  * •

MCTest (?): This task requires machines to answer questions about fictional
stories, directly tackling the high-level goal of open-domain machine
comprehension.

  * •

COIN Task 1 (?): The data for the task is short narrations about everyday
scenarios with multiple-choice questions.

### Implementation Details

Our model is evaluated based on the pre-trained language model BERT (?) and
XLNet (?) which both have small and large versions. The basic version
BERT${}_{\text{base}}$ has 12-layer transformer blocks, 768 hidden-size, and
12 self-attention heads, totally 110M parameters. The large version
BERT${}_{\text{large}}$ has 24-layer transformer blocks, 1024 hidden-size, and
16 self-attention heads, totally 340M parameters. Two versions of XLNet have
the similar sizes as BERT.

In our experiments, the max input sequence length is set to 512. A dropout
rate of 0.1 is applied to every BERT layer. We optimize the model using
BertAdam (?) optimizer with a learning rate 2e-5. We train for 10 epochs with
batch size 8 using eight 1080Ti GPUs when BERT${}_{\text{large}}$ and
XLNet${}_{\text{large}}$ are used as the encoder. Batch size is set to 16 when
using BERT${}_{\text{base}}$ and XLNet${}_{\text{base}}$ as the encoder111Our
code is at https://github.com/Qzsl123/dcmn..

Model | RACE-M/H | RACE
---|---|---
HAF (?) | 45.0/46.4 | 46.0
MRU (?) | 57.7/47.4 | 50.4
HCM (?) | 55.8/48.2 | 50.4
MMN (?) | 61.1/52.2 | 54.7
GPT (?) | 62.9/57.4 | 59.0
RSM (?) | 69.2/61.5 | 63.8
OCN (?) | 76.7/69.6 | 71.7
XLNet (?) | 85.5/80.2 | 81.8
BERT${}_{\text{base}}$∗ | 71.1/62.3 | 65.0
BERT${}_{\text{large}}$∗ | 76.6/70.1 | 72.0
XLNet${}_{\text{large}}$∗ | 83.7/78.6 | 80.1
Our Models |  |
BERT${}_{\text{base}}$∗ \+ DCMN | 73.2/64.2 | 67.0
BERT${}_{\text{large}}$∗ \+ DCMN | 79.2/72.1 | 74.1
BERT${}_{\text{large}}$∗ \+ DCMN + P${}_{\text{SS}}$ \+ A${}_{\text{OI}}$ | 79.3/74.4 | 75.8
XLNet${}_{\text{large}}$∗ \+ DCMN + P${}_{\text{SS}}$ \+ A${}_{\text{OI}}$ | 86.5/81.3 | 82.8
Human Performance |  |
Turkers | 85.1/69.4 | 73.3
Ceiling | 95.4/94.2 | 94.5

Table 4:  Experiment results on RACE test set. All the results are from single
models. P${}_{\text{SS}}$: Passage Sentence Selection; A${}_{\text{OI}}$:
Answer Option Interaction. ∗ indicates our implementation.

### Evaluation and Ablation Study on RACE

Table 4 reports the experimental results on RACE and its two subtasks: RACE-M
and RACE-H. In the table, Turkers is the performance of Amazon Turkers on a
randomly sampled subset of the RACE test set and Ceiling is the percentage of
the unambiguous questions with a correct answer in a subset of the test set.
Here we give the results of directly fine-tuned BERT${}_{\text{base}}$,
BERT${}_{\text{large}}$ and XLNet${}_{\text{large}}$ on RACE and get the
accuracy of 65.0%, 72.0% and 80.1%, respectively. Because of the limited
computing resources, the largest batch size can only be set to 8 in our
experiments which leads to 1.7% decrease (80.1% vs. 81.8%) on XLNet compared
to the result reported in (?)222The implementation is very close to the result
80.3% in (?) when using batch size 8 on RACE..

The comparison indicates that our proposed method obtains significant
improvement over pre-trained language models (75.8% vs. 72.0% on
BERT${}_{\text{large}}$ and 82.8% vs. 80.1% on XLNet${}_{\text{large}}$) and
achieves the state-of-the-art result on RACE.

In Table 5, we focus on the contribution of main components (DCMN, passage
sentence selection and answer option interaction) in our model. From the
results, the bidirectional matching strategy (DCMN) gives the main
contribution and achieves further improvement by integrating with the two
reading strategies. Finally, we have the best performance by combining all
components.

| BERT${}_{\text{base}}$ | BERT${}_{\text{large}}$ | XLNet${}_{\text{large}}$
---|---|---|---
base encoder | 64.6 | 71.8 | 80.1
\+ DCMN | 66.0 (+1.4) | 73.8 (+2.0) | 81.5 (+1.4)
\+ DCMN + P$\\_{SS}$ | 66.6 (+2.0) | 74.6 (+2.8) | 82.1 (+2.0)
\+ DCMN + P$\\_{OI}$ | 66.8 (+2.2) | 74.4 (+2.6) | 82.2 (+2.1)
\+ DCMN + ALL (DCMN+) | 67.4 (+2.8) | 75.4 (+3.6) | 82.6 (+2.5)

Table 5:  Ablation study on RACE dev set. P${}_{\text{SS}}$: Passage Sentence Selection. A${}_{\text{OI}}$: Answer Option Interaction. DCMN+: DCMN + P${}_{\text{SS}}$ \+ A${}_{\text{OI}}$. Task | Previous STOA |  | BERT | DCMN_BERT | XLNet | DCMN_XLNet
---|---|---|---|---|---|---
SemEval Task 11 | (?) | 89.5 | 90.5 | 91.8 (+1.3) | 92.0 | 93.4 (+1.4)
ROCStories | (?) | 91.8 | 90.8 | 92.4 (+1.6) | 93.8 | 95.8 (+2.0)
MCTest-MC160 | (?) | 81.7 | 73.8 | 85.0 (+11.2) | 80.6 | 86.2 (+5.6)
MCTest-MC500 | (?) | 82.0 | 80.4 | 86.5 (+6.1) | 83.4 | 86.6 (+3.2)
COIN Task 1 | (?) | 84.2 | 84.3 | 88.8 (+4.5) | 89.1 | 91.1 (+2.0)
Average |  | 85.8 | 84.0 | 88.9 (+4.9) | 87.8 | 90.6 (+2.8)
Table 6:  Results on the test set of SemEval Task 11, ROCStories, MCTest and
the development set of COIN Task 1. The test set of COIN is not public.
DCMN_BERT: BERT + DCMN + P${}_{\text{SS}}$ \+ A${}_{\text{OI}}$. Previous
SOTA: previous state-of-the-art model. All the results are from single models.

### Evaluation on Other Multi-choice Datasets

The results on four other multi-choice MRC challenges are shown in Table 6.
When adapting our method to the non-conventional MRC dataset ROCStories which
requires to choose the correct ending to a four-sentence incomplete story from
two answer options (?), the question context is left empty as no explicit
questions are provided. Passage sentence selection is not used in this dataset
because there are only four sentences as the passage. Since the test set of
COIN is not publicly available, we report the performance of the model on its
development set.

Model | RACE | Model | RACE | Model | RACE
---|---|---|---|---|---
BERT${}_{\text{base}}$ | 64.6 |  |  |  |
\+ Unidirectional |  |  |  |  |
$[S^{P\\_O};S^{P\\_Q};S^{O\\_Q}]$ | 65.0 | $[S^{P\\_Q};S^{Q\\_O}]$ | 63.4 | $[S^{P\\_Q};S^{O\\_Q}]$ | 64.5
$[S^{P\\_O};S^{Q\\_P};S^{O\\_Q}]$ | 65.2 | $[S^{P\\_O};S^{Q\\_O}]$ | 63.6 | $[S^{Q\\_P};S^{O\\_Q}]$ | 65.2
$[S^{P\\_Q};S^{P\\_O}]$ (HCM) | 64.4 | $[S^{P\\_O};S^{O\\_Q}]$ | 64.2 | $[S^{P\\_Q};S^{O\\_P}]$ | 64.7
$[S^{P\\_O};S^{P\\_Q};S^{Q\\_O}]$ (HAF) | 64.2 | $[S^{P\\_O};S^{Q\\_P};S^{Q\\_O}]$ | 64.4 | $[S^{Q\\_P};S^{Q\\_O}]$ | 64.3
$[S^{Q\\_O};S^{O\\_Q};S^{P\\_Q};S^{P\\_O}]$ (MMN) | 63.2 |  |  |  |
\+ Bidirectional |  |  |  |  |
$[M^{P\\_Q};M^{P\\_O}]$ | 66.4 | $[M^{P\\_O};M^{Q\\_O}]$ | 66.0 | $[M^{P\\_Q};M^{Q\\_O}]$ | 65.5
$[M^{P\\_Q};M^{P\\_O};M^{Q\\_O}]$ (DCMN) | 67.1 |  |  |  |

Table 7:  Performance comparison with different combination methods on the
RACE dev set. (HCM) (?), (HAF) (?), (MMN) (?) are previous methods. We use
BERT${}_{\text{base}}$ as our encoder here. $[;]$ indicates the concatenation
operation. $S^{P\\_O}$ and $M^{P\\_O}$ are the unidirectional and
bidirectional representation referred in Eq. 8.

As shown in Table 6, we achieve state-of-the-art (SOTA) results on all
datasets and obtain 3.1% absolute improvement in average accuracy over the
previous average SOTA (88.9% vs. 85.8%) by using BERT as encoder and 4.8%
(90.6% vs. 85.8%) by using XLNet as encoder. To further investigate the
contribution of our model, we also report the results of directly fine-tuned
BERT/XLNet on the target datasets. From the comparison, we can see that our
model obtains 4.9% and 2.8% absolute improvement in average accuracy over the
baseline of directly fine-tuned BERT (88.9% vs. 84.0%) and XLNet (90.6% vs.
87.8%), respectively. These results indicate our proposed model has a heavy
impact on the performance no matter how strong the adopted pre-trained
language model itself is.

### Comparison with Unidirectional Methods

Here we focus on whether the bidirectional matching method works better than
previous unidirectional methods. In Table 7, we enumerate all the combinations
of unidirectional matching strategies333Here we omit the combinations with
$S^{O\\_P}$ because we find the combinations with $S^{P\\_O}$ works better
than $S^{O\\_P}$. which only use passage-aware question representation
$S^{Q\\_P}$ or question-aware passage representation $S^{P\\_Q}$ when modeling
the relationship between the passage and question. Specially, we roughly
summarize the matching methods in previous work (i.e. HCM, HAF, MMN) using our
model notations which meet their general ideas except some calculation
details.

From the comparison, we observe that previous matching strategies (HCM 64.4%,
HAF 64.2%, MMN 63.2%) fail to give further performance improvement over the
strong encoder (64.6%). In contrast, all bidirectional combinations work
better than the encoder. All three pairwise matching representations
($M^{P\\_Q}$, $M^{P\\_O}$, $M^{Q\\_O}$) are necessary and by concatenating
them together, we achieve the highest performance (67.1%).

Top K | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
RACE-cos | 58.4 | 60.1 | 63.3 | 65.8 | 66.5 | 66
RACE-bi | 59.5 | 60.5 | 63.4 | 66.8 | 66.4 | 66.2
COIN-cos | 81.0 | 82.0 | 83.5 | 83.0 | 82.5 | 82.4
COIN-bi | 81.7 | 82.0 | 82.6 | 82.8 | 82.4 | 82.2
Table 8:  Results on RACE and COIN dev set with cosine and bilinear score in
P${}_{\text{SS}}$. We use BERT${}_{\text{base}}$ as encoder here.

### Results with Different Settings in P${}_{\text{SS}}$

Table 8 shows the performance comparison with different scoring methods, and
we observe that both methods have their advantages and disadvantages. Cosine
score method works better on COIN dataset (83.5% vs. 82.8%) and bilinear score
works better on RACE dataset (66.8% vs. 66.5%).

Figure 2 shows the results of passage sentence selection (Pss) on COIN and
RACE dev set with different numbers of selected sentences (Top $K$). The
results without Pss module are also shown in the figure (RACE-w and COIN-w).
We observe that Pss mechanism consistently shows a positive impact on both
datasets when more than four sentences are selected compared to the model
without Pss (RACE-w and COIN-w). The highest performance is achieved when top
3 sentences are selected on COIN and top 5 sentences on RACE where the main
reason is that the questions in RACE are designed by human experts and require
more complex reasoning.

Figure 2: Results of sentence selection on dev sets of RACE and COIN when
selecting different numbers of sentences (Top $K$). We use
BERT${}_{\text{base}}$ as encoder and cosine score method here. RACE/COIN-w
indicates the results on RACE/COIN without passage sentence selection module.

### Why Previous Methods Break Down?

As shown in Table 7, applying previous models to a strong BERT encoder fails
to give performance increase over directly fine-tuned BERT. The contrast is
clear that our proposed model achieves more than 3.8% absolute increase over
the BERT baseline. We summarize the reasons resulting in such contrast as
follows: (i) the unidirectional representations cannot well capture the
relationship between two sequences, (ii) previous methods use elementwise
subtraction and multiplication to fuse $\textbf{E}^{q}$ and $\textbf{H}^{o}$
in Eq. 7 (i.e.,
$[\textbf{E}^{q}\ominus\textbf{H}^{o};\textbf{E}^{q}\otimes\textbf{H}^{o}]$)
which is shown suboptimal as such processing breaks the symmetry of equation.
Symmetric representations from both directions show essentially helpful for
bidirectional architecture.

### Evaluation on Different Types of Questions

Inspired by (?), we further analyze the performance of the main components on
different question types. Questions are roughly divided into five categories:
_detail_ , _inference_ , _main_ , _attitude_ and _vocabulary_ (?; ?). We
annotate all the instances of the RACE development set. As shown in Figure 3,
all the combinations of components work better than the BERT baseline in most
question types. Bidirectional matching strategy (DCMN) consistently improves
the results across all categories. DCMN+P${}_{\text{SS}}$ works best on the
inference and attitude categories which indicates PSS module may effectively
improve the reasoning ability of the model. DCMN+A${}_{\text{OI}}$ works
better than DCMN on detail and main categories which indicates that the model
achieves better distinguish ability with answer option interaction
information.

Figure 3: Results on different question types, tested on the RACE dev set.
BERT${}_{\text{large}}$ is used as encoder here. OI: Answer Option
Interaction. SS: Passage Sentence Selection.

## Related Work

Neural network based methods have been applied to several natural language
processing tasks, especially to MRC (?; ?).

The task of selecting sentences to answer the question has been studied across
several question-answering (QA) datasets, by modeling the relevance between a
sentence and the question (?; ?; ?; ?; ?). (?) apply distant supervision to
generate imperfect labels and then use them to train a neural evidence
extractor. (?) propose a simple sentence selector to select the minimal set of
sentences then feed into the QA model. They are different from our work in
that (i) we select the sentences by modeling the relevance among sentence-
question-option triplet, not sentence-question pair. (ii) Our model uses the
output of language model as the sentence embedding and computes the relevance
score using these sentence vectors directly, without the need of manually
defined labels. (iii) We achieve a generally positive impact by selecting
sentences while previous sentence selection methods usually bring performance
decrease in most cases.

Most recent works attempting to integrate answer option interaction
information focus on building attention mechanism at word-level (?; ?; ?)
whose performance increase is very limited. Our answer option interaction
module is different from previous works in that: (i) we encode the comparison
information by modeling the bilinear representation among the options at
sentence-level which is similar to modeling passage-question sequence
relationship, other than attention mechanism. (ii) We use gated mechanism to
fuse the comparison information into the original answer option
representations.

## Conclusion

This paper proposes dual co-matching network integrated with two reading
strategies (passage sentence selection and answer option interaction) to
enhance multi-choice machine reading comprehension. In terms of strong pre-
trained language models such as BERT and XLNet as encoder, our proposed method
achieves state-of-the-art results on five representative multi-choice MRC
datasets including RACE. The experiment results consistently indicate the
general effectiveness and applicability of our model.

## References

  * [Chen et al. 2019] Chen, Z.; Cui, Y.; Ma, W.; and Wang, S.  2019\.  Convolutional Spatial Attention Model for Reading Comprehension with Multiple-Choice Questions.  In AAAI 2019.
  * [Choi et al. 2017] Choi, E.; Hewlett, D.; Uszkoreit, J.; Polosukhin, I.; Lacoste, A.; and Berant, J.  2017\.  Coarse-to-Fine Question Answering for Long Documents.  In ACL 2017, 209–220.
  * [Devlin et al. 2019] Devlin, J.; Chang, M.-W.; Lee, K.; and Toutanova, K.  2019\.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  In NAACL 2019, 4171–4186.
  * [Hermann et al. 2015] Hermann, K. M.; Kocisky, T.; Grefenstette, E.; Espeholt, L.; Kay, W.; Suleyman, M.; and Blunsom, P.  2015\.  Teaching Machines to Read and Comprehend.  In NIPS 2015.
  * [Lai et al. 2017] Lai, G.; Xie, Q.; Liu, H.; Yang, Y.; and Hovy, E.  2017\.  RACE: Large-scale ReAding Comprehension Dataset From Examinations.  In EMNLP 2017, 785–794.
  * [Li, Ding, and Liu 2019] Li, Z.; Ding, X.; and Liu, T.  2019\.  Story Ending Prediction by Transferable BERT.  In IJCAI-19, 1800–1806.
  * [Luo, Xiao, and Zhao 2019] Luo, Y.; Xiao, F.; and Zhao, H.  2019\.  Hierarchical Contextualized Representation for Named Entity Recognition.  In AAAI 2020.
  * [Min et al. 2018a] Min, S.; Zhong, V.; Socher, R.; and Xiong, C.  2018a.  Efficient and Robust Question Answering from Minimal Context over Documents.  In ACL 2018, 1725–1735.
  * [Min et al. 2018b] Min, S.; Zhong, V.; Socher, R.; and Xiong, C.  2018b.  Efficient and Robust Question Answering from Minimal Context over Documents.  In ACL 2018, 1725–1735.
  * [Mostafazadeh et al. 2016] Mostafazadeh, N.; Chambers, N.; He, X.; Parikh, D.; Batra, D.; Vanderwende, L.; Kohli, P.; and Allen, J.  2016\.  A Corpus and Cloze Evaluation for Deeper Understanding of Commonsense Stories.  In NAACL 2016, 839–849.
  * [Nguyen et al. 2016] Nguyen, T.; Rosenberg, M.; Song, X.; Gao, J.; Tiwary, S.; Majumder, R.; and Deng, L.  2016\.  MS MARCO: A Human Generated MAchine Reading COmprehension Dataset.  CoRR abs/1611.09268.
  * [Ostermann et al. 2018] Ostermann, S.; Roth, M.; Modi, A.; Thater, S.; and Pinkal, M.  2018\.  SemEval-2018 Task 11: Machine Comprehension Using Commonsense Knowledge.  In Proceedings of The 12th International Workshop on Semantic Evaluation.
  * [Pujari and Goldwasser 2019] Pujari, R., and Goldwasser, D.  2019\.  Using Natural Language Relations between Answer Choices for Machine Comprehension.  In NAACL-HLT 2019, 4010–4015.
  * [Qian and Schedl 2004] Qian, D., and Schedl, M.  2004\.  Evaluation of an In-Depth Vocabulary Knowledge Measure for Assessing Reading Performance.  Language Testing - LANG TEST 21:28–52.
  * [Radford 2018] Radford, A.  2018\.  Improving Language Understanding by Generative Pre-Training.  In OpenAI preprint.
  * [Raiman and Miller 2017] Raiman, J., and Miller, J.  2017\.  Globally Normalized Reader.  In EMNLP 2017, 1059–1069.
  * [Rajpurkar et al. 2016] Rajpurkar, P.; Zhang, J.; Lopyrev, K.; and Liang, P.  2016\.  SQuAD: 100,000+ Questions for Machine Comprehension of Text.  In EMNLP 2016, 2383–2392.
  * [Ran et al. 2019] Ran, Q.; Li, P.; Hu, W.; and Zhou, J.  2019\.  Option Comparison Network for Multiple-choice Reading Comprehension.  CoRR abs/1903.03033.
  * [Richardson, Burges, and Renshaw 2013] Richardson, M.; Burges, C. J.; and Renshaw, E.  2013\.  MCTest: A Challenge Dataset for the Open-Domain Machine Comprehension of Text.  In EMNLP 2013, 193–203.
  * [Srivastava, Greff, and Schmidhuber 2015] Srivastava, R. K.; Greff, K.; and Schmidhuber, J.  2015\.  Highway networks.  In ICML 2015.
  * [Sun et al. 2019] Sun, K.; Yu, D.; Yu, D.; and Cardie, C.  2019\.  Improving Machine Reading Comprehension with General Reading Strategies.  In NAACL 2019.
  * [Tang, Cai, and Zhuo 2019] Tang, M.; Cai, J.; and Zhuo, H. H.  2019\.  Multi-Matching Network for Multiple Choice Reading Comprehension.  In AAAI 2019.
  * [Tay, Tuan, and Hui 2018] Tay, Y.; Tuan, L. A.; and Hui, S. C.  2018\.  Multi-range Reasoning for Machine Comprehension.  CoRR abs/1803.09074.
  * [Trischler et al. 2017] Trischler, A.; Wang, T.; Yuan, X.; Harris, J.; Sordoni, A.; Bachman, P.; and Suleman, K.  2017\.  NewsQA: A machine comprehension dataset.  In Proceedings of the 2nd Workshop on Representation Learning for NLP, 191–200.
  * [Wang et al. 2018a] Wang, L.; Sun, M.; Zhao, W.; Shen, K.; and Liu, J.  2018a.  Yuanfudao at SemEval-2018 Task 11: Three-way Attention and Relational Knowledge for Commonsense Machine Comprehension.  In Proceedings of The 12th International Workshop on Semantic Evaluation, 758–762.
  * [Wang et al. 2018b] Wang, S.; Yu, M.; Jiang, J.; and Chang, S.  2018b.  A Co-Matching Model for Multi-choice Reading Comprehension.  In ACL 2018, 746–751.
  * [Wang et al. 2019] Wang, H.; Yu, D.; Sun, K.; Chen, J.; Yu, D.; Roth, D.; and McAllester, D. A.  2019\.  Evidence Sentence Extraction for Machine Reading Comprehension.  CoRR abs/1902.08852.
  * [Xiao et al. 2019] Xiao, F.; Li, J.; Zhao, H.; Wang, R.; and Chen, K.  2019\.  Lattice-Based Transformer Encoder for Neural Machine Translation.  In ACL 2019, 3090–3097.
  * [Yang et al. 2019] Yang, Z.; Dai, Z.; Yang, Y.; Carbonell, J. G.; Salakhutdinov, R.; and Le, Q. V.  2019\.  XLNet: Generalized Autoregressive Pretraining for Language Understanding.  In NIPS 2019.
  * [Zhang et al. 2018] Zhang, Z.; Huang, Y.; Zhu, P.; and Zhao, H.  2018\.  Effective Character-augmented Word Embedding for Machine Reading Comprehension.  In NLPCC 2018.
  * [Zhang et al. 2019a] Zhang, Z.; Wu, Y.; Zhao, H.; Li, Z.; Zhang, S.; Zhou, X.; and Zhou, X.  2019a.  Semantics-aware BERT for Language Understanding.  In AAAI 2020.
  * [Zhang et al. 2019b] Zhang, Z.; Wu, Y.; Zhou, J.; Duan, S.; Zhao, H.; and Wang, R.  2019b.  SG-Net: Syntax-Guided Machine Reading Comprehension.  In AAAI 2020.
  * [Zhang et al. 2019c] Zhang, Z.; Zhao, H.; Ling, K.; Li, J.; Li, Z.; and He, S.  2019c.  Effective Subword Segmentation for Text Comprehension.  In IEEE/ACM Transactions on Audio, Speech, and Language Processing, 2019.
  * [Zhang, Huang, and Zhao 2018] Zhang, Z.; Huang, Y.; and Zhao, H.  2018\.  Subword-augmented Embedding for Cloze Reading Comprehension.  In COLING 2018, 1802–1814.
  * [Zhou and Zhao ] Zhou, J., and Zhao, H.  Head-Driven Phrase Structure Grammar Parsing on Penn Treebank.  In ACL 2019, 2396––2408.
  * [Zhu et al. 2018] Zhu, H.; Wei, F.; Qin, B.; and Liu, T.  2018\.  Hierarchical Attention Flow for Multiple-choice Reading Comprehension.  In AAAI 2018.


# Effects of turbulent diffusion and back-reaction on the dust distribution
around two resonant planets

Francesco Marzari1 and Gennaro D’Angelo,2
1Department of Physics and Astronomy, University of Padova, via Marzolo 8,
I-35131, Padova, Italy
2Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545,
USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted . Received ; in original form )

###### Abstract

In evolved and dusty circumstellar discs, two planets with masses comparable
to Jupiter and Saturn that migrate outwards while maintaining an orbital
resonance can produce distinctive features in the dust distribution. Dust
accumulates at the outer edge of the common gas gap, which behaves as a dust
trap, where the local dust concentration is significantly enhanced by the
planets’ outward motion. Concurrently, an expanding cavity forms in the dust
distribution inside the planets’ orbits, because dust does not filter through
the common gaseous gap and grain depletion in the region continues via inward
drifting. There is no cavity in the gas distribution because gas can filter
through the gap, although ongoing gas accretion on the planets can reduce the
gas density in the inner disc. Such behaviour was demonstrated by means of
simulations neglecting the effects of dust diffusion due to turbulence and of
dust backreaction on the gas. Both effects may alter the formation of the dust
peak at the gap outer edge and of the inner dust cavity, by letting grains
filter through the dust trap. We performed high resolution hydrodynamical
simulations of the coupled evolution of gas and dust species, the latter
treated as pressureless fluids, in the presence of two giant planets. We show
that diffusion and backreaction can change some morphological aspects of the
dust distribution but do not alter some main features, such as the outer peak
and the expanding inner cavity. These findings are confirmed for different
parametrizations of gas viscosity.

###### keywords:

accretion, accretion discs — methods: numerical — planets and satellites:
gaseous planets — planet–disc interactions

††pubyear: 2024††pagerange: Effects of turbulent diffusion and back-reaction
on the dust distribution around two resonant planets–References

##  1 Introduction

In a recent study, Marzari et al. (2019) examined the distribution of dust
particles around two resonant planets embedded in a circumstellar disc. The
two planets, with masses comparable to those of Jupiter and Saturn, had orbits
in a resonant configuration, with ratios of the mean motion equal to either
2:1 or 3:2. Because of the type of resonance and of the applied disc
conditions, the planets tend to migrate outwards and dust particles tend to
accumulate outside of the orbit of the exterior planet. Concurrently, the
inward migration of dust grains that move inside of the orbit of the interior
planet leads to an enlargement of the dust gap compared to the gap in the gas
and to a dynamical decoupling between the gaps in the gas and dust
distributions. The build-up of the dust density at the outer edge of the gap
surrounding the planets is markedly higher in the case of the 2:1 mean-motion
resonance and may appear as a bright ring (at appropriate wavelengths) in
resolved observations of discs. A similar phenomenon was also found for lower-
mass planets (in the Super-Earth mass range, Marzari & D’Angelo, 2020),
although less pronounced. All those simulations were performed without
including the effects of possible dust diffusion due to gas turbulence and of
dust back-reaction on the gas. Here we explore the consequences of these two
mechanism on the accumulation of dust grains at the outer edge of the gap and
on the formation of a wider gap in the dust distribution compared to the
density gap in the gas. This latter feature may lead to the formation of a
transition disc, if the planets are close enough to the star (low gas surface
density inside the planets’ orbits can be caused by ongoing accretion of gas
on star and planets).

Circumstellar discs are likely turbulent (e.g., Hughes et al., 2011). Various
mechanisms have been proposed as drivers of turbulence, such as convective
instability (e.g., Klahr & Hubbard, 2014; Lyra, 2014), vertical shear
instability (e.g., Urpin, 2003; Nelson et al., 2013; Stoll et al., 2017), and
magneto-rotational instability (e.g., Balbus & Hawley, 1991, 1998, 2003). Gas
turbulence may force dust grains to diffuse (over a length-scale dictated by
the type of turbulence), a process that not only may affect dust accumulation
but can also hinder the efficiency of dust entrapment at radial location of
gas pressure maxima. In fact, according to Sierra et al. (2019) and Pinilla et
al. (2020), dust diffusion may reduce, or even prevent, significant
concentrations of grains at locations of gas “bumps”, since the grains can
disperse out of the region. This process might affect the conclusions of our
previous results on the accumulation of dust at the outer edge of the gaseous
gap of two planets in resonance by letting dust seep through the gap and reach
the inner disc regions. If the effect is large enough, the concentration of
dust at those radial locations, obtained in the numerical simulations by
Marzari et al. (2019) and Marzari & D’Angelo (2020), may be severely depleted
or even largely absent.

In addition to dust diffusion, the back-reaction of dust on gas can also
impact the formation of grain concentration at a local pressure maximum.
According to Taki et al. (2016), dust back-reaction can deform the pressure
gradient of the gas when high-enough values of the dust-to-gas mass ratio are
reached. This may be the case of the dust concentration attained at the outer
edge of gaseous gaps, observed in the simulations of two planets in resonance
migrating away from the star.

To test the relevance of these two mechanisms, diffusion and back-reaction, on
the formation of dust over-dense regions caused by the outward migration of
two planets in resonance, we performed simulations of the evolution of two
planets in resonance in which both these two mechanisms are taken into
account. In Section 2, we describe the numerical model exploited to study the
coupled evolution of dust and gas in presence of the two resonant planets. In
Section 3, we outline the dust behaviour when the planets evolve in the 3:2
mean-motion resonance whereas, in Section 4, we analyse the case of the 2:1
mean-motion resonance. In Section 5, we test the robustness of our results by
changing the viscosity parameterization, including a constant kinematic
viscosity the one that applies a constant value of the $\alpha$ parameter.
Finally, in Section 6, we discuss our results.

##  2 Methods and Algorithms

In past work on the coupled evolution of dust and gas in protoplanetary discs,
we adopted a Lagrangian description of the solid component. Instead, an
Eulerian formalism is applied in the present study, since solids are treated
as pressureless fluids. Some details on the involved equations are provided
below to highlight the differences between the two approaches. Marzari et al.
(2019) and Marzari & D’Angelo (2020) used the two-dimensional (2D) FARGO
hydrodynamics code (Masset, 2000), modified to include the dynamical evolution
of dust particles embedded in the gaseous disc in a Lagrangian fashion
(Picogna & Kley, 2015; Picogna et al., 2018).

The drag force on the particles was computed from the local gas density
according to the equation (Woitke & Helling, 2003)

$\mathbf{F}=\left(\frac{3K}{3K+1}\right)^{2}\mathbf{F}_{\mathrm{E}}+\left(\frac{1}{3K+1}\right)^{2}\mathbf{F}_{\mathrm{S}},$
(1)

where $\mathbf{F}_{\mathrm{E}}$ is the Epstein drag contribution, given by

$\mathbf{F}_{\mathrm{E}}=-\frac{4}{3}\left(1+\frac{9\pi}{128}M^{2}\right)^{1/2}s^{2}\rho_{g}v_{\mathrm{th}}\mathbf{v}_{\mathrm{rel}},$
(2)

and $\mathbf{F}_{\mathrm{S}}$ is the Stokes drag component, given by

$\mathbf{F}_{\mathrm{S}}=-\frac{1}{2}C_{D}\pi
s^{2}\rho_{g}v_{\mathrm{rel}}\mathbf{v}_{\mathrm{rel}}.$ (3)

In the above equations, $\rho_{g}$ is the gas density, $s$ is the radius of
the particle, $v_{\mathrm{th}}$ is the local thermal velocity of the gas and
$\mathbf{v}_{\mathrm{rel}}$ is the velocity of the dust particle relative to
the gas. The quantity $K$ is the Knudsen number and $M$ is the Mach number
(computed from the particle’s relative velocity $v_{\mathrm{rel}}$). Quantity
$C_{D}$ is the drag coefficient for the Stokes regime.

In this paper, to test the relevance of diffusion and back-reaction on the
formation of dust over-dense regions caused by the outward migration of two
planets in resonance, we carry out simulations with the code FARGO3D (Benítez-
Llambay et al., 2019). In this version of the code, the dust particles are
treated as additional pressureless fluids where momentum is transferred
between the gas and each of the dust species (but not among dust species). The
dust fluid is affected by Epstein drag, which imparts a force _per unit
volume_ to a dust species given by

$\mathbf{F}_{d}=-\rho_{d}\frac{\Omega}{\tau_{s}}(\mathbf{v}_{d}-\mathbf{v}_{g}),$
(4)

where $\rho_{d}$ is the dust density, $\tau_{s}$ is the Stokes number of the
particle and $\Omega$ is the Keplerian frequency of the gas. An equal and
opposite force is imparted to the gas. Note that, in Equation (4), information
regarding the drag coefficient are incorporated into $\tau_{s}$.

A term is added to the continuity equation to model the diffusion of the dust
species within the gas (Morfill & Voelk, 1984)

$\frac{\partial\rho_{d}}{\partial
t}=\nabla\cdot\left(D_{d}\rho\nabla\frac{\rho_{d}}{\rho}\right),$ (5)

where $\rho=\rho_{g}+\rho_{d}$ and $\rho_{d}$ is the density of individual
dust species. Equation (5) is only applied to the pressureless fluids and it
assumes the same diffusion coefficient for all dust species, which is taken
equal to the value of the gas kinematic viscosity (Brauer et al., 2008)

$D_{d}=\nu.$ (6)

The effects of this choice are not tested.

The original version of the code applies Stokes numbers, $\tau_{s}$, in
Equation (4) that are constants. The code was modified so that we can select
each pressureless fluid (dust species) not according to a Stokes number, which
varies as a function of the local properties of the gas (density, temperature
and velocity), but rather according to the particle size.

##  3 Dust distribution near planets in the 3:2 resonance

We first investigate the case of a pair of planets that become captured in the
3:2 mean-motion resonance and migrate outward thereafter. The interior planet
has the mass of Jupiter whereas the exterior planet has the mass of Saturn. A
more massive inner planet is a condition conducive to outward migration. The
planets orbit in a cold, local-isothermal disc with a fixed aspect ratio
$H/r=h=0.02$. (Calculations with a larger ratio, $h=0.05$, are also
presented.)

The disc extends from $0.4$ to $12\,\mathrm{au}$ and is discretised over a
grid of $512\times 1024$ area elements (in the radius and azimuth,
respectively). The initial surface density of the gas declines as

$\Sigma(r)=\Sigma_{0}\left(\frac{r_{0}}{r}\right),$ (7)

where $\Sigma_{0}=200\,\mathrm{g\,cm^{-2}}$ is the density at the reference
radius $r_{0}=1\,\mathrm{au}$.

Three different populations of icy grains (bulk density
$1\,\mathrm{g\,cm^{-3}}$) are included in the simulations, whose sizes are
$100\,\mu\mathrm{m}$, $1\,\mathrm{mm}$ and $1\,\mathrm{cm}$. For the applied
disc conditions, these particles have Stokes numbers less than $\approx 0.1$.

The initial dust-to-gas mass ratio for each of the three dust species is
$0.0033$, so that the overall dust-to-gas mass ratio adds up to $0.01$, which
is a typical value adopted for circumstellar discs and is based on values
found in the interstellar medium. However, dust needs not be primordial in
origin, that is, part of the inventory of solids from which the planets
formed. In fact, the dust populations may represent, or contain, second
generation dust produced by collisions among leftover planetesimals, after the
planets became massive enough (Turrini et al., 2019; D’Angelo & Marzari, 2022;
Marta Bernabò et al., 2022). The equations describing the evolution of the
system are solved in a non-inertial reference frame centered on the star,
including the indirect terms arising from the planets’ and disc’s gravity.

Various values of the gas kinematic viscosity, $\nu$, are considered because
this parameter affects the tidal interactions between the planets and the gas,
and also determines the amount of dust diffusion through Equation (6).
Additionally, it can also affect the speed of the planets’ outward migration.
In these models, we adopt a constant value of kinematic viscosity. The impact
of $\nu$ on the efficiency of the outward migration is illustrated in Figure
1, for a given value of the initial gas density at the reference radius
$r_{0}$. In these models, the planets start to migrate at the beginning of the
simulations, when the distributions of gas and dust are unperturbed (hence the
initial rapid inward migration of the planets). As the outer planet approaches
the inner planet and their orbits become caught in resonance, the tidal
perturbations exerted by the outer planet on the circumstellar material alter
the torque balance on the inner planet. Consequently, the inner planet first
slows down and then migrates away from the star, pushing the outer planet
outward through the resonance forcing.

Figure 1: The top panel shows the semi-major axis of pairs of planets during
their migration, for three values of a constant kinematic viscosity, $\nu$, of
the gas (in units of $r^{2}_{0}\Omega_{0}$, see text). In these cases, the
planet pair is locked in the 3:2 mean-motion resonance. In the bottom panel
illustrates the evolution of the orbital eccentricity. The inner, more massive
planet has a lower eccentricity in all cases. Data are averaged over a $250$
yr window; see text for further details.

The evolution of the semi-major axis of the outer planet, $a_{2}$, is shown
for three different values of the gas kinematic viscosity: $\nu=10^{-6}$,
$5\times 10^{-6}$, and $10^{-5},$ in units of $r^{2}_{0}\Omega_{0}$
($\Omega_{0}$ is the Keplerian frequency at $r_{0}$). Note that a constant
kinematic viscosity corresponds to a variable $\alpha$ parameter,
$\alpha\propto\nu/(h^{2}\sqrt{r})$. With our choice of parameters, at
$r=r_{0}$, said values of $\nu$ would correspond to $\alpha_{0}=0.0025$,
$0.0125$, and $0.025$, respectively. However, in the disc regions where the
planets orbit, $\alpha$ would be smaller by a factor of $2$ or more.

After a different behaviour at the beginning of the evolution, prior to or
shortly after the capture into the 3:2 orbital resonance (see Figure 1), the
planets undergo sustained outward migration locked in mean-motion resonance.
The migration speed of the pair is related to $\nu$ and determined by the
shape of the common (or overlapping) gaseous gap of the two planets. A similar
outcome is obtained for the cases involving the 2:1 resonance, as shown in the
next section. Notice that the outer planet is subject to a negative torque
exerted by the disc material exterior to its orbit, hence it would tend to
move inward, whereas the resonance forcing pushes it outward. These two
opposing torques allow the resonance to be maintained during the outward
migration of the pair.

In the calculations presented herein, all disc material exerts torques on the
planets, including material within the planets’ Hill spheres. Since the
numerical resolution is limited and density variations close to the planets
may not be properly resolved, some spurious effects may arise that affect
outward migration. To quantify possible differences, one model was also
performed by removing the torques exerted from within the planets’ Hill
spheres. The orbital evolution is comparable to that of the calculation with
default setup, although at a somewhat reduced outward migration speed. The
amount of kinematic viscosity can also alter the local distribution of
material around the planets and unresolved density gradients can possibly
impact the resulting migration velocity. Nonetheless, it must be pointed out
that for the purposes of this study the details of the outward migration
process are not important, henceforth tests on the response of the system to
numerical parameters are unnecessary. The only requirement is that the planet
pair becomes locked in resonance and moves away from the star for a prolonged
amount of time.

Figure 2: Surface density (averaged in azimuth around the star) of dust
particles of different sizes, ranging from 100$~{}\mu$m to $1$ cm. These
profiles are compared to the re-scaled gas density (i.e., multiplied by
$0.0033$; the total dust-to-gas mass ratio is $0.01$). The top panel refers to
a kinematic viscosity equal to $\nu=10^{-6}$, the middle panel to $\nu=5\times
10^{-6}$ and the bottom panel to $\nu=10^{-5}$, in units of
$r^{2}_{0}\Omega_{0}$.

When the outer planet reaches $6.5$ au from the star, we compare the dust
distributions in the three cases with different viscosity. This comparison is
shown in Figure 2. The dust density profiles show a peak at the outer border
of the gas gap carved by the planets’ tidal perturbations. This peak is more
marked in the distributions of the largest grains, which are less coupled to
the gas (i.e., they have a larger Stokes number and therefore a longer
coupling timescale). Dust-to-gas mass ratios at the peak location range from
$0.04$ to $0.08$, increasing as $\nu$ decreases. In the region inside the
inner edge of the gap there is a significant depletion of dust due to drifting
motion towards the star. Re-supply of dust to this region is reduced, or
halted, by the dust trap at the outer edge of the gas gap. After some time,
the disc would appear as a transition disc with an inner cavity in the dust
density, which expands outward over time due to the outward migration of the
planets. These effects appear more evident at lower gas viscosity, which may
be due to the lower level of diffusion but it may also be related to the
different migration velocity of the planet pair (see Figure 1), as discussed
below. Beyond $\approx 10$ au the disc is depleted of mm- and cm-grains (but
not of the smallest grains), an effect associated to the inward drift of the
particles via gas drag (which does not affect as much the smallest grains).
This is a boundary effect related to the fact that particles are not flowing
inward from greater distances (i.e., there is not re-supply of solids at the
outer boundary). Test simulations not reported here confirm this issue.

The accumulation efficiency of particles at the outer edge of the gas gap
depends on the stopping time of the particles ($\tau_{s}/\Omega$) and the
timescale over which the radial pressure gradient of the gas moves (in our
case, due to planet migration). For a given stopping time, the shorter the
outward migration timescale, the less time dust grains have to accumulate.
Stated differently, for a similar orbital configuration of the planets, i.e.,
similar orbital frequencies $\Omega$, a more rapid outward migration can
facilitate the filtering process of dust toward the star. The overall outcome
is a less depleted inner disc in the cases with more vigorous outward
migration.

For the same reason, the different migration speed also controls the build up
of the dust at the outer edge of the gas gap. As a consequence, the reduction
in the peak density at the outer edge of the gap observed at higher viscosity
values, shown in Figure 2, is likely due to the combination of the two
effects: a higher diffusion rate and a decrease in the trapping efficiency of
dust grains caused by the faster outward migration rate of the planets.

To better characterize the role of diffusion and back-reaction in the
formation, size and shape of the outer peak in dust density, we ran a
simulation in which both these two processes were neglected, adopting a
kinematic viscosity $\nu=10^{-5}$ $r^{2}_{0}\Omega_{0}$ (see Figure 3). In
this model, the migration speed is broadly consistent with that of the model
illustrated in the bottom panel of Figure 2, but the peak in the dust density
distributions of the largest particles tend to be higher and sharper compared
to corresponding one in the the bottom panel of Figure 2. The implication is
that diffusion is indeed affecting the concentration of the dust by spreading
large grains over a broader radial region. Nonetheless, the reduced efficiency
in collecting dust caused by diffusion is not strong enough to prevent the
formation of a prominent density peak outside the planets’ orbits, which may
be detected by high resolution observations.

Figure 3: Dust and gas profiles (as in Figure 2) with $\nu=10^{-5}$
$r^{2}_{0}\Omega_{0}$, but without the inclusion of diffusion and back-
reaction.

In Figure 4 we plot the surface density distribution of the gas and of dust
grains of different sizes for the case with lowest viscosity, $\nu=10^{-6}$
$r^{2}_{0}\Omega_{0}$. The gas gap is significantly narrower than that of the
dust and the width of the latter is larger for larger grain size. For $1$ cm
size dust grains, the density distribution is confined in an over-dense ring
at the outer edge of the gas gap. Therefore, according to these results, dust
diffusion and back-reaction are not able to prevent the formation of narrow
rings in the dust distribution at the outer edge of the gas gap. In our
simulations there is no re-supply of dust at the outer boundary of the
simulated disc, but it is expected that a continuous distribution of solids
beyond the grid boundary would supply dust to the inner disc regions. In this
case, at the outer edge of the gap, we would observe an enhanced density, as
predicted by our simulations. Beyond this density peak, however, there would
be a continuous distribution of dust originating from more distant regions.
Over timescales much longer than those simulated here, dust drifting from
larger distances may also affect the density peaks at the outer edge of the
gas gap. In fact, an additional simulation with a wider radial boundary (not
reported here) does show enhanced peaks in large grains, due to solids
drifting from farther distances. Also the population of small grains would
increase at the peaks over longer times, but at a slower rate, dictated by the
drift velocity. Nonetheless, the depletion of dust within the inner edge of
the gas gap would not be affected by this process because the outer dust trap
appears to be efficient enough to halt (or severely impede) refilling of
grains. In fact, if refilling of the inner disc was sustained, it would occur
within the time of our simulations but it is not observed.

Continued supply of dust from farther out in the disc may, at some point,
raise the density in the peak regions beyond some threshold value to make the
peaks unstable. For example, the dust-to-gas mass ratio may become large
enough to induce a back-reaction response on the gas that redistributes the
particles over some radial region via collisions and/or enhanced dust
coagulation may ensue (coagulation into larger particles would reduce the
back-reaction force exerted on the gas, because of the lower surface-to-mass
ratio). The back-reaction response may also smooth out the gas density
gradient in the radial direction, altering the radial pressure gradient and
reducing its ability to retain grains. Such processes are not considered in
this study.

Figure 4: Surface density maps illustrating gas and dust distributions around
two planets locked in the 3:2 resonance ($\nu=10^{-6}$ $r^{2}_{0}\Omega_{0}$).
The top panel shows the gas density, whereas the second, third and fourth
panels represent $100$ $\mu$m, $1$ mm and $1$ cm particles, respectively.

Figure 5:  Dust and gas density profiles (as in Figure 2) in a disc with
$\nu=10^{-5}$ $r^{2}_{0}\Omega_{0}$ and aspect ratio $H/r=h=0.05$. In the top
panel $\Sigma_{0}=800\,\mathrm{g\,cm^{-2}}$ whereas
$\Sigma_{0}=200\,\mathrm{g\,cm^{-2}}$ in the bottom panel.

To test the impact of a larger aspect ratio (i.e., warmer disc) on the
formation of the inner dust cavity and of the peak external to the gas gap, we
performed two additional simulations adopting $h=0.05$. In the first model we
used a higher gas surface density ($\Sigma_{0}=800\,\mathrm{g\,cm^{-2}}$), in
order to increase the speed of planet, migration while in the second we used
the same density as in the previous cases
($\Sigma_{0}=200\,\mathrm{g\,cm^{-2}}$). In the latter simulation, since the
rate of outward migration is smaller, the planets are located closer to the
star in the plot. In the top panel of Figure 5, the density profiles are shown
for the different grain sizes in the high density case. The bottom panel
illustrates those profiles for the low density case. In both simulations the
density patterns are similar to those in Figure 2, suggesting that a higher
aspect ratio does not impair the ability of two planets in resonance to carve
an inner gap in the distribution of the larger grains and to trap dust at the
gap’s outer edge (dust-to-gas mass ratios around the peak region are
$0.02$–$0.03$). Thus, dust filtering through the outer edge of the gas gap is
not increased by the different morphology of the dust trap in these warmer
discs.

##  4 Dust distribution near planets in the 2:1 resonance

To test the evolution of the dust when the planets are captured in a 2:1 mean-
motion resonance, we decreased the gas density to
$\Sigma_{0}=40\,\mathrm{g/cm^{2}}$ in order to induce orbital locking in this
resonance. It is known that capture in this mean-motion resonance is a more
delicate process than capture in the 3:2 resonance. If the relative migration
velocity with which the pair of planets approach each other is above a certain
threshold, the resonance forcing is overcome and convergent migration
continues. Additionally, once the 2:1 orbital resonance is established,
migration of the two planets typically proceeds inward because of the way the
two gas gaps overlap. In this configuration, outward migration may be obtained
by choosing an appropriately low kinematic viscosity so that a wide, common
gaseous gap forms around the orbits of the two planets.

As in the the models of the previous section, the planets start to migrate in
unperturbed distributions of gas and therefore the initial inward migration of
the planets is artificially rapid. This choice affects the radius at which the
planets are trapped in resonance but is otherwise not much relevant for the
purposes of this study. As for the models discussed in Section 3, the model
with a larger value of the kinematic viscosity results in a more vigorous
outward migration of the planets once they become trapped in resonance. This
is illustrated in Figure 6 for cases with $\nu=10^{-6}$ and $\nu=10^{-7}$
$r^{2}_{0}\Omega_{0}$ ($\alpha_{0}=2.5\times 10^{-3}$ and $2.5\times 10^{-4}$,
respectively) . Note, however, that the orbits also become significantly
eccentric, which also affects outward migration (D’Angelo et al., 2006). We
also tested a larger value, $\nu=5\times 10^{-6}$ $r^{2}_{0}\Omega_{0}$, but
the outer planet crosses the 2:1 resonance with the inner planet. The pair is
temporarily trapped in the 5:3 resonance at which point it begins migrating
outwards. That resonance is then broken and the pair becomes finally captured
in the 3:2 resonance, continuing to migrate outward. For even higher values of
the kinematic viscosity, there is no trapping in the 2:1 resonance, which
might be attained by further reducing $\Sigma_{0}$. However, this possibility
was not tested since it could lead to initial dust-to-gas mass ratios too
dissimilar from the other simulations presented herein.

Figure 6:  The top panel shows the migration of the exterior planet of a pair,
for two values of a constant kinematic viscosity of the gas, $\nu$, as
indicated in the legend. The planet pair is locked in the 2:1 mean-motion
resonance. Only the semi-major axis of the exterior planet is plotted because
of the larger mutual distance compared, for example, to the 3:2 case. The
orbital eccentricity of both interior and exterior planets is illustrated in
the bottom panel. Data are averaged over a $250$ yr window as in the case of
the 3:2 resonance (Figure 1).

The dust distributions are shown in Figure 7 for the two different values of
kinematic viscosity and, in the bottom panel, for a model without the
inclusion of diffusion and back-reaction. Because of the significantly
different migration speed, we only compare the dust density distributions at
the end of the simulation (when the planet are located at different orbital
radii). For this resonance too a significant dust enhancement develops at the
outer edge of the gas gap, where dust-to-gas mass ratios can achieve values of
order unity. The inner cavity is larger at all grain sizes, compared to that
produced by the 3:2 resonant configuration. This is possibly related to the
slower outward migration rate of the planets in the 2:1 resonance, which can
reduce filtering of the dust through the orbits of the planets and toward the
star. Additionally, the lower gas density increases the Stokes numbers of the
grains, reducing their drift timescale and facilitating grain removal from the
inner disc.

Comparing the top and bottom panels of Figure 7, one can notice significant
differences in the dusty features exterior to the planets orbits. Since the
gas density is also different in the two models (at those times), it is
unclear how much of the difference is caused by the action of diffusion and
back-reaction of the solids.

The orbital eccentricity of the planets can drive asymmetries in the
distribution of the disc material. Consequently, both the gas gap outer edge
and dusty rings can become asymmetric (see Figure 8). We did not perform a
detailed analysis of ring asymmetries. However, one possible explanation for
the more asymmetric features arising from the planets in the 2:1 resonance
(compared to those in the 3:2 resonance, see Figure 4), may be the larger
eccentric perturbation driven by the inner (more massive) planet, which tends
to have a larger orbital eccentricity in the 2:1 resonance than in the 3:2
resonance (the outer planets can have comparable eccentricities).

In Figure 8, the dust density distributions are shown for grains of various
sizes for the case with lowest viscosity, $\nu=10^{-7}$ $r^{2}_{0}\Omega_{0}$.
As mentioned, the width of the gap is significantly larger compared to the
case of the 3:2 resonance and the dusty ring at the outer edge of the gas gap
is evident for all particles, more than in the 3:2 resonance. For the 2:1
resonance, the effects of the dust trap appear more marked, both at the inner
and at the outer edge of the gas gap. It is expected that after sufficient
time from the beginning of the outward migration, the disc reduces to a single
overdense ring at the outer edge of the gap, close to the outer planet orbit.
As discussed in the previous section, this feature would result from the lack
of dust supply from larger orbital radii, beyond the outer boundary of the
grid.

Figure 7: Density profile of dust particles of different sizes ranging from
$100$ $\mu$m to $1$ cm with the planets in a 2:1 resonance. These profiles are
compared to the re-scaled gas density (multiplied by $0.0033$). The top panel
is for a kinematic viscosity equal to $\nu=10^{-6}$ $r^{2}_{0}\Omega_{0}$ and
the middle panel for $\nu=10^{-7}$ $r^{2}_{0}\Omega_{0}$. The bottom panel
refers to a model with $\nu=10^{-6}$ $r^{2}_{0}\Omega_{0}$, but without
diffusion and back-reaction.

Figure 8: Density maps illustrating gas and dust distributions around planets
locked in the 2:1 orbital resonance ($\nu=10^{-7}$ $r^{2}_{0}\Omega_{0}$) .
The top panel refers to the surface density of the gas. The other panels
refer, respectively, to the distributions of $100$ $\mu$m, $1$ mm and $1$ cm
dust particles.

##  5 Comparison with Alpha Viscosity

To test the robustness of our results, we performed two additional
simulations, for the 2:1 resonance, in which the kinematic viscosity is set as
$\nu=\alpha H^{2}\Omega$, where $H$ is the local pressure scale height of the
disc and the parameter $\alpha$ is assumed to be a constant. Since we work
with a local-isothermal disc and a constant aspect ratio, $H\propto r$ and
$\nu\propto\alpha\sqrt{r}$. In the first model, we set $\alpha=10^{-3}$
whereas we set $\alpha=10^{-4}$ in the second. These parameters result in
$\nu\approx 10^{-6}$ and $\approx 10^{-7}$ $r^{2}_{0}\Omega_{0}$ around the
middle radius of the computational domain.

The outcomes of these simulations are illustrated in Figure 9, after $4\times
10^{4}$ yrs of evolution. The density peak at the outer edge of the gas gap is
clearly visible in both cases, although there are differences in over-density
morphology. In the top panel, at $\alpha=10^{-3}$, the peak appears similar at
all sizes while, in the bottom panel ($\alpha=10^{-4}$), the peak is split in
two for the largest particles ($s=1$ mm and $s=1$ cm) and composed of three
separate maxima for $s=100\,\mu$m. This behaviour was not observed for the 3:2
resonance but it was already present in the simulation involving the 2:1 mean-
motion resonance obtained with a constant kinematic viscosity (see Figure 7,
top panel). A possible interpretation is that for the 2:1 resonance multiple
dust traps develop at the outer border of the gap. By comparing the gas
density distribution in the two cases shown in Figure 9, this hypothesis
appears to be confirmed since, in the case with $\alpha=10^{-4}$, $\Sigma$ at
the outer border of the gas gap appears more variable compared to that of the
case with $\alpha=10^{-3}$.

The splitting of the peak in various maxima can be also observed in the
simulations of Marzari et al. (2019), for the case involving a pair of planet
locked in the 2:1 resonance, even if their physical model (radiative disc),
the initial parameters for the gas density and the viscosity values are
different.

Figure 9: Dust and gas density profiles (as in Figure 7) around a pair of
planets locked in a 2:1 orbital resonance, after $4\times 10^{4}$ yrs, in an
$\alpha$-viscosity disc. In the top panel $\alpha=10^{-3}$ and in the bottom
panel $\alpha=10^{-4}$.

##  6 Conclusions

According to the models presented by Marzari et al. (2019), two planets in the
mass range of Jupiter and Saturn, embedded in a circumstellar disc and
migrating outwards, can strongly affect the dust distribution in their
surroundings. They can create a transition disc with an inner cavity in the
dust distribution that expands following the outward migration of the planets.
They can also build a strong peak in the dust density at the outer edge of the
gas gap carved by the resonant planets, which would appear as a bright dust
ring. However, these findings may be altered by the dust diffusion due to gas
turbulence and by a change in gap morphology due to the back-reaction of the
dust on the gas. It is, therefore, important to test if these two processes
weaken the dust trap at the outer edge of the gap, allowing the dust to filter
through the gap and preventing the formation of an inner dust cavity and of
the outer over-density ring.

We used the code FARGO3D (Benítez-Llambay et al., 2019), in which dust species
are treated as additional pressureless fluids, and performed a set of local-
isothermal, high resolution simulations involving giant planets locked in the
3:2 and 2:1 orbital resonances and including the effects of diffusion and dust
back-reaction. For the 3:2 resonance, the diffusion and back-reaction slightly
affect the dust distribution by reducing the height of the overdense region at
the outer edge of the gas gap but without removing it. This has been tested
for different values of the disc viscosity, which determines the amount of
diffusion. This outcome proves that the dust trap is still strong enough to
halt the inward drift of the dust and to lead to the formation of an inner
cavity in the dust distribution and an overdense ring at the outer edge of the
gap. The details of peak heights can also be affected by continued supply of
solids from large radial distances and, therefore, depend on boundary effects
and evolution timescales. Given the limited radial extent of the models
presented herein, the possible feedback of very dense rings of solids on the
gas distribution was not investigated.

The 2:1 resonance also provides a robust mechanism capable of creating an
efficient dust trap. The height of the outer peak (outside the planets’
orbits) is not much affected by diffusion. Its morphology, however, appears
more complex than it does in the 3:2 resonance situation, since in some cases
two dust peaks form at the outer edge of the gas gap, possibly due to the
development of multiple dust traps. The robustness of these results for this
resonance were tested by performing two additional simulations in which we
switched from a constant kinematic viscosity, $\nu$, to a constant $\alpha$
viscosity parameter. In terms of large-scale features, the outcome of these
simulations do not significantly differ from those with constant $\nu$,
showing that the formation of the inner dust cavity and the outer peak(s) are
not due to the viscosity parametrization (although details can depend on the
type of viscosity).

In the models presented herein, the gas distribution interior to the planets’
orbits is not significantly depleted. However, ongoing accretion on the
planets, neglected here, is expected to reduce the gas mass flux toward the
inner disc (Lubow & D’Angelo, 2006), possibly leading to the formation of an
inner cavity in the gas distribution as well.

## Acknowledgements

We thank the reviewer, Clément Baruteau, whose comments helped us improve this
paper. GD acknowledges support provided by NASA’s Research Opportunities in
Space and Earth Science. Computational resources supporting this work were
provided by the NASA High-End Computing (HEC) Program through the NASA
Advanced Supercomputing (NAS) Division at Ames Research Center.

## Data Availability

The data underlying the research results described in the article will be
shared upon reasonable request to the authors.

## References

  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Reviews of Modern Physics, 70, 1
  * Balbus & Hawley (2003) Balbus S. A., Hawley J. F., 2003, in Falgarone E., Passot T., eds, , Vol. 614, Turbulence and Magnetic Fields in Astrophysics. pp 329–348, doi:10.48550/arXiv.astro-ph/0203353
  * Benítez-Llambay et al. (2019) Benítez-Llambay P., Krapp L., Pessah M. E., 2019, ApJS, 241, 25
  * Brauer et al. (2008) Brauer F., Dullemond C. P., Henning T., 2008, A&A, 480, 859
  * D’Angelo & Marzari (2022) D’Angelo G., Marzari F., 2022, MNRAS, 509, 3181
  * D’Angelo et al. (2006) D’Angelo G., Lubow S. H., Bate M. R., 2006, ApJ, 652, 1698
  * Hughes et al. (2011) Hughes A. M., Wilner D. J., Andrews S. M., Qi C., Hogerheijde M. R., 2011, ApJ, 727, 85
  * Klahr & Hubbard (2014) Klahr H., Hubbard A., 2014, ApJ, 788, 21
  * Lubow & D’Angelo (2006) Lubow S. H., D’Angelo G., 2006, ApJ, 641, 526
  * Lyra (2014) Lyra W., 2014, ApJ, 789, 77
  * Marta Bernabò et al. (2022) Marta Bernabò L., Turrini D., Testi L., Marzari F., Polychroni D., 2022, arXiv e-prints, p. arXiv:2202.13891
  * Marzari & D’Angelo (2020) Marzari F., D’Angelo G., 2020, A&A, 641, A125
  * Marzari et al. (2019) Marzari F., D’Angelo G., Picogna G., 2019, AJ, 157, 45
  * Masset (2000) Masset F., 2000, A&AS, 141, 165
  * Morfill & Voelk (1984) Morfill G. E., Voelk H. J., 1984, ApJ, 287, 371
  * Nelson et al. (2013) Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610
  * Picogna & Kley (2015) Picogna G., Kley W., 2015, A&A, 584, A110
  * Picogna et al. (2018) Picogna G., Stoll M. H. R., Kley W., 2018, A&A, 616, A116
  * Pinilla et al. (2020) Pinilla P., Pascucci I., Marino S., 2020, A&A, 635, A105
  * Sierra et al. (2019) Sierra A., Lizano S., Macías E., Carrasco-González C., Osorio M., Flock M., 2019, ApJ, 876, 7
  * Stoll et al. (2017) Stoll M. H. R., Kley W., Picogna G., 2017, A&A, 599, L6
  * Taki et al. (2016) Taki T., Fujimoto M., Ida S., 2016, A&A, 591, A86
  * Turrini et al. (2019) Turrini D., Marzari F., Polychroni D., Testi L., 2019, ApJ, 877, 50
  * Urpin (2003) Urpin V., 2003, A&A, 404, 397
  * Woitke & Helling (2003) Woitke P., Helling C., 2003, A&A, 399, 297