11institutetext: Indian Institute of Technology Gandhinagar, India
22institutetext: Bosch Research and Technology Centre, Bangalore, India
22email<EMAIL_ADDRESS><EMAIL_ADDRESS>

# Sparse Graph Representations for Procedural Instructional Documents

Shruti Singh∗ 11 Rishabh Gupta 22

###### Abstract

Computation of document similarity is a critical task in various NLP domains
that has applications in deduplication, matching, and recommendation. Graph
representations such as Joint Concept Interaction Graph (JCIG) represent a
pair of documents as a joint undirected weighted graph. JCIGs facilitate an
interpretable representation of document pairs as a graph. However, JCIGs are
undirected, and don’t consider the sequential flow of sentences in documents.
We propose two approaches to model document similarity by representing
document pairs as a directed and sparse JCIG that incorporates sequential
information. We propose two algorithms inspired by Supergenome Sorting and
Hamiltonian Path that replace the undirected edges with directed edges. Our
approach also sparsifies the graph to $O(n)$ edges from JCIG’s worst case of
$O(n^{2})$. We show that our sparse directed graph model architecture
consisting of a Siamese encoder and GCN achieves comparable results to the
baseline on datasets not containing sequential information and beats the
baseline by ten points on an instructional documents dataset containing
sequential information.

###### Keywords:

Document Similarity Concept Interaction Graphs GCN

**footnotetext: work done during internship at Bosch

##  1 Introduction and Related Works

Document similarity has applications in detecting duplicate documents
collected from diverse sources, paraphrase detection, and recommendation of
similar articles. Several datasets such as quora question pairs [22],
microsoft research paraphrase corpus [8], STS benchmark [4] for sentence
similarity exist that address the task of semantic similarity computation.
Recently, all these dataset leaderboards are dominated by transformer-based
models [21]. However, all these datasets contain short documents (few
sentences only) and transformer-based models don’t scale quadratically to long
documents.

Joint concept interaction graphs (JCIG) proposed by [12] represents a pair of
documents as a joint undirected weighted graph, where nodes of the graph
represent the concepts in the documents, and the edges represent the
interaction between the concepts. JCIGs show state-of-the-art performance in
detecting similar chinese news articles, also beating a fine-tuned BERT
baseline. However, JCIGs don’t encode sequential information in documents
which is crucial for certain domains such as recipes, instruction manuals, and
compliance documents. In such domains, sequence of steps is crucial, and
ignoring the order and following steps randomly can lead to damage of devices
or undesirable end products. Hence the need to encode the sequential
information while matching documents for such domains is important. Several
procedural text datasets have been curated such as RecipeQA [23] consisting of
recipe documents, EMQAP [15] of emanual documents, and ProPara [6] designed to
track changing entity states in procedural text. However, all these works
model the task in question-answering format, while our work focuses on the
document similarity task. In this work, we focus on JCIG [12] and propose
improvements over it to encode sequential information.

(a) Directed JCIG construction from document pair.

(b) Siamese encoder for node representations.

(c) GCN encoder and MLP classifier.

Figure 1: Overall pipeline of the document similarity architecture that has
three major steps: (a) Directed JCIG construction, (b) Siamese encoder, and
(c) GCN encoder and MLP classifier.

We propose to add directions to JCIGs that encode the sequential information
and then sparsify the graph by deducing the dominant direction ($D^{3}$) of
information flow between a pair of nodes. We propose two algorithms for
$D^{3}$: Supergenome Sorting (SGS) and Hamiltonian Path (HP) and showcase that
our model beats the JCIG baseline by 10 points on the iFixit dataset which
consists of instruction manuals for appliances. For the Chinese news articles
datasets (CNSE and CNSS), the sequential information is not crucial, and hence
we don’t observe significant improvements but also don’t observe any drop in
performance. The paper is organized as follows. We discuss JCIG construction
in Section 2. We describe the importance of encoding sequential information in
Section 3, direction pseudograph in Section 4, and SGS and HP algorithms in
Section 5. This is followed by a description of datasets and experiments in
Section 6 and results in Section 7. A schematic representation of our
architecture is presented in Figure 1.

##  2 Joint Concept Interaction Graph Representation

We briefly discuss the Joint Concept Interaction Graph (JCIG) [12] for
representing a pair of documents along with our ablation studies, as it forms
the fundamental basis on which we extend in our work. We begin by discussing
the construction of a concept interaction graph for a single document in
Section 2.1, followed by the construction of a joint concept interaction graph
(JCIG) for a pair of documents in Section 2.2.

###  2.1 Concept Interaction Graph

A Concept Interaction Graph (CIG) represents a document $\mathcal{D}_{A}$ as
an undirected weighted graph $\mathcal{G}_{A}(V,E)$, where $V$ is the vertex
set and $E$ is the edge set. The vertex set consists of keywords
$\mathcal{K}_{A}$ extracted from document $\mathcal{D}_{A}$. Next, we describe
the primary steps in the construction of a CIG.

---

Summa |  RAKE |  TextRank |  KeyBERT
---|---|---|---
shell gently |  camera using tweezers |  difficulty |  cable gently
opener lever |  gently pry back cover away |  spudger or plastic opener |  motherboard away
orange |  probably need tweezers |  gold |  screwdriver gently

Table 1: For a sample document on ‘camera replacement’ from the iFixit
dataset, we present flawed keywords extracted by Summa, RAKE, TextRank, and
KeyBERT. Based on manual inspection of keywords extracted by each algorithm,
we went ahead with the TextRank implementation by [16]. The document on the
left highlights the keywords extracted by Summa that are spread across
consecutive sentences (marked with underline), which are not optimal keywords.

#### Keyword Extraction

Keywords are extracted from document $\mathcal{D}_{A}$, using popular keyword
extraction algorithms: RAKE [19], Summa [1], TextRank [13, 16], KeyBERT [11].
Let the keywords extracted in this stage for a document $\mathcal{D}_{A}$ be
denoted by $\mathcal{K}_{A}=\\{kw_{1},kw_{2},...,kw_{n}\\}$.

TextRank [13, 16] constructs an undirected unweighted graph whose nodes are
tokens in the text (tokens filtered via syntactic rules to remove stopwords),
and edges are added between two nodes if they co-occur within a window of N
tokens. PageRank [3] is used to compute the score for each vertex, and top-k
vertices are selected as keywords. KeyBERT [11] computes the document
representation and n-gram phrase representations of n-grams in the document
via BERT [7]. It computes the cosine similarity between the document and
n-gram phrase representations and selects the most similar n-grams as
keywords. KeyBERT requires the range of n-grams to be considered and stopwords
to be ignored. We manually inspected the keywords extracted with RAKE, Summa,
TextRank, and KeyBERT for different configurations and observed that the most
optimal keywords are extracted via Summa. KeyBERT selects verb phrases, RAKE
has a preference for long phrases (usually 3 or 4-grams). Summa performs
similarly to TextRank with SpaCy backend, but Summa has some unsuitable
phrases which are spread across consecutive sentences. We present a sample of
flawed keywords extracted via each algorithm in Table 1. Based on manual
analysis of keywords from multiple documents, we employed TextRank
implementation by [16].

#### Keyword Graph ($\mathcal{KG}$)

After extraction of keywords in a document is Keyword Graph ($\mathcal{KG}$)
construction. The vertices of the Keyword Graph for the document
$\mathcal{D}_{A}$ are initialized from the keyword set $\mathcal{K}_{A}$. An
edge is added between keywords $kw_{i}$ and $kw_{j}$ if they co-occur in the
same sentence. Let $S_{A}$ be the set of sentences in the document
$\mathcal{D_{A}}$, and each sentence $s_{i}\in S_{A}$ is a set of words and
keywords. $\mathcal{KG}_{A}(V)$ and $\mathcal{KG}_{A}(E)$ denote the vertex
set and edge set of $\mathcal{KG}_{A}$.

$\displaystyle S_{A}$ $\displaystyle=$ $\displaystyle[s_{1},s_{2},...,s_{m}]$
$\displaystyle s_{i}$ $\displaystyle=$
$\displaystyle\\{w_{1},w_{2},...,kw_{j},kw_{k},..\\},$ $\displaystyle s_{i}\in
S_{A}$ $\displaystyle\mathcal{KG}_{A}(V)$ $\displaystyle=$
$\displaystyle\mathcal{K}_{A}=\\{kw_{1},kw_{2},...,kw_{n}\\}$
$\displaystyle\mathcal{KG}_{A}(E)$ $\displaystyle=$
$\displaystyle\\{(v_{i},v_{j}):v_{i},v_{j}\in\mathcal{K}_{A}\text{ and }$
$\displaystyle v_{i},v_{j}\in s_{p},s_{p}\in S_{A}\\}$

#### Concept Detection for CIG Node Construction

The $\mathcal{KG}$ captures the co-occurrence of keywords in sentences.
However, since there could be multiple keywords in a long document, the
vertices are highly specific and don’t represent the primary themes/concepts
in the document. The $\mathcal{KG}$ nodes are constructed at a fine
granularity, while the CIG envisons nodes to represent major concepts in the
document. Additionally, the $\mathcal{KG}$ is densely connected, and cannot be
directly used as the concept interaction representation graph. To reduce the
graph size, highly similar nodes in the $\mathcal{KG}$ can be clustered, and
each cluster can represent a single node. To cluster highly correlated
keywords, i.e., keywords that discuss similar concepts, community detection on
the Keyword Graph $\mathcal{KG}$ is employed. $\mathcal{KG}$ contains densely
connected components as keywords discussing a similar idea frequently co-occur
in the same sentence. The community detection algorithm segregates keywords
into clusters (or communities) by finding densely connected subgraphs. The
individual cluster or community can be intuitively understood as a _‘Concept’_
, such that concept $C^{i}=\\{kw_{m},kw_{n},kw_{o},...,kw_{p}\\}$ is a bag of
keywords. The union of all concept communities of a document $\mathcal{D}_{A}$
is equal to the set of keywords $\mathcal{K}_{A}$. However, a keyword could
belong to multiple communities or exclusively to just one community, depending
on the community detection algorithm being used. Let us denote the set of
concepts (or communities) for $\mathcal{D}_{A}$ as $\mathcal{C}_{A}$.

$\displaystyle\mathcal{C}_{A}=\bigcup_{i=1}^{p}C_{A}^{i}=\mathcal{K}_{A}=\mathcal{KG}_{A}(V)$

As part of the ablation study, we experiment with various community detection
algorithms, namely Girvan Newmann [10], Greedy Modularity [5], Louvain [2],
and Leiden [20]. The results indicate that the community detection algorithm
does not impact the performance on the iFixit dataset, i.e. all community
detection algorithms perform comparably.

Each community is a bag of keywords. To construct the Concept Interaction
Graph (CIG) of document $\mathcal{D}_{A}$, the vertex set $CIG_{A}(V)$ is
initialized with the communities, i.e., CIG vertices are the communities
detected in the keyword graph $\mathcal{KG}_{A}$. In addition, a dummy vertex
$C^{\text{D}}$ is also added (details in the upcoming paragraph on Sentence
Assignment).

