# Minimally tough chordal graphs with toughness at most $1/2$

Humara Khan  Gyula Y. Katona  Department of Computer Science and Information
Theory
Budapest University of technology and Economics
Budapest, Hungary  Humara Khan  Department of Computer Science and Information
Theory
Budapest University of technology and Economics
Budapest, Hungary

(June 2022)

###### Abstract

Let $t$ be a positive real number. A graph is called _$t$ -tough_ if the
removal of any vertex set $S$ that disconnects the graph leaves at most
$|S|/t$ components. The toughness of a graph is the largest $t$ for which the
graph is $t$-tough. A graph is minimally $t$-tough if the toughness of the
graph is $t$ and the deletion of any edge from the graph decreases the
toughness. A graph is _chordal_ if it does not contain an induced cycle of
length at least $4$. We characterize the minimally $t$-tough, chordal graphs
for all $t\leq 1/2$. As a corollary, a characterization of minimally
$t$-tough, interval graphs is obtained for $t\leq 1/2$.

##  1 Introduction

All graphs considered in this paper are finite, simple and undirected. Let
$\omega(G)$ denote the number of components and $\alpha(G)$ denote the
independence number and $\kappa(G)$ denote the connectivity number of the
graph $G$. (Using $\omega(G)$ to denote the number of components may be
confusing, however, most of the literature on toughness uses this notation.)
For a connected graph $G$, a vertex set $S\subseteq V(G)$ is called a cutset
if $\omega(G-S)>1$.

The notion of toughness was introduced by Chvátal in [1].

######  Definition 1.

Let $t$ be a real number. A graph $G$ is called $t$-tough if $|S|\geq
t\cdot\omega(G-S)$ for any cutset $S\subseteq V(G)$. The toughness of $G$,
denoted by $\tau(G)$, is the largest $t$ for which G is $t$-tough, taking
$\tau(K_{n})=\infty$ for all $n\geq 1$. We say that a cutset $S\subseteq V(G)$
is a tough set if $\omega(G-S)=|S|/\tau(G)$.

Note that a graph is disconnected if and only if its toughness is $0$.

######  Definition 2.

A graph $G$ is said to be minimally $t$-tough if $\tau(G)=t$ and $\tau(G-e)<t$
for all $e\in E(G)$.

It follows directly from the definition that every $t$-tough noncomplete graph
is $2t$-connected, implying $\kappa(G)\geq 2\tau(G)$ for noncomplete graphs.
Therefore, the minimum degree of any $t$-tough noncomplete graph is at least
$\lceil 2t\rceil$.

The following conjecture is motivated by a theorem of Mader [2], which states
that every minimally $k$-connected graph has a vertex of degree $k$.

######  Conjecture 3 (Kriesell [3]).

Every minimally $1$-tough graph has a vertex of degree $2$.

This conjecture can be naturally generalized.

######  Conjecture 4 (Generalized Kriesell’s Conjecture).

Every minimally $t$-tough graph has a vertex of degree $\lceil 2t\rceil$.

######  Definition 5.

A graph is chordal if it does not contain an induced cycle of length at least
$4$.

A vertex $v$ of a graph $G$ is simplicial if its closed neighborhood $N[v]$
forms a clique in $G$.

It was already known that the conjecture is true for chordal graphs when
$t\leq 1$.

######  Theorem 6 ([4]).

Let $t\leq 1$ be a positive rational number. If $t\leq\frac{1}{2}$, then every
simplicial vertex of any minimally $t$-tough, chordal graph has degree $1$. If
$\frac{1}{2}<t\leq 1$, then there exist no minimally $t$-tough, chordal
graphs.

Our goal is to characterize minimally tough, chordal graphs with toughness at
most $1/2$. This also gives a new proof of the Generalized Kriesell’s
Conjecture for these graphs. In the proof, we will need the notion of clique
trees, which we introduce now.

Let ${\cal K}_{G}$ be the set of maximal cliques of a graph $G$. A clique tree
of $G$ is a tree $T$ with vertex set ${\cal K}_{G}$ such that it satisfies the
following clique-intersection property: For every pair of distinct cliques
$K,K^{\prime}\in{\cal K}_{G}$, the set $K\cap K^{\prime}$ is contained in
every clique on the path connecting $K$ and $K^{\prime}$ in the tree [7, 8].
We can also assign weights to the edges of a clique tree. If
$K,K^{\prime}\in{\cal K}_{G}$ are adjacent in the clique tree, then let
$|K\cap K^{\prime}|$ be the weight of this edge. One can show that each weight
is at least 1 (see [7]).

It is also proved in [7], that clique trees also satisfy the induced-subtree
property: For every vertex $v\in V(G)$ the set of cliques containing $v$
induces a subtree of the clique tree $T$.

######  Theorem 7 ([7]).

  * •

A connected graph $G$ is chordal if and only if there exists a tree $T=({\cal
K}_{G},{\cal E}_{T})$ for which the clique intersection-property holds.

  * •

A connected graph $G$ is chordal if and only if there exists a tree $T=({\cal
K}_{G},{\cal E}_{T})$ for which the induced-subtree property holds.

In other words, a connected graph $G$ is chordal if and only if it has a
clique tree.

######  Definition 8.

We call a graph a _TT-graph_ if it can be obtained from a tree of maximum
degree $\Delta\geq 3$ in the following way. Let $Y$ be a subset of vertices
satisfying one of the following conditions:

  1. (a)

if $\Delta=3$, then $Y$ is the set of all degree 3 vertices so that $Y$ is an
independent set and every neighbor of a vertex in $Y$ has degree $2$, or

  2. (b)

if $\Delta\geq 3$, then $Y$ is a subset of some (or all, or none) of the
degree 3 vertices so that $Y$ is an independent set and every neighbor of a
vertex in $Y$ has degree $\Delta$.

Now remove all vertices of $Y$ (in one step) and for each removed vertex join
their three neighbors with a triangle.

Figure 1: Creating a TT-graph from a tree by Definition 8

Notice that a TT-graph $G$ may contain triangles, but there is no larger
clique in it. Also, any two triangles may only have one common vertex. Trees
are also TT-graphs. It is also obvious that TT-graphs are chordal graphs. For
TT-graphs we define the _modified degree_ for each vertex: $md(v)$ is the
number of components in $G-v$. This is clearly the same as $d_{G}(v)$ minus
the number of triangles containing $v$. If a TT-graph $G$ is constructed from
a tree $T$, then $md_{G}(v)=d_{T}(v)$. Furthermore, it follows from the
construction method that all vertices contained in a triangle must have
maximal modified degree. Let $\mu(G)$ denote the maximum modified degree in
$G$.

During the proof we will use the following general theorem to characterize
those graphs that are not minimally $t$-tough.

######  Theorem 9 ([6]).

Let $G$ be a connected graph that is not complete and let $t=\tau(G)$. Then
$G$ is not minimally $t$-tough if and only if $G$ contains an edge $e=uv$ such
that the following conditions are met.

  1. (a)

There exist at least $2t+1$ internally vertex-disjoint $u$-$v$ path in $G$
(including $uv$).

  2. (b)

Every cutset $S$ in $G$ that is also a $u$-$v$ cutset in $G-e$ satisfies

$|S|\geq(\omega(G-S)+1)t.$

Our main result is a characterization of chordal graphs with toughness $t\leq
1/2$.

######  Theorem 10.

For any positive rational number $t\leq 1/2$, a chordal graph is minimally
$t$-tough if and only if it is a TT-graph with toughness $t=1/\mu(G)$.

The proof is divided to several lemmas. Each lemma will cover a subset of the
chordal graphs depending on some properties of their clique trees.

##  2 Proof

In the rest of the paper we assume that $t\leq 1/2$ is a positive rational
number.

######  Lemma 11.

If $G$ is a chordal graph such that in every clique tree, all the edges have
weight $1$, then every vertex of $G$ is either a simplicial vertex or a cut-
vertex.

###### Proof.

Notice that the clique intersection property implies that if there are two
maximal cliques with intersection of size at least $2$, then the clique tree
will have an edge of weight at least $2$ even if the corresponding vertices
are not adjacent in the clique tree. This would contradict the assumption of
the lemma.

Clearly, every vertex of the graph must be contained in a maximal clique. If a
vertex $v$ is contained in only one, then it is a simplicial vertex. So assume
that $v$ is contained in at least two maximal cliques, say $Q_{1}$ and
$Q_{2}$. The assumption implies that $|Q_{1}\cap Q_{2}|=1$. Now we show that
$v$ is a cut-vertex.

Suppose to the contrary that it is not, then for any vertex $x\in Q_{1}-Q_{2}$
and any vertex $y\in Q_{2}-Q_{1}$ there is a path $P$ connecting $x$ and $y$
and avoiding $v$, choose a shortest such path. If $xy\in E(G)$, then the
intersection of the maximal clique containing the triangle $vxy$ and $Q_{1}$
has size at least $2$, a contradiction. If $xy\notin E(G)$, then the path $P$
with the path $xvy$ forms a cycle of length at least $4$. Since $G$ is
chordal, this cycle must be triangulated. Moreover, since we chose a shortest
path, $v$ must be connected to all the other vertices on the path $P$. Now if
$x^{\prime}y^{\prime}$ is the first edge of the path leaving $Q_{1}$, then the
intersection of the maximal clique containing the triangle
$vx^{\prime}y^{\prime}$ and $Q_{1}$ has size at least $2$, a contradiction. ∎

######  Lemma 12.

If $G$ is a chordal graph with $\tau(G)=t\leq 1/2$ which has a clique tree
containing an edge of weight at least $2$, then $G$ is not minimally
$t$-tough.

###### Proof.

Since the toughness of $G$ is $t$ and $t\leq 1/2$, it is clear that the graph
is not complete, thus we can use Theorem 9 to prove that the graph is not
minimally $t$-tough. By the assumption, a clique tree of $G$ contains an edge
that has weight at least $2$. This means that $G$ contains two cliques,
$Q_{1}$ and $Q_{2}$, that share at least 2 vertices. If there is more than one
edge in the clique tree having weight at least $2$, then we have to choose a
suitable one. Clearly, removing any edge of the clique tree will separate the
clique tree into two components. Pick an edge with weight at least $2$ so that
one of the components contains only edges of weight 1, and let $Q_{1}$ be the
clique that corresponds to the vertex in this component. Notice that this
implies that if another clique contains a vertex of $Q_{1}-Q_{2}$, then it
cannot contain any other vertex of $Q_{1}$.

Now let $u$ and $v$ be two vertices in $Q_{1}\cap Q_{2}$. Since $u$ and $v$
belong to $Q_{1}$, they are connected by an edge, denote it by $e$. It will be
shown that the conditions of Theorem 9 are satisfied for $e=uv$.

All vertices of the two cliques $Q_{1}$ and $Q_{2}$ are adjacent to both $u$
and $v$. Since these cliques are different, there must be a vertex in
$Q_{1}-Q_{2}$ which is a common neighbor of $u$ and $v$. This gives a path
between $u$ and $v$, so together with the $uv$ edge there are two paths
between them. Since $2t+1\leq 2$ follows from our assumption, condition a of
Theorem 9 is satisfied. Notice that there must also be a vertex also in
$Q_{2}-Q_{1}$, this will be used later.

To show that condition b holds as well, consider an arbitrary cutset $S$ in
$G$ that is also a $u$-$v$ cutset in $G-e$. Now that all common neighbors of
$u$ and $v$ must belong to $S$. Let $x\in Q_{1}-Q_{2}\subseteq N(u)\cup N(v)$
be an arbitrary element. Thus $x\in S$ holds.

We claim that $x$ cannot have a neighbor $y$ outside of $Q_{1}$ that is
connected to $u$ by a path avoiding $x$. Suppose that $yy_{1}y_{2}\dots y_{m}$
is such a path with $y_{m}=u$. Let $k$ be the smallest index for which
$y_{k}\in Q_{1}$; such an index exists since $y_{m}=u\in Q_{1}$. Since $x\in
Q_{1}$ holds, $yy_{1}y_{2}\ldots y_{k}x$ is a cycle, which must be
triangulated because $G$ is chordal. One of these triangles must contain the
edge $xy_{k}$, let the third vertex be $y_{k^{\prime}}$ with $k^{\prime}<k$.
This vertex is not in $Q_{1}$ by the choice of $k$. There is a maximal clique
containing the triangle $xy_{k}y_{k^{\prime}}$, however, this contradicts the
choice of $Q_{1}$.

The above claim implies that either all neighbors of $x$ belong to $Q_{1}$ or
$x$ is a cut-vertex (that separates $u$ from some neighbors of $x$). In the
first case $\omega(G-x)t\leq\omega(G-x)=1$ holds since $t\leq 1/2$. In the
second case, $\omega(G-x)t\leq 1$ holds since $G$ is $t$-tough. So in both
cases we have $\omega(G-x)\leq 1/t$.

Let $S^{\prime}=S-\\{x\\}$.

Case 1: If $S^{\prime}$ is a cutset in $G$, then
$|S^{\prime}|\geq\omega(G-S^{\prime})t$ must hold since $G$ is $t$-tough.
Consider the components of $G-S^{\prime}$. The vertex $x$ is in one of these
components, so when we delete $x$ as well, then this component will fall apart
into at most $\omega(G-x)$ components. This follows from the observation that
by the above claims the neighborhood of $x$ spans disjoint cliques, so it is
impossible to have a path between two such cliques that avoids $x$. Therefore
we have

$\omega(G-S)\leq\omega(G-S^{\prime})+\omega(G-x)-1\leq\omega(G-S^{\prime})+\frac{1}{t}-1,$
$\omega(G-S)t+t\leq\omega(G-S^{\prime})t+1.$

This implies that

$|S|=|S^{\prime}|+1\geq\omega(G-S^{\prime})t+1\geq\omega(G-S)t+t=(\omega(G-S)+1)t.$

Thus condition b of Theorem 9 is satisfied in this case.

Case 2: If $S-{x}$ is not a cutset, then we prove that
$\omega(G-S)\leq\omega(G-x)$. It is easy to see that it is implied by the
following claim.

######  Claim 13.

If $S-x$ is not a cutset and two vertices of $G-S$ are in the same component
of $G-x$, then they are also in the same component of $G-S$.

###### Proof.

Let $C_{0},C_{1},\ldots$ denote the components of $G-x$. Since $u\notin S$,
one of the components must contain $u$, let $C_{0}$ be this component. We have
shown before that the neighbors of $x$ in $C_{0}$ must belong to $Q_{1}$. If
there are no other neighbors of $x$ then the only component is $C_{0}$,
otherwise there are at least two components. Now let $p$ and $q$ be two
vertices of a component of $G-x$ and suppose to the contrary that they belong
to different components in $G-S$.

Subcase (i): $p,q\in C_{0}$.

Since $S-x$ is not a cutset, this means that there exists a $p$-$q$ path
avoiding $S-x$, thus it must contain $x$. Let
$p,\ldots,x^{-},x,x^{+},\ldots,q$ be such a path. Clearly all vertices of this
path must be in $C_{0}$ because of the properties of $x$. This implies that
$x^{-},x^{+}\in Q_{1}$. Thus $x^{-},x^{+}$ are connected with an edge, so
$p,\ldots,x^{-},x^{+},\ldots,q$ is a $p$-$q$ path that does not contain $x$
and avoids $S-x$, so it avoids $S$, a contradiction.

Subcase (ii): $p,q\in C_{i}$ for some $i>0$.

The choice of $Q_{1}$ implies that $C_{i}$ is a chordal graph that satisfies
the conditions of Lemma 11. If $p$ and $q$ belong to the same clique, then
they belong to the same component of both $G-S$, a contradiction. If they
belong to different cliques, then any path connecting them must contain a cut-
vertex which belongs to $S-\\{x\\}$. This implies that $S-\\{x\\}$ is a
cutset, a contradiction. ∎

Now we complete the proof of the lemma. We noticed before that $S$ must
contain all common neighbors of $u$ and $v$, and that there exists a common
neighbor other than $x$. Thus $S-x$ is not empty. Hence, using that
$\omega(G-x)t\leq 1$ and $\omega(G-S)\leq\omega(G-x)$, we have

$|S|=1+|S-\\{x\\}|\geq\omega(G-x)t+1\geq\omega(G-S)t+t=(\omega(G-S)+1)t.$

So condition b of Theorem 9 is satisfied in this case, too. ∎

######  Lemma 14.

Let $G$ be a chordal graph so that in every clique tree all the edges have
weight $1$ and there exist a clique of size at least $4$ in $G$. This implies
that $G$ is not minimally $t$-tough.

###### Proof.

Let $e=uv$ be an arbitrary edge of a clique $Q$ of size at least $4$.

Condition a of Theorem 9 is clearly satisfied since $2t+1\leq 2$ and there are
two paths of length 2 between $u$ and $v$ through two other vertices of $Q$
(and there is a third path, the $uv$ edge). To show that condition b holds as
well, consider an arbitrary cutset $S$ in $G$ that is also a $u$-$v$ cutset in
$G-e$. Clearly, $Q-\\{u,v\\}\subseteq S$.

Let $x$ be an arbitrary vertex in $Q-\\{u,v\\}$ and $S^{\prime}=S-{x}$. Now
$S^{\prime}$ is not empty since $|Q|\geq 4$. Now we can repeat the same
argument as in the proof of Lemma 12, which implies that Condition a of
Theorem 9 holds in this case, too. ∎

In the introduction, the modified degree was defined for TT-graphs. However,
it can be easily extended for any graph in which each vertex is either
simplicial vertex or a cut-vertex:

######  Definition 15.

If $G$ is graph in which each vertex is either simplicial vertex or a cut-
vertex, then the modified degree of a vertex $v$, denoted by $md(v)$, is the
number of components in $G-v$. Let $\mu(G)$ denote the maximum modified degree
in $G$.

For these graphs it is easy to determine the toughness.

######  Lemma 16.

If $G$ is a connected, noncomplete graph in which each vertex is either a
simplicial vertex or a cut-vertex, then $\tau(G)=1/\mu(G)$.

###### Proof.

If $md(v)=\mu(G)$, then clearly the number of components of $G-v$ is $\mu(G)$,
so $\tau(G)\leq 1/\mu(G)$.

Let $S$ be a tough set. We claim that $|S|=1$. If $|S|>1$, then let $x\in S$
be a vertex having the smallest modified degree, and let $S^{\prime}=S-x$. If
$x$ is a simplicial vertex, then $\omega(G-S^{\prime})=\omega(G-S)$, which
contradicts the assumption that $S$ is a tough set. Thus $x$ is a cut-vertex.
Consider the components of $G-S^{\prime}$. The vertex $x$ belongs to one of
these components. Since $x$ is a cut-vertex, we have
$\omega(G-S)\leq\omega(G-S^{\prime})+\mu(G)-1$, because deleting $x$ besides
$S^{\prime}$ would give new $\mu(G)$ components instead of this old component.
Now

$\tau(G)\omega(G-S^{\prime})\geq\tau(G)\omega(G-S)-\tau(G)(\mu(G)-1)=\\\
|S|-\mu(G)\tau(G)+\tau(G)>|S|-1=|S^{\prime}|,$

since $\mu(G)\tau(G)<1$ by our first claim, and $\tau(G)>0$. On the other
hand, since $\tau(G)$ is the toughness, we have
$\omega(G-S^{\prime})\tau(G)\leq|S^{\prime}|$, a contradiction.

Hence, it is enough to consider the cutsets of size one to determine the
toughness, and it is clear that the most number of components is obtained by
deleting a vertex with maximum modified degree. ∎

######  Lemma 17.

Let $G$ be a connected, noncomplete chordal graph so that in every clique tree
all the edges have weight $1$ and all cliques have size at most $3$ in $G$. If
$G$ is minimally $t$-tough, then it is a TT-graph and $t=1/\mu(G)$.

###### Proof.

In Claim 11 we showed that each vertex of $G$ is either a simplicial vertex or
a cut-vertex, so we can use Definition 15.

We show that if $G$ is minimally $t$-tough and $v$ is contained in a triangle
$Q$, then $md(v)=\mu(G)$. Suppose on the contrary that $md(v)\leq\mu(G)-1$.
Let $e=uw$ be the edge of $Q$ that is not incident to $v$. Consider now
$G^{\prime}=G-e$. If $u$ or $w$ was a cut-vertex in $G$, then the same will
hold in $G^{\prime}$. If either one of them was a simplicial vertex in $G$,
then it means that it had degree 2 in $G$. Thus in $G^{\prime}$ it will have
degree 1, implying that it is still a simplicial vertex. So all vertices of
$G^{\prime}$ are still either simplicial or cut-vertices since for other
vertices nothing will change.

Notice that $md_{G^{\prime}}(v)=md_{G}(v)+1\leq\mu(G)$. Also,
$\mu(G^{\prime})=\mu(G)$, since the modified degree does not change for any
other vertex besides $v$. Thus $\tau(G^{\prime})=\tau(G)$ holds by Lemma 16,
implying a contradiction.

When $\mu(G)\geq 3$, then if we replace each triangle with a claw $K_{1,3}$ so
that the leaves of the claw correspond to the vertices of the triangle and the
center is a news vertex, then we clearly obtain a tree with maximum degree
$\mu(G)$. Also, the leaves of each claw have maximum degree in the tree. So
reversing the replacements show that $G$ can be obtained from this tree as
described in case (b) of the definition of TT-graphs.

In the special case when $\mu(G)=2$, replacing the triangles with a claw will
result in a tree of maximum degree $3$, but all the neighbors of the degree
$3$ vertices will have degree $2$. So reversing the replacements show that $G$
can be obtained from this tree as described in case (a) of the definition of
TT-graphs. Therefore we proved that $G$ is a TT-graph with toughness
$1/\mu(G)$. ∎

###### Proof of Theorem 10.

Let $G$ be a minimally $t$-tough, chordal graph with $t\leq 1/2$. Since $G$ is
chordal and connected, it has a set of clique trees. By Lemma 12, none of
these clique trees contain an edge of weight at least $2$, and by Lemma 14,
there is no clique of size at least $4$ in $G$. Thus the conditions of Lemma
17 are satisfied, so applying the lemma, we conclude that $G$ is a TT-graph.

To complete the proof, we show that any TT-graph $G$ is minimally $t$-tough.
It is clear that in a TT-graph every vertex is either a leaf, and therefore
simplicial, or a cut-vertex. Thus Lemma 16 can be applied, implying that
$\tau(G)=1/\mu(G)$. It is also clear that TT-graphs have two kinds of edges:
every edge is either a bridge or it is contained in a triangle.

If an edge $e$ is a bridge, then $G-e$ is disconnected, thus its toughness is
$0$. So the toughness of the graph clearly decreases by removing such an edge.

Now let $e$ be an edge of a triangle $Q$. By the definition of TT-graphs, we
know that all vertices of $Q$ have maximum modified degree $\mu(G)$. This
implies that in $G-e$, the vertex of $Q$ not incident to $e$ has modified
degree $\mu(G)+1$. Since every vertex of $G-e$ is still either simplicial or a
cut-vertex, Lemma 16 shows that the toughness of $G-e$ is
$1/(\mu(G)+1)<1/\mu(G)$. So the toughness of the graph decreases by removing
this kind of edges, too. ∎

##  3 Implication for interval graphs

######  Definition 18.

A graph is called an interval graph if it is obtained from intervals on the
real line in such a way that each interval represents a vertex and the
intersection of any two intervals represent an edge. Such graphs are also
known as the intersection graph of intervals.

An independent triple of vertices $x,y,z$ is called an _asteroidal triple_
(AT, for short) if between any pair in the triple there exists a path that
avoids the neighborhood of the third vertex. _Caterpillars_ are those trees
for which removing the leaves produces a path.

######  Theorem 19 ([5]).

A graph is an interval graph if and only if it is chordal and asteroidal
triple-free.

######  Corollary 20.

If $G$ is a minimally $t$-tough, interval graph with $t\leq 1/2$, then $G$ is
a caterpillar graph.

###### Proof.