$\displaystyle
CIG_{A}(V)=\\{C_{A}^{i}:C_{A}^{i}\in\mathcal{C}_{A}\\}\cup\\{C_{A}^{\text{D}}\\}$

#### Sentence Assignment

Each node in the CIG is a concept, which is a bag of keywords. The bag of
keywords is represented using a TF-IDF vector. Based on the TF-IDF vector
similarity, each sentence in the document is assigned to concept nodes if the
similarity is greater than a threshold $\text{sim}_{\text{thresh}}$.
$\text{sim}_{\text{thresh}}$ is a predefined threshold used for filtering out
low-quality sentence-concept matches. Sentences that are not assigned to any
concept are matched to a dummy vertex $C^{\text{D}}$ in the graph. The
sentence assignment function $\mathcal{F}$ for $CIG_{A}(V,E)$ of
$\mathcal{D}_{A}$, assigns sentences $s_{i}\in S_{A}$ to a set of concept
vertices $SA_{i}\subseteq\mathcal{C}_{A}$ or the dummy vertex
$C_{A}^{\text{D}}$.

$\displaystyle\mathcal{F}$ $\displaystyle:S_{A}\to\mathcal{C}_{A}\cup
C_{A}^{\text{D}}$ $\displaystyle SA_{i}$ $\displaystyle=\\{C_{A}^{j}:1\leq
j\leq|\mathcal{C}_{A}|,$
$\displaystyle\text{TFIDFSimilarity}(s_{i},C_{A}^{j})\geq\text{sim}_{\text{thresh}}\\}$

$SA_{i}$ denotes the set of concepts that a sentence $s_{i}$ is assigned to
based on the threshold. A sentence can be assigned to more than one concept.
If the similarity of a sentence is less than the threshold similarity
$\text{sim}_{\text{thresh}}$ with all the concept nodes, then it is assigned
to the dummy vertex $C_{A}^{\text{D}}$.

$\mathcal{F}(s_{i})=\begin{cases}SA_{i},&SA_{i}\neq\phi\\\
C_{A}^{\text{D}},&SA_{i}=\phi\end{cases}$

#### Weighted Edge Computation

The Sentence Assignment stage assigns sentences to the concept vertices. Each
vertex is represented by the TF-IDF vector of all the sentences assigned to
it. Weighted edges of format (source, target, weight) are added to the CIG of
document $\mathcal{D_{A}}$:

$\displaystyle CIG_{A}(E)=\\{(v_{i},v_{j},w):v_{i},v_{j}\in CIG(V),$
$\displaystyle i\neq j,w=\text{TFIDFSim}(v_{i},v_{j})\\}$

To summarize, keywords are extracted and clustered using community detection
to form concepts, which form the vertices in the CIG. Based on the occurrence
of keywords in sentences, each sentence is assigned to a concept node.
Sentences assigned to a concept node are utilized to form node representations
using TF-IDF vectors. Weighted edges are added to the graph by computing the
similarity between TF-IDF node representations.

###  2.2 Joint Concept Interaction Graph

CIGs represent a document as an undirected weighted graph. A joint concept
interaction graph, motivated by CIGs, represents a pair of documents by
capturing the similarities between the documents on different topics. The
vertex set of JCIG for a pair of documents $\mathcal{D}_{A}$ and
$\mathcal{D}_{B}$, represented as $JCIG_{AB}(V)$, is the union of the vertices
of $CIG_{A}$ and $CIG_{B}$. Similar to sentence assignment $\mathcal{F}$ in
CIG construction, the sentences in documents $\mathcal{D}_{A}$ and
$\mathcal{D}_{B}$ are assigned to the vertices of the $JCIG_{AB}$. Weighted
edge construction also follows the process of similarity computation between
TF-IDF vertex representations and edges with weights greater than a threshold
$w_{\text{thesh}}$ are added to the JCIG. However, it should be noted that the
TF-IDF node representations of vertices in the JCIG consist of sentences from
both documents, unlike CIG, which has sentences only from one document. For
JCIG, the formulations of the vertex set $JCIG_{AB}(V)$, edge set
$JCIG_{AB}(E)$, and the sentence assignment function $\mathcal{F}$ are as
follows:

$\displaystyle\mathcal{F}:S_{A}\cup S_{B}\to$
$\displaystyle\;\mathcal{C}_{A}\cup\mathcal{C}_{B}\cup C_{AB}^{D}$
$\displaystyle JCIG_{AB}(V)=$ $\displaystyle CIG_{A}(V)\cup CIG_{B}(V)$
$\displaystyle JCIG_{AB}(E)=$ $\displaystyle\\{(v_{i},v_{j},w):$
$\displaystyle v_{i},v_{j}\in JCIG_{AB}(V),i\neq j,$ $\displaystyle
w=\text{TFIDFSimilarity}(v_{i},v_{j}),$ $\displaystyle w\geq
w_{\text{thesh}}\\}$

JCIG vertices capture the important topics in the documents. Sentences
assigned to concepts capture the sentences from each document that discuss the
topic. However, the sequential nature of the document is destroyed in the
graph. For certain domains, such as instruction manuals, recipes, legal
documents, and compliance documents, the sequence of steps is crucial. To
retain the sequential information, we propose adding directed edges to the
graph. To this end, we propose our framework for encoding sequential
information in JCIGs.

##  3 Encoding Sequential Information

Sequential information is crucial in certain domains of text such as recipe
documents, instruction manuals, and compliance documents. Often these
documents contain a sequence of steps/instructions/actions that necessarily
need to be performed in the specified sequence. For example, an instruction
manual on replacing keys in a keyboard specifically mentions that the right
hook should be unlocked first, followed by lifting of the left pin hook. The
sequence is important, and the reverse procedure won’t help, and in some cases
could also lead to damage to the device. Hence, similarity computation of a
pair of documents should consider the sequence of information flow in each
document while reasoning if they are same or not. As we have observed, JCIG
does not consider sequential flow of information and the information is lost
in JCIG representation. To overcome this limitation, we propose our algorithm
of adding directed edges that encode the sequential information. The process
involves two components: (i) Direction Pseudograph (details in Section 4), and
(ii) Dominant Direction Deduction (details in Section 5).

##  4 Direction Pseudograph

Sentences are the most granular component of a document that can be
interpreted individually and convey information about the document. Hence, to
encode the directionality of information flow among sentences in JCIG, we
propose a simple solution that starts by constructing a pseudograph, which is
a multigraph i.e. multiple edges exist between vertices. In this section, we
will consider the construction of a pseudograph for a document that will add
directed edges as a representative of the sequential flow. Consider that we
initialize a pseudo concept graph, by adding all extracted concepts
(communities found via community detection on the keyword graph) as vertices.
Next, we add edges to the graph that encode the sequential information of the
document. Our motivation is to capture the sequential flow of information
between concepts. We start with the first and the second sentence in the
document and find the concepts they have been assigned to. Suppose $CS(s_{1})$
and $CS(s_{2})$ denote the concepts sets that sentence 1 and sentence 2 have
been assigned. Next, we add a directed unweighted edge from each concept
$cs_{i}\in CS(s_{1})$ to $cs_{j}\in CS(s_{2})$ to denote the direction of
information flow from $cs_{i}$ to $cs_{j}$. Then the same process is repeated
for sentence 2 and sentence 3. Similarly, we add directed unweighted edges
between concepts contained in each consecutive pair of sentences, which lead
to multiple edges between a pair of nodes. The edge set of the psuedograph is
defined as follows, where edge $e_{mn}=(m,n)$ denotes a directed unweighted
edge from node $m$ to $n$, and $|S_{A}|$ is the set of sentences in
$\mathcal{D}_{A}$:

$\displaystyle Pseudograph(V)$ $\displaystyle=JCIG(V)$ $\displaystyle
Pseudograph(E)$ $\displaystyle=\\{e_{mn}=(m,n):$ $\displaystyle
m\in\text{CS}(s_{i}),n\in\text{CS}(s_{i+1}),$ $\displaystyle m\neq n,1\leq
i\leq|S_{A}|-1\\}$

##  5 Dominant Direction Deduction for Sparse Directed Graphs

In this section, we discuss two algorithms for constructing sparse directed
graphs. To this end, we employ the pseudograph and CIGs or JCIGs. In the
previous step, we encoded the pairwise sequential information between concept
nodes based on the sequence of sentences in the document. However, the
pseudograph can be densely connected and edges between a pair of nodes could
exist in both directions, rendering the addition of sequential information
useless. Multiple sentences assigned to multiple concepts will lead to
multiple directed edges between concept vertices, which necessitates the
deduction of a dominant direction between a pair of nodes. Intuitively, the
dominant direction between a node pair captures the information flow from one
concept to another. We formulate the concept of dominant direction deduction
($D^{3}$), which finds the dominant flow of sequential information given
multiple parallel and opposite edges exist between a pair of nodes in the
direction pseudograph.

To reiterate, CIG is an undirected weighted graph for a single document and
JCIG is an undirected weighted graph for a pair of documents. As we want to
add directed edges from the psuedograph to CIG or JCIG to denote direction, we
design two variants of determining dominant edges. The dominant direction
deduction algorithm also leads to graph sparsification as multiple edges in
the psuedograph between a pair of nodes in both direction, are reduced to just
one edge in majority pairs. We represent the two variants of applying dominant
direction deduction ($D^{3}$) algorithm in Figure 2, (i) I-$D^{3}$: applying
individually to the CIG, and (ii) C-$D^{3}$: applying on the combined JCIG of
the document pair. Applying the variants with the two $D^{3}$ algorithms
Supergenome Sorting (SGS) and Hamiltonian Path (HP), we discuss all four
variants: Individual and Combined Supergenome Sorting (I-SGS and C-SGS
respectively, in Section 5.1) and Individual and Combined Hamiltonian Path
(I-HP and C-HP respectively, in Section 5.2). Our experiments reveal that C-HP
performs the best.

(a) Individual - Dominant Direction Deduction (I-$D^{3}$)

(b) Combined - Dominant Direction Deduction (C-$D^{3}$)

Figure 2: Variants of the dominant direction deduction algorithm. Applying
$D^{3}$ to a graph G, leads to a Seq G, a graph that encodes sequential
information. There are two $D^{3}$ algorithms, SGS and HP, so overall, we
propose four variants of $D^{3}$: I-SGS, C-SGS, I-HP, and C-HP.

###  5.1 Supergenome Sorting (SGS)