Since every interval graph is chordal, by Theorem 10, $G$ must be a TT-graph.
Now we show that $G$ is a tree. Suppose to the contrary that $G$ is not a
tree, i.e. it contains a triangle. Each vertex in a triangle must have a
distinct neighbor outside the triangle. These neighbors form an AT in the
graph, a contradiction. So $G$ is a tree. However, any tree that is not a
caterpillar contains the AT subgraph which is obtained from the claw by
subdividing each edge with a vertex. Thus $G$ must be a caterpillar. ∎

## Acknowledgment

The authors are grateful to Kitti Varga for her invaluable help in improving
the manuscript. This research was supported by the Ministry of Innovation and
Technology and the National Research, Development and Innovation Office within
the Artificial Intelligence National Laboratory of Hungary.

## References

  * [1] V. Chvátal, Tough graphs and hamiltonian circuits, Discrete Mathematics, 5, pp. 215–228 (1973)
  * [2] W. Mader, Eine Eigenschaft der Atome endlicher Graphen, Archiv der Mathematik, 22, pp. 333–336 (1971)
  * [3] T. Kaiser, Problems from the workshop on dominating cycles, http://iti.zcu.cz/history/2003/Hajek/problems/hajek-problems.ps
  * [4] G.Y. Katona, K. Varga, Minimally toughness in special graph classes, arXiv:1802.00055 [math.CO] (2018)
  * [5] C.G. Lekkerkerker and J.Ch. Boland, Representation of a finite graph by a set of intervals on the real line,  Fund. Math., 51 (1962), pp. 45–64.
  * [6] C. Dallard, B. Fernández, G. Y. Katona, M. Milanič, Conditions for minimally tough graphs, manuscript (2021)
  * [7] J.R.S. Blair and B. Peyton, An introduction to chordal graphs and clique trees. In Graph theory and sparse matrix computation, volume 56 of IMA Vol. Math. Appl., pp. 1–29. Springer, New York, (1993) doi:10.1007/978-1-4613-8369-7-1
  * [8] F. Gavril, The intersection graphs of subtrees in trees are exactly the chordal graphs. J. Combinatorial Theory Ser. B, 16 pp. 47–56, (1974) doi:10.1016/0095-8956(74)90094-x


11institutetext: EPCC at the University of Edinburgh, 47 Potterrow, Edinburgh,
UK

# Performance characterisation of the 64-core SG2042 RISC-V CPU for HPC

Nick Brown 11 0000-0003-2925-7275 Maurice Jamieson 11 0000-0003-1626-4871

###### Abstract

Whilst RISC-V has grown phenomenally quickly in embedded computing, it is yet
to gain significant traction in High Performance Computing (HPC). However, as
we move further into the exascale era, the flexibility offered by RISC-V has
the potential to be very beneficial in future supercomputers especially as the
community places an increased emphasis on decarbonising its workloads.
Sophon’s SG2042 is the first mass produced, commodity available, high-core
count RISC-V CPU designed for high performance workloads. First released in
summer 2023, and at the time of writing now becoming widely available, a key
question is whether this is a realistic proposition for HPC applications.

In this paper we use NASA’s NAS Parallel Benchmark (NPB) suite to characterise
performance of the SG2042 against other CPUs implementing the RISC-V, x86-64,
and AArch64 ISAs. We find that the SG2042 consistently outperforms all other
RISC-V solutions, delivering between a 2.6 and 16.7 performance improvement at
the single core level. When compared against the x86-64 and AArch64 CPUs,
which are commonplace for high performance workloads, we find that the SG2042
performs comparatively well with computationally bound algorithms but
decreases in relative performance when the algorithms are memory bandwidth or
latency bound. Based on this work, we identify that performance of the
SG2042’s memory subsystem is the greatest bottleneck.

###### Keywords:

RISC-V Sophon SG2042 NAS Parallel Benchmark suite (NPB) High Performance
Computing (HPC)

##  1 Introduction

RISC-V is an open Instruction Set Architecture (ISA) that, since it was first
released over a decade ago, has gained significant traction. At the time of
writing it was recently announced that over 13 billion RISC-V CPU cores have
been manufactured, but many of these are in embedded computing such as
automotive, space, and micro-controllers. RISC-V has yet to become commonplace
in High Performance Computing (HPC), but as the HPC community moves further
into the exascale era and there is an increased emphasis on decarbonisation of
workloads, we need to consider how to best deliver both increased performance
and greater energy efficiency. To this end, there is a renewed interest in new
hardware solutions and technologies built atop RISC-V have a strong potential
here as they can offer specialisation whilst still providing a common software
ecosystem.

Sophon’s SG2042 is the first high core count commodity available RISC-V CPU
designed for high performance workloads. First released in summer 2023, this
mass produced, 64-core RISC-V CPU is aimed at high performance workloads. Not
only does this processor provide significantly more cores that existing, SoC
based, mass produced RISC-V CPUs, but furthermore the T-Head XuanTie C920
cores themselves have been designed for high performance. Consequently this
new RISC-V CPU is very interesting to the HPC community and previous work [2]
found that, for the RAJAPerf suite [3], it delivers a considerable performance
uplift compared to existing commodity available RISC-V CPUs, but struggled to
match a set of x86-based CPUs that are commonplace in HPC machines. In this
paper we leverage NASA’s NAS Parallel Benchmark (NPB) suite to undertake more
in depth performance characterisation of the SG2042. Running this suite across
CPUs that implement the RISC-V, x86-86 and AArch64 ISAs, and in the later two
categories because we have selected CPUs that are used in production
supercomputers, we are able to better understand the types of workloads that
the SG2042 suits and where it might fall short.

##  2 Background

###  2.1 The Sophon SG2042

The Sophon SG2042 CPU is a 64-core processor running at 2GHz and organised in
clusters of four XuanTie C920 cores. Each 64-bit core, designed by T-Head, is
designed for high performance workloads and adopts a 12-stage out-of-order
multiple issue superscalar pipeline design [7]. Implementing the RV64GCV
instruction set, the C920 has three decode, four rename/dispatch, eight
issue/execute and two load/store execution units. Version 0.7.1 of the
vectorisation standard extension (RVV v0.7.1) is supported [11], with a vector
width of 128 bits. Each C920 core contains 64KB of L1 instruction (I) and data
(D) cache, 1MB of L2 cache which is shared between the cluster of four cores,
and 64MB of L3 system cache which is shared by all cores in the package. The
SG2042 also provides four DDR4-3200 memory controllers, and 32 lanes of PCI-E
Gen4. The CPU we use for the benchmarking in this paper is contained in a
Pioneer Box by Milk-V which has 128GB of DDR4 RAM.

The SG2042’s C920 core only provides RVV v0.7.1 which is not supported by
mainline GCC or LLVM. To this end, T-Head have provided their own fork of the
GNU compiler (XuanTie GCC) which has been optimised for their processors and
supports RVV v0.7.1. It has been found [5] that GCC8.4, which is part of their
20210618 release, provides the best auto-vectorisation capability and-so this
is the version we use for the benchmarking experiments undertaken in this
paper. Their version of the compiler generates Vector Length Specific (VLS)
RVV assembly which specifically targets the 128-bit vector width of the C920.
All codes are compiled at optimisation level three, and all reported results
are averaged over five runs. At the time of execution each benchmark run
reported in this paper was making exclusive use of the machine.

In [2] the authors benchmarked the SG2042 using the RAJAPerf suite, however
this was across a large number of individual kernels and from the results it
was difficult to isolate and identify individual performance patterns. By
contrast, in this paper we characterise and explore each individual benchmark
of the NPB suite to better classify the performance properties of the SG2042.

###  2.2 NAS Parallel Benchmarks (NPB) suite

The NAS Parallel Benchmark (NPB) suite [1] is a collection of benchmarks
developed by NASA’s Advanced Supercomputing (NAS) division to characterise HPC
systems, especially for Computational Fluid Dynamics (CFD) applications. First
released in the mid 1990s, in this paper we leverage the original eight
benchmarks in the suite, which comprises five kernels and three pseudo
applications. The kernels capture key algorithmic patterns that are ubiquitous
throughout HPC codes and test key performance characteristics that are
important across many workloads. The pseudo applications combine multiple
kernels to provide more complicated workloads. All these benchmarks are
configured using a variety of problem sizes known as classes. There are a
variety of implementations of the suite provided by NAS, including the OpenMP
and MPI versions that we use here, and throughout this paper use the official
code without any modifications.

Table 1: Summary of memory behaviour for NPB benchmarks on a Xeon Platinum 8170 Benchmark |  Clock ticks cache stall |  Clock ticks DDR stall |  Time DDR bandwidth bound
---|---|---|---
Integer Sort (IS) | 35% | 0% | 16%
Multi Grid (MG) | 34% | 20% | 88%
Embarrassingly Parallel (EP) | 11% | 0% | 0%
Conjugate Gradient (CG) | 19% | 18% | 0%
Fast Fourier Transform (FT) | 13% | 9% | 18%
Block Tridiagonal (BT) | 8% | 9% | 0%
LU Gauss Seidel (LU) | 12% | 11% | 0%
Scalar Pentadiagonal (SP) | 20% | 21% | 0%

Table 1 summarises, for each benchmark in the suite, the memory behaviour when
run using OpenMP on all 26 physical cores of a Xeon Platinum 8170. The _Clock
ticks cache stall_ and _Clock ticks DDR stall_ columns report how often the
CPU was stalled on cache and main memory accesses respectively, and the _Time
DDR bandwidth bound_ column reports the percentage of execution time that
there was a high DDR bandwidth utilisation.

The IS kernel tests indirect, random, memory accesses which it can be seen
stalls a significant fraction of the CPU due to cache accesses. It can be seen
that the MG kernel is heavily memory bound both in terms of time stalled on
cache and main memory accesses, and also the percentage of execution time
where DDR is under high utilisation. By contrast, the EP benchmark is designed
to test compute performance and there are far fewer cycles stalled on memory
access, and no time spent with high DDR bandwidth utilisation. CG comprises
irregular memory access and nearest neighbour communication, which results in
around 37% of clock ticks stalled on cache or DDR accesses, and the FT
benchmark requires all-to-all communications between ranks to undertake a
parallel transposition of data. For FT it can be seen that whilst there is
only 22% of clock ticks stalled, which is lower than the five other kernels
apart from EP, the kernel is utilising a high DDR bandwidth for 18% of the
time.

The BT, LU and SP pseudo application benchmarks are more complicated than the
five NPB kernels, and represent common, real-world, HPC use-cases. All three
of these pseudo applications compute a finite difference solution to the 3D
compressible Navier Stokes equations, where the LU benchmark solves this via a
block-lower block-upper triangular approximation based upon Gauss Seidel
iterative method [8]. The BT and SP benchmarks solve the same problem as LU,
but base their solution on a Beam-Warming approximation. In BT the resulting
equations are block-tridiagonal whereas in SP are fully diagonalised [8]. Both
these systems are solved using Gaussian elimination. It can be seen from Table
1 that, out of these three pseudo applications, BT stalls the least on memory
accesses and SP the most.

##  3 RISC-V core comparison

In this section we compare performance of existing commodity RISC-V solutions.
Due to the difference in core counts between RISC-V CPUs, we focus here on
single cores performance to understand how the XuanTie C920 core of the Sophon
SG2042 performs against other widely available RISC-V cores. We compare
against the U74 core [9] which is contained in the JH7200 and JH7100 SoCs of
the VisionFive V2 and V1 respectively, and both of these boards contain 8GB of
DRAM. We also compare against the SiFive Freedom U740 SoC, also containing the
U74 core and 16GB of DDR, and the T-Head XuanTie C906 [7] in the AllWinner D1
SoC with 1GB of memory.

Table 2: Single core comparison between RISC-V technologies with performance reported in Mops/s (Higher is better) using NPB kernels running at class B. In red is the percentage performance delivered compared to the C920 core of the SG2042. Benchmark | SG2042 | VisionFive V2 | VisionFive V1 | SiFive U740 | All Winner D1
---|---|---|---|---|---
IS | 60.6 |  17.84 _(29%)_ |  6.36 _(10%)_ |  9.09 _(15%)_ |  5.41 _(9%)_
MG | 1210.05 |  288.65 _(24%)_ |  72.31 _(6%)_ |  90.28 _(7%)_ |  163.19 _(13%)_
EP | 31.35 |  12.01 _(38%)_ |  7.55 _(24%)_ |  9.08 _(29%)_ |  9.23 _(29%)_
CG | 205.25 |  43.61 _(21%)_ |  21.96 _(11%)_ |  20.09 _(10%)_ |  12.99 _(6%)_
FT | 857.64 |  245.99 _(29%)_ |  88.35 _(10%)_ |  116.59 _(14%)_ | DNR

Table 2 reports a single core performance comparison between these RISC-V
technologies, for the five NPB kernels at class B, with performance reported
in million operations per second (Mop/s) and a higher number is better. In
italicised red is the percentage performance that a single core of this CPU
delivers compared to a single C920 core found in the SG2042. It can be seen
that, irrespective of the kernel, the C920 significantly out performs all
other RISC-V technologies. Consistently, the U74 of the VisionFive V2 performs
closest to the C920, but is still only delivering between 21% and 38% the
performance of the C920. Whilst the VisionFive V1 and SiFive U740 both contain
the same U74 core as the VisionFive V2, they are significantly slower and this
is broadly in agreement with [2].

The C906 of the All Winner D1 is out performed by the C920 and the U74 of the
V2 quite considerably. However, this is the cheapest of the SoCs considered
here, and the C906 outperforms the V1 and U740 for the EP and MG benchmarks.
Given the performance profile of the benchmarks reported in Table 1, this
suggests that the raw compute power of the C906 is similar to that of the U74
and the memory bandwidth is greater on the All Winner D1 than the VisionFive
V1 and SiFive U740. However, for those benchmarks with more complex,
irregular, memory patterns such as IS and CG, the C906 seems to struggle
compared to the other RISC-V cores. Incidentally, it was not possible to run
the FT benchmark on the All Winner D1 due to the limited 1GB of memory
becoming exhausted.

In this section we therefore conclude that the C920 core of the SG2042
significantly outperforms all other commodity available RISC-V CPU cores.
Whilst this is in agreement with [2], in this section we have compared against
a wider range of RISC-V CPUs than [2] and for specific algorithmic patterns
that are very commonly found in HPC codes, especially for CFD. Considering
that they use the same U74 core, it is surprising that the VisionFive V2
outperforms the V1 and U740 by such a large margin, but again this is in
agreement with [2] and [6], and one of the reasons for this is that the V2 is
running at 1.5GHz compared to 1.2GHz for both the V1 and U740.

##  4 Comparing the SG2042 against other architectures

In Section 3 we compared the performance of the SG2042’s C920 core against
other RISC-V commodity available CPU cores. Whilst it is interesting to
explore performance against RISC-V CPUs, and indeed the C920 delivers
impressive performance compared to other RISC-V hardware, to understand
whether the SG2042 is a contender for HPC it is far more instructive to
benchmark against CPUs of other architectures that are commonly used for HPC
workloads.

Table 3: Summary of CPUs that are benchmarked in this section CPU | ISA | Part | Base clock |  Number of cores | Vector
---|---|---|---|---|---
AMD EPYC | x86-64 | EPYC 7742 | 2.25GHz | 64 | AVX2
Intel Skylake | x86-64 | Xeon Platinum 8170 | 2.1 GHz | 26 | AVX512
Marvell ThunderX2 | ARMv8.1 | CN9980 | 2 GHz | 32 | NEON
Sophon SG2042 | RV64GCV | SG2042 | 2 GHz | 64 | RVV v0.7.1

In this section we compare against CPUs of other architectures that are
commonplace in HPC, and these are summarised in Table 3. The AMD EPYC is the
Rome series of AMD CPUs, containing the Zen-2 micro architecture and we run
this on ARCHER2, a Cray EX and the UK national supercomputer. Similarly to the
SG2042, the AMD EPYC contains 64 physical cores across four NUMA regions, each
with 16 cores, but has eight instead of four memory controllers and memory
channels. Each core in the AMD EPYC contains 32KB of I and D L1 cache, 512 KB
of L2 cache, and there is 16MB of L3 cache shared between four cores.
Providing AVX2, the EPYC 7742 has 256-bit wide vector registers, which is
double that of the SG2042, but is capable of processing two AVX-256
instructions per cycle. ARCHER2 contains 256GB of DDR memory. We use GCC
version 11.2 when compiling on ARCHER2. Simultaneous Multithreading (SMT) is
disabled for our runs, which is the default configuration on ARCHER2.

We also compare against an Intel Skylake Xeon Platinum 8170, which is the same
CPU used to profile the NPB benchmarks in Table 1. This Skylake-SP CPU
contains 26 cores, each with 32KB of I and D L1 cache, 1MB of L2 cache and
1.375MB of L3 cache (the later is shared across all cores). The Skylake
supports AVX512, double and quadruple the width of the EPYC 7742 and SG2042
respectively, and each Skylake core has two FPUs. The machine we run on has
192GB of DDR4 memory, and we use GCC version 8.4.

Lastly, we compare against the CN9980 Marvell ThunderX2 which contains 32
cores implementing the ARMv8.1 (AArch64) ISA via the Vulcan micro
architecture. Each core contains 32KB of I and D L1 cache, as well as 256KB of
L2. There is a total of 32MB L3 cache, 1MB per core, shared by the entire
chip. NEON is supported, which provides 128 bit wide vector registers and this
is interesting because it matches the vector width of the C920 core in the
SG2042. Similarly to the Skylake, the Marvell ThunderX2 has two FPUs per core.
This is the CPU used in Fulhame, an HPE Apollo 70 system, with 128GB of DDR
per node, we use GCC version 9.2 and SMT is also disabled in our runs.

For the performance comparison undertaken in this section, we run class C of
the NASA Parallel Benchmarks and run over multiple cores of the CPUs by using
the OpenMP implementations of the benchmarks [4]. Each thread is mapped to an
individual physical CPU core, all reported results are averaged over five runs
and all codes are built at optimisation level three.

###  4.1 Integer Sort (IS)

As described in Section 2.2, the Integer Sort (IS) benchmark is concerned with
integer comparison and indirect, random, memory access performance. Figure 1
illustrates the performance results for this benchmarks across our CPUs of
interested, reported in Mops/s (higher is better). It can be seen that the
SG2042 performs considerably worse than all other CPUs with performance
plateauing at 16 cores. By contrast, the ThunderX2 and AMD EPYC delivering
similar performance until the 32 cores of the ThunderX2 are exhausted. The
Skylake performs better than all the other CPUs, but is limited by its lower
core count, where the ThunderX2 catches up to the Skylake at 32 cores and the
EPYC outperforms it at 64 cores.

Figure 1: IS benchmark performance (higher is better) parallelised via OpenMP

It can be seen from Figure 1 that the SG2042 struggles significantly with this
benchmark where as was seen in Table 1 the irregular, random, memory accesses
result in a comparatively large number of time stalled due to cache access and
DDR bandwidth utilisation is high for a small fraction of the runtime. A
hypothesis is that could be due to the cache hierarchy, where the Skylake
which performs the best has the largest L2 cache, 1MB per core, compared to
256KB (per core, 1MB shared between four cores) for the SG2042, 256KB for the
ThunderX2 and 512KB for the AMD EPYC. The surprise here is in the performance
difference between the SG2042 and the ThunderX2, as per core they both have
the same amount of L2 and L3 cache.

###  4.2 Multi Grid (MG)

It was illustrated in Section 2.2 that the Multi Grid (MG) benchmark is
heavily memory bandwidth bound, and results of executing this benchmark kernel
on the CPUs of interest is illustrated in Figure 2. It can be seen that the
AMD EPYC provides considerably best performance, with the Skylake and
ThunderX2 delivering similar performance and both plateauing at 16 cores where
memory bandwidth is likely saturated. By contrast, the SG2042 lags the other
CPUs considerably, also plateauing at 16 and 32 cores but then with a
performance increase at 64 cores.

Figure 2: MG benchmark performance (higher is better) parallelised via OpenMP

The memory configuration of the CPUs partially helps to explain the relative
performance reported in Figure 2. The AMD EPYC has 8 memory controllers and 8
memory channels, connected to DDR4-3200 memory. By contrast, the Skylake and
ThunderX2 both only have 2 memory controllers and are both connected to
DDR4-2666 albeit with the ThunderX2 having 8 memory channels compared to 6
memory channels in the Skylake. The SG2042 has four memory controllers and
only four memory channels, connected to DDR4-3200. Whilst there are fewer
memory channels on the SG2042 than the other CPUs, it also has double the
memory controllers than the Skylake and ThunderX2 CPUs, and also faster
memory, but lags performance compared with those CPUs. Details around the
memory subsystem on the SG2042 are difficult to come by, but it is our
hypothesis that the memory controllers on the SG2042 are considerably less
advanced than the other CPUs considered in this section.

The behaviour of the MG benchmark also helps explain one of the anomalies of
the IS benchmark performance. It was our hypothesis that the L2 and L3 cache
design was in part governing performance of the SG2042 compared to other CPUs,
however the ThunderX2 also has the same cache design but was faster than the
SG2042 for that benchmark. However, as seen for the MG benchmark, the SG2042
is also severely memory bandwidth bound and this likely explains the gap in
performance between the SG2042 and ThunderX2 for the IS benchmark.

###  4.3 Embarrassingly Parallel (EP)

The Embarrassingly Parallel (EP) benchmark is compute bound, and results of
this on our CPUs are illustrated in Figure 3. It can be seen that across the
CPUs being tested, there are two groups; the SG2042 and ThunderX2 share very
similar performance but with the SG2042’s 64 cores then making a significant
difference compared to the 32 cores of the ThunderX2. The EPYC and Skylake
both deliver similar performance, which is greater than the SG2042 and
ThunderX2, but the 26 cores of the Skylake are a disadvantage against the
SG2042 which then significantly out performs the Skylake at 64 cores. The AMD
EPYC performs best out of all the CPUs, especially at the larger core counts.

Figure 3: EP benchmark performance (higher is better) parallelised via OpenMP

This performance behaviour is in stark contrast to the IS and MG benchmarks,
and demonstrates that the SG2042 and ThunderX2 deliver very similar compute
performance at the same number of cores. This makes some sense given that they
both provide 128-bit vectorisation, albeit with the ThunderX2 having two FPUs
per core compared to one on the SG2042. By contrast, the Skylake and AMD EPYC
CPUs provides wider vectorisation, 512-bit and 256-bit respectively and this
in part helps explain the performance difference between these two groups.
However, given the Skylake provides AVX512, and the AMD EPYC only AVX2, and
that the ThunderX2 has two FPUs per core and the SG2042 only one, clearly the
GCC compiler is not able to fully vectorise the code of this benchmark and
make full use of the FPUs.

Given that, at the largest number of cores each technology provides, the
SG2042 performs second best for this benchmark out of all CPUs in our
comparison, this demonstrates that for compute bound problems the large core
count of the SG2042 is beneficial.

###  4.4 Conjugate Gradient (CG)

As can be seen from Table 1, the Conjugate Gradient (CG) benchmark also spends
considerable time stalled on cache and DDR memory accesses, and this is
because it comprises of irregular memory access and nearest neighbour
communications. Figure 4 illustrates the performance of this benchmark kernel
across our CPUs. Given the performance of the IS and MG benchmarks it is no
surprise that the SG2042 falls short of the other technologies, but it is
closer to the ThunderX2 than we had expected delivering around 50% the
performance of the ThunderX2 at 32 cores.

Figure 4: CG benchmark performance (higher is better) parallelised via OpenMP

Potentially, what is making the difference here is in the size of the L3
cache, where the AMD EPYC has 16MB L3 shared between four cores (4MB per
core), whereas the Skylake has 1.3MB L3 cache per core shared across all cores
(potentially helped by the larger 1MB L2 cache). By contrast, both the
ThunderX2 and SG2042 have the same size of 256MB L2 and 1MB L3 cache per core.
This would help explain the performance differences, with the additional
memory bandwidth limitations of the SG2042 causing additional the overhead
which reduces performance further.

###  4.5 Fast Fourier Transform (FT)