The graph resembles the supergenome graph presented in [9], and we use a
similar method for dominant direction deduction. The process starts by finding
the pair of concept nodes in the graph (CIG in case individual $D^{3}$ is
applied, or JCIG in case combined $D^{3}$) with the maximum number of edges,
say concept $C^{i}$ and $C^{j}$. The number of edges from $C^{i}$ to $C^{j}$
and from $C^{j}$ to $C^{i}$ are counted. The final dominant direction between
$C^{i}$ and $C^{j}$ corresponds to the direction with a greater number of
edges. The same procedure is repeated for each pair of nodes in the
pseudograph until the number of edges between a pair of nodes is reduced to
one. In rare cases when the count of edges between a pair of nodes is same in
both directions, a bidirectional edge is added. To summarize, supergenome
sorting algorithm takes the multigraph psuedograph and replaces multiple edges
between a pair of nodes with weighted directed edges, with edge weight equal
to the number of edges. Finally, the dominant direction between a pair of
nodes corresponds to the direction which has the higher weight edge. The
pseudograph with domainant directions received after this stage is then merged
with the original CIG or JCIG by taking an intersection of the edges. The
vertex set of pseudograph and CIG/JCIG is the same, i.e. concepts set. For the
intersecting edges, the direction is adopted from the pseudograph and the
weight of the edge is adopted from the CIG/JCIG, resulting in a sequential
weighted directed graph (Seq $\text{CIG}_{A}$ or Seq $\text{JCIG}_{AB}$). In
the worst case, $n^{2}$ edges in the original graph are reduced to
$\frac{n^{2}}{2}$ after applying supergenome sorting.

Individual SGS (I-SGS) is applied to CIG of two documents $\mathcal{D}_{A}$
and $\mathcal{D}_{B}$ individually. I-SGS constructs the pseudographs
(Pseudograph(A) and Pseudograph(B)) individually for the documents
$\mathcal{D}_{A}$ and $\mathcal{D}_{B}$ from their CIGs, $\text{CIG}_{A}$ and
$\text{CIG}_{B}$. I-SGS updates the edges in the pseudographs, by adding the
dominant edges and removing other edges. After merging $\text{CIG}_{A}$ with
Pseudograph(A) and $\text{CIG}_{B}$ with Pseudograph(B) (where pseudograph is
updated with dominant directions), resulting in Seq $\text{CIG}_{A}$ and Seq
$\text{CIG}_{B}$, which are finally combined to construct the Seq
$\text{JCIG}_{AB}$ as discussed in Section 2.2. On the other hand, Combined
SGS (C-SGS) constructs the Pseudograph(AB) from the $\text{JCIG}_{A}B$,
applies the supergenome sorting on the joint psuedograph of document pairs,
i.e. Pseudograph(AB), to get the sparse directed Pseudograph(AB). The
resultant Pseudograph(AB) is intersected with $\text{JCIG}_{A}B$ resulting in
a Seq $\text{JCIG}_{A}B$. It should be noted that the steps of supergenome
sorting can be applied without any modification to the pseudograph of a CIG or
JCIG as the steps are specific to a multigraph. The resulting pseudograph
contains the dominant directions (in the form of a single directed edge
between a pair of nodes), which is used to add directions to the final JCIG.

###  5.2 Hamiltonian Path (HP)

A Hamiltonian Path in a graph is a path that visits each vertex only once.
Given the pseudograph, which is a directed unweighted multigraph, each edge is
assigned a weight. For an edge e = (m, n), denoting an edge from concept node
m to concept node n, the weight is computed as the cosine similarity between
the TF-IDF representations of concept m and concept n. Similar to the sentence
assignment phase, each concept is represented by TF-IDF vectors based on the
sentences that are assigned to the concept. Then the directed weighted
multigraph is converted to a directed weighted simple graph by adding the
weights of parallel edges between the same pair of nodes. Note that this is
similar to constructing an edge with weight equivalent to the product of
number of edges (in the same direction) between concept node pair and cosine
similarity of concept nodes. It should be noted that the original pseudograph
is densely connected and the resultant pseudograph after applying HP is a
tournament graph in majority cases; and a hamiltonian path exists in a
tournament graph. However, in cases the pseudograph is not a tournament graph,
a simple addition of cosine similarity weighted edge between concept node
pairs that don’t have any edge, can convert the graph into a tournament graph.
The resultant directed weighted simple pseudograph contains dominant
directions, and is then merged with the CIG or JCIG in a similar way discussed
in SGS in Section 5.1.

In the worst case, $n^{2}$ edges in the original graph are reduced to $n$
edges after applying hamiltonian path algorithm. This sparisifies the graph
while adding directed edges which encode sequential information.
Sparsification leads to an advantage in the downstream architecture, where a
GCN classifier operates on the graph, as we observe speedup in message passing
as the number of edges is significantly reduced.

###  5.3 Similarity Classification Architecture

The JCIG encodes the information content from a pair of documents on concepts,
which are major themes in the documents. Each concept node consists of
sentences from both the documents if the concept occurs in both the documents.
In case a concept occurs only in one document, then the node contains
information about that concept from a single document only. We employ the
architecture for document pair matching by [12]. The architecture consists of
a Siamese encoder [17] which learns match vectors for each vertex, followed by
a graph convolution layer, node vector aggregation layer, and a multi-layer
perceptron for binary classification. The schematic representation of the
entire architecture is presented in Figure 1. The graph nodes contain
sentences for both documents, encoded via Word2Vec [14] and passed to the
Siamese encoder. In case a concept does not contain sentences for a document,
an average vector (average of all words in the vocabulary) is used to
represent the document. Siamese encoder learns match vectors for each node by
computing similarities and differences among sentences from both documents
belonging to the concept. The learned match vector for each node is the node
feature and the directed edges from our $D^{3}$ algorithm gives the adjacency
matrix of the graph. We use a three-layer GCN, and aggregate the hidden node
representations using an aggregate layer that computes mean of hidden node
vectors in the last layer. Finally, the mean of node vectors is passed to an
MLP for binary classification, predicting if the document pairs are similar or
not. Let $c_{A}(v)$ and $c_{B}(v)$ denote the context vectors obtained via the
siamese encoder for concept node $v$, for documents $\mathcal{D}_{A}$ and
$\mathcal{D}_{B}$. The siamese encoder takes the Word2Vec representations of
sentences from both documents, encodes the sentences assigned to node $v$ from
each of the documents, and a match vector is constructed as follows, where
$\circ$ denotes Hadamard product:

$\displaystyle m_{AB}(v)=(|c_{A}(v)-c{B}(v)|,c_{A}(v)\circ c_{B}(v))$

##  6 Datasets and Experimental Setup

We discuss the dataset details and the experimental setup in the upcoming
subsections.

###  6.1 CNSE, CNSS, and iFixit Dataset

CNSE and CNSS Dataset: The Chinese News Same Event dataset (CNSE) and Chinese
News Same Story dataset (CNSS) are curated by [12] to evaluate the performance
of JCIG in document matching. The dataset contains news stories collected from
news providers such as Sina, Sohu, Tencent, WeChat, etc., and the task is to
identify stories discussing the same event or story. The dataset has a pair of
news articles with binary labels denoting if the two articles are discussing
the same event/story. CNSE and CNSS datasets are not sensitive to sequential
information, and hence we did not initially propose to use our methodology for
such documents. However, we evaluate our methodology on CNSE and CNSS dataset
and show that performance does not deteriorate on these datasets in comparison
to the JCIG baseline.

iFixit Dataset: We curate instruction manuals from the iFixit website, which
contains step-by-step guides in English for repair and assembly of appliances
and devices across fifteen categories. In this work, we only curate manuals
for the category ‘Appliances’, consisting of 1333 manuals across twelve
categories such as refrigerator, washer, water dispenser, kitchen appliances,
etc. The dataset contains manuals that discuss the same procedure (such as
assembly, disassembly, and replacement of a specific part) for appliances from
different brands. For example, there are ‘water filter replacement’ guides for
Whirpool and Samsung refrigerators. Similarly, there are guides for other
tasks such as screen replacement, thermal fuse replacement, condenser
cleaning, etc. which are organized by different subjects. Each guide also has
a brand tag, which indicates the brand of the appliance or device and tags
about the components focused in the guide. To construct positive document
pairs from the iFixit Appliances, we filter guides with the same subject and
having at least two common tags. We construct negative document pairs by
grouping documents that don’t have the same subject.

###  6.2 Experimental Setup

Our model is implemented using PyTorch and we follow standard 60-20-20 data
splits for training, validation, and test set. We utilized 3 Quadro RTX 5000
GPU (16G memory) for running our experiments. We train all models for 50
epochs with early stopping. We run ablations with learning rates: 0.01, 0.005,
0.001, 0.0005, and 0.0001. We achieve the best results with a learning rate of
0.0005 across all configurations. We also experimented with 2, 3, and 4 GCN
layers, and 3-layer GCN performed the best. We train word2vec using gensim
library [18], and experiment with word vectors of sizes 100 and 200. All
documents were tokenized with NLTK library before word2vec training. We use
community detection algorithms available in NetworkX during concept detection
stage.

Method$\to$ | JCIG | SGS | HP
---|---|---|---
Dataset$\downarrow$ | Acc | F1 | Acc | F1 | Acc | F1
CNSE | 80.8 | 78.5 | 81.1 | 78.7 | 81.2 | 79.0
CNSS | 86.4 | 86.7 | 86.3 | 86.6 | 86.6 | 87.1
iFixit | 60.9 | 37.8 | 64.4 | 48.1 | 70.7 | 58.5
Table 2: Results on the CNSE, CNSS, and iFixit dataset. HP shows significant
improvement in accuracy and F1 scores over the JCIG baseline for the iFixit
dataset. We observe that HP performs superior to SGS in all cases. SGS and HP
perform equally well on the CNSE and CNSS datasets, but significant
improvements are not observed as sequential information is not crucial in news
articles.

##  7 Results and Observations

Even though CNSE and CNSS datasets don’t contain any crucial sequential
information, our proposed SGS and HP based dominant direction deduction
methods achieve performance comparable to the JCIG baseline. The iFixit
dataset contains instruction manuals for appliances and sequential information
is critical. HP algorithm shows significant improvement in accuracy (Table 2)
and achieves an improvement in accuracy of ten points in comparison to JCIG
baseline. SGS also outperforms JCIG on the iFixit dataset by four points.
During ablation, we observed that HP outperforms SGS and C-$D^{3}$ outperforms
I-$D^{3}$.

##  8 Conclusion

Encoding sequence information in graph representations of document pairs is
crucial for recipe documents, instruction manuals, and compliance documents.
While JCIGs perform well for detecting similar news articles, we showcase that
they perform poorly for instruction manuals in the iFixit dataset. We propose
Supergenome sorting and Hamiltonian path algorithms for deducing dominant
direction between graph nodes and also sparsify the graph by adding directed
edges. Our proposed HP algorithm beats the JCIG baseline by ten points on the
iFixit dataset.

##  9 Limitations

The CNSE and CNSS datasets are Chinese news articles, while the iFixit dataset
is in English. We posit that the performance of our methods will carry forward
similarly to other languages. However, our method is limited by keyword
detection algorithms and the availability of text-processing libraries in
other languages which can significantly impact the downstream performance.

##  10 Ethics Statement

The curated dataset of iFixit manuals does not contain any information that is
related to individuals, and neither does it contain any offensive content. The
proposed method for document similarity has no ethical considerations.