The fast Fourier Transform (FT) benchmark requires all-to-all communication
between ranks, and the performance of this benchmark can be seen in Figure 5.
As described in Table 1, for this benchmark on the Skylake there was some
stalling due to cache and DDR access (13% and 9% respectively) but also for
18% of the time DDR was under high utilisation. Once again, the SG2042 is
significantly slower than the other CPUs, with the ThunderX2 sitting around
half way between the performance of the x86 CPUs and the SG2042 and this is
likely for the same reasons explored for the CG benchmark.

Figure 5: FT benchmark performance (higher is better) parallelised via OpenMP

###  4.6 Pseudo Applications

Table 4 reports performance for the BT, LU and SP benchmarks and this is
expressed as how many times faster this CPU is than the SG2042. Given the
findings that the SG2042 struggles to perform when there is increased pressure
on the memory subsystem, and based upon the stall numbers reported in Table 1,
it was our expectation that the SG2042 would perform best for the BT benchmark
and worst for the SP benchmark with the LU benchmark in between. This is
broadly the case based upon the figures in Table 4, where each number reports
the number of times faster than each of the other CPUs are for each pseudo
application at a specific core count. It can be seen that three other CPUs
significantly outperform the SG2042 for the three pseudo applications.

Table 4: For each pseudo application, the number of times faster a specific CPU is than the SG2042 at the given number of cores Number | BT benchmark | LU benchmark | SP benchmark
---|---|---|---
cores | EPYC | Skylake | ThunderX2 | EPYC | Skylake | ThunderX2 | EPYC | Skylake | ThunderX2
16 | 3.23 | 3.28 | 2.43 | 3.65 | 4.15 | 2.86 | 5.01 | 3.91 | 3.65
26 | 3.57 | 2.97 | 2.69 | 3.20 | 3.16 | 2.62 | 6.25 | 3.48 | 3.57
32 | 3.68 | - | 2.64 | 3.40 | - | 2.94 | 5.26 | - | 3.22
64 | 4.19 | - | - | 2.95 | - | - | 4.22 | - | -

##  5 MPI vs OpenMP on the Sophon SG2042

In Section 4, the NPB benchmarks were all run using the official NAS OpenMP
implementation. This is sensible given that execution is occurring within a
single memory space, however a question is whether, for best performance, one
should write their parallel code using OpenMP or MPI within a node. These two
models are very different, where OpenMP follows a thread based approach and by
default all threads share the same memory area and explicit marshalling and
protection of shared memory is required. In contrast, when using MPI tasks run
as independent processes which share no data and instead communicate with each
other via explicit messages.

As there are both OpenMP and MPI implementations of the NPB benchmark suite,
we are able to undertake a direct comparison. This is illustrated in Figure 6
which reports the percentage performance delivered by the MPI implementation
compared to the OpenMP version, for each benchmark at different numbers of
cores. It can be seen that there are significant differences in performance
between OpenMP and MPI for some configurations.

Figure 6: Percentage performance of MPI based NPB benchmark compared to the
OpenMP benchmark implementation

It can be observed, in Figure 6, that over one or two cores the benchmarks do
not benefit from MPI compared to OpenMP. However, from four cores upwards the
CG and FT benchmarks run faster when using MPI than OpenMP. At 64 cores, all
the MPI implementations run faster than their OpenMP counterparts. When
exploring MPI against OpenMP on other architectures, MPI was always slower on
the AMD EPYC and this was also true on the Skylake apart for the CG benchmark
where MPI was on average two times faster than the OpenMP implementation. When
profiling the CG MPI benchmark on the Skylake, it was found that clock ticks
stalls due to cache reduced from 19% in the OpenMP implementation to 5.5% with
the MPI version. The MPI implementation experienced no clock ticks stalled due
to DDR accesses, down from 18% for the OpenMP version.

We surmise that the way these benchmarks are implemented means that the
structure imposed by the MPI implementation tends to put less pressure on the
memory subsystem when undertaking communications. For codes with heavy inter-
core communication such as CG (nearest neighbour point to point) and FT (all
to all collective) then this can be beneficial on the SG2042.

##  6 Conclusions

The Sophon SG2042 is an impressive RISC-V CPU and, using NASA’s NAS Parallel
Benchmark (NPB) suite, we have demonstrated that for benchmarks that closely
represent ubiquitous HPC algorithms, especially in CFD, it significantly out
performs existing RISC-V solutions. When compared against CPUs which implement
other ISAs and whose use is widespread for high performance workloads, the
Sophon SG2042 is outperformed by the AMD EPYC by between 1.77 and 15.06 times,
the Intel Skylake between 0.59 and 5.98 times, and the Marvel ThunderX2
between 0.59 and 5.91 times. The SG2042 is most competitive for
computationally bound codes, and its high core count ultimately out performed
the Skylake and ThunderX2 for the EP benchmark. However the SG2042 CPU
struggled with algorithms that are memory bandwidth or latency bound.

From this work we conclude that the memory subsystem of the SG2042 is a
bottleneck, and Sophon recently announced the SG2044 which is reported to have
three times the DDR memory bandwidth [10], as well as implementing RVV v1.0.
This has the potential to provide a very significant performance improvement
over the SG2042 for many of the benchmarks explored in this paper, and the
improved memory performance will likely ameliorate the bottlenecks we have
observed with the high core count continuing to deliver good computational
performance. It is therefore our conclusion that the Sophon SG family of
RISC-V CPUs has strong potential in HPC, and whilst the current generation
SG2042 is an impressive first generation, based on announcements made by
Sophon it looks likely that the SG2044 will address the key performance
challenges that we have observed in this paper and make for a compelling
RISC-V product family.

##  7 Acknowledgements

This work has been funded by the ExCALIBUR H&ES RISC-V testbed. This work used
the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk).
The Fulhame HPE Apollo 70 system is supplied to EPCC as part of the Catalyst
UK programme. For the purpose of open access, the author has applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising from this submission.

## References

  * [1] Benchmarks, N.P.: Nas parallel benchmarks. CG and IS (2006)
  * [2] Brown, N., Jamieson, M., Lee, J., Wang, P.: Is risc-v ready for hpc prime-time: Evaluating the 64-core sophon sg2042 risc-v cpu. In: Proceedings of the SC’23 Workshops of The International Conference on High Performance Computing, Network, Storage, and Analysis. pp. 1566–1574 (2023)
  * [3] Hornung, R.D., Hones, H.E.: Raja performance suite. Tech. rep., Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States) (2017)
  * [4] Jin, H.Q., Frumkin, M., Yan, J.: The openmp implementation of nas parallel benchmarks and its performance (1999)
  * [5] Lee, J.K., Jamieson, M., Brown, N.: Backporting risc-v vector assembly. In: International Conference on High Performance Computing. pp. 433–443. Springer (2023)
  * [6] Lee, J.K., Jamieson, M., Brown, N., Jesus, R.: Test-driving risc-v vector hardware for hpc. In: International Conference on High Performance Computing. pp. 419–432. Springer (2023)
  * [7] Open xuantie c906 (2023), https://xrvm.com/cpu-details?id=4056751997003636736
  * [8] Saphir, W., Van der Wijngaart, R.F., Woo, A., Yarrow, M.: New implementations and results for the nas parallel benchmarks 2. In: PPSC. Citeseer (1997)
  * [9] Sifive u74-mc core complex manual (2021), https://starfivetech.com/uploads/u74mc_core_complex_manual_21G1.pdf
  * [10] Sophgo risc-v roadmap (2024), https://github.com/RISCVtestbed/riscvtestbed.github.io/blob/main/assets/files/hpcasia24/hpc_asia_wang.pdf
  * [11] C920: Specifications (2023), https://xuantie.t-head.cn/product/xuantie/4082464366237126656


# Envelope Ejection and the Transition to Homologous Expansion in Common-
Envelope Events

Vinaya Valsan1, Sarah V. Borges1, Logan Prust1,2, Philip Chang1
1 Department of Physics, University of Wisconsin-Milwaukee, 3135 North
Maryland Avenue, Milwaukee, WI 53211, USA
2 Kavli Institute for Theoretical Physics, University of California, Santa
Barbara, CA 93106, USA<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We conduct a long-timescale ($5000\,$d) 3-D simulation of a common-envelope
event with a $2\,M_{\odot}$ red giant and a $1\,M_{\odot}$ main sequence
companion, using the moving-mesh hydrodynamic solver MANGA. Starting with an
orbital radius of $52\,R_{\odot}$, our binary shrinks to an orbital radius of
$5\,R_{\odot}$ in $200\,$d. We show that over a timescale of about $1500\,$d,
the envelope is completely ejected while $80$ per cent is ejected in about
$400\,$d. The complete ejection of the envelope is solely powered by the
orbital energy of the binary, without the need for late-time reheating from
recombination or jets. Motivated by recent theoretical and observational
results, we also find that the envelope enters a phase of homologous expansion
about $550\,\rm d$ after the start of our simulation. We also run a simplified
1-D model to show that heating from the central binary in the envelope at late
times does not influence the ejection. This homologous expansion of the
envelope would likely simplify calculations of the observational implications
such as light curves.

###### keywords:

binaries: close – hydrodynamics – stars: winds, outflows – methods: numerical

††pubyear: 2022††pagerange: Envelope Ejection and the Transition to Homologous
Expansion in Common-Envelope Events–Acknowledgements

##  1 Introduction

Common-envelope evolution (CEE; Paczynski 1976) is believed to be responsible
for the production of many close binary systems such as X-ray binaries, binary
neutron stars, binary black holes, and white dwarf binaries including
cataclysmic variables (for a review, see Ivanova et al., 2013). The physics
and astrophysics of these systems have been a subject of continuous study for
the last $50$ yr. However, the complex physics of this process, which includes
gravity, hydrodynamics, nuclear burning, recombination, and radiation, has
precluded much analytic progress. In the last decade, advances in computing
and algorithms have made high-resolution, long-timescale simulations of CEE
possible.

These long-timescale simulations have begun to unravel the relevant physics of
CEE and its astrophysical impact. In particular, studies carried out by
several groups on low-mass binaries such as an asymptotic giant branch (AGB)
primary with a main-sequence (MS) companion star (Sand et al., 2020; Chamandy
et al., 2020; Ondratschek et al., 2022) and a red-giant branch (RGB) primary
with an MS star (Iaconi et al., 2019) have demonstrated that the envelope is
completely ejected on sufficiently long timescales (years to decades). These
groups differ in their conclusions of what physics is important. For instance,
Iaconi et al. (2019), Sand et al. (2020), and Ondratschek et al. (2022) argue
that recombination energy is crucial in envelope ejection, but Chamandy et al.
(2020) argue otherwise. This difference is likely due to the limited number of
cases studied and future studies will likely bring this physics into sharper
focus.

In addition, recent work by Iaconi et al. (2019) argues that the expansion of
the envelope leads to homologous expansion on long timescales. This should not
come as a surprise because any expansion leads to homologous expansion so long
as the trajectories remain ballistic. However, if the envelope expands
homologously, this helps to simplify the theory of CEE and would provide an
easier means to compute observables.

Recent observations have also hinted at the presence of homologously expanding
ejecta in CEE and stellar merger events. For instance, the properties of
observed CO emission in V4332 Sgr is well reproduced with the homologous
expansion model (Kamiński et al., 2018). Additionally, Kamiński et al. (2020)
showed that the observed properties of the molecular remnant of Nova 1670 (CK
Vulpeculae) are satisfactorily reproduced by linear velocity fields.

In this paper, we study the physics of the ejection of the envelope using
long-timescale simulations. Starting with our recent work (Prust & Chang,
2019; Prust, 2020), we optimize our numerical techniques to allow for an order
of magnitude increase in integration time. We show that for the case of a
$2\,M_{\odot}$ RGB and a $1\,M_{\odot}$ MS companion, we achieve complete
envelope ejection in $1500\,$d. We also show that the envelope enters a
homologous phase early on (about $550\,$d) and that the morphology of the
ejected material is roughly spherical.