## References

  * [1] Barrios, F., López, F., Argerich, L., Wachenchauzer, R.: Variations of the similarity function of textrank for automated summarization. CoRR abs/1602.03606 (2016), http://arxiv.org/abs/1602.03606
  * [2] Blondel, V.D., Guillaume, J.L., Lambiotte, R., Lefebvre, E.: Fast unfolding of communities in large networks. Journal of statistical mechanics: theory and experiment 2008(10), P10008 (2008)
  * [3] Brin, S., Page, L.: The anatomy of a large-scale hypertextual web search engine. Computer Networks and ISDN Systems 30(1), 107–117 (1998). https://doi.org/https://doi.org/10.1016/S0169-7552(98)00110-X, https://www.sciencedirect.com/science/article/pii/S016975529800110X, proceedings of the Seventh International World Wide Web Conference
  * [4] Cer, D., Diab, M., Agirre, E., Lopez-Gazpio, I., Specia, L.: SemEval-2017 task 1: Semantic textual similarity multilingual and crosslingual focused evaluation. In: Proceedings of the 11th International Workshop on Semantic Evaluation (SemEval-2017). pp. 1–14. Association for Computational Linguistics, Vancouver, Canada (Aug 2017). https://doi.org/10.18653/v1/S17-2001, https://aclanthology.org/S17-2001
  * [5] Clauset, A., Newman, M.E., Moore, C.: Finding community structure in very large networks. Physical review E 70(6), 066111 (2004)
  * [6] Dalvi, B., Huang, L., Tandon, N., Yih, W.t., Clark, P.: Tracking state changes in procedural text: a challenge dataset and models for process paragraph comprehension. In: Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers). pp. 1595–1604. Association for Computational Linguistics, New Orleans, Louisiana (Jun 2018). https://doi.org/10.18653/v1/N18-1144, https://aclanthology.org/N18-1144
  * [7] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: 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). pp. 4171–4186. Association for Computational Linguistics, Minneapolis, Minnesota (Jun 2019). https://doi.org/10.18653/v1/N19-1423, https://aclanthology.org/N19-1423
  * [8] Dolan, W.B., Brockett, C.: Automatically constructing a corpus of sentential paraphrases. In: Proceedings of the Third International Workshop on Paraphrasing (IWP2005) (2005), https://aclanthology.org/I05-5002
  * [9] Gärtner, F., Höner zu Siederdissen, C., Müller, L., Stadler, P.F.: Coordinate systems for supergenomes. Algorithms for Molecular Biology 13(1), 1–19 (2018)
  * [10] Girvan, M., Newman, M.E.: Community structure in social and biological networks. Proceedings of the national academy of sciences 99(12), 7821–7826 (2002)
  * [11] Grootendorst, M.: Keybert: Minimal keyword extraction with bert. (2020). https://doi.org/10.5281/zenodo.4461265, https://doi.org/10.5281/zenodo.4461265
  * [12] Liu, B., Niu, D., Wei, H., Lin, J., He, Y., Lai, K., Xu, Y.: Matching article pairs with graphical decomposition and convolutions. In: Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics. pp. 6284–6294. Association for Computational Linguistics, Florence, Italy (Jul 2019). https://doi.org/10.18653/v1/P19-1632, https://aclanthology.org/P19-1632
  * [13] Mihalcea, R., Tarau, P.: TextRank: Bringing order into text. In: Proceedings of the 2004 Conference on Empirical Methods in Natural Language Processing. pp. 404–411. Association for Computational Linguistics, Barcelona, Spain (Jul 2004), https://aclanthology.org/W04-3252
  * [14] Mikolov, T., Chen, K., Corrado, G., Dean, J.: Efficient estimation of word representations in vector space. arXiv preprint arXiv:1301.3781 (2013)
  * [15] Nandy, A., Sharma, S., Maddhashiya, S., Sachdeva, K., Goyal, P., Ganguly, N.: Question answering over electronic devices: A new benchmark dataset and a multi-task learning based QA framework. In: Findings of the Association for Computational Linguistics: EMNLP 2021. pp. 4600–4609. Association for Computational Linguistics, Punta Cana, Dominican Republic (Nov 2021). https://doi.org/10.18653/v1/2021.findings-emnlp.392, https://aclanthology.org/2021.findings-emnlp.392
  * [16] Nathan, P.: PyTextRank, a Python implementation of TextRank for phrase extraction and summarization of text documents (2016). https://doi.org/10.5281/zenodo.4637885, https://github.com/DerwenAI/pytextrank
  * [17] Neculoiu, P., Versteegh, M., Rotaru, M.: Learning text similarity with Siamese recurrent networks. In: Proceedings of the 1st Workshop on Representation Learning for NLP. pp. 148–157. Association for Computational Linguistics, Berlin, Germany (Aug 2016). https://doi.org/10.18653/v1/W16-1617, https://aclanthology.org/W16-1617
  * [18] Řehůřek, R., Sojka, P.: Software Framework for Topic Modelling with Large Corpora. In: Proceedings of the LREC 2010 Workshop on New Challenges for NLP Frameworks. pp. 45–50. ELRA, Valletta, Malta (May 2010), http://is.muni.cz/publication/884893/en
  * [19] Rose, S., Engel, D., Cramer, N., Cowley, W.: Automatic Keyword Extraction from Individual Documents, chap. 1, pp. 1–20. John Wiley and Sons, Ltd (2010). https://doi.org/https://doi.org/10.1002/9780470689646.ch1, https://onlinelibrary.wiley.com/doi/abs/10.1002/9780470689646.ch1
  * [20] Traag, V.A., Waltman, L., Van Eck, N.J.: From louvain to leiden: guaranteeing well-connected communities. Scientific reports 9(1), 5233 (2019)
  * [21] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017)
  * [22] Wang, Z., Hamza, W., Florian, R.: Bilateral multi-perspective matching for natural language sentences. In: Proceedings of the 26th International Joint Conference on Artificial Intelligence. p. 4144–4150. IJCAI’17, AAAI Press (2017)
  * [23] Yagcioglu, S., Erdem, A., Erdem, E., Ikizler-Cinbis, N.: RecipeQA: A challenge dataset for multimodal comprehension of cooking recipes. In: Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing. pp. 1358–1368. Association for Computational Linguistics, Brussels, Belgium (Oct-Nov 2018). https://doi.org/10.18653/v1/D18-1166, https://aclanthology.org/D18-1166


# Rotation reduces convective mixing in Jupiter and other gas giants

J. R. Fuentes  Department of Applied Mathematics, University of Colorado
Boulder, Boulder, CO 80309-0526, USA  Evan H. Anders  Center for
Interdisciplinary Exploration and Research in Astrophysics, Northwestern
University, Evanston, Illinois 60201, USA  Andrew Cumming  Department of
Physics and Trottier Space Institute, McGill University, Montreal, QC H3A 2T8,
Canada  Bradley W. Hindman  Department of Applied Mathematics, University of
Colorado Boulder, Boulder, CO 80309-0526, USA  JILA, University of Colorado
Boulder, Boulder, CO 80309-0440, USA

###### Abstract

Recent measurements of Jupiter’s gravitational moments by the Juno spacecraft
and seismology of Saturn’s rings suggest that the primordial composition
gradients in the deep interior of these planets have persisted since their
formation. One possible explanation is the presence of a double-diffusive
staircase below the planet’s outer convection zone, which inhibits mixing
across the deeper layers. However, hydrodynamic simulations have shown that
these staircases are not long-lasting and can be disrupted by overshooting
convection. In this letter, we suggests that planetary rotation could be
another factor for the longevity of primordial composition gradients. Using
rotational mixing-length theory and 3D hydrodynamic simulations, we
demonstrate that rotation significantly reduces both the convective velocity
and the mixing of primordial composition gradients. In particular, for Jovian
conditions at $t\sim 10^{8}~{}\mathrm{yrs}$ after formation, rotation reduces
the convective velocity by a factor of 6, and in turn, the kinetic energy flux
available for mixing gets reduced by a factor of $6^{3}\sim 200$. This leads
to an entrainment timescale that is more than two orders of magnitude longer
than without rotation. We encourage future hydrodynamic models of Jupiter and
other gas giants to include rapid rotation, because the decrease in the mixing
efficiency could explain why Jupiter and Saturn are not fully mixed.

UAT Keywords: Jupiter (873), Saturn (1426), Solar system gas giant planets
(1191), Planetary interior (1248), Hydrodynamics (1963), Hydrodynamical
simulations (767), Convective zones (301)

##  1 Introduction

The history of the solar system is written in the interiors of the giant
planets. In particular, since Jupiter accounts for 75% of the Solar System’s
planetary mass, understanding its interior (or more precisely, determining the
distribution and total mass of heavy elements within the planet) can provide
important clues on the first stages of the formation of the solar system (see
the recent review by Miguel & Vazan, 2023). Traditionally, Jupiter has been
considered a fully convective planet, with a well-defined interface separating
a dense core made of heavy elements and the hydrogen-helium envelope (Pollack
et al., 1996). However, new data from the Juno mission indicates that Jupiter
does not have a well-defined central core of heavy elements, nor is it
homogeneously mixed (e.g., Bolton et al., 2017b, a; Wahl et al., 2017; Debras
& Chabrier, 2019; Howard et al., 2023). Interestingly, ring seismology
suggests that Saturn is likewise not homogeneously mixed (Mankovich & Fuller,
2021). These new observational constraints on the internal structure of gas
giants require updates to formation and evolution models (Helled et al.,
2022). Composition gradients decrease the efficiency of the heat transport
throughout the planet’s interior and thus affect its cooling, and mass-radius
relationship at a given age (e.g., Chabrier & Baraffe, 2007; Leconte &
Chabrier, 2012, 2013).

Primordial composition gradients are expected to naturally result from planet
formation (e.g., Helled & Stevenson, 2017; Müller et al., 2020; Stevenson et
al., 2022) or collisional events during the planet’s evolution (Liu et al.,
2019), but the persistence of composition gradients against strong convective-
mixing over evolutionary timescales is poorly understood. One mechanism for
maintaining a composition gradient appears in one-dimensional evolution
models, where convective staircases (i.e., multiple convective layers
separated by thin, stably-stratified, diffusive interfaces) form underneath
the outer convective envelope, preventing mixing deeper in the interior (e.g.,
Vazan et al., 2015, 2018; Stevenson et al., 2022). However, 3D simulations of
layered convection in incompressible flows have shown that convective
staircases fully mix on short timescales (e.g., Wood et al., 2013; Garaud,
2018). Recently, Fuentes et al. (2022) explored the formation of a convective
staircase beneath a growing convection zone, and found that compositional
mixing across the interface between the convection zone and the stable region
prevents the formation of the staircase. Similar results were found by Anders
et al. (2022) in the context of stellar convection, concluding that
compositional gradients are ineffective barriers against convective mixing
over evolutionary timescales.