The paper is organized as follows. In Section 2, we discuss the numerical
setup of our calculations. We follow Prust & Chang (2019) and Prust (2020),
but describe a significant improvement in how we generate initial conditions,
giving substantial speedups. In Section 3, we discuss the results of these
simulations and show complete envelope ejection after $1500\,$d. We also
demonstrate that the envelope enters a homologous phase early on and that it
can be approximated as spherical. Motivated by these results, we develop a 1-D
numerical model in Section 4 and compare these simplified calculations with
the full 3-D calculation. In Section 5, we discuss the major theoretical and
observational implications of our results and close in Section 6.

##  2 Numerical Setup

We use the moving-mesh hydrodynamic solver for CHANGA, which we call MANGA
(Chang et al., 2017; Prust & Chang, 2019), to study CEE. MANGA is a moving-
mesh hydrodynamic simulation code based on the algorithms described by
Springel (2010). A detailed description of MANGA is presented by Chang et al.
(2017). Improvements such as the multiple time-stepping algorithms and
integration of stellar equations of state are presented by Prust & Chang
(2019). A discussion of its use for simulations of main sequence tidal
disruption events is presented by Spaulding & Chang (2021). Finally, recent
code improvements are discussed by Chang et al. (2020) for radiation
hydrodynamics, Chang & Etienne (2020) for general relativistic hydrodynamics,
Prust (2020) for moving-boundary conditions, and Prust & Chang (in
preparation) for magnetohydrodynamics. We refer the interested reader to this
literature for a detailed description of MANGA.

###  2.1 Initial Conditions

We use the same procedure as Prust & Chang (2019) to construct initial
conditions. We evolve a $2\,M_{\odot}$ star with MESA (Paxton et al., 2011;
Paxton et al., 2013, 2015, 2018, 2019; Jermyn et al., 2022) from the pre-main
sequence to the red giant phase and stop when it reaches a radius of around
$52\,R_{\odot}$. As done by Prust & Chang (2019), we construct a star of mass
$2\,M_{\odot}$, whose entropy profile matches that of the original star. The
core of the newly-constructed star is modelled as a dark matter particle with
gravitational softening for computational expediency. We then map the radial
profile of density and temperature to an unstructured particle (mesh-
generating point) mesh. The simulation consists of 430K mesh-generating
points, of which 80K model the star. The companion star, which is also
modelled as a softened dark matter particle, is then placed at the surface of
the red giant.

We should note that placing the companion on the surface of the red giant as
an initial condition is unrealistic for a couple of reasons. First, the red
giant model used at the beginning of the simulation is only in hydrostatic
equilibrium in isolation and thus does not account for the effects of the
companion. Second, a more realistic scenario would involve the binary system
evolving, allowing the red giant to slowly expand on a nuclear timescale until
it fills its Roche lobe. At this point, it undergoes unstable mass transfer,
driving the system into CEE. However, this more realistic scenario is not
easily realizable in numerical simulation as it involves the slow evolution of
the red giant in the thermal timescale and would demand a realistic treatment
of nuclear burning and radiation. The computational costs of such a simulation
would be prohibitive. Thus, we have simplified the initial conditions to the
ones stated and anticipate that for the long-term evolution of the CEE event,
the initial conditions do not play a large role.

In this paper, we have made a number of modifications in an effort to reduce
the computational cost. First, we use an adiabatic equation of state
($\gamma=5/3$) instead of a MESA equation of state. The simulation using the
MESA equation of state is significantly more computationally expensive than
the adiabatic case. The primary advantage of using the MESA equation of state
is that it encodes additional information in regard to recombination energy,
which may be important in ejecting the envelope for low-mass systems.
Nevertheless, as we will demonstrate below, orbital energy alone achieves
complete envelope ejection without the need for recombination energy. We also
improved the grid generation for the tenuous atmosphere surrounding the stars.
In particular, we increase the spacing between grid points exponentially
outside the stars up to the final coarsest resolution instead of the power-law
increase used previously by Prust & Chang (2019). This reduces the number of
nearest-neighbour searches for mesh-generating points near the boundary of the
star and atmosphere. We find that this improves the performance of the code by
around a factor of 2.

These reductions in computational costs and improved computing power allow us
to run our simulation for $5000$ d, which is longer than our previous
simulations (Prust & Chang, 2019) by a factor of around $20$. It is also in
line with other recent long-timescale simulations (Iaconi et al., 2019; Sand
et al., 2020).

##  3 Results

Figure 1: Projection of density on to the $x-z$ plane (left panel) and $x-y$
plane (right panel) at different time slices ($503\,$d, $1006\,$d, $2013\,$d)
from the simulation. The ‘$+$’ sign marks the centre of mass of the system.
Figure 2: Smoothed separation between the centre of mass of the red giant and
the main-sequence star as a function of $t$. There is an initial rapid plunge
of the two centres toward each other, but this plateaus to about
$5\,R_{\odot}$ after $200\,$d.

We show a series of projected density plots in Fig.1, projected along the axis
of the orbital plane (z-axis, right column) and projected along a direction in
the orbital plane (x-axis, left column) at $503\,$d, $1006\,$d and $2013\,$d.
The ‘+’ sign marks the centre of mass of the system. Initially, these
projections demonstrate that the ejected matter is axisymmetric, but not
spherical. The overall shape remains fairly constant, though it does become
more spherical as it evolves. This constancy of the overall shape and
evolution to greater spherical symmetry will be important in our discussion of
homologous expansion and the spherical approximation.

In Fig. 2, we plot the orbital separation, $a$, between the centres of the
stars. Each $a$ in the plot is taken to be the average of $a$ in a time
interval of 14 d. In the beginning, the two centres plunge toward one another
in a period of rapid orbital decay. This starts to slow down and at $200$ d,
the orbital decay plateaus to an orbital separation of $a\approx
5\,R_{\odot}$, which is similar to our previous result of $a\approx
3.6\,R_{\odot}$ (Prust & Chang, 2019).

Figure 3: Fraction of unbound mass, $f_{\rm unb}$, as a function of time, $t$.
We only consider the mechanical energy in this case.

For each $\rm i$-th particle/mesh-generating point in the system, the total
mechanical energy is defined as

$\displaystyle E_{\rm mech,i}=m_{\rm i}\left(\frac{1}{2}v_{\rm
i}^{2}+\phi_{\rm i}\right),$ (1)

where $m_{\rm i}$ is the total mass of the mesh-generating point, $v_{\rm i}$
is the fluid velocity of the mesh-generating point relative to the bound
centre of mass of the system, and $\phi_{\rm i}$ is the gravitational
potential. Particles with $E_{\rm mech,i}<0$ are bound to the binary while
those with $E_{\rm mech,i}>0$ are unbound. As discussed by Prust & Chang
(2019), we must carefully define the velocities relative to that of the centre
of mass of the bound material. This involves an iterative computation to find
the bound mass and the centre of mass velocity. The unbound mass fraction is
then defined as the fractional mass of the material with positive total
energy.

In Fig. 3 we plot the unbound mass fraction, $f_{\rm unb}$, as a function of
time. Nearly all gas from the red giant is unbound after $1500\,$d and over
$80$ per cent at $400\,$d. At $250\,$d, the fraction of mass that is unbound
compared to our previous result (Prust & Chang, 2019) is substantially larger
(40 per cent vs. 10 per cent), considering just the mechanical energy.
However, the equation of state is different between the two simulations (ideal
gas vs MESA). In addition, envelope ejection occurs in the absence of
additional late-time energy injection via hydrogen recombination and/or jets.
Other work has also recently demonstrated complete envelope ejection on a time
scale of about $1000\,$d, but these results can rely on additional late-time
energy injection.

###  3.1 Homologous Expansion

Recently, Iaconi et al. (2019) showed in their long-timescale CEE simulations
that the envelope ejection follows a homologous expansion approximation. In
their work, they simulate the evolution of a $0.88\,M_{\odot},83\,R_{\odot}$
RGB and a $0.6\,M_{\odot}$ companion star and follow the system for about $15$
yr. They show that the external layers of the envelope become homologous as
soon as they are ejected, but that it takes about $14$ yr for the bulk of the
unbound gas to enter homologous expansion. Motivated by this result, we
investigate the onset of homologous expansion in our simulations.

To begin, we recall that the distinguishing characteristic of homologous
expansion is that the velocity follows a radial profile $v\propto r$. In
essence, this means that fluid elements are on ballistic trajectories with
little or no interaction between fluid elements or external forces. As such,
the radial position of a fluid element can then be written as

$\displaystyle r(t)=v_{\rm r}t_{\rm h},$ (2)

where $v_{\rm r}$ is the radial velocity of the fluid element and $t_{\rm h}$
is the homologous expansion time. We note that while the formalism of
homologous expansion is relatively simple and discussed widely in the
literature, we discuss it here again to define it in terms of computational
quantities like time, which is defined as zero at the beginning of a
simulation and has no relation to the zero time in homologous expansion. Thus,
we define $t_{\rm h}=t-t_{0}$, where $t$ is the time since the start of the
simulation and $t_{0}$ is some fitted time that defines the $t=0$ point of
homologous expansion. Indeed, $t_{\rm h}$ is mapped exactly to the time in
most discussions of homologous expansion. We can rewrite equation (2) as

$\displaystyle v_{\rm r}=\frac{r}{t-t_{0}}=\frac{r}{t_{\rm h}},$ (3)

and by differentiating and integrating the above equation, we can write the
position of a fluid particle at time $t$ with respect to an initial time
$t_{\rm i}$ as

$\displaystyle r(t)=r_{\rm i}\frac{t-t_{0}}{t_{\rm i}-t_{0}}=r_{\rm
i}\frac{t_{\rm h}}{t_{\rm{h,i}}},$ (4)

where $r_{\rm i}$ is the radial position at a time $t_{\rm i}>t_{0}$. From
equation (4), we now define a scaled radius with respect to the initial time
given the current radius for any fluid element in the simulation

$\displaystyle r_{\rm s}(r,t)=r\frac{t_{\rm h,i}}{t_{\rm h}}=r\frac{t_{\rm
i}-t_{0}}{t-t_{0}}.$ (5)

In other words, $r_{\rm s}$ maps the position of a fluid element, $r$, at a
particular time $t$ to the position of a fluid element at the initial time
$t_{\rm i}$. Having defined the scaled radius, $r_{\rm s}$, we also define the
scaled density and velocity as

$\displaystyle\rho_{\rm s}(r,t)$ $\displaystyle=$ $\displaystyle\rho_{\rm
i}(r_{\rm s}(r,t))\left(\frac{t_{\rm h}}{t_{\rm h,i}}\right)^{3}$ (6)
$\displaystyle v_{\rm r,s}(r,t)$ $\displaystyle=$ $\displaystyle v_{\rm
r,i}(r_{\rm s}(r,t)).$ (7)

The simple rescaling given in equations (6) and (7) is insufficient to
describe the entire system. While it works for the expanding envelope, it does
not describe the tenuous atmosphere. Toward that end, we define the radius of
the envelope, $R(t)$, which is the position of the outer boundary of the
homologously expanding region inside which equations (5), (6) and (7) are
valid. We define a dimensionless radius $\eta$ as

$\displaystyle\eta=\frac{r}{R},$ (8)

so that $\eta=0$ at the centre of the CE and $\eta=1$ at the envelope’s outer
boundary. From equation (4), we can estimate $R$ for any given time with
respect to the initial time $t_{\rm i}$: $R=R_{\rm i}({t_{\rm h}}/{t_{\rm
h,i}})$. We can then redefine scaled density, $\rho_{\rm s}$, and scaled
velocity, $v_{\rm s}$, for regions inside and outside of the envelope. This
gives

$\displaystyle\rho_{\rm s}$ $\displaystyle=$
$\displaystyle\begin{cases}\rho_{\rm i}(r_{\rm s}(r,t))\left(\frac{t_{\rm
h}}{t_{\rm h,i}}\right)^{3}&\text{if $\eta\leq 1$}\\\ \rho_{\rm b}&\text{if
$\eta$ > 1}\end{cases}$ $\displaystyle v_{\rm r,s}$ $\displaystyle=$
$\displaystyle\begin{cases}v_{\rm r,i}(r_{\rm s}(r,t))&\text{if $\eta\leq
1$}\\\ 0&\text{if $\eta$ > 1}\end{cases}$ (9)

Figure 4: Radial velocity (top) and scaled density (bottom) as a function of
scaled radius for different times ($2000$ d, $2500$ d and $3000$ d). We also
show a power-law fit (solid black line) for $v_{\rm r}\propto r_{\rm
s}^{0.86}$ (top). Figure 5: The fractional change in the velocity of all
particles compared to the final average velocity as a function of time. The
dashed line represents the mean and the shaded area represents the standard
deviation of the fractional change in the particle velocity. Figure 6: Ratio
of thermal energy to kinetic energy as a function of time.

To show that our simulations follow the scaling defined by equations (5) and
(9), we select $t_{\rm i}\approx 800\,$d and fit $t_{0}\approx 550\,$d so that
the re-scaled radial velocities $v_{\rm r,s}$ match one another between a few
$\times 10^{2}\,R_{\odot}$ and a few $\times 10^{3}\,R_{\odot}$. We show this
result in the top plot of Fig. 4. As this plot shows, the velocities for
$r_{\rm s}>$ a few $\times 10^{2}\,R_{\odot}$ match each other for a number of
different time steps. We also fit a power law between $6\times
10^{2}\,R_{\odot}$ and $5\times 10^{3}\,R_{\odot}$, and the resulting fit is

$\displaystyle v_{\rm r}=4.8\left(\frac{r_{\rm
s}}{100\,R_{\odot}}\right)^{0.86}{\rm km\,s^{-1}}.$ (10)

The radial power law exponent is about $1$, which is consistent with
homologous expansion. Thus we find that by $2000\,$d since the start of the
simulation, homologous expansion is definitively reached.

We also plot the scaled density $\rho_{\rm s}$ in the bottom plot of Fig. 4.
Similar to the behaviour of $v_{\rm r}$, we observe that the $\rho_{\rm s}$
match one another for a few different epochs when rescaled by $r_{\rm s}$.
This is expected in the case of homologous expansion when the (scaled) density
structure is frozen.

In Fig. 5, we show the fractional change in the absolute 3-D velocity of all
fluid elements relative to their asymptotic (late-time) velocities approach
zero as the envelope evolves. The fractional change is computed based on the
late-time velocity of each fluid element, defined as the average of the
velocity between $2000\,$d to $2500\,$d. The dotted line in Fig. 5 represents
the mean of the fractional change in the particle velocity. The shaded region
represents the standard deviation from the mean. Thus, the velocities of fluid
elements do not change by more than 5 per cent either in the magnitude or
direction after about $1500\,$d on average. This implies that the fluid
elements are on ballistic trajectories.

In addition, we plot the ratio of total thermal energy in the system to the
total kinetic energy as a function of time in Fig. 6. The thermal energy is
smaller than 5 per cent of the kinetic energy after 500 d. The fact that this
thermal energy does not continue to drop due to adiabatic expansion is because
we use a temperature floor in our simulations to maintain numerical stability.
In any case, thermal energy is a negligible fraction of the energy budget of
the system.

Finally, for a homologously expanding system under adiabatic conditions, the
average density scales like

$\displaystyle\bar{\rho}\propto t_{\rm h}^{-3}.$ (11)

This has also been previously demonstrated numerically by Iaconi et al. (2019)
for their SPH simulations. Here we confirm the same result by plotting the
average density of unbound particles (solid line) in the envelope as a
function of $t$ in Fig. 7. We also plot a $t^{-3}$ power law fit (dashed line)
that is fitted for $t\in[500,2000]$ d. The average density from our simulation
follows the $t^{-3}$ power law.

Figure 7: Mean density of the envelope as a function of $t$ (blue solid
curve). The black dotted curve is the $t^{-3}$ fit and is fitted for
$t\in[500,2000]$ d . Note that the average density profile follows a $t^{-3}$
power-law closely.

##  4 1-D Model

Figure 8: Scaled density, $\rho_{\rm s}$ as a function of scaled radius,
$r_{\rm s}$ from 1-D simulations, for $10$ (red), $20$ (blue), $50$ (magenta),
and $100$ (yellow) yr and a heating parameter of $\lambda=1$.  Figure 9:
Radial velocity, $v_{\rm r}$, as a function of scaled radius from 1-D
simulations. The times and $\lambda$ are the same as in Fig.8. The solid black
line is the best fit of $v_{\rm r}\propto r_{\rm s}^{0.95}$.

Motivated by the results of the previous section, we now study a simplified
1-D spherically-symmetric model of the ejected envelope. We have developed a
1-D finite-volume spherically-symmetric hydrodynamics code in Python that uses
an HLLE Riemann solver (Harten et al., 1983; Einfeldt, 1988) with piecewise-
constant (first-order) reconstruction to study this ejected envelope111This
1-D code will be shared on reasonable request to the corresponding author..
Our models consist of 350 grid points that are logarithmically spaced,
starting from $r=3\,R_{\odot}$ to $r=3\times 10^{7}\,R_{\odot}$, giving 50
grid points per decade. We set the origin to the centre of mass of the binary
and use free (von Neumann) boundary conditions on the inner and outer
boundaries. While the discussion of the hydrodynamic equations can be widely
found in the literature, we will briefly recap them here for completeness.
These equations can be written in compact notation by introducing a state
vector $\boldsymbol{\mathcal{U}}=(\rho,\rho v_{\rm r},\rho e)$

$\displaystyle\frac{\partial\boldsymbol{\mathcal{U}}}{\partial
t}+\frac{1}{r^{2}}\cdot\frac{\partial\boldsymbol{\mathcal{F}}}{\partial
r}=\boldsymbol{\mathcal{S}},$ (12)

where $\boldsymbol{\mathcal{F}}=(r^{2}\rho v_{\rm r},r^{2}(\rho v_{\rm
r}v_{\rm r}+P),r^{2}(\rho e+P)v_{\rm r})$ is the flux function,
$\boldsymbol{\mathcal{S}}=(0,-\rho\frac{GM_{\rm r}}{r^{2}},-\rho v_{\rm
r}\frac{GM_{\rm r}}{r^{2}}+\mathcal{S}_{\rm h})$ is the source function, $e$
is the specific energy, $G$ is the gravitational constant and $M_{\rm r}$
accounts for the mass of the central binary and the envelope enclosed within
radius $r$. The extra term, $\mathcal{S}_{\rm h}$, is added to study the
effect of the heating supplied to the envelope from the central binary. We
discuss this below. For the initial conditions of the 1-D model, we take the
fitted results from the 3-D numerical simulation at $t=800$ d.

The 1-D model computes $1000$ outputs over $100$ yr, each one separated by
$0.1$ yr. We present the results in Figs. 8 and 9 for $10$, $20$, $50$, and
$100$ yr. These times are relative to the start of the 3-D simulation so that
the same time between the two simulations can be directly compared. The scaled
radius and density follow equations (5) and (6), respectively. We calculate
the best fit of the linear part to get the power-law relation between $v_{\rm
r}$ and $r_{\rm s}$. We fit each time step separately and average them to
produce a best fit of $v_{\rm r}\propto r_{\rm s}^{0.95}$.

We note that the 1-D model is intentionally not constrained to adhere to the
homologous expansion model, though we did use initial conditions that
correspond to the beginning of the homologous phase. One feature that is
observed in these 1-D models, but not in the full 3-D models, is evident in
Fig. 9. Here at a radius below about a few $\times 10\,R_{\odot}$, the radial
velocity becomes negative. This manifests as the vertical rise in Fig. 9. This
is due to the gravitational potential from the inner binary being much larger
than the total energy of the envelope in this region. However, these negative
radial velocities are not seen in full 3-D models though the data is quite
noisy in this region (Fig. 4).

This difference may be attributed to the periodic forcing of the orbiting
binary on the gas in this inner region. To examine this, we develop a simple
model of binary heating for this 1-D model. The heating from the central
binary can be thought of as a periodic forcing from the binary driving a
damped simple harmonic oscillator with a frequency equal to the epicyclic
frequency, $\kappa=\Omega$, where $\Omega$ is the Keplerian orbital frequency.
A discussion of this heating term is in Appendix A, but the resulting
parameterized heating term is

$\displaystyle\mathcal{S}_{\rm h}=\lambda\Omega\frac{GM_{\rm
bin}\rho}{r}\left(\frac{a_{\rm bin}}{r}\right)^{2},$ (13)

where $\lambda$ is a free parameter, $M_{\rm bin}$ is the binary mass and
$a_{\rm bin}$ is the binary separation. Here, we set $M_{\rm
bin}=1.36\,M_{\odot}$, and $a_{\rm bin}=5\,R_{\odot}$.

Figure 10: Scaled density ($\rho_{\rm s}$) and radial velocity ($v_{\rm r}$)
as a function of scaled radius ($r_{\rm s}$), for 10 yr and different heating
parameters $\lambda$.

Fig. 10 shows $v_{\rm r}$ and $\rho_{\rm s}$ at $10$ yr for different heating
parameters $\lambda$. The heating term $\mathcal{S}_{\rm h}$ impacts the inner
envelope and has little impact in the outer regions ($r_{\rm
s}>10^{4}\,R_{\odot}$). For the case of no heating ($\lambda=0$, red line),
the inner regions follow a free-fall inflow solution. For larger heating
rates, this inflow is suppressed but not completely eliminated, though the
infall region occurs at substantially smaller radii than is effectively probed
by our 3-D simulations. In either case, the heating rate makes no impact on
either the density or velocity profile at large radii. We thus conclude that
binary heating or any other late-time heating has little effect on the
ultimate expansion and ejection of the envelope, but is required to prevent
the infall of the innermost envelope. In the 3-D simulations, angular momentum
or turbulence of the gas close to the binary may play a similar role, but this
is not well modelled in the 1-D simulation. In any case, it is evident that
there is extra physics that is not entirely accounted for in the 1-D model
that leads to homologous expansion (at smaller radii) in the 3-D model.

##  5 Discussion

In this work, we study the long-timescale evolution of CEE and the homologous
nature of envelope expansion. We simulate a common-envelope event using MANGA
for a $2\,M_{\odot}$ RGB and a $1\,M_{\odot}$ MS star binary system. We show
that nearly all gas from the red giant is unbound in 1500 d using an adiabatic
equation of state and relying only on the release of orbital energy.

This is in agreement with and in contrast to other work. For instance,
Chamandy et al. (2020) evolved a binary system of an AGB + white dwarf or MS
star system through 20 orbits using an adiabatic equation of state. In
agreement with our findings, they show that the envelope unbinds at a constant
rate and would become unbound in less than $10$ yr.

On the other hand, a number of others suggest that additional energy injection
is necessary. Ondratschek et al. (2022) studied a binary system consisting of
an AGB primary similar to Sand et al. (2020) with a white dwarf or a MS star
such that the mass ratio is 0.25 using the OPAL equation of state (Rogers &
Nayfonov, 2002). They find a complete envelope ejection in about $1000$ d when
considering thermal and ionization energy along with mechanical energy, and in
about $3400$ d when ignoring the thermal and ionization energy. Sand et al.
(2020) studied the fraction of unbound mass in two different simulations, one
using the ideal gas equation of state and the other using the OPAL equation of
state (Rogers & Nayfonov, 2002) for a binary system with an AGB primary and a
white dwarf or an MS companion star such that the mass ratio is 0.5. In the
case of the ideal gas equation of state, only $20$ per cent of the mass
becomes unbound and the rate of mass ejection is slower if the internal energy
is ignored. In their OPAL runs, $80$ per cent of the mass is unbound in about
$2500$ d considering only the mechanical energy and $100$ per cent is ejected
in about $1000$ d considering mechanical energy along with thermal and
recombination energy. Finally, Iaconi et al. (2019) find that the envelope of
a binary system with $0.88\,M_{\odot}$ RGB and a $0.6\,M_{\odot}$ MS companion
star is completely ejected in about $500$ d when considering mechanical and
recombination energies.

One crucial difference between our work and that of Sand et al. (2020) and
Ondratschek et al. (2022) is that the cores of the stars end up in a much
tighter binary in our case. In particular, Sand et al. (2020) starts the
binary at $236~{}R_{\odot}$ but ends up at $41~{}R_{\odot}$. In our case, we
start at around $50~{}R_{\odot}$, but end up at $5~{}R_{\odot}$. Hence, our
orbit shrinks by a factor of nearly $10$ whereas Sand et al. (2020)’s orbit
shrinks by a factor of about $5$. The corresponding relative gravitational
energy release is hence a factor of 2 greater in our case.