Rotation is another factor that could prevent convective mixing in Jupiter and
other gas giants. Experiments and numerical simulations of thermal convection
have shown that rotation reduces the convective velocities when the flow has
small Rossby number, $\mathrm{Ro}$, the ratio of the rotation period to the
convective turnover time (e.g., Stevenson, 1979; Fernando et al., 1991; Barker
et al., 2014; Dietrich & Wicht, 2018; Aurnou et al., 2020; Hindman et al.,
2020). Further, numerical experiments wherein rapidly-rotating convection
zones grow via entrainment have shown that the vertical transport of buoyancy
and kinetic energy across the convective layer, and thus the entrainment rate,
are reduced by rotation (Julien et al., 1996a). Convective flows in Jupiter
and other gas giants can be significantly slow compared to rotation
($\mathrm{Ro}\lesssim 10^{-6}$, e.g. Guillot et al., 2004), so the effects of
rotation on the erosion of primordial composition gradients are likely to be
significant. Indeed, Leconte & Chabrier (2012) suggested that rotation would
hamper convection and impede the mixing of heavy elements.

In this Letter, we examine how rotation modifies the rate of mixing of a
stable fluid layer by a neighboring convection zone. In Section 2, we use
mixing-length theory to argue that rotation reduces the convective velocity
and the kinetic energy flux available for mixing. Then, in Section 3 we use
the scalings from rotating mixing length theory to show that rotation
decreases the rate of entrainment and mixing at the boundary between the
stable and unstable layers. In Section 4, we use 3D numerical simulations to
support the predictions from the previous sections. Finally, in Section 5 we
discuss the implications of our results for Jupiter and other gas giants.

##  2 Mixing-length theory for rotating convection

In this section, we summarize the mixing-length scalings for the convective
velocity and kinetic energy flux in non-rotating and rotating flows.

In the absence of rotation, convective motions are nearly isotropic with a
characteristic length scale, $\ell$, that is roughly the depth of the
convection zone (e.g., in a laboratory experiment) or a pressure scale height
(e.g., in stellar and planetary interiors). Further, the velocity scale of the
convective flows is approximately the overturning buoyant “free-fall”
velocity, $U_{\mathrm{NR}}=\sqrt{g\alpha\delta T_{\mathrm{NR}}\ell}$, where
the subscript $\mathrm{NR}$ means _no rotation_ , $g$ is the acceleration due
to gravity, and $\alpha=-(\partial\ln\rho/\partial T)|_{C}$ is the coefficient
of thermal expansion at constant composition $C$, and $\delta T_{\mathrm{NR}}$
is the characteristic temperature perturbation of convective parcels. The heat
flux carried by the convective motions is $F_{H}\sim\rho c_{P}\delta
T_{\mathrm{NR}}U_{\mathrm{NR}}$, where $\rho$ is the mass density, and $c_{P}$
is the specific heat capacity at constant pressure. Then, the convective
velocity can be estimated as

$U_{\mathrm{NR}}\sim\left(\dfrac{\alpha gF_{H}\ell}{\rho
c_{P}}\right)^{1/3}\,.$ (1)

Typical values for Jupiter’s deep interior at $t\sim 10^{8}~{}\mathrm{yrs}$
after formation, when the convective luminosity is high enough so that a large
fraction of the planet is expected to be mixed, $\ell\sim H_{P}\sim
10^{9}~{}\mathrm{cm}$ (a pressure scale height), $c_{P}/\alpha g\sim
10^{10}~{}\mathrm{cm}$, $\rho\sim 1~{}\mathrm{g~{}cm^{-3}}$, and $F_{H}\sim
6\times 10^{5}~{}\mathrm{cgs}$ (e.g., Cumming et al., 2018; Helled et al.,
2022), give $U_{\mathrm{NR}}\sim 40~{}\mathrm{cm\,s^{-1}}$ (and accordingly,
the convective turnover time is $\ell/U_{\mathrm{NR}}\sim 0.8~{}\mathrm{yr}$).

In the limit of rapid rotation, convection is quasi-geostrophic, and
convective flows are set by a balance between the Coriolis, inertial, and
buoyancy (Archimedean) forces (called CIA balance, e.g., Stevenson, 1979). In
this limit, a new length scale, $\ell_{\perp}$, associated with motions
perpendicular to the rotation vector, $\bm{\Omega}$, emerges (e.g., Barker et
al., 2014; Aurnou et al., 2020; Vasil et al., 2021). From CIA balance, one
obtains

$\dfrac{2\Omega
U_{\mathrm{R}}}{\ell}\sim\dfrac{U^{2}_{\mathrm{R}}}{\ell^{2}_{\perp}}\sim\dfrac{\alpha
g\delta T_{\mathrm{R}}}{\ell_{\perp}}\,,$ (2)

where the subscript R means _rotating_ and $\Omega$ is the angular frequency
of the rotation vector. Using the definition of the heat flux $F_{H}\sim\rho
c_{P}\delta T_{\mathrm{R}}U_{\mathrm{R}}$, and eliminating $\delta
T_{\mathrm{R}}$ and $\ell_{\perp}$ from Equations (2), it follows that

$\displaystyle U_{\mathrm{R}}\sim\left(\dfrac{g\alpha F_{H}}{\rho
c_{P}}\right)^{2/5}\left(\dfrac{\ell}{2\Omega}\right)^{1/5}\,\sim U_{\rm
NR}\left(\frac{U_{\rm NR}}{2\Omega\ell}\right)^{1/5},$ (3)

Plugging in the previous values for Jupiter, together with $2\pi/\Omega\sim
10$ hours, we obtain $U_{\mathrm{R}}\sim 6~{}\mathrm{cm\,s^{-1}}$. Rotation
therefore reduces the convective velocity by a factor of $\sim 6$; then, the
convective turnover time increases by a factor of $\sim 6.6$, giving
$\ell/U_{\mathrm{R}}\sim 5~{}\mathrm{yrs}$.

Although the effect of rotation on the convective velocity could be considered
minor in the sense that the resulting turnover time is still many orders of
magnitude smaller than the age of the planet, the important quantity for
entrainment is the kinetic energy flux, which is reduced by a much larger
factor. This can be seen from the ratio between the kinetic energy flux
$F_{\rm K}\sim(1/2)\rho U_{\mathrm{conv}}^{3}$ associated with the convective
motions and the thermal heat flux carried by convection. Without rotation,
Equation (1) gives $F_{\mathrm{K,NR}}/F_{H}=(1/2)\rho
U_{\mathrm{NR}}^{3}/F_{H}\sim(\alpha g\ell/c_{P})$. In the limit of rapid
rotation, Equation (3) gives

$\frac{F_{\mathrm{K,R}}}{F_{H}}\sim\left(\frac{\alpha
g\ell}{c_{P}}\right)\left(\frac{U_{\mathrm{R}}}{2\Omega\ell}\right)^{1/2}\sim\mathrm{Ro}^{1/2}\frac{F_{\rm
K,NR}}{F_{H}}\,,$ (4)

where $\mathrm{Ro}=U_{\mathrm{R}}/2\Omega\ell\ll 1$ is the Rossby number of
the flow based on the rotationally-constrained convective velocity. Then, _for
a given heat flux_ , rotation reduces the kinetic energy flux by a factor of
$\mathrm{Ro}^{1/2}\ll 1$. Using the values above for Jupiter, we estimate
$\mathrm{Ro}\sim 10^{-5}$, and consequently, the kinetic energy flux available
to mix the fluid is reduced by a factor of $\sim 4\times 10^{-3}$. This
reduction is important because the energy source for the upward transport and
mixing of heavier material is the kinetic energy of the convective motions. We
investigate this in more detail in Section 3.

##  3 Entrainment in rapidly rotating flows

Mixing of material (entrainment) across convective boundaries (e.g., a core-
envelope interface) has been investigated with great detail in laboratory
experiments and numerical simulations (e.g., Turner, 1968; Fernando, 1987;
Molemaker & Dijkstra, 1997; Fuentes & Cumming, 2020). In those studies, a
convection zone advances into a stable region, and the returning flows from
overshooting convection entrain fluid with the chemical composition of the
stable layer. This material rapidly mixes, so the convection zone grows over
time (see Anders & Pedersen, 2023, for a review of boundary mixing processes).
By reducing the convective velocities, rotation reduces the entrainment
(mixing) efficiency, and consequently, decreases the growth rate of the
convective layer.

Consider a plane-parallel fluid layer of depth $H$, stably-stratified with
composition decreasing linearly with height, $C_{0}(z)=|dC_{0}/dz|(H-z)$. Note
that $C$ could represent the spatial variation in mass fraction of heavy
elements, or the mean molecular weight of the fluid, as would be induced by
salt in water or heavy elements in hydrogen. Then, a heat flux $F_{0}$ escapes
through the top of the stable layer, such that convection is driven from above
and mixes the composition gradient from top to bottom. For simplicity, we
adopt the Boussinesq approximation (Spiegel & Veronis, 1960), wherein the
fluid density is constant everywhere except for buoyant perturbations, which
is a decent approximation in the vicinity of a convective boundary.

We now estimate the rate of entrainment $dh/dt$, where $h(t)$ is the depth of
the convection zone at time $t$. Lifting and mixing a quantity $\Delta C$ of
heavy fluid requires an energy expenditure $\Delta E=(1/6)\rho_{0}\beta\Delta
Cgh^{2}$ (see, e.g., Fuentes & Cumming, 2020), where mass conservation
dictates that $\Delta C=0.5|dC_{0}/dz|h$. Here, $\rho_{0}$ is the density of
the layer, and $\beta=(\partial\ln\rho/\partial C)_{T}$ is the coefficient of
compositional expansion (both quantities assumed constant). Differentiating
$\Delta E$ with respect to time, we write the rate of change of potential
energy due to mixing as $d(\Delta
E)/dt=(1/4)\rho_{0}\beta|dC_{0}/dz|gh^{2}\dot{h}$. We adopt the entrainment
hypothesis, which states that the rate of change of potential energy across
the layer is set by the kinetic energy flux of the convective motions
$F_{\mathrm{KE}}=(1/2)\rho_{0}U^{3}_{\mathrm{conv}}$ (Linden, 1975). It
follows that

$\dfrac{1}{4}\rho_{0}\beta
g\bigg{|}\dfrac{dC_{0}}{dz}\bigg{|}h^{2}\dfrac{dh}{dt}\approx\dfrac{1}{2}\rho_{0}U^{3}_{\mathrm{conv}}\,.$
(5)

Equation (5) assumes that all the kinetic energy is used to do work and lift
material from below, and that thermal effects do not significantly affect the
rate of change in the potential energy across the layer (a good approximation
since $\alpha\Delta T\ll\beta\Delta C$).

We see that given $U_{\mathrm{conv}}$, Equation (5) can be integrated to
obtain $h(t)$. Without rotation, we use $\rho=\rho_{0}$, $F_{H}=F_{0}$, and
$\ell=h$ and assume $U_{\mathrm{conv}}\approx
U_{\mathrm{NR}}\propto\ell^{1/3}$ (Equation 1), to obtain (within a factor of
order unity due to the assumptions above)111A full derivation including the
mixing efficiency, interfacial heat transport, and thermal effects on the
density can be found in Fuentes & Cumming (2020).