In addition, we also show that the envelope reaches homologous expansion
starting around a few hundred days ($550$ d). In comparison, Iaconi et al.
(2019) showed that it takes about $5000$ d for the bulk of the unbound gas to
become homologously expanding, even though the external layers of the envelope
become homologous as soon as they are ejected. This difference may be due to
the analysis methodology. Namely, Iaconi et al. (2019) traced the ballistic
trajectories of SPH particles whereas we looked at the velocities of fluid
elements and fit the velocities and associated radii to a $t_{\rm h}=0$ point.
We are also simulating these events at substantially higher resolution through
a combination of greater particle count and the use of a moving-mesh
methodology.

Inspired by the homologous expansion in our 3-D simulation, we also study a
1-D model for $t>800$ d. Unsurprisingly, we find that initializing the 1-D
simulation with the spherical approximation of the 3-D simulation data
produces homologous expansion in the bulk of the envelope. However, the inner
regions require some additional physics not modelled in the 1-D simulation to
preclude fallback and, hence, the breaking of homology. Here, we attribute
this additional physics to heating from the periodic forcing of the binary but
note that turbulence or angular momentum may play the same role.

The fact that simulated common-envelope events follow both a (roughly)
spherical and homologous expansion approximation is likely useful for their
theoretical and observational studies. First, CEE need not be numerically
simulated for extremely long timescales. Instead, they only need to be
simulated to the point where they begin homologous expansion, which occurs on
a timescale of years as opposed to decades. This will result in significant
computational savings and an associated expansion in the parameter space that
can be explored.

Secondly, the fact that they obey both the spherical and homologous expansion
approximations implies that radiation transfer codes that are used to
calculate supernova light curves and spectra can be adapted to compute light
curves and spectra from CEE events. Homologous expansion kinematic models are
widely accepted as a good first-order approximation to model radiation
transfer in supernova explosions (Röpke, 2005; Kerzendorf & Sim, 2014; Liu et
al., 2018). These radiation transfer codes typically assume a spherical and
homologous expansion profile for the ejecta to greatly speed up the
calculation. This fact has already been utilized by Kamiński et al. (2018) to
model the CO emission in V4332 Sgr as a homologously expanding bipolar flow.
Similarly Kamiński et al. (2020) show that the observed properties of the
molecular remnant of Nova 1670 (CK Vulpeculae) are reproduced by assuming
linear velocity fields.

Finally, we note that we have neglected effects, such as magnetic fields
(Ondratschek et al., 2022) and jets (see for instance Soker & Kaplan, 2021)
that could cause the outflow to become more non-spherical at late times. The
effect of jets in CEE is still unclear, as Zou et al. (2022) found that jets
are quickly choked within the envelope. On the other hand, Ondratschek et al.
(2022) showed that late-time jet-like outflows produce the bipolar morphology
seen in many planetary nebular systems. These effects will be a topic of
future work.

##  6 Conclusions

In this work, we have analyzed the nature of expansion and the timescale of
complete envelope ejection in common-envelope evolution. We used the moving-
mesh hydrodynamic solver MANGA to perform a long-timescale simulation of a CE
system involving a $2\,M_{\odot}$ red giant and a $1\,M_{\odot}$ main sequence
star. We let the system evolve for around $13$ yr. Starting at an orbital
radius of $52\,R_{\odot}$, the binary plateaus to an orbital radius of
$5\,R_{\odot}$ in $200$ d. We observe that nearly all envelope material is
unbound after 1500 d, and $80$ per cent is unbound in 400 d. This timescale is
similar to that found by others who also studied envelope ejection in low-mass
binary systems. However, we find that there is no need for additional energy
injection from recombination.

Motivated by previous work by Iaconi et al. (2019), we also show that the
envelope enters a phase of homologous expansion after 550 d. This is likely
important for theoretical and observational work on CEE. First, one can have
significant computational savings by doing numerical simulations up to the
homologous start time and using the homologous expansion model afterwards.
Secondly, the radiative transfer codes used for finding light curves and
spectra for supernovae can be adapted for use in CEE simulations.

Finally, we study the homologous expansion model in 1-D simulations using the
power-law fits of the homologous phase as initial conditions. From this study,
we found that periodic heating from the binary star at late times can affect
the inner regions of the envelope but does not impact the homologous
expansion.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## Acknowledgements

VV acknowledges support from the NSF through grants PHY-1912649 and
PHY-2207728. SvB, LP, and PC acknowledge support from the NASA ATP program
through NASA grant NNH17ZDA001N-ATP, the NSF through grant AST-2108269, and
the UWM Research Assistance Fund. LP also acknowledges support from the UWM R1
Distinguished Dissertator Award, NSF through grant PHY-1748958 and a grant
from the Simons Foundation (216179, LB). MANGA simulations were completed on
the Mortimer HPC System at UWM, which was funded by the NSF Campus
Cyberinfrastructure Award OAC-2126229 and UWM. We also use the yt software
platform for the analysis of the data and generation of plots in this work
(Turk et al., 2011).

## References

  * Chamandy et al. (2020) Chamandy L., Blackman E. G., Frank A., Carroll-Nellenback J., Tu Y., 2020, MNRAS, 495, 4028
  * Chang & Etienne (2020) Chang P., Etienne Z. B., 2020, MNRAS, 496, 206
  * Chang et al. (2017) Chang P., Wadsley J., Quinn T. R., 2017, MNRAS, 471, 3577
  * Chang et al. (2020) Chang P., Davis S. W., Jiang Y.-F., 2020, MNRAS, 493, 5397
  * Einfeldt (1988) Einfeldt B., 1988, SIAM Journal on Numerical Analysis, 25, 294
  * Harten et al. (1983) Harten A., Lax P. D., Leer B. v., 1983, SIAM Review, 25, 35
  * Iaconi et al. (2019) Iaconi R., Maeda K., De Marco O., Nozawa T., Reichardt T., 2019, MNRAS, 489, 3334
  * Ivanova et al. (2013) Ivanova N., et al., 2013, A&ARv, 21, 59
  * Jermyn et al. (2022) Jermyn A. S., et al., 2022, arXiv e-prints, p. arXiv:2208.03651
  * Kamiński et al. (2018) Kamiński T., Steffen W., Tylenda R., Young K. H., Patel N. A., Menten K. M., 2018, A&A, 617, A129
  * Kamiński et al. (2020) Kamiński T., Menten K. M., Tylenda R., Wong K. T., Belloche A., Mehner A., Schmidt M. R., Patel N. A., 2020, A&A, 644, A59
  * Kerzendorf & Sim (2014) Kerzendorf W. E., Sim S. A., 2014, MNRAS, 440, 387
  * Liu et al. (2018) Liu L.-D., Zhang B., Wang L.-J., Dai Z.-G., 2018, ApJ, 868, L24
  * Ondratschek et al. (2022) Ondratschek P. A., Röpke F. K., Schneider F. R. N., Fendt C., Sand C., Ohlmann S. T., Pakmor R., Springel V., 2022, A&A, 660, L8
  * Paczynski (1976) Paczynski B., 1976, in Eggleton P., Mitton S., Whelan J., eds, Vol. 73, Structure and Evolution of Close Binary Systems. p. 75
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton B., et al., 2019, ApJS, 243, 10
  * Prust (2020) Prust L. J., 2020, MNRAS, 494, 4616
  * Prust & Chang (2019) Prust L. J., Chang P., 2019, MNRAS, 486, 5809
  * Rogers & Nayfonov (2002) Rogers F. J., Nayfonov A., 2002, ApJ, 576, 1064
  * Röpke (2005) Röpke F. K., 2005, A&A, 432, 969
  * Sand et al. (2020) Sand C., Ohlmann S. T., Schneider F. R. N., Pakmor R., Röpke F. K., 2020, A&A, 644, A60
  * Soker & Kaplan (2021) Soker N., Kaplan N., 2021, Research in Astronomy and Astrophysics, 21, 090
  * Spaulding & Chang (2021) Spaulding A., Chang P., 2021, MNRAS, 501, 1748
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, The Astrophysical Journal Supplement Series, 192, 9
  * Zou et al. (2022) Zou Y., Chamandy L., Carroll-Nellenback J., Blackman E. G., Frank A., 2022, MNRAS, 514, 3041

##  Appendix A Simple Harmonic oscillator model for heating from central
binary

For a fluid element in a circumbinary orbit around two point masses, we model
the periodic forcing from the inner binary as a forced simple harmonic
oscillator with a natural resonant frequency of $\kappa$, which is the
epicyclic frequency. The periodic forcing term has a frequency of
$m(\Omega_{\rm b}-\Omega)$, where $m$ is some natural number, $\Omega_{\rm b}$
is the orbital frequency of the binary and $\Omega$ is the orbital frequency
of the fluid element.

A general driven damped harmonic oscillator equation is then

$\displaystyle\frac{d^{2}\delta r}{dt^{2}}+\Gamma\frac{d\delta
r}{dt}+\kappa^{2}\delta r=f_{0}\cos{m(\Omega_{\rm b}-\Omega)t},$ (14)

where $\Gamma$ is a damping term that arises from fluid dissipation and
$f_{0}$ is the overall forcing amplitude. Letting $\omega=m(\Omega_{\rm
b}-\Omega)$, the complex form of this equation is

$\displaystyle\frac{d^{2}\delta r}{dt^{2}}+\Gamma\frac{d\delta
r}{dt}+\kappa^{2}\delta r=\frac{F_{0}}{m}\exp(i\omega t),$ (15)

Starting with an ansatz, $\delta=\mathcal{A}\exp{i\omega t}$, we can show

$\displaystyle\mathcal{A}=f_{0}\frac{\kappa^{2}-\omega^{2}-i\omega\Gamma}{(\kappa^{2}-\omega^{2})^{2}+(\Gamma\omega)^{2}},$
(16)

where we can write

$\displaystyle\mathcal{A}=A\exp{i\delta}$ (17)

where

$\displaystyle
A=\frac{f_{0}}{\sqrt{(\kappa^{2}-\omega^{2})^{2}+(\Gamma\omega)^{2}}}$ (18)
$\displaystyle\delta=\tan^{-1}\left(-\frac{\omega\Gamma}{\kappa^{2}-\omega^{2}}\right)$
(19)

The specific energy is then

$\displaystyle\epsilon=\frac{1}{2}\dot{\delta
r}^{2}+\frac{1}{2}\kappa^{2}{\delta r}^{2}$ (20)

The maximum amplitude, $A$, occurs when $\omega=\kappa$ and is:

$\displaystyle
A=\frac{f_{0}}{\Gamma\kappa},\qquad\textrm{and}\qquad\delta=\frac{\pi}{2},$
(21)

and the maximum specific energy is:

$\displaystyle\epsilon_{\rm max}=\frac{1}{2}\frac{(f_{0})^{2}}{\Gamma^{2}}$
(22)

The specific energy dissipation is thus,

$\displaystyle\dot{\epsilon}=2\epsilon_{\rm
max}\Gamma=\frac{f_{0}^{2}}{\Gamma}$ (23)

For our case, we can model the forcing from the inner binary as a tidal
forcing

$\displaystyle f_{0}=\beta\frac{GM_{\rm bin}a_{\rm bin}}{r^{3}},$ (24)

where $\beta$ is some constant, $M_{\rm bin}$ is the mass of the binary, $r$
is the orbital radius of the fluid element, and $a_{\rm bin}$ is the
separation of the two stars. Likewise, we can set the damping rate as some
fraction of the epicyclic frequency $\Gamma=\alpha\kappa$. The volumetric
heating rate is then

$\displaystyle\mathcal{S}_{h}=\frac{2\beta^{2}}{\alpha}\frac{\kappa GM_{\rm
bin}\rho}{r}\left(\frac{a_{\rm bin}}{r}\right)^{2}.$ (25)

Here we define the relevant constant as $\lambda\equiv 2\beta^{2}/\alpha$.
Combining these terms together and setting $\kappa=\Omega$ in Keplerian
potentials, we arrive at equation (13).