$h_{\mathrm{NR}}(t)\approx 2\left(\dfrac{\alpha
gF_{0}}{\rho_{0}c_{P}N^{3}}\right)^{1/2}\left(Nt\right)^{1/2}\,,$ (6)

where $N=\sqrt{g\beta|dC_{0}/dz|}$ is the buoyancy frequency. From Equation
(6), the rate at which the convection zone expands inward depends on how
quickly it cools down (which sets the convective velocity magnitude) and on
the magnitude of the initial composition gradient.

On the other hand, in the limit of rapid rotation we should use
$U_{\mathrm{conv}}\approx U_{\mathrm{R}}\propto\ell^{1/5}$ (Equation 3), and
integration of Equation (5) gives

$h_{\mathrm{R}}(t)\approx\left(\frac{24}{5}\right)^{5/12}\left(\dfrac{\alpha
gF_{0}}{\rho_{0}c_{P}N^{3}}\right)^{1/2}\left(\dfrac{N}{2\Omega}\right)^{1/4}\left(Nt\right)^{5/12}\,,$
(7)

or

$h_{\rm
R}(t)\approx\left[0.96\left(\dfrac{N}{2\Omega}\right)^{1/4}(N\,t)^{-1/12}\right]h_{\rm
NR}(t)\,.$ (8)

We can see that rotation slows down the propagation of the convection zone, it
reduces the power-law index from 0.5 to approximately 0.41, but the most
significant effect comes from the prefactor ($N\,t)^{-1/12}$, which can be
very small over long timescales.

It is instructive to compare the mixing timescales in both cases. From
Equations (6) and (7), the time required for convection to mix a stable layer
of vertical size $L$ (i.e., $h(t)=L)$ is

$\tau_{\mathrm{NR}}\approx N^{-1}\left(\dfrac{\rho_{0}c_{P}N^{3}L^{2}}{4\alpha
gF_{0}}\right)\,,$ (9)

in the non-rotating case, and

$\displaystyle\tau_{\mathrm{R}}$
$\displaystyle\approx\dfrac{5}{24}N^{-1}\left(\dfrac{\rho_{0}c_{P}N^{3}L^{2}}{\alpha
gF_{0}}\right)^{6/5}\left(\dfrac{2\Omega}{N}\right)^{3/5}\approx\dfrac{5}{6}\mathrm{Ro}^{-1/2}~{}\tau_{\mathrm{NR}}\,,$
(10)

in the rotating case. Therefore, the mixing timescale increases approximately
by a factor of $(5/6)\mathrm{Ro}^{-1/2}$, a commensurate change to the
decrease in the kinetic energy flux—compare Equations (4) and (10).

In a system like Jupiter where $\mathrm{Ro}\sim 10^{-5}$ at $t\sim
10^{8}~{}\mathrm{yrs}$ after formation, this corresponds to a mixing timescale
that is more than two orders of magnitude longer. This difference is
substantial, e.g., for $L\sim 10^{9}~{}\mathrm{cm}$, $N^{2}\sim g/H_{P}\sim
10^{-6}~{}\mathrm{s^{-2}}$, Equations (9) and (10) give
$\tau_{\mathrm{NR}}\sim 10^{8}~{}\mathrm{yrs}$ and $\tau_{\mathrm{R}}\sim
10^{10}~{}\mathrm{yrs}$, i.e., for the rotating case the mixing timescale
would be more than the age of the Solar System.

##  4 Numerical experiments and results

Figure 1: Volume renderings of the $z$ component of the vorticity
$\bm{\omega}=\nabla\times\bm{u}$ for the non-rotating case (left) and rotating
case (right). For both cases, the snapshots are shown at times when the bottom
of the convection zone has expanded beyond the initial location (from $z=1$ to
0.75). The Rossby number in the rotating case is $\mathrm{Ro_{sim}}\approx
0.065$, computed using the averaged flow velocity in the convection zone, and
thickness of the convection zone.
Figure 2: Horizontally-averaged profiles of composition $\overline{C}$,
convective composition flux $\overline{wC}$, and the kinetic energy flux in
the vertical direction $(1/2)\overline{w\bm{u}^{2}}$ (left, middle, and right
panel, respectively). The overline denotes that variables are averaged over
the horizontal directions. Note that $\bm{u}^{2}=u^{2}+v^{2}+w^{3}$, where
$u$, $v$, and $w$ are the $x$, $y$, and $z$ components of the velocity. The
results are shown at at a times where the convection boundary is located at
$z=0.75$ (i.e., $h=1.25$) for both runs. This time corresponds to $t=3750$ for
the non rotating case, and $t=11500$ for the rotating case. The background
colors distinguish between the convection zone and the stable region.

In the following, we use two hydrodynamic simulations (with and without
rotation) to test the theory above. We model 3D thermal convection in a plane-
parallel binary mixture of fluid (heavy and light elements). The domain
initially consists of two layers: a convection zone of uniform composition, on
top of a region stabilized by a composition gradient. The top of the
convection zone is constantly cooled so that convection propagates downwards.
We focus on how rotation affects the mixing of heavy material at the lower
convective boundary.

In our simulations, as in our derivation above, we assume that the length
scale of the fluid motions is much smaller than a density scale height, so
that the effects of curvature can be ignored and the Boussinesq approximation
is valid (Spiegel & Veronis, 1960). Although the density in giant planets can
vary by many orders of magnitude, this approximation fully captures nonlinear
mixing near the convective boundary, which is the focus of our study. In the
rotating case, we consider uniform rotation throughout the fluid, and assume
that gravity and rotation are aligned and point in the z direction,
$\bm{\Omega}=\Omega\hat{\bm{z}}$, and $\bm{g}=-g\hat{\bm{z}}$, i.e., the
simulation is at a polar latitude. The density of the fluid depends on both
composition and temperature fluctuations. In what follows, all results are
presented in dimensionless form. For details on the numerical setup,
nondimensionalization, and code, we refer the reader to Appendices A and B.

Figure 1 shows 3D snapshots of the $z$ component of the vorticity for both the
non-rotating and rotating simulations. Initially, in both cases, the
convection zone spans the upper half of the domain, $z\in[1,2]$. After
significant evolution ($t\sim 3750$ for the non-rotating case, and $t\sim
11500$ for the rotating case), the bottom of the convection zone has advanced
to $z=0.75$, mixing the primordial composition gradient. Note that in the
rotating case, it takes approximately 3 times longer for the zone to increase
its vertical extent by $\Delta h=0.25$. We show later that this difference
agrees with Equation (10). Rotation visibly alters the flow morphologies. In
the non-rotating case, the characteristic flow length scale is similar in all
directions, while in the rotating case, the horizontal length scale of the
flow is much smaller than the vertical scale (as expected from the Taylor-
Proudman theorem, e.g., Proudman, 1916; Taylor, 1917).

Figure 2 shows horizontally-averaged profiles of the composition (left),
vertical convective composition flux (middle), and vertical kinetic energy
flux (right); fluxes are depicted at the same instant in time as displayed in
Figure 1. Rotation does not affect the composition profile, which is
homogeneous in the convection zone with a transition to the primordial linear
distribution in the stable region. This is expected because mass conservation
requires that the amount of material transported and mixed in the convection
zone depends only on the thickness of the convective layer, $\Delta
C=0.5|dC_{0}/dz|h$. The initial composition in the convection zone is $C=0.5$,
and with $|dC_{0}/dz|=1$ and $\Delta h=0.25$ giving $\Delta C=0.125$, the
expected evolved convection zone composition is $C=0.625$, which is what we
observe (see left panel). In the rotating case, both the convective
composition flux and the kinetic energy flux are smaller by a factor of
$\approx 3$ when compared with the values of the non-rotating case (see middle
and right panels), supporting both our hypotheses that rotation decreases the
vertical transport of dense material and that this reduction is related to the
kinetic energy flux. Although we only focus on the compositional transport and
kinetic energy flux, it is worth mentioning that the thermal structure of the
fluid across the box is similar to that of the composition. The only
difference is that in the non-rotating case, the convection zone tends to be
isothermalized by strong turbulence, while in the rotating case, a small
temperature gradient tends to be sustained in the convection zone (this has
been reported in previous studies of thermal convection in rotating flows,
see, e.g., Julien et al., 1996b; Cheng et al., 2020).

Figure 3: Panel (a): Thickness of the convection zone $h$ as a function of
time $t$. The results are shown for both the rotating and non-rotating cases.
The solid lines correspond to the best fits of the form
$h(t)=[h_{0}^{1/n}+K(t-t_{0})]^{n}$ for each case. We obtain $n\approx 0.49$
for the non-rotating case (compared to the analytic prediction $n=0.5$) and
$n\approx 0.40$ for the rotating case (compared to the analytic prediction
$n=5/12\approx 0.41$). See text for more details. Panel (b):
$\mathrm{Ro_{sim}}$, $\mathrm{Ro_{sim}}^{1/2}$, and
$F_{K,\mathrm{R}}/F_{K,\mathrm{NR}}$, as a function of the thickness of the
convection zone, $h$. Note that $\mathrm{Ro_{sim}}^{1/2}$ has similar values
to $F_{K,\mathrm{R}}/F_{K,\mathrm{NR}}$, giving support to Equation (4).

We expect from the entrainment rates calculated in Section 3 that the
convection zone will grow more rapidly in the nonrotating case than the
rotating case. We quantify this difference by measuring the size of the
convection zone $h$ as a function of time. We define the bottom of the
convection zone as the location where the convective composition flux is
maximized (see Figure 2, middle), so that $h$ is the distance between that
location and the top of the box. Figure 3 (a) shows $h$ vs. time for the two
cases; as expected, $h$ in the non-rotating case grows more rapidly than in
the rotating case. We fit a power-law function of the form
$h(t)=[h_{0}^{1/n}+K(t-t_{0})]^{n}$, where $h_{0}$ is the size of the
convection zone at initial time $t_{0}$, and $K$ is a constant. In both
simulations, we demand that $h_{0}=1$, since initially the upper half of the
box is mixed. The time $t_{0}$ is formally the time when the convection zone
starts to grow. This only occurs after the convective turbulence is fully
established. The initial state in our simulations is one of quiescence. In the
non-rotating case, the convection continues to form and equilibrate up to
$t_{0}=225$ and only then does entrainment of the stable fluid begin. The
convection in the rotating case takes longer to be established, $t_{0}=325$.We
find $K\approx 1.61\times 10^{-4}\,(6.62\times 10^{-5})$ and $n\approx
0.49\,(0.4)$ in the non-rotating (rotating) case. These values are close to
what is predicted by the simplified model in Section 3, where the predicted
powerlaw index is $n=1/2=0.5$ for the non-rotating case, and $n=5/12\approx
0.41$ for the rotating case. In dimensionless form, the predicted values of
the prefactor are $K=4(R_{\rho}/\mathrm{Pe_{ff}})\sim 10^{-4}$ for the non-
rotating case, and
$K=(24/5)R_{\rho}(\mathrm{Ro_{ff}}/\mathrm{Pe_{ff}}^{2})^{3/5}\sim 2\times
10^{-5}$ for the rotating case, where the input dimensionless parameters
$R_{\rho}$, $\mathrm{Pe_{ff}}$, and $\mathrm{Ro_{ff}}$, are the density ratio,
and the Peclet and Rossby numbers based on the free-fall velocity,
respectively. The reader is refered to Appendix A for more details on the
nondimensionalization.

Finally, we measure the Rossby number of the flow, $\mathrm{Ro_{sim}}=\lVert
u\rVert/2\Omega h$, where $\lVert u\rVert$ and $h$ are the magnitude of the
flow velocity averaged within the convection zone and the thickness of the
convection zone, respectively. We find $\mathrm{Ro_{sim}}\approx 0.065$, which
decreases slightly as the convection zone evolves and deepens (see Figure 3b).
When compared with the ratio between the kinetic energy flux in the rotating
and non-rotating case (as predicted by mixing-length theory, see Equation 4),
we find $\mathrm{Ro_{sim}}^{1/2}/(F_{K,\mathrm{R}}/F_{K,\mathrm{NR}})\approx
0.87$. Finally, our theory (Equation 10) predicts that mixing a depth $\Delta
h$ should take longer in the rotating case by a factor of
$(5/6)\mathrm{Ro}^{-1/2}$. In our simulations this difference is approximately
a factor of 3, close to $(5/6)\mathrm{Ro_{sim}}^{-1/2}\approx 3.3$.

##  5 Summary and Discussion

In this work, we studied the effect of rotation on the compositional transport
across the boundary between a convection zone and a stable region (e.g., a
core-envelope interface in a gas giant). Mixing length theory predicts that
rotation reduces the convective velocity (compare Equations 1 and 3). For
Jovian conditions at $t\sim 10^{8}~{}\mathrm{yrs}$ since formation, this
corresponds to a velocity reduction of approximately a factor of 6.
Consequently, the kinetic energy flux available for the entrainment and mixing
of heavy material gets reduced by a factor of
$(U_{\mathrm{NR}}/U_{\mathrm{R}})^{3}\sim\mathrm{Ro}^{-1/2}\sim 6^{3}\sim 200$
(see Equation 4). Using the prescriptions for the convective velocity from
mixing-length theory, we derived a simple analytical model for the rate of
entrainment and the propagation of a convection zone. The thickness of the
outer convection zone follows $h(t)\propto t^{1/2}$ in the non-rotating case
(Equation 6), and $h(t)\propto\Omega^{1/4}t^{5/12}$ in the rotating case
(Equation 7), supporting the hypothesis that rotation decreases compositional
mixing. From this model, the timescale to mix a stable layer of vertical size
$L$ is larger by a factor of $(5/6)\mathrm{Ro}^{-1/2}$, where $\mathrm{Ro}$ is
the Rossby number, when rotation is taken into account in the convective
velocity (Equation 10). We ran two numerical experiments (including and not
including rotation) of convective entrainment and found a good agreement
between the analytic predictions and the simulations (e.g., see Figures 2–3).

The agreement between our numerical simulations and the mixing-length theory
predictions is encouraging. The mixing time can increase by two orders of
magnitude for Jupiter-like values at $t\sim 10^{8}$ years after formation.
This could significantly reduce the rate of core erosion, preventing strong
mixing without the need to invoke a double-diffusive staircase. For example,
Moll et al. (2017) derived that the rate at which core mass is transported
into the envelope is given by $\dot{M}_{\mathrm{core}}=\gamma^{-1}\alpha
L_{\mathrm{core}}/c_{P}$, where $L_{\mathrm{core}}$ is the luminosity at the
core-envelope interface, and $\gamma^{-1}$ is ratio of the buoyancy flux
associated with heavy-element transport to the buoyancy flux associated with
heat transport. Using 3D simulations of non-rotating incompressible flows,
Moll et al. (2017) constrained $\gamma^{-1}\gtrsim 0.5$ for Jupiter
conditions. In the entrainment model investigated in this work, the rate at
which mass is transported into the outer convection zone is
$\dot{M}_{\mathrm{core}}\propto hdh/dt$. Then, it follows that rotation
reduces the erosion rate by a factor of
$\dot{M}_{\mathrm{core,R}}/\dot{M}_{\mathrm{core,NR}}=\gamma^{-1}_{\mathrm{R}}/\gamma^{-1}_{\mathrm{NR}}\approx\left(N/2\Omega\right)^{1/2}\left(24Nt/5\right)^{-1/6}$
(where we have used Equations 6 and 7). Using values for Jupiter at $t\sim
10^{8}$ years after formation, we estimate
$\gamma^{-1}_{\mathrm{R}}/\gamma^{-1}_{\mathrm{NR}}\approx 0.01$. Although we
have derived this ratio with a rather simple model, the fraction is
significantly smaller than unity and illustrates the importance of rotation
for convective boundary mixing.

Although our numerical experiments are done in a Cartesian box and assume that
gravity and rotation are parallel to each other (as expected at the poles of
the planet), the reduction in the convective velocity does not depend on the
exact geometry or shape of the composition profile. However, numerical
simulations with misaligned gravity and rotation have shown that thermal
mixing is not uniform across latitude (e.g., Dietrich & Wicht, 2018; Currie et
al., 2020). A logical next step would be to investigate this problem using
spherical geometry and a primordial composition gradient that is more
consistent with formation models (e.g., Stevenson et al., 2022).

Our simulations do not include density stratification, as we adopt the
Bousinessq approximation which formally limits the vertical scale to be much
less than a density scale height. Including density stratification and
compressibility will change the mixing rate by a combination of two effects.
First, the potential energy required to mix a region should decrease; while
composition is uniformly mixed in the convection zone in both cases, the
compressible case mixes the density to the adiabatic gradient, whereas
Boussinesq convection mixes the density to a constant value, which requires a
greater energy expenditure. Second, stratification increases the asymmetry
between upflows and downflows and in turn increases the net kinetic energy
flux, so it should matter whether a stable entraining region is at the top or
bottom of a convection zone (see e.g., Fig. 4 Meakin & Arnett, 2007). At this
time we cannot speculate on the net result of those effects. To our knowledge,
there are no simulations of convective entrainment accounting for both
composition gradients and rotation in a compressible flow. We will investigate
this problem in a future paper.

We encourage future work on hydrodynamic models of gas giants, as improved
mixing prescriptions at convective boundaries provide a promising avenue for
better aligning 1D evolution models with observational inferences of gas giant
interiors from NASA’s missions Juno and Cassini.

We thank Nick Featherstone and Keith Julien for useful conversations on
rotating convection. A.C. acknowledges support by NSERC Discovery Grant
RGPIN-2023-03620. J.R.F. and B.W.H. are supported by NASA through grants
80NSSC18K1125, 80NSSC19K0267 and 80NSSC20K0193. E.H.A. is supported by a CIERA
Postdoctoral Fellowship. A.C. is a member of the Centre de Recherche en
Astrophysique du Québec (CRAQ) and the Trottier Institute for Research on
Exoplanets (iREx). Computations were conducted with support from the NASA High
End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS)
Division at Ames Research Center on Pleiades with allocation GID s2276.

## References

  * Anders et al. (2023) Anders, E., Fuentes, R., Cumming, A., & Hindman, B. W. 2023, evanhanders/rotation_reduces_entrainment: Code Release, v.1.0, Zenodo, doi:10.5281/zenodo.7950969
  * Anders et al. (2022) Anders, E. H., Jermyn, A. S., Lecoanet, D., et al. 2022, ApJ, 928, L10
  * Anders & Pedersen (2023) Anders, E. H., & Pedersen, M. G. 2023, Galaxies, 11, 56
  * Aurnou et al. (2020) Aurnou, J. M., Horn, S., & Julien, K. 2020, PhRvR, 2, 043115
  * Barker et al. (2014) Barker, A. J., Dempsey, A. M., & Lithwick, Y. 2014, ApJ, 791, 13
  * Bolton et al. (2017a) Bolton, S. J., Adriani, A., Adumitroaie, V., et al. 2017a, Sci, 356, 821
  * Bolton et al. (2017b) Bolton, S. J., Lunine, J., Stevenson, D., et al. 2017b, Space Sci. Rev., 213, 5
  * Burns et al. (2020) Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., & Brown, B. P. 2020, PhRvR, 2, 023068
  * Chabrier & Baraffe (2007) Chabrier, G., & Baraffe, I. 2007, ApJ, 661, L81
  * Cheng et al. (2020) Cheng, J. S., Madonia, M., Aguirre Guzmán, A. J., & Kunnen, R. P. J. 2020, PhRvF, 5, 113501
  * Cumming et al. (2018) Cumming, A., Helled, R., & Venturini, J. 2018, MNRAS, 477, 4817
  * Currie et al. (2020) Currie, L. K., Barker, A. J., Lithwick, Y., & Browning, M. K. 2020, MNRAS, 493, 5233
  * Debras & Chabrier (2019) Debras, F., & Chabrier, G. 2019, ApJ, 872, 100
  * Dietrich & Wicht (2018) Dietrich, W., & Wicht, J. 2018, FrEaS, 6, 189
  * Fernando (1987) Fernando, H. J. S. 1987, JFM, 182, 525
  * Fernando et al. (1991) Fernando, H. J. S., Chen, R.-R., & Boyer, D. L. 1991, JFM, 228, 513
  * French et al. (2012) French, M., Becker, A., Lorenzen, W., & et al. 2012, ApJS, 202, 5
  * Fuentes & Cumming (2020) Fuentes, J. R., & Cumming, A. 2020, PhRvF, 5, 124501
  * Fuentes et al. (2022) Fuentes, J. R., Cumming, A., & Anders, E. H. 2022, PhRvF, 7, 124501
  * Garaud (2018) Garaud, P. 2018, AnRFM, 50, 275
  * Guillot et al. (2004) Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, in Jupiter. The Planet, Satellites and Magnetosphere, ed. F. Bagenal, T. E. Dowling, & W. B. McKinnon, Vol. 1, 35–57
  * Helled & Stevenson (2017) Helled, R., & Stevenson, D. 2017, ApJ, 840, L4
  * Helled et al. (2022) Helled, R., Stevenson, D. J., Lunine, J. I., et al. 2022, Icarus, 378, 114937
  * Hindman et al. (2020) Hindman, B. W., Featherstone, N. A., & Julien, K. 2020, ApJ, 898, 120
  * Howard et al. (2023) Howard, S., Guillot, T., Bazot, M., et al. 2023, A&A, 672, A33
  * Julien et al. (1996a) Julien, K., Legg, S., McWilliams, J., & Werne, J. 1996a, DyAtO, 24, 237
  * Julien et al. (1996b) —. 1996b, JFM, 322, 243
  * Leconte & Chabrier (2012) Leconte, J., & Chabrier, G. 2012, A&A, 540, A20
  * Leconte & Chabrier (2013) —. 2013, NatGe, 6, 347
  * Linden (1975) Linden, P. F. 1975, JFM, 71, 385
  * Liu et al. (2019) Liu, S.-F., Hori, Y., Müller, S., et al. 2019, Nature, 572, 355
  * Mankovich & Fuller (2021) Mankovich, C. R., & Fuller, J. 2021, Nature Astronomy, 5, 1103
  * Meakin & Arnett (2007) Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448
  * Miguel & Vazan (2023) Miguel, Y., & Vazan, A. 2023, Remote Sensing, 15, 681
  * Molemaker & Dijkstra (1997) Molemaker, M. J., & Dijkstra, H. A. 1997, JFM, 331, 199–229
  * Moll et al. (2017) Moll, R., Garaud, P., Mankovich, C., & et al. 2017, ApJ, 849, 24
  * Müller et al. (2020) Müller, S., Helled, R., & Cumming, A. 2020, A&A, 638, A121
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  * Proudman (1916) Proudman, J. 1916, RSPSA, 92, 408
  * Spiegel & Veronis (1960) Spiegel, E. A., & Veronis, G. 1960, ApJ, 131, 442
  * Stevenson (1979) Stevenson, D. J. 1979, GApFD, 12, 139
  * Stevenson et al. (2022) Stevenson, D. J., Bodenheimer, P., Lissauer, J. J., & D’Angelo, G. 2022, PSJ, 3, 74
  * Taylor (1917) Taylor, G. I. 1917, RSPSA, 93, 99
  * Turner (1968) Turner, J. S. 1968, JFM, 33, 183
  * Vasil et al. (2021) Vasil, G. M., Julien, K., & Featherstone, N. A. 2021, PNAS, 118, e2022518118
  * Vazan et al. (2018) Vazan, A., Helled, R., & Guillot, T. 2018, A&A, 610, L14
  * Vazan et al. (2015) Vazan, A., Helled, R., Kovetz, A., & Podolak, M. 2015, ApJ, 803, 32
  * Wahl et al. (2017) Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, GeoRL, 44, 4649
  * Wang & Ruuth (2008) Wang, D., & Ruuth, S. J. 2008, JCM, 26, 838
  * Wood et al. (2013) Wood, T. S., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 157

##  Appendix A Fluid equations & Initial Conditions

We study convective entrainment in rapidly rotating flows with composition and
temperature gradients. To do so, we adopt the Boussinesq approximation of the
fluid equations and evolve both temperature $T$ and composition $C$. We
normalize the fluid equations using as units of length and time the initial
size of the convection zone $h_{0}$ and the free fall time across it
$t_{\mathrm{ff}}=h_{0}/u_{\mathrm{ff}}$ (where $u_{\mathrm{ff}}=\sqrt{\alpha
gh_{0}\Delta T_{F}}$, $g$ is the acceleration due to gravity, and
$\alpha\equiv-(\partial\ln\rho/\partial T)|_{C}$ is the coefficient of thermal
expansion). Further, for the characteristic temperature scale we utilize
$\Delta T_{F}=F_{0}h_{0}/\rho_{0}\kappa_{T}c_{P}$, with $\rho_{0}$ being the
background density, $F_{0}$ the cooling energy flux at the surface of the
convection zone, $\kappa_{T}$ the thermal diffusivity, and $c_{P}$ the
specific heat at constant pressure. The units of temperature and composition
are respectively $[T]=\Delta T_{F}$, and $[C]=\Delta C_{0}$ (i.e., the initial
composition change across the stable region). The corresponding normalizations
for the kinetic, heat and composition fluxes are
$[F_{K}]=(1/2)\rho_{0}u_{\mathrm{ff}}^{3}$, $[F_{H}]=\rho
c_{P}u_{\mathrm{ff}}\Delta T_{F}$, and $[F_{X}]=\rho_{0}u_{\mathrm{ff}}\Delta
C_{0}$, respectively.

The nondimensional equations are

$\displaystyle\bm{\nabla}\cdot\bm{u}=0$ (A1)
$\displaystyle\partial_{t}\bm{u}+\bm{u}\cdot\bm{\nabla}\bm{u}+\bm{\nabla}\varpi+\mathrm{Ro_{ff}}^{-1}\hat{\bm{z}}\times\bm{u}=\left(T-\frac{C}{\mathrm{R}_{\rho}}\right)\hat{\bm{z}}+\frac{\mathrm{Pr}}{\rm{Pe_{ff}}}\bm{\nabla}^{2}\bm{u},$
(A2)
$\displaystyle\partial_{t}T+\bm{u}\cdot\bm{\nabla}T=\frac{1}{\rm{Pe_{ff}}}\bm{\nabla}^{2}T,$
(A3)
$\displaystyle\partial_{t}C+\bm{u}\cdot\bm{\nabla}C=\frac{\tau_{0}}{\rm{Pe_{ff}}}\bm{\nabla}^{2}\overline{C}+\frac{\tau}{\rm{Pe_{ff}}}\bm{\nabla}^{2}C^{\prime},$
(A4)

where $\bm{u}=(u,v,w)$ is the velocity field, and $\varpi$ is the reduced
pressure. Overbars denote horizontal averages and primes denote fluctuations
around that average, such that e.g., $C=\overline{C}+C^{\prime}$. Note that we
reduce diffusion on $\overline{C}$ by choosing $\tau_{0}<\tau$ to ensure its
evolution is due to advection.

In the equations above, there are 5 dimensionless numbers that characterize
the evolution of the flow. These are the freefall Rossby number, Peclet
number, density ratio, Prandtl number, and inverse Lewis number, which are
defined respectively as

$\displaystyle\begin{split}&\mathrm{Ro_{ff}}=\frac{u_{\rm{ff}}}{2\Omega
h_{0}},\hskip
2.84544pt\mathrm{Pe_{ff}}=\frac{u_{\rm{ff}}h_{0}}{\kappa_{T}},\hskip
2.84544ptR_{\rho}=\frac{\alpha\Delta T_{F}}{\beta\Delta C_{0}},\hskip
2.84544pt\mathrm{Pr}=\frac{\nu}{\kappa_{T}},\hskip
2.84544pt\tau=\frac{\kappa_{C}}{\kappa_{T}},\end{split}$ (A5)

where $\Omega$ is the rotation rate, $\alpha\equiv(\partial\ln\rho/\partial
T)|_{C}$ and $\beta\equiv(\partial\ln\rho/\partial C)|_{T}$ are respectively
the coefficient of thermal and compositional expansion/contraction, $\nu$ is
the kinematic viscosity, and $\kappa_{T}$ and $\kappa_{C}$ are the thermal and
compositional diffusivity, respectively. Note that $R_{\rho}$ is the ratio
between the destabilizing and stabilizing effect of the thermal and
compositional buoyancy.

The initial structure of the fluid is composed of two layers, defined across
$z\in[0,2]$, where the temperature and composition are respectively

$\displaystyle T(z)=2-z\,,\qquad C_{0}(z)=\begin{cases}0.5&z>1\\\ 2-z&z<1\\\
\end{cases}.$ (A6)

Step functions are not well represented in pseudospectral codes, so we use a
smoothed Heaviside function centered at $z=z_{0}$,

$H(z;z_{0},d_{w})=\frac{1}{2}\left(1+\mathrm{erf}\left[\frac{z-z_{0}}{d_{w}}\right]\right).$
(A7)

where erf is the error function; to construct $C_{0}(z)$ in Equation. A6, we
set

$C_{0}(z)=(2-z)[1-H(z;1,0.05)]+\frac{1}{2}H(z;1,0.05).$ (A8)

We fix the flux of temperature and composition at the top of the domain and we
fix their values at the bottom of the domain, so the boundary conditions for
temperature and composition are $\partial_{z}T(z=2)=-1$, $T(z=0)=2$,
$\partial_{z}C(z=2)=0$, and $C(z=0)=2$. The boundary conditions for the
velocity are stress free at the upper boundary
($\hat{z}\cdot\bm{u}=\hat{x}\cdot\partial_{z}\bm{u}=\hat{y}\cdot\partial_{z}\bm{u}=0$
at $z=2$) and no-slip at the lower boundary in the stable fluid
($\bm{u}=\bm{0}$ at $z=0$).

The simulations in this work use $\mathrm{Pe_{ff}}=4.4\times 10^{3}$,
$R_{\rho}^{-1}=10$, $\mathrm{Pr}=\tau=0.5$, and $\tau_{0}=10^{-3}$. We set the
Coriolis force to be identically zero in the non-rotating case
($\mathrm{Ro_{ff}}=\infty$) , while in the rotating case
$\mathrm{Ro_{ff}}=0.81$. Under Jovian conditions,
$\mathrm{Ro_{ff}}=\mathcal{O}(10^{-5})$,
$\mathrm{Pe_{ff}}=\mathcal{O}(10^{10})$
$\tau\approx\mathrm{Pr}=\mathcal{O}(10^{-2})$, and
$R^{-1}_{\rho}=\mathcal{O}(10)$ (see e.g., Guillot et al., 2004; French et
al., 2012), but we make clear that the density ratio is very uncertain given
our ignorance on the compositional and thermal stratification in the deep
interiors of the gas giants.) Our simulations are as turbulent as possible and
are qualitatively in the same regime as gas giants ($\mathrm{Ro_{ff}}<1$,
$\mathrm{Pe_{ff}}\gg 1$, $\mathrm{Pr}<1$, $\tau<1$, $R_{\rho}^{-1}>1$). It is
worth mentioning that although the input Rossby number $\mathrm{Ro_{ff}}$ is
close to unity, the Rossby number computed from the simulated flows is much
smaller ($\approx 0.065$). It would be possible to decrease the Rossby number,
but the time that would be required for the convection zone to mix material
would be significantly longer and beyond our computational resources.

##  Appendix B Simulation Details

We time-evolve equations A1-A4 using the Dedalus pseudospectral solver (Burns
et al., 2020) version 3222We use the code in the master branch of the
https://github.com/DedalusProject/dedalus repository, and use the commit with
short-sha 29f3a59. using timestepper SBDF2 (Wang & Ruuth, 2008) and CFL safety
factor 0.2. All variables are represented using a Chebyshev series with 512
terms for $z\in[0,L_{z}]$, with $L_{z}=2$, and Fourier series in the periodic
$x$ and $y$ directions. In the non-rotating case, $x\in[0,L_{x}]$,
$y\in[0,L_{y}]$, with $L_{x}=L_{y}=3$, and each direction has 192 terms. In
the rotating case, $L_{x}=L_{y}=2$, and each direction has 128 terms. To avoid
aliasing errors, we use the 3/2-dealiasing rule in all directions. To start
our simulations, we add random noise temperature perturbations sampled from a
normal distribution with a magnitude of $10^{-5}$ to the initial temperature
field. The Python scripts used to run the simulations in this work are
available online as a Git
repository333https://github.com/evanhanders/rotation_reduces_entrainment,
v1.0., and backed up in a Zenodo repository (Anders et al., 2023).