# A Theoretical Analysis of Granulometry-based Roughness Measures on Cartosat
DEMs

K Nagajothi, Sravan Danda, Aditya Challa, and B S Daya Sagar K. Nagajothi is
with Regional Remote Sensing Centre, Indian Space Research Organisation,
Bangalore 560037, India (e-mail: nagajothi_k@nrsc.gov.in).Sravan Danda is with
Computer Science and Information Systems, APPCAIR, BITS-Pilani Goa Campus,
Zuari Nagar, 403726, India (sravan8809@gmail.com)Aditya Challa is with
Computer Science and Automation, Indian Institute of Science, Bangalore,
560012, India (aditya.challa.20@gmail.com)B. S. D. Sagar is with Systems
Science and Informatics Unit, Indian Statistical Institute-Bangalore Centre,
8th Mile, Mysore Road, RV College PO, Bangalore-560059, India
<EMAIL_ADDRESS>
###### Abstract

The study of water bodies such as rivers is an important problem in the remote
sensing community. A meaningful set of quantitative features reflecting the
geophysical properties help us better understand the formation and evolution
of rivers. Typically, river sub-basins are analysed using Cartosat Digital
Elevation Models (DEMs), obtained at regular time epochs. One of the useful
geophysical features of a river sub-basin is that of a roughness measure on
DEMs. However, to the best of our knowledge, there is not much literature
available on theoretical analysis of roughness measures. In this article, we
revisit the roughness measure on DEM data adapted from multiscale
granulometries in mathematical morphology, namely multiscale directional
granulometric index (MDGI). This measure was classically used to obtain shape-
size analysis in greyscale images. In earlier works, MDGIs were introduced to
capture the characteristic surficial roughness of a river sub-basin along
specific directions. Also, MDGIs can be efficiently computed and are known to
be useful features for classification of river sub-basins. In this article, we
provide a theoretical analysis of a MDGI. In particular, we characterize non-
trivial sufficient conditions on the structure of DEMs under which MDGIs are
invariant. These properties are illustrated with some fictitious DEMs. We also
provide connections to a discrete derivative of volume of a DEM. Based on
these connections, we provide intuition as to why a MDGI is considered a
roughness measure. Further, we experimentally illustrate on Lower-Indus,
Wardha, and Barmer river sub-basins that the proposed features capture the
characteristics of the river sub-basin.

###### Index Terms:

Digital Elevation Model, Cartosat, Granulometric Index, Mathematical
Morphology.

##  I Introduction

The study of geophysical properties of rivers is an important problem in the
remote sensing community. A study of a river sub-basin at regular time epochs
over a large span of time helps understand the evolution of the river. The
evolution of river sub-basins provides information required to prioritize the
rivers that need immediate attention for conservation/ identify natural
calamities etc. However, such a study is highly dependent on extracting
meaningful geophysical features of the river sub-basins. For example, the
complexity of the surficial roughness of a river sub-basin provides
information as to what the dominant wind directions are, in that region.

Recall that Cartosat Digital Elevation Models (DEMs) are typically used to
compute geophysical features of river sub-basins. In literature, several
studies indicate that surficial roughness is an important characteristic of a
river sub-basin [13, 2, 8, 7, 10, 9, 11, 1, 12]. However, these studies lack a
detailed theoretical analysis of the roughness measures proposed. In this
article, we analyse in detail, a surficial roughness measure that was proposed
in [5] i.e. multiscale directional granulometric index (MDGI), a special case
of a more general measure namely a multiscale granulometric index.

Multiscale granulometric index was originally proposed in [3] to obtain a
shape-size analysis of objects in greyscale images. As greyscale images can be
viewed as digital surfaces with the greyscale intensity at each pixel
representing the height of the surface, these measures have been adapted to
DEMs [12]. It was shown experimentally that such an adaptation is indeed
useful from an application point-of-view i.e. to classify river sub-basins
[5]. A natural question would then be to ask: Can we characterize the
equivalence classes of DEMs obtained by the equivalence relation - two DEMs
are equivalent if their MDGI are identical? In other words, can we find
necessary and sufficient conditions on the structure of a DEM under which a
MDGI is invariant? In this article, we partially answer this question and
provide theoretical insights on how a MDGI varies with the structure of a DEM.

In particular, the contributions of this article are as follows:

  1. 1.

We provide an alternate visualization of the definition of a MDGI proposed in
[5].

  2. 2.

Using the alternate visualization, we characterize non-trivial sufficient
conditions on DEMs under which a MDGI is invariant. The invariance properties
are intuitively explained and illustrated on fictitious DEMs.

  3. 3.

We analyse the relation between a MDGI and a discrete derivative of volume of
a DEM. Using this analysis, we provide an intuition as to why a MDGI is
considered a roughness measure.

  4. 4.

Further, a preliminary application of MDGI is shown on real data i.e.
Cartosat-1 DEM data of Lower-Indus, Wardha, and Barmer river sub-basins to
show that these measures capture characteristics of the sub-basin.

The rest of the article is organized as follows: In section II, we provide the
definitions of basic morphological operators and multiscale granulometric
index. Also, the existing literature on the usage of directional granulometric
indices is briefly described. Section III contains the core contributions of
the article i.e. an alternate visualization of a MDGI, the invariance
properties, relation to discrete derivative on the volume of a DEM, and an
intuition as to why a MDGI is considered a roughness measure. Section IV
contains experiments on real data i.e. on watersheds of Indus, Wardha and
Barmer river sub-basins.

##  II Multiscale Granulometric Index

In this section, we recall the formal definitions of a multiscale
granulometric index and briefly describe the existing literature. First, we
start with the basic definitions.

###  II-A Elementary Morphological Operators

###### Definition 1

Let $A\subset\mathbb{Z}^{2}$ be a finite set. A digital elevation model (DEM)
of a river basin/sub-basin is represented as a function $f:A\rightarrow H$
where $H\subset\mathbb{Z}^{+}$ is a finite set.

Figure 1: A non-rectangular fictitious DEM. Each circle represents a physical
square area. The values inside each of the circles are the elevations.

Each point $a\in A$ represents a small physical, square area and $f(a)$
represents the discretized average height of the physical area. Also, observe
that $A$ is possibly non-rectangular i.e. $A=\cup_{i=n_{1}}^{n_{2}}A_{i}$
where $A_{i}=\\{(i,j):m_{1,i}\leq j\leq m_{2,i}\\}$). See Fig 1 for an
illustration. This is in contrast to greyscale images where $A$ is always
rectangular i.e. $m_{1,i}$ and $m_{2,i}$ are independent of $i$. The number of
distinct elements in $H$ are comparable to the number of grey levels in a
greyscale image. A higher cardinality of $H$ indicates a finer resolution in
the elevations and is analogous to a finer spectral resolution in greyscale
images.

Next, we need the definition of a structuring element. Using the notion of a
structuring element, dilation and erosion, the fundamental blocks of roughness
measures based on multiscale granulometric index are then defined. Note that
we restrict the definition of structuring element i.e. assume that the
structuring element contains its origin and is symmetric. This definition
suffices for the purposes of this article.

###### Definition 2

A structuring element $SE\subset\mathbb{Z}^{2}$ is a finite set such that: (1)
$(0,0)\in SE$, (2) $(i,j)\in SE\Rightarrow(-i,-j)\in SE$.

The different types of structuring elements used in this article are given by
$B_{1}=\\{(-1,1),(0,0),(1,-1)\\}$, $B_{2}=\\{(0,1),(0,0),(0,-1)\\}$,
$B_{3}=\\{(-1,-1),(0,0),(1,1)\\}$, $B_{4}=\\{(-1,0),(0,0),(1,0)\\}$, and
$B=\\{(x,y)\in\mathbb{Z}^{2}:-1\leq x,y\leq 1\\}$. Fig 2 provides a pictorial
representation of these structuring elements. Observe that each of the
structuring elements $B_{1},B_{2},B_{3},B_{4}$ are $3$ units long and are
effectively one-dimensional.

Figure 2: Left to right: Types of structuring elements used in this article.
The centre element refers to $(0,0)$. A value of $1$ indicates that the point
corresponding to the coordinates are present in the structuring element and a
value of $0$ otherwise.

Recall that a greyscale dilation and a greyscale erosion are defined as
follows:

###### Definition 3

Let $f:A\rightarrow H$ be a DEM and let $SE$ be a structuring element, then a
dilation of $f$ by $SE$ is given by

$\left[f\oplus{SE}\right](x,y)=max_{(s,t)\in{SE}}\\{f(x+s,y+t)\\}$ (1)

where $SE$ is a structuring element

###### Definition 4

Let $f:A\rightarrow H$ be a DEM and let $SE$ be a structuring element, then an
erosion of $f$ by $SE$ is given by

$\left[f\ominus{SE}\right](x,y)=min_{(s,t)\in{SE}}\\{f(x+s,y+t)\\}$ (2)

where $SE$ is a structuring element

Next, we need the definition of a morphological opening and a multiscale
morphological opening.

###### Definition 5

Let $f:A\rightarrow H$ be a DEM and let $SE$ be a structuring element, then an
opening of $f$ is given by

$\left[f\circ{SE}\right](x,y)=\left[\left[f\ominus{SE}\right]\oplus{SE}\right](x,y)$
(3)

###### Definition 6

Let $f:A\rightarrow H$ be a DEM and let $SE$ be a structuring element, then a
multiscale opening of $f$ is given by

$\left[f\circ n{SE}\right](x,y)$ (4)

where $nSE=SE\oplus SE\oplus\cdots\oplus SE$ with the number of dilations in
the telescoping expression being $n-1$.

###  II-B Directional Multiscale Granulometric Index

Before we provide a formal definition of a multiscale granulometric index, we
need to define the notion of volume of a DEM.

###### Definition 7

Let $f:A\rightarrow H$ be a DEM. The volume of $f$, $V(f)$ is defined as
follows:

$V(f)=\sum_{a\in A}f(a)$ (5)

Intuitively, the volume of a DEM captures the physical volume of a DEM on and
above the altitude chosen to be zero. For example, the volume of DEM shown in
Fig 1 is $46$. It is easy to see that an application of a multiscale opening
results in a DEM with lower volume as $n$ increases. Also, it is easy to see
that there exists $N_{0}\in\mathbb{N}$ such that $V(f\circ
n{SE})=V(f\circ(n+1){SE})\ \forall n\geq N_{0}$. Recall that the definition of
multiscale granulometric index is given by

###### Definition 8

$GI_{SE}(f)=-\sum_{n=0}^{\infty}p_{n}log(p_{n})$ (6)

where

$p_{n}=\frac{V(f\circ n{SE})-V(f\circ(n+1){SE})}{V(f)}$ (7)

Note that the existence of $N_{0}\in\mathbb{N}$ such that $V(f\circ
n{SE})=V(f\circ(n+1){SE})\ \forall n\geq N_{0}$ ensures that the summation is
finite. The terms inside the summation for $n\geq N_{0}$ have to be
interpreted as zero. When the structuring element $SE$ is chosen to be one of
$B_{1},B_{2},B_{3},B_{4}$, the obtained multiscale granulometric index is said
to be a directional multiscale granulometric index or MDGI. Intuitively, this
makes sense as each of $B_{1},B_{2},B_{3},B_{4}$ are linear and indicate four
primary directions.

###  II-C Existing Literature

Multiscale granulometric index was first introduced in [3] to perform a shape-
size analysis of objects in greyscale images. Then it was used to analyse
textures in greyscale images [6]. Later, these ideas were generalized to
analyse soil section image analysis [14]. Multiscale granulometric index was
theoretically analyzed from the perspective of identifying shapes and sizes of
objects in greyscale images. The utility of granulometries in greyscale images
led to the development of efficient algorithms for specialized classes of
structuring elements [15, 4].

Very recently, these ideas were adapted to DEMs. It was experimentally shown
in [5] that multiscale granulometric indices obtained using specific
structuring elements retain characteristic information of the river basins. A
natural question would then be to ask: Can we find necessary and sufficient
conditions on the structure of a DEM under which the directional granulometric
index is invariant? In the next section, we partially answer this question by
identifying non-trivial sufficient conditions on DEMs such that all DEMs
satisfying such conditions have the same directional granulometric index.

##  III Theoretical Analysis of Directional Granulometric Indices

In this section, we analyse the MDGIs from a theoretical perspective. Firstly,
we recall some modified definitions from graph theory to suit the purposes of
subsequent analysis on a MDGI. Secondly, we provide an alternate way to view a
MDGI using graphs. Then, by building on this visualization of MDGI, we provide
intuition on sufficient conditions under which DEMs have the MDGI. Then, we
prove the main result of this article formalizing the intuition i.e.
characterization of non-trivial sufficient conditions on structure of a DEM
such that the MDGI is invariant. This is followed by a short subsection
analysing the relation between MDGI and a discrete derivative of volume of a
DEM. Using this analysis, we provide intuition as to why a MDGI is considered
a roughness measure.

###  III-A Some Modified Graph Definitions

###### Definition 9

$\mathcal{G}=(V,E,W)$ is said to be a node-weighted graph if $V$ denotes the
set of nodes is a finite set, $E\subset\\{\\{v_{i},v_{j}\\}:v_{i}\neq v_{j},\
v_{i},v_{j}\in V\\}$ denotes the set of edges, and $W:V\rightarrow H$ is a
non-negative integer-valued function on $V$ such that $H\subset\mathbb{Z}^{+}$
is a finite set.

Figure 3: A node-weighted graph constructed on the DEM given by Fig 1. A
$4$-adjacency relation is used to construct the set of edges.

A node-weighted graph can be used to model DEMs taking into account the
spatial relations of neighbouring physical areas on which the elevations are
stored. For example, Fig 3 shows a node-weighted graph constructed on a
fictitious DEM illustrated by Fig 1.

###### Definition 10

Let $\mathcal{G}=(V,E,W)$ be a node-weighted graph. Let $W:V\rightarrow H$ and
$h\in H$. $\mathcal{G}_{\geq h}=(V_{\geq h},E_{\geq h},W|_{V_{\geq h}})$ is
said to be an upper-thresholded subgraph of $\mathcal{G}=(V,E,W)$ at elevation
$h$, where $V_{\geq h}=\\{v\in V:W(v)\geq h\\}$, $E_{\geq
h}=\\{\\{v_{i},v_{j}\\}:\\{v_{i},v_{j}\\}\in E\ and\ W(v_{i})\geq
h,W(v_{j})\geq h\\}$, and $W|_{V_{\geq h}}$ denotes the restriction of the
function $W$ to $V_{\geq h}$

Figure 4: Upper-thresholded graph constructed on the node-weighted graph given
by Fig 3. The threshold is set at elevation $3$.

Intuitively, an upper-thresholded subgraph of a node-weighted graph
constructed on a DEM provides an abstraction of the sub-structure of the the
DEM that is above an elevation level. For example, Fig 4 illustrates the
upper-thresholded subgraph at elevation $3$ on the node-weighted graph given
by Fig 3.

###### Definition 11

Let $\mathcal{G}=(V,E,W)$ be a node-weighted graph. A subset of nodes
$V_{1}\subset V$ is said to be connected if for every pair of nodes
$v_{s},v_{t}\in V_{1}$, there exists a sequence of nodes
$<v_{s}=v_{0},v_{1},\cdots,v_{r-1},v_{r}=v_{t}>$ such that
$\\{v_{i},v_{i+1}\\}\in E$ for every $0\leq i\leq r-1$. A subset of nodes
$V_{1}\subset V$ is said to be maximally connected if (1) $V_{1}$ is
connected, and (2) $V_{1}\subset V_{2}\subset V$ and $V_{2}$ is connected
implies $V_{2}=V_{1}$.

We remark that given any node-weighted graph $\mathcal{G}=(V,E,W)$, the set
$V$ can be decomposed uniquely as a disjoint union of maximally connected
subsets of $V$. For example, the node-weighted graph in Fig 3 has one
maximally connected subset which is the vertex set itself. Similarly, the
upper-thresholded graph in Fig 4 which is also a node-weighted graph has two
maximally connected subsets of the vertex set.

###  III-B Another Interpretation of a Directional Granulometric Index

Recall from subsection II-B that a multiscale granulometric index is given by
Def 8 (see Eq 6 and Eq 7). Intuitively, a multiscale granulometric index
measures the entropy of the volume loss on the series of morphological
openings with increasing sizes of the structuring element.

Assume that the structuring element $SE$ is given by one of
$B_{1},B_{2},B_{3},B_{4}$ as defined in subsection II-A. Each of the four
structuring elements are effectively one-dimensional. Hence, a directional
granulometric index effectively measures volume loss on linear scans (but in
different directions). Fig 5 shows an illustration of the scans obtained for
SEs $B_{4}$ and $B_{3}$.

Figure 5: One dimensional scans on the DEM in Fig 1 required to obtain a
theoretical analysis of $GI_{B_{4}}(f)$ and $GI_{B_{3}}(f)$.

Thus, same theoretical analysis on MDGI holds for each of
$B_{1},B_{2},B_{3},B_{4}$. Let $f:A\rightarrow H$ be a DEM. In order to
understand the directional granulometric index better, we first try to analyze
a one-dimensional DEM i.e. working with DEMs restricting the domain to
horizontal scans. The analysis for a generic two-dimensional set $A$ would be
a straightforward extension with slightly involved notation. Mathematically,
such a restriction would be equivalent to working with sets of type:

$A_{i_{0}}=\\{(i_{0},j):m_{1}\leq j\leq m_{2}\\}\subset A$ (8)

for a fixed $i_{0}\in\mathbb{Z}$ and $m_{1}<m_{2}\in\mathbb{Z}$. In general,
$m_{1},m_{2}$ depend on $i_{0}$ as $A$ is not necessarily rectangular.
However, we blur this detail to work with a simplified notation.

We now analyse a MDGI obtained by multiscale openings using horizontal linear
structuring elements $\\{L_{n}:n\in\mathbb{Z}^{+}\\}$, where $L_{n}$ denotes a
horizontal structuring element with $n$ consecutive $1$s. A similar analysis
holds for $\\{nB_{4}:n\in\mathbb{Z}^{+}\\}$ because
$\\{nB_{4}:n\in\mathbb{Z}^{+}\\}\subsetneq\\{L_{n}:n\in\mathbb{Z}^{+}\\}$ (in
particular $nB_{4}=L_{2n+1}$ for each $n\in\mathbb{Z}^{+}$). We are now ready
to examine $GI_{L_{1}}(f|_{A_{i_{0}}})$ given by Eq 6 and Eq 7 where
$f:A\rightarrow H$ is a DEM. Let
$\mathcal{G}^{f}=(A_{i_{0}},E_{chain},f|_{A_{i_{0}}})$ denote a node-weighted
graph with

$E_{chain}=\\{\\{(i_{0},j),(i_{0},j+1)\\}:m_{1}\leq j\leq m_{2}-1\\}$ (9)

Consider the sequence of upper-thresholded subgraphs of the node-weighted
graph $\mathcal{G}^{f}=(A_{i_{0}},E_{chain},f|_{A_{i_{0}}})$ at all possible
elevations i.e.

$\\{\mathcal{G}^{f}_{\geq h}\\}_{h=min(H)}^{max(H)}=\\{(V^{f}_{\geq
h},E^{f}_{\geq h},f|_{V^{f}_{\geq h}})\\}_{h=min(H)}^{max(H)}$ (10)

Let $V^{f}_{\geq h}=\cup_{r=1}^{n_{h}}V^{f,r}_{\geq h}$ denote the disjoint
union of maximally connected subsets for each $h\in[min(H),max(H)]$. Denote
$n_{t,h}$ as

$n_{t,h}=|\\{V^{f,r}_{\geq h}:|V^{f,r}_{\geq h}|=t\\}|$ (11)

Here $n_{t,h}$ denotes the number of maximally connected subsets of
$V^{f}_{\geq h}$ which are exactly $t$ units long. It is easy to see that the
probabilities given by Eq 7 satisfy

$p_{k}\varpropto k\sum_{h=min(H)}^{max(H)}n_{k,h},$ (12)

for each $k\in\mathbb{Z}^{+}$. This is because $L_{k}=kL_{1}$ is $k$ units
long for each $k\in\mathbb{Z}^{+}$. An opening with $kL_{1}$ removes any
maximally connected subset of length less than $k$ units. Hence, probability
$p_{k}$ is proportional to the volume obtained by slices of rectangular blocks
that are $k$ units long on the DEM $f$.

We are now ready to extend these ideas to a generic two-dimensional set $A$.
In the two-dimensional case $n_{t,h}$ for each horizontal scan given by Eq 11
would be dependent on $i_{0}$ i.e. the choice of row, denoted by
$n_{t,h}^{(i)}$ . Assuming $f:A\rightarrow H$ is the DEM on which we wish to
compute the MDGI, Eq 12 would transform to:

$p_{k}\varpropto\sum_{i=n_{1}}^{n_{2}}k\sum_{h=min(H)}^{max(H)}n_{k,h}^{(i)}$
(13)

###  III-C Invariances of Directional Granulometric Indices

Recall from Sec II that a multiscale granulometric index is given by Def 8
(see Eq 6 and Eq 7). We are interested to characterize sufficient conditions
on the structure of a DEM such that all DEMs that satisfy those conditions
have the same MDGI. Mathematically, we need to find non-trivial collections of
DEMs $\mathcal{F}_{c}=\\{f|GI_{B_{4}}(f)=c\\}$ where $c>0$ is a positive
constant. A sufficient condition for a MDGI to be invariant is that the
probabilities given by Eq 7 remain the same. On a closer look at Eq 12, it is
easy to see that if $n_{t,h}$ given by Eq 11 remains constant for each $t,h$
then the MDGI for each such DEM is the same. We now state the result formally:

###### Theorem 1

Let $f_{1}:A\rightarrow H$ and $f_{2}:A\rightarrow H$ be distinct DEMs. If the
number of maximally connected subsets $n_{t,h}^{(i)}$ (given by Eq 11) of
upper-thresholded subgraphs (given by Eq 10) for every row $i$, every length
$t$, and every elevation $h$ are identical for both $f_{1}$ and $f_{2}$ then
$GI_{L_{1}}(f_{1})=GI_{L_{1}}(f_{2})$.

The proof follows from Eq 13, Eq 6 and Eq 7. To see that the sufficient
conditions imposed in Theorem 1 are non-trivial, consider $f:A\rightarrow H$.
We will now construct a ‘large’ collection of DEMs different from $f$ whose
MDGI w.r.t. $B_{4}$ is identical to that of $f$. Let
$A=\cup_{i=n_{1}}^{n_{2}}A_{i}$ where $A_{i}=\\{(i,j):m_{1,i}\leq j\leq
m_{2,i}\\}$. Define $\hat{f}_{i}:A_{i}\rightarrow H$ as

$\hat{f}_{i}(i,j)=f(i,m_{2,i}-j+m_{1,i})$ (14)

for each $m_{1,i}\leq j\leq m_{2,i}$. Intuitively, $\hat{f}_{i}$ is the
mirror-reflection of $f|_{A_{i}}$. Now, consider the collection
$\mathcal{F}_{Reflection}(f)=\\{g:A\rightarrow H:g|_{A_{i}}=f|_{A_{i}}\ or\
\hat{f}_{i}\\}$. It is easy to see that
$|\mathcal{F}_{Reflection}(f)|=2^{n_{2}-n_{1}+1}$. One of the elements of
$\mathcal{F}_{Reflection}(f)$ is $f$. Hence, we could construct
$2^{n_{2}-n_{1}+1}-1$ different DEMs with the same MDGI. Fig 6 provides an
illustration of this construction i.e. two DEMs different from the DEM
provided by Fig 1 with same MDGI $GI_{B_{4}}(.)$.

Figure 6: Two DEMs different from the DEM provided by Fig 1 with the same MDGI
$GI_{B_{4}}(.)$.

In general, it is possible to construct much larger set of DEMs with same
MDGI. This is because the domain of DEMs are usually different i.e. arbitrary
shaped DEMs are encountered in practice. Also, the conditions characterized by
Eq 14 are relatively more restrictive sufficient conditions compared to the
sufficient conditions provided by Theorem 1. For example, among the DEMs
illustrated in Fig 7, one DEM cannot be constructed from the other using the
construction provided by Eq 14. However, both these DEMs have identical
$n_{t,h}$ for each $t,h$ and hence have the same MDGI $GI_{B_{4}}(.)$.
Further, the conditions provided by Theorem 1 are sufficient but not
necessary.

Figure 7: The top row corresponds to one DEM and the second row corresponds to
another DEM. These DEMs have identical $GI_{B_{4}}(.)$ but they cannot be
constructed from each other using the construction provided by Eq 14.

Next, we discuss another type of invariance which involves transformations on
the elevation level i.e. comparing DEMs $f_{1}:A\rightarrow H_{1}$ and
$f_{2}:A\rightarrow H_{2}$ where $H_{1}$ and $H_{2}$ are finite subsets of
$\mathbb{Z}^{+}$. Formally, we have the following result:

###### Theorem 2

Let $A\subset\mathbb{Z}^{+}$ and $H\subset\mathbb{Z}^{+}$ be finite sets.
Assume that $T:H\rightarrow\mathbb{Z}^{+}$ is a scaling transformation i.e.
$T(x)=kx$ for some arbitrary but fixed $k\in\mathbb{Z}^{+}$. The DEMs
$f:A\rightarrow H$ and $(T\cdot f):A\rightarrow T(H)$ have identical MDGIs
i.e. $GI_{B_{l}}(f)=GI_{B_{l}}(T\cdot f)$ for each $1\leq l\leq 4$. Here
$\cdot$ denotes composition of functions.

The proof of Theorem 2 follows from the fact that each of the terms on right
side of Eq 13 is scaled up by the same constant and hence has no effect on the
LHS of the same equation.

###  III-D Relation to Discrete Derivative of Volume of a DEM

In this subsection, we relate a discrete derivative of the volume of a DEM to
the MDGI. Firstly, for the sake of simplified notation, assume that we are
working on horizontal slices $A_{i}$ of the domain $A$ of the DEM
$f:A\rightarrow H$, i.e. subsets of the type Eq 8 of $A$. Define

$\Phi_{h_{0}}(f|_{A_{i}})=\sum_{h=h_{0}}^{max(H)}\sum_{v\in V^{f,r}_{\geq
h}\subset V^{f}_{\geq h}}W(v)$ (15)

The quantity $\Phi_{h_{0}}(f|_{A_{i}})$ denotes the volume of DEM $f$ on and
above elevation $h_{0}$ on the horizontal slice $A_{i}$. In particular,
$\sum_{i=n_{1}}^{n_{2}}\Phi_{min(H)}(f|_{A_{i}})$ denotes the total volume of
DEM $f$, $V(f)$ given by Eq 5. $\Phi_{h_{0}}(f|_{A_{i}})$ can be rewritten as

$\Phi_{h_{0}}(f)=\sum_{h=h_{0}}^{max(H)}\sum_{t}tn_{t,h}$ (16)

The discrete derivative of the volume of a DEM w.r.t. elevation is hence given
by

$\Phi_{h_{0}}(f)-\Phi_{h_{0}+1}(f)=\sum_{t}tn_{t,h_{0}}$ (17)

Notice that this expression can be interpreted as the sum of areas of
maximally connected subsets of $V^{f}_{\geq h_{0}}$ on the slice $A_{i}$. In
general, when we consider a two-dimensional domain $A$ of the DEM $f$, Eq 17
transforms to:

$\Phi_{h_{0}}(f)-\Phi_{h_{0}+1}(f)=\sum_{i=n_{1}}^{n_{2}}\sum_{t}tn_{t,h_{0}}^{(i)}$
(18)

###  III-E Why is a Multiscale Directional Granulometric Index considered a
Roughness Measure?

In this subsection, we provide an intuition as to why a MDGI is regarded as a
roughness measure. To accomplish this, we consider a special class of DEMs
given by:

$\mathcal{F}_{UniPeak}=\\{f:A\rightarrow H:\sum_{t}n_{t,h}=1\forall\ h\in
f(A)\\},$ (19)

where $A$ denotes a one-dimensional set i.e. $\exists m_{1}\leq
m_{2}\in\mathbb{Z}$ such that $A=\\{(i_{0},j)|m_{1}\leq j\leq m_{2}\\}$. It is
easy to see that in such a case the set of discrete derivatives given by Eq 17
would effectively be a permutation of the probabilities given by Eq 12. This
means that the entropy calculated on the successive discrete derivatives of
volume of a DEM is identical to the MDGI when computed on a one-dimensional
uni-peak DEM.

##  IV Experiments

In this section, we provide empirical evidence to show that MDGIs capture the
characteristic features of a river sub-basin.

###  IV-A Study Area and Data Used

We consider lower Indus sub-basin (fluvial), Wardha sub-basin (floodplain),
and Barmer sub-basin (desert) for the experiments. Lower Indus sub-basin, one
of the $14$ sub-basins lies in between the geographical coordinates of
$73^{\circ}11’$ to $76^{\circ}44’$ East longitudes and $34^{\circ}42^{\prime}$
to $36^{\circ}9^{\prime}$ North latitudes, is divided into $31$ watersheds of
sizes ranging between $319$ sq.km and $1270$ sq.km. Wardha sub-basin, one of
the principal tributaries of Godavari river, is situated in between the
geographical coordinates of $19^{\circ}18’N$ and $21^{\circ}58’N$ latitudes,
and $77^{\circ}20’E$ and $79^{\circ}45’E$ longitudes. This sub-basin has $69$
watersheds. Barmer is another sub-basin of Indus Basin situated between
$69^{\circ}48’$ and $71^{\circ}43’$ East longitudes, and between
$25^{\circ}28’$ to $27^{\circ}69’$ North latitudes. It is fully under Thar
Desert and is divided into $38$ watersheds. Cartosat DEMs of the Lower Indus
sub-basins, Wardha sub-basins, and Barmer sub-basins are illustrated in Fig 8.

Figure 8: (a), (e), and (i) Cartosat DEMs of Lower Indus, Wardha, and Barmer
sub-basins. The delineations highlight distinct watersheds within each of the
sub-basins. (b), (f), and (j) Watersheds classified based on normalized
multiscale directional granulometric indices. The MDGIs are all scaled down
within each sub-basin and are color-coded as per the ranges mentioned in the
legend. The ranges are arbitrarily chosen. (c), (g), and (k) high directional
granulometric indices with the colors and the texture highlighting the
corresponding SE for which the MDGI is highest, and (d), (h), and (l) low
directional granulometric indices with the colors and the texture highlighting
the corresponding SE for which the MDGI is lowest, of $31$ watersheds of the
Lower Indus, $69$ watersheds of Wardha and $38$ watersheds of Barmer sub-
basins respectively.

###  IV-B Some Preliminary Observations

Figure 9: (a), (c), and (e): The histograms of the highest order-statistic
among normalized HDGIs obtained with $B_{1},B_{2},B_{3},B_{4}$ as the
structuring elements on Lower-Indus, Wardha, and Barmer sub-basins
respectively. (b), (d), and (f): The histograms of the lowest order-statistic
among normalized HDGIs obtained with $B_{1},B_{2},B_{3},B_{4}$ as the
structuring elements on Lower-Indus, Wardha, and Barmer sub-basins
respectively.

We compute the MDGIs of all watersheds in each of the river sub-basins using
the structuring elements $B_{1},B_{2},B_{3},B_{4}$. As there is a variability
in the size of each of the watersheds i.e. the domain of each DEM is of
different cardinalities, we normalize these MDGIs with the multiscale
granulometric index obtained by using $B$ as the structuring element. As a
first observation, for each river, we check the order-statistics of the
normalized MDGIs. Fig 9 shows a plot of the histograms of the maximum and
minimum of the normalized MDGIs among the four primary directions across the
river sub-basins. These histograms can be viewed as empirical probability
distributions. The maximum (respectively minimum) among the normalized MDGIs
is denoted as high directional (respectively low directional) in Figs 8 and 9.
It is easy to see that the Barmer sub-basins show a different pattern in the
order-statistics as illustrated in Fig 9. In particular, it is often the case
that the maximum is along the direction of $B_{3}$ and the lowest is given by
$B_{2}$ which is not the case with the other two river sub-basins. However,
these order-statistics do not help in identifying the differences between
Wardha and Lower-Indus sub-basins.

###  IV-C More Observations

To identify the differences between the watersheds of all three river sub-
basins, we construct features $X[0],\cdots,X[15]$ based on the normalized
MDGIs. The details of the construction are as follows:

Let $GI_{B_{i}}(f)/GI_{B}(f)=Z_{i}$ for $i=1,2,3,4$

$X[{(i-1)4+(j-1)}]=\begin{cases}&Z_{i}\ if\ Z_{(j)}=Z_{i}\\\ &0\
otherwise\end{cases}$ (20)

for each $i=1,2,3,4$ and $j=1,2,3,4$, where $Z_{(j)}$ denotes the $j^{th}$
order-statistic among $Z_{1},Z_{2},Z_{3},Z_{4}$ i.e. $Z_{(1)}\leq Z_{(2)}\leq
Z_{(3)}\leq Z_{(4)}$ form a permutation of $Z_{1},Z_{2},Z_{3},Z_{4}$. As the
number of watersheds ($138$ altogether) is small, we do not split the data
into training and test sets. Instead, we try to obtain interpretable rules
that can classify the watersheds into appropriate sub-basins based on the
constructed features. A simple decision tree with a depth of $2$ is
constructed. The depth of the decision tree is restricted so as to ensure that
we do not over fit the data. Also, the reason for choosing $2$ as the depth is
that the minimum depth of a decision tree required to classify $3$ classes is
$2$. We observe that such a decision tree (see Fig 10) is capable of obtaining
$\approx 71\%$ accuracy. Note that a random classifier on the other hand can
obtain a maximum accuracy of $69/138=50\%$. This is because the data is class-
unbalanced with $31$ Indus, $69$ Wardha and $38$ Barmer watersheds.

Figure 10: A decision tree with a depth of $2$ is built using the features
$X[0],\cdots,X[15]$ based on the MDGIs as described in the text. This decision
tree is the smallest possible classification tree that is capable of
separating three classes. This classification tree obtains a classification
accuracy of $\approx 71\%$ while a random classifier can obtain a maximum
classification accuracy of $50\%$ on classifying $138$ watersheds.

This experiment indicates that the features computed from MDGIs carry
characteristic information of the sub-basins. To explore the discernibility of
MDGI-based features, we built decision trees of larger depths until a maximum
depth of $9$. We observed that the accuracies of decision trees of depths $5$,
$6$, and $9$ are $\approx 86\%$, $\approx 89\%$, and $\approx 94\%$
respectively. This indicates that MDGI-based features are useful for river
sub-basin classification. However, these features cannot be used as standalone
features to classify sub-basins. In order to build better classifiers, one
needs to incorporate domain knowledge on the river sub-basins through some
form of remotely-sensed data or otherwise.

##  V Conclusions and Perspectives

In this article, we revisit the roughness measure on DEM data adapted from
multiscale granulometries in mathematical morphology, namely multiscale
directional granulometric index (MDGI). In earlier works, MDGIs were
introduced to capture the characteristic surficial roughness of a river sub-
basin along specific directions. They are known to be useful features for
classification of river sub-basins. In this article, we provided a theoretical
analysis of a MDGI. In particular, we characterized non-trivial sufficient
conditions on the structure of DEMs under which MDGIs are invariant. These
properties are illustrated with some fictitious DEMs. We also provided
connections to a discrete derivative of volume of a DEM. Based on these
connections, we provided intuition as to why a MDGI is considered a roughness
measure. Further, we experimentally illustrated on Lower-Indus, Wardha, and
Barmer river sub-basins that the proposed features capture the characteristics
of the river sub-basin.

Building on the ideas from this article, one can explore at least two
directions: 1) building on the main theorem, one can investigate more
sufficient conditions ultimately trying to characterize sufficient and
necessary conditions on the structure of a DEM such that MDGI is invariant, 2)
on the experimental side, use the features proposed in the article alongside
other features on river sub-basins to build better classifiers.

## Acknowledgment

Nagajothi Kannan would like to thank Regional Remote Sensing Centre, Indian
Space Research Organisation. Sravan Danda would like to acknowledge the
funding received from BPGC/RIG/2020-21/11-2020/01 (Research Initiation Grant
provided by BITS-Pilani K K Birla Goa Campus) and thank APPCAIR, and Computer
Science and Information Systems, BITS-Pilani Goa. Aditya Challa would like to
thank Indian Institute of Science (IISc) for the Raman Post Doctoral
fellowship. The work of B. S. D. Sagar was supported by the DST-ITPAR-Phase-IV
project under the Grant number INT/Italy/ITPAR-IV/Telecommunication/2018.

## References

  * [1] L Chockalingam and BS Daya Sagar.  Morphometry of network and nonnetwork space of basins.  Journal of Geophysical Research: Solid Earth, 110(B8), 2005.
  * [2] Igor V Florinsky.  Derivation of topographic variables from a digital elevation model given by a spheroidal trapezoidal grid.  International Journal of Geographical Information Science, 12(8):829–852, 1998.
  * [3] Petros Maragos.  Pattern spectrum and multiscale shape representation.  IEEE Transactions on pattern analysis and machine intelligence, 11(7):701–716, 1989.
  * [4] Vincent Morard, Petr Dokládal, and Etienne Decenciere.  One-dimensional openings, granulometries and component trees in $\mathcal{O}$(1) in per pixel.  IEEE Journal of selected topics in signal processing, 6(7):840–848, 2012.
  * [5] Kannan Nagajothi and BS Daya Sagar.  Classification of geophysical basins derived from srtm and cartosat dems via directional granulometries.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12(12):5259–5267, 2019.
  * [6] Mark Rzadca.  Multivariate granulometry and its application to texture segementation.  1994\.
  * [7] Behara Seshadri Daya Sagar.  Mathematical morphology in geomorphology and GISci.  CRC Press, 2013.
  * [8] BS Daya Sagar.  Universal scaling laws in surface water bodies and their zones of influence.  Water resources research, 43(2), 2007.
  * [9] BS Daya Sagar and L Chockalingam.  Fractal dimension of non-network space of a catchment basin.  Geophysical research letters, 31(12), 2004.
  * [10] BS Daya Sagar and Tay Lea Tien.  Allometric power-law relationships in a hortonian fractal digital elevation model.  Geophysical research letters, 31(6), 2004.
  * [11] Lea Tien Tay, BS Daya Sagar, and Hean Teik Chuah.  Allometric relationships between traveltime channel networks, convex hulls, and convexity measures.  Water resources research, 42(6), 2006.
  * [12] LT Tay, BS Daya Sagar, and HT Chuah.  Granulometric analyses of basin-wise dems: a comparative study.  International Journal of Remote Sensing, 28(15):3363–3378, 2007\.
  * [13] Donald L Turcotte.  Fractals in geology and geophysics.  Pure and applied Geophysics, 131(1):171–196, 1989.
  * [14] Costas S Tzafestas and Petros Maragos.  Shape connectivity: multiscale analysis and application to generalized granulometries.  Journal of Mathematical Imaging and Vision, 17(2):109–129, 2002\.
  * [15] Luc Vincent.  Granulometries and opening trees.  Fundamenta Informaticae, 41(1, 2):57–90, 2000.

|  Nagajothi Kannan (M’18) received the B.Sc. degree in physics from
Bharathiar University, Coimbatore, India, in 1986, and the master’s degree in
computer applications from Indira Gandhi National Open University, New Delhi,
India, in 2002. He is currently a Senior Scientist with the Indian Space
Research Organisation, Bangalore, India. He has more than 20 years of work
experience in the area of digital image processing, geographical information
system, software development, and system management. He has worked on
geospatial application projects of national missions and large-scale user
projects that involved design and development of large database, system
integration, and customized solutions. His research interests include digital
image processing, mathematical morphology, GISci and geospatial applications,
and decision support systems. Mr. Nagajothi is a member of IEEE GRSS.
---|---
|  Sravan Danda received the B.Math.(Hons.) degree in Mathematics from the
Indian Statistical Institute - Bangalore, and the M.Stat. degree in
Mathematical Statistics from the Indian Statistical Institute - Kolkata, in
2009, and 2011, respectively. From 2011 to 2013, he worked as a Business
Analyst at Genpact - Retail Analytics, Bangalore. He completed his PhD in
computer science from Systems Science and Informatics Unit, Indian Statistical
Institute - Bangalore in 2019 under the joint supervision of B.S.Daya Sagar
and Laurent Najman. He is currently working as an Assistant Professor at
Department of Computer Science and Information Systems, BITS Pilani K K Birla
Goa Campus. His current research interests are Discrete Mathematical
Morphology and Discrete Optimization in Machine Learning.
---|---
|  Aditya Challa received the B.Math.(Hons.) degree in Mathematics from the
Indian Statistical Institute - Bangalore, and Masters in Complex Systems from
University of Warwick, UK - in 2010, and 2012, respectively. From 2012 to
2014, he worked as a Business Analyst at Tata Consultancy Services, Bangalore.
He completed his PhD in computer science from Systems Science and Informatics
Unit, Indian Statistical Institute - Bangalore in 2019. He is currently a
Raman PostDoc Fellow at Indian Institute of Science, Bangalore. His current
research interests focus on using techniques from Mathematical Morphology in
Machine Learning.
---|---
|  B S Daya Sagar (M’03-SM’03) is a Full Professor of the Systems Science and
Informatics Unit (SSIU) at the Indian Statistical Institute. Sagar received
his MSc and Ph.D. degrees in Geoengineering and Remote Sensing from the
Faculty of Engineering, Andhra University, Visakhapatnam, India, in 1991 and
1994 respectively. He is also the first Head of the SSIU. Earlier, he worked
in the College of Engineering, Andhra University, and Centre for Remote
Imaging Sensing and Processing (CRISP), The National University of Singapore
in various positions during 1992-2001. He served as Associate Professor and
Researcher in the Faculty of Engineering & Technology (FET), Multimedia
University, Malaysia, during 2001-2007. Sagar has made significant
contributions to the field of geosciences, with special emphasis on the
development of spatial algorithms meant for geo-pattern retrieval, analysis,
reasoning, modeling, and visualization by using concepts of mathematical
morphology and fractal geometry. He has published over 85 papers in journals
and has authored and/or guest-edited 11 books and/or special theme issues for
journals. He recently authored a book entitled ‘Mathematical Morphology in
Geomorphology and GISci,’ CRC Press: Boca Raton, 2013, p. 546. He recently co-
edited two special issues on ”Filtering and Segmentation with Mathematical
Morphology” for IEEE Journal of Selected Topics in Signal Processing (v. 6,
no. 7, p. 737-886, 2012), and ”Applied Earth Observation and Remote Sensing in
India” for IEEE Journal of Selected Topics in Applied Earth Observation and
Remote Sensing (v. 10, no. 12, p. 5149-5328, 2017). His recent book “Handbook
of Mathematical Geosciences”, Springer Publishers, p. 942, 2018 reached 750000
downloads. He was elected as a member of the New York Academy of Sciences in
1995, as a Fellow of the Royal Geographical Society in 2000, as a Senior
Member of the IEEE Geoscience and Remote Sensing Society in 2003, as a Fellow
of the Indian Geophysical Union in 2011. He is also a member of the American
Geophysical Union since 2004, and a life member of the International
Association for Mathematical Geosciences (IAMG). He delivered the ”Curzon & Co
- Seshachalam Lecture - 2009” at Sarada Ranganathan Endowment Lectures
(SRELS), Bangalore, and the ”Frank Harary Endowment Lecture - 2019” at
International Conference on Discrete Mathematics - 2019 (ICDM - 2019). He was
awarded the ’Dr. Balakrishna Memorial Award’ of the Andhra Pradesh Academy of
Sciences in 1995, the Krishnan Medal of the Indian Geophysical Union in 2002,
the ’Georges Matheron Award - 2011 with Lectureship’ of the IAMG, and the
Award of IAMG Certificate of Appreciation - 2018. He is the Founding Chairman
of the Bangalore Section IEEE GRSS Chapter. He has been recently appointed as
an IEEE Geoscience and Remote Sensing Society (GRSS) Distinguished Lecturer
(DL) for a two-year period from 2020 to 2022. He is on the Editorial Boards of
Computers & Geosciences, Frontiers: Environmental Informatics, and
Mathematical Geosciences. He is also the Editor-In-Chief of the Springer
Publishers’ Encyclopedia of Mathematical Geosciences.
---|---


# Efficient Search in Graph Edit Distance: Metric Search Trees vs. Brute Force
Verification

Wenqi Marshall Guo
<EMAIL_ADDRESS>Jeffrey Uhlmann
<EMAIL_ADDRESS>
(March 2024)

Technical Report for Multidisciplinary Undergraduate Research Conference
(MURC) 2024 111For more information, please access
https://github.com/weathon/GED-CMT.

##  1 Introduction

Graph similarity search (GSS) plays a key role in domains such as molecular
and protein similarity search[1]. Specifically, Graph Edit Distance (GED)–the
measure of similarity based on the steps required to edit one graph into
another–is one of the most commonly used metrics in GSS. However, the
computation of exact GED traditionally requires exponential time, which can
slow down the search process considerably[1].

Previous researchers have proposed various methods to speed up the GED-based
GSS process. One of them is a brute-force-based method that first uses lower
bounds to filter the graphs in the database (only the graphs with lower bounds
less than the threshold will be kept) and then verification in which the
computation of exact GED is replaced by a procedure to verify if the GED of
two graphs is lower than a given threshold[1].

Another promising approach involves augmenting similarity searches with metric
search trees, such as the Cascading Metric Tree (CMT)[3]. CMT has been shown
to enhance similarity searches using different metrics, like Euclidean
distance[3] and Kendall-Tau distances[2]. In this study, we aim to investigate
whether adapting the CMT for use with GED using upper-and-lower bounds (UBLB)
can outperform the simple brute force verification relayed in [1]. Our results
suggest that, contrary to our initial hypothesis, CMT does not routinely
outperform brute force verification strategy in most contexts.

##  2 Methodology

As mentioned previously, exact GED is extremely hard to compute, this leads to
the use of UBLB in the querying process of the algorithm. However, since the
tree will be pre-built only once, to enable more effective pruning, the exact
distance calculation is still employed during the tree’s construction phase.

Our query algorithm could be described as follows:

  * •

Two empty answer sets will be prepared, ’confirmed set’ and ’suspected set’

  * •

The UBLB information is used to prune the search tree with a similar process
in [3].

  * •

If the upper bonds of the largest distance in a subtree are less than the
threshold, all the graphs of the subtree will be added to the confirmed set.

  * •

When it reaches the leaf nodes:

    * –

If the upper bounds are lower than the threshold, the corresponding graph is
added to the confirmed set.

    * –

If the lower bounds are lower than the threshold, the graph is added to the
suspected set.

  * •

A brute force verification step then examines each member of the suspected set
to filter out any false positives from the UBLB search.

We tested our implementation on Euclidean distance to debug and ensure it is
working as intended. Then, the algorithm was tested on graph data derived from
PubChem and the performance of CMT was compared to that of brute force
verification.

##  3 Results and Discussion

Results indicate that CMT does not routinely outperform brute force
verification. In fact, in many cases, CMT significantly lags behind brute
force in terms of speed. One possible explanation for this could be even for
UBLB of the GED, the computation is still expensive. On the other hand, the
verification could be done very fast with a small query radius (which is
usually the case for practical applications). This surprising finding
underlines the need for continuous innovation in methodologies to speed up GED
computation.

##  4 Future Work and Limitations

Our study is a preliminary investigation of the relative performance between
tree search strategies like CMT and brute force verification and it has
several limitations which indicate possible directions for future work. These
include:

There are different hyper-parameters that will affect the result such as the
iterations of the upper and lower bounds algorithm. Different hyper-parameters
might yield different results.

We only tested the algorithm in a very limited search space. The dataset size
could largely limit the performance of CMT.

We only used the graph data from PubChem, there might be a bias in the dataset
and the result could be different on other datasets.

The structure of an efficient search tree typically leads to better relative
performance with larger datasets, a trend partially confirmed by our study. At
the same time, our research also suggests that smaller graphs may deliver
better relative performance. However, this is a “blank space” in our dataset:
we do not have many small molecules. Examining this space might be valuable to
find potential space where CMT is more efficient.

We did not optimize our implementation of the CMT, by optimizing it the
performance could be improved.

##  5 Conclusion

This study investigated the performance of Cascading Metric Trees in GED
computation for graph similarity search in comparison to brute force
verification. Surprisingly, our results showed that despite its proven
efficacy in other metrics, CMT did not routinely outperform brute force
verification in computing GED. These findings highlight the need for continued
exploration of methods for speeding up GED computation and verification, with
a particular focus on optimizing the computation of GED upper and lower bounds
and investigating the performance of these methods on a wider range of
datasets and search spaces.

##  6 ACKNOWLEDGEMENT

Part of the computation for this work was performed on the high-performance
computing infrastructure provided by Research Support Solutions and in part by
the National Science Foundation under grant number CNS-1429294 at the
University of Missouri, Columbia MO. DOI:
https://doi.org/10.32469/10355/69802https://doi.org/10.32469/10355/69802

OpenAI ChatGPT is used as an aid in this project such as brainstorming,
language refining, etc. This tech report is produced by ChatGPT based on an
author-written outline and proofread and largely edited by the author.

##  7 References

[1] L. Chang, X. Feng, X. Lin, L. Qin, W. Zhang, and D. Ouyang, ‘Speeding up
ged verification for graph similarity search’, in 2020 IEEE 36th International
Conference on Data Engineering (ICDE), Apr. 2020, pp. 793–804. doi:
10.1109/ICDE48307.2020.00074.

[2] W. Guo and J. Uhlmann, ‘Metric search for rank list compatibility matching
with applications’. arXiv, Aug. 10, 2023. doi: 10.48550/arXiv.2303.11174.

[3] J. Uhlmann and M. R. Zuniga, ‘The cascading metric tree’. arXiv, Dec. 20,
2021. doi: 10.48550/arXiv.2112.10900.

[4] D. B. Blumenthal, S. Bougleux, J. Gamper, and L. Brun, ‘Gedlib: a c++
library for graph edit distance computation’, in Graph-Based Representations
in Pattern Recognition, D. Conte, J.-Y. Ramel, and P. Foggia, Eds., in Lecture
Notes in Computer Science. Cham: Springer International Publishing, 2019, pp.
14–24. doi: 10.1007/978-3-030-20081-7_2.

[5] Z. Abu-Aisheh, R. Raveaux, and J.-Y. Ramel, ‘A graph database repository
and performance evaluation metrics for graph edit distance’, in Graph-Based
Representations in Pattern Recognition, C.-L. Liu, B. Luo, W. G. Kropatsch,
and J. Cheng, Eds., in Lecture Notes in Computer Science. Cham: Springer
International Publishing, 2015, pp. 138–147. doi:
10.1007/978-3-319-18224-7_14.

[6] K. Riesen, A. Fischer, and H. Bunke, ‘Computing upper and lower bounds of
graph edit distance in cubic time’, in Artificial Neural Networks in Pattern
Recognition, N. El Gayar, F. Schwenker, and C. Suen, Eds., Cham: Springer
International Publishing, 2014, pp. 129–140. doi:
10.1007/978-3-319-11656-3_12.


# More Evidence for Variable Helium Absorption from HD 189733b

Michael Zhang Department of Astronomy, California Institute of Technology,
Pasadena, CA 91125, USA  P. Wilson Cauley Laboratory of Atmospheric and Space
Physics, University of Colorado Boulder, Boulder, CO 80309  Heather A. Knutson
Division of Geological and Planetary Sciences, California Institute of
Technology  Kevin France Laboratory of Atmospheric and Space Physics,
University of Colorado Boulder, Boulder, CO 80309  Laura Kreidberg Max Planck
Institute for Astronomy, Heidelberg, Germany  Antonija Oklopčić Anton
Pannekoek Institute of Astronomy, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, Netherlands  Seth Redfield Astronomy Department and Van
Vleck Observatory, Wesleyan University, Middletown, CT 06459, USA  Evgenya L.
Shkolnik School of Earth and Space Exploration, Arizona State University,
Tempe, AZ 85287, USA

###### Abstract

We present a new Keck/NIRSPEC observation of metastable helium absorption from
the upper atmosphere of HD 189733b, a hot Jupiter orbiting a nearby moderately
active star. We measure an average helium transit depth of $0.420\pm 0.013$%
integrated over the [-20, 20] km/s velocity range. Comparing this measurement
to eight previously published transit observations with different instruments,
we find that our depth is 32% (9$\sigma$) lower than the average of the three
CARMENES transits, but only 16% (4.4$\sigma$) lower than the average of the
five GIANO transits. We perform 1D hydrodynamical simulations of the outflow,
and find that XUV variability on the order of 33%–common for this star–can
change the helium absorption depth by up to 60%, although a more typical
change is 15%. We conclude that changes in stellar XUV flux can explain the
observational variability in helium absorption, but that variability in the
stellar He line cannot be excluded. 3D models are necessary to explore other
sources of variability, such as shear instability and changing stellar wind
conditions.

##  1 Introduction

HD 189733b is currently one of the most extensively studied exoplanets in the
literature. Its exceptional observational favorability has made it a popular
target for atmospheric characterization studies. It is a large (1.1 $R_{J}$)
planet orbiting a somewhat small K dwarf (0.78 $R_{\sun}$) on a relatively
tight orbit (P=2.2 d), and is located just 20 pc from Earth (Addison et al.,
2019; Rosenthal et al., 2021). Since its discovery in 2005 (Bouchy et al.,
2005), it has been observed in wavelengths ranging from the X-ray to the mid-
infrared, using a wide variety of observational approaches. HD 189733b is
currently the only exoplanet with a possible X-ray transit detection
(Poppenhaeger et al., 2013), and was the first planet with a measured infrared
phase curve (Knutson et al., 2007). It has an extensively characterized
transmission and emission spectrum (e.g. Sing et al. 2011; McCullough et al.
2014; Morello et al. 2014), which provide us with some of the most precise
constraints on its atmospheric composition (e.g. Zhang et al. 2020; Harrington
et al. 2022). This planet has been the focus of a diverse array of theoretical
studies seeking to explain its observed properties, including its atmospheric
circulation patterns (e.g. Showman et al. 2009) and the nature of its
atmospheric aerosols (e.g. Steinrueck et al. 2021).

HD 189733b is also well suited for mass loss studies, for many of the same
reasons that make it favorable for observations of the lower atmosphere.
Atmospheric mass loss is an important phenomenon that could reshape exoplanet
demographics, potentially turning mini-Neptunes into smaller and denser super-
Earths, as well as creating the Fulton gap (Fulton et al., 2017; Fulton &
Petigura, 2018) and the Neptunian desert (e.g. Szabó & Kiss 2011; Fulton &
Petigura 2018). However, there are many theoretical uncertainties in mass loss
models that make exoplanet atmospheric evolution difficult to model, including
the effects of magnetic fields (Adams, 2011) and of interaction with the
stellar wind. The 10833 Å line has long been proposed as a possible probe of
exoplanet atmospheres (Seager & Sasselov, 2000; Oklopčić & Hirata, 2018), but
observations have only recently succeeded in detecting exoplanetary absorption
in this line (Spake et al., 2018; Nortmann et al., 2018; Allart et al., 2019).
It has been suggested that K-type stars have the most suitable stellar
spectrum for populating the metastable helium state (Oklopčić, 2019;
Poppenhaeger, 2022). HD 189733b, aside from being large, nearby, and highly
irradiated, orbits a K-type star.

In 2016–2017, Salz et al. (2018) observed three transits of HD 189733b with
CARMENES and detected clear He I absorption in all three, with consistent
depths. Guilluy et al. (2020) observed five transits with GIANO and obtained
inconsistent depths indicating that HD 189733b’s outflow may be varying in
time. In this paper, we present a new transit observation with Keck/NIRSPEC.
Section 2 describes the data acquisition, Section 3 the data analysis, Section
4 compares our results to prior work, Section 5 uses a 1D hydrodynamic model
to explore how much the varying high-energy stellar flux should change the
signal, and Section 6 summarizes our conclusions.

##  2 Observations

We observed a transit of HD 189733b with Keck/NIRSPEC on July 14, 2020, from
08:28 UTC to 14:41 UTC, using the NIRSPEC-1 filter and the 0.288x12$"$ slit.
These observations cover a wavelength range of 0.947–1.121 µm at a resolution
of R=37,500. Over the course of these 6 hours, we took 252 exposures in a ABBA
nod pattern, with each exposure taking slightly over a minute and each ABBA
pattern taking roughly 5 minutes. These data cover 1.0 h of pre-transit
baseline, the 1.8 h transit, and 3.4 h of post-transit baseline. The airmass
decreased from 1.2 at the beginning of the night to 1.0 around 0.5 h after
mid-transit, before rising again to 1.66 at the end of the observations. The
resulting SNR per spectral pixel in the continuum surrounding the He I line
ranged from 150–280, with a SNR peak 2 hours after transit, and an unfortunate
SNR trough coinciding with mid-transit.

Two of the raw spectra have multiple traces due to telescope nodding errors:
one during egress, and one an hour after the end of egress. We exclude these
spectra and their nod companions from the analysis, for a total of four
excluded reduced spectra.

The NIRSPEC observation was intended to be simultaneous with a Ly$\alpha$
observation by the Hubble Space Telescope’s COS spectrograph (GO 15710).
Unfortunately, HST failed to acquire the guide star, so we did not collect any
Ly$\alpha$ data.

##  3 Analysis

To analyze these observations, we use the methods outlined in Zhang et al.
(2022a). We make a master dark by stacking 5 darks with an exposure time of
1.5 s and 5 coadds. We make a master flat by subtracting the master dark from
each flat and stacking the 29 flats. We create four difference images out of
each $A_{1}B_{1}B_{2}A_{2}$ nod group: $A_{1}-B_{1}$, $B_{1}-A_{1}$,
$B_{2}-A_{2}$, $A_{2}-B_{2}$. This subtracts off bias, dark current and
skyglow. We then use a custom variant of optimal extraction to extract the 1D
spectrum, while masking both bad pixels identified in the dark and flat
processing stages, and cosmic rays identified by the optimal extraction
algorithm itself. Afterwards, we generate a template spectrum by combining a
PHOENIX model (Husser et al., 2013) with a telluric model, accounting for the
radial velocity of the star relative to Earth. We use this template spectrum
to obtain a wavelength solution for order 70 (containing the helium line) in
each individual exposure and then run molecfit to correct for tellurics. This
last step is not strictly necessary for these data because–due to a fortuitous
combination of the star’s radial velocity and Earth’s orbital velocity–there
are no telluric lines (with the possible exception of microtellurics) that
overlap with the helium line.

Having obtained wavelength calibrated and telluric corrected spectra, we
interpolate all spectra onto a common wavelength grid with a uniform
logarithmic spacing of $\lambda/110,000$ and a range of 10,810–10,850 Å. The
resolution of 110,000 is chosen so that the spacing of the grid is equal to
that of the spectral pixels at 10833 Å. We remove fringing by applying a notch
filter twice, using the same parameters as in Zhang et al. (2021). We divide
each spectrum by the continuum, take the logarithm of the entire spectral
grid, and subtract the mean of every row and column from that row and column.
The end result is a $N_{obs}\times N_{wav}$ grid of numbers representing the
relative deviation of a pixel from the mean for that row and column.

After obtaining this residuals grid, for every column (wavelength), we
subtract the mean of the out-of-transit part of the residuals image for that
column; we then invert the residuals image. The resulting residuals image now
shows the excess absorption relative to the out-of-transit baseline. However,
there are still continuum variations that contribute structure in this image.
We correct for these variations by masking out the strong lines (including the
helium line), fitting a 3rd order polynomial to each row (epoch) with respect
to wavelength, and subtracting off the polynomial. The result is shown in
Figure 1. Finally, we shift all of the in-transit residuals spectra to the
planetary rest frame and stack them to obtain the 1D excess absorption
spectrum shown in Figure 2.

Figure 1: Excess absorption in percent, as a function of time and wavelength.
The diagonal red lines represent the wavelengths of the helium lines. The
horizontal white light represent the beginning of ingress and end of egress in
white light. Figure 2: Average in-transit excess absorption spectrum from
beginning of ingress to end of egress, in the planetary frame. The red lines
represent the wavelengths of the three helium lines. Figure 3: Light curve of
the region [-20, 15] km/s from the main helium peak, in the stellar frame. The
beginning of ingress and end of egress are marked by the vertical black dotted
and solid lines. Note the stellar activity that follows the transit.

##  4 Comparison to Published Observations

Figure 4: Comparison of the average depth of the helium line that we find (in
green) to literature values, measured by CARMENES and GIANO. All averages are
computed between the end of ingress and beginning of egress in time, and
between -20 and 20 km/s in planet-centric wavelength space. Figure 5: Excess
absorption spectra obtained by the three instruments averaged over 3 transits
(CARMENES), 5 transits (GIANO-B), and 1 transit (NIRSPEC). CARMENES and
GIANO-B spectra are convolved to the NIRSPEC resolution of 37,500, resulting
in the artificial smoothness of the CARMENES spectrum. These spectra are
computed from the end of ingress to beginning of ingress, unlike Figure 2,
which is computed from the beginning of ingress to end of egress. CARMENES and
GIANO-B spectra were kindly provided by the authors of the respective papers.

We compare our results to the eight literature values for HD 189733b’s helium
absorption signal. Salz et al. (2018) observed 3 transits with the fibre-fed
Echelle spectrograph CARMENES (R $\sim$ 80,400), situated on the 3.5m Calar
Alto telescope. Guilluy et al. (2020) obtained 5 transits with the slit-fed
Echelle spectrograph GIANO-B (R $\sim$ 50,000) on the 3.6m Telescopio
Nazionale Galileo telescope. The CARMENES observations have a SNR of 160–240
per spectral pixel in the continuum surrounding the helium line (Salz et al.,
2018), but their exposure time is closer to five minutes than one minute, and
their wavelength dispersion is 70% higher. Binning our observations to their
time resolution and binning their observations to our spectral resolution, we
find that our SNR is 1.6–2.0x times their single-night SNRs. The same
comparison, performed for the GIANO data, reveals that our SNR is 3.6–4.4x
their single-night SNRs. The GIANO comparison is not perfect, as Guilluy et
al. (2020) report their SNRs over the whole order containing the helium line,
not the immediate region surrounding the line.

Guilluy et al. (2020) report the depths for each individual transit, which
they calculate by averaging from the end of ingress to the beginning of egress
for a wavelength range of -20 to 20 km/s in planet-centric velocity space. We
use their definition to calculate the analogous average transit depths for the
CARMENES transits (using values taken from their Fig. A.3) and for our NIRSPEC
transit. We did not take into account the differing resolutions of the
instruments when defining the equivalent pixel range for our wavelength
bounds. We quantified our sensitivity to the instrumental resolution by
convolving the CARMENES spectrum down to the NIRSPEC resolution of 37,500 and
the NIRSPEC pixel sampling before calculating the average transit depth, and
found that it did not change the result by more than 0.01% for any of the
three transits.

We compare the resulting transit depths for each observation in Figure 4. This
figure shows that the consistency of the three CARMENES transits over a year
appears to have been a coincidence, and is not representative of the broader
sample. The GIANO transits show widely varying depths, while our NIRSPEC depth
is lower than all but one of the previous transit observations. Combining the
depths for every instrument with a weighted average, we find an average
transit depth of $0.617\pm 0.017$% for CARMENES, $0.508\pm 0.015$% for GIANO,
and $0.420\pm 0.013$% for NIRSPEC. The error bars are the standard deviation
of the mean, not the sample. The CARMENES average differs from the NIRSPEC
result by more than 9$\sigma$, and from the GIANO result by 4.4$\sigma$.

###  4.1 Effect of Stellar Variability on Observed He Transit Depths

We consider two ways in which stellar variability can affect the observed
transit depths. First, if the stellar disk is inhomogenous and the planet does
not transit a perfectly representative chord, the transit depth will be
biased. Second, the stellar He line may change during the observations, making
the transit appear deeper or shallower than it really is.

Guilluy et al. (2020) used the five GIANO transit observations, which also
included simultaneous H$\alpha$ measurements, to extensively explore the first
type of variability. H$\alpha$ is an excellent tracer of stellar activity,
making it a useful proxy for the spot coverage fraction during each visit. In
active regions on the star, H$\alpha$ is seen in emission while the stellar He
I line becomes deeper. Guilluy et al. (2020) detected changes in the H$\alpha$
line during the two transit observations with the highest and lowest He I
depths, suggesting that the planet was transiting an inhomogenous star. The
remaining three transits had consistent He I depths ($\sim 0.5$%) and no
detectable variations in the stellar H$\alpha$ line. The deep He I transit
could be explained if the planet transited an especially quiet region of the
star (thus preferentially leaving the active regions unocculted). However, it
is difficult to explain why the other transit that took place during a period
of enhanced activity is shallower than expected. Guilluy et al. (2020) argue
that this observation is unlikely to be caused by a starspot occultation, as
they do not observe a concurrent weakening in the Si 10830 Å line, which is
sensitive to changes in temperature. They tentatively suggest that the planet
might have occulted a filament instead. Our new transit depth is also
shallower than expected, and we also do not see any evidence for a concurrent
change in the Si 10830 Å line during the transit. However, unlike Guilluy et
al. (2020), we do not see a mid-transit spike in the He band-integrated light
curve, nor any other abnormalities that would indicate the occultation of a
filament or starspot.

We can quantify the effects of stellar active region contamination on the
helium transmission spectrum using the techniques outlined in Cauley et al.
(2018). In order to reproduce the strength and shape of the stellar He I
10833Å triplet we found that an active region filling factor of $\approx 65\%$
and an optical depth at line center of $\tau\approx 0.7$ are required. We note
that this is the facular or plage filling factor and not a spot filling
factor. The large filling factor is in line with previous estimates for helium
absorption on active K-dwarfs (Andretta et al., 2017) and is similar to the
value of $75\%$ estimated by Guilluy et al. (2020). After estimating the
filling factor, we simulated the transit of HD 189733 b across the active
stellar surface and tuned the atmospheric parameters to match the observed
transmission spectrum. We then tested various filling factors and geometries
(e.g., uniform distribution, active latitudes) to see how the transmission
spectrum changes while the atmospheric parameters are held constant. We find
that even drastic changes in the filling factor ($\approx 30\%-40\%$) are not
enough to reproduce the difference between the atmospheric absorption observed
by Salz et al. (2018) and our NIRSPEC transmission spectrum.

The strength of the stellar helium line does not change much from epoch to
epoch: Guilluy et al. (2020) found that the line core flux varied by only
$\approx 4\%$ across a 1.5 year baseline. The NIRSPEC spectrum displays a
similar line depth ($\approx 0.7$) compared to Guilluy et al. (2020) and Salz
et al. (2018), suggesting that the surface activity level, at least as
measured by absorbing metastable helium regions, is fairly constant. This
implies that the filling factor cannot change dramatically from epoch to
epoch, otherwise we would have observed large variations in the stellar line
strength (of order the percent change in filling factor). The relative
stability of the star’s filling factor and the large size of the changes in
the filling factor required to account for the epoch to epoch variability in
the observed helium absorption signal make it unlikely that this variability
is due to inhomogeneities on the stellar surface.

On the other hand, variability in the stellar He line is a plausible
contributing factor to the differing transit depths. Figure 3 shows the light
curve of the helium line, centered slightly blueward of the main peak. The
planetary absorption is evident, but so is substantial variability after
transit, including a sudden drop 2 h after mid-transit and a slow 0.6% rise in
the 2 hours following the drop. This variability can also be seen from a close
examination of Figure 1, but the band-integrated light curve makes it clearer.
If similar stellar variability occurred during the transit or in the immediate
out-of-transit baseline, the planetary excess absorption could appear
$\sim$0.2% larger or smaller than it actually is.

###  4.2 Potential Evidence for Variability in the Planetary Outflow

It is plausible to think that the observed variability might instead be due to
changes in the planet’s atmospheric outflow, which might be caused by
variations in the star’s high-energy radiation or stellar wind environments.
In this case, we might also expect to see visit-to-visit variations in the
shape of the wavelength-dependent planetary absorption signal. Figure 5
compares the excess absorption spectrum obtained by the three instruments,
which tells the same story as the average depths: NIRSPEC sees a similar
signal as GIANO-B (all transits combined), which in turn sees a weaker signal
than CARMENES.

One interesting difference between the GIANO and NIRSPEC excess absorption
spectra is that despite having a very similar primary peak, the height of the
secondary peak at 10832Å is lower for GIANO by $\sim$0.1%. Adopting the formal
errors on the two spectra (0.07% for GIANO and 0.05% for NIRSPEC) and assuming
statistical independence, we calculate that the difference is significant at
the $\sim 3.5\sigma$ level. However, the significance of this feature is
likely overestimated, as noise in high resolution spectra is often correlated
across adjacent wavelengths. If the difference is real, we can use it to
constrain the optical thickness of the outflow. A completely optically thin
outflow would have a primary-to-secondary peak ratio of 8:1 (derived from the
$g_{i}f_{ik}$ of the three lines), whereas a completely optically thick
outflow would have a peak ratio of 1:1, as the gas would absorb all light in
both peaks. The higher peak ratio in the GIANO data could indicate that the
absorption is coming from a more optically thick region than in our NIRSPEC
observations. However, the fact that we see significant absorption extending
between the two peaks makes it difficult to determine how much of the
absorption at 10832Å is due to the line there and how much is due to
absorption from strongly blueshifted gas absorbing at a rest wavelength of
10833.3Å. Given the right kinematic structure–namely a long tail of gas being
accelerated away from the star by the stellar wind–it is conceivable that even
the CARMENES observations, for which the apparent peak ratio is 2.8 (Salz et
al., 2018), might be consistent with an optically thin outflow.

##  5 Modeling Variability in the Planetary Outflow

We can use atmosphere models to quantify the effect of stellar variability on
the helium absorption signal. Here we focus on the effect of variations in the
stellar XUV flux, which can be captured using relatively simple 1D
hydrodynamic models. The code we use is The PLUTO-CLOUDY Interface (TPCI)
(Salz et al., 2015a), a combination of two sophisticated and widely used
codes: the hydrodynamic solver PLUTO (Mignone et al., 2007), and the plasma
simulation and spectral synthesis code CLOUDY (Ferland et al., 2013). TPCI has
been extensively used to study photoevaporation (e.g. Salz et al. 2015b, 2016;
Kasper et al. 2020; Zhang et al. 2022a).

We obtain the nominal stellar spectrum in two ways. In the first method, we
use the XUV spectrum derived by Lampón et al. (2021b) (their Fig. 2, data
provided by M. Lamp’on) for the wavelength range of 5–1450 Å. We chose a
cutoff of 1450 Å because their spectrum redward of that wavelength is not
based on observations. In the second method, the X-ray spectrum is taken from
the thermal plasma model of Poppenhaeger et al. (2013), the Ly$\alpha$ flux at
1 AU is taken to be 11.8 erg s-1 cm-2 (France et al., 2013), and the EUV
spectrum is derived from the Ly$\alpha$ flux using the scaling relations of
Linsky et al. (2014). For both methods, fluxes at longer wavelengths are taken
from a PHOENIX model (Husser et al., 2013) ($T_{\rm eff}=5000$, log(g)=4.5,
solar metallicity). Of particular importance is the spectrum at 1230–2588 Å,
which we call mid ultraviolet, because this radiation ionizes metastable
helium, but cannot contribute to producing metastable helium because it cannot
ionize ground state hydrogen. Fortunately, we have observational constraints
for the flux at these wavelengths from XMM-Newton’s Optical Monitor, a
photometer which operates simultaneously with the X-ray instruments. We
reviewed all archival observations taken by the OM in the UVW2 (212 nm, width
50 nm) and UVM2 (231 nm, width 48 nm) filters, and found that they do not
differ by more than a few percent. (For a summary of all XMM-Newton
observations of HD 189733, see Pillitteri et al. (2022).) To roughly match the
observations, we boosted the model spectrum between 1230–2588 Å by 30% and
used it to predict the count rates that OM should see. The prediction was only
6% too low for UVW2 and 3% too high for UVM2, indicating our model MUV flux is
accurate. The final spectra constructed using the two methods is shown in
Figure 6, while Table 1 shows the band-integrated fluxes.

Comparing the two methods, the X-ray flux is in excellent agreement: Lampón et
al. (2021b) reports a X-ray luminosity that corresponds to a flux of 7.8 erg
s-1 cm-2 at 1 AU, while we obtain 7.6 erg s-1 cm-2. However, the EUV (100–912
Å) flux is dramatically discrepant. Lampón et al. (2021b) imply a flux of 57
erg s-1 cm-2 at 1 AU, while we obtain 6.5 erg s-1 cm-2, a difference of 9x.
Other authors have obtained 11 (Sanz-Forcada et al., 2011), 22 (Bourrier et
al., 2020), and 36 erg s-1 cm-2 (Poppenhaeger et al., 2013). This large
uncertainty in the EUV flux is ultimately because the stellar EUV is only
measurable from space and no space telescopes currently have EUV capabilities
(France et al., 2022). The uncertainty has been noted by many previous
publications. For example, Oklopčić (2019) compared the results of two
different methods of EUV reconstruction for the same star that differed by an
order of magnitude, which produced helium excess absorption depths that
differed by a factor of 3. France et al. (2022) compared four different
methods of reconstructing the EUV spectrum for Proxima Centauri and found that
they were discrepant by 3–100x, depending on the wavelength.

Figure 6: Reconstructed stellar spectrum using two methods. The two methods
differ blueward of 1450 Å, where we use Lampón et al. (2021b) for method 1 and
our own reconstruction for method 2. Note the large discrepancy in the EUV.
Table 1: Band-integrated fluxes from both stellar spectrum reconstruction
methods

Band | Wavelengths | Flux (method 1) | Flux (method 2)
---|---|---|---
| Å |  erg/s/cm^2 |  erg/s/cm^2
X-ray | 5-100 | 7.6 ±2.2 | 7.8 ±2.2
EUV | 100-912 | 57 | 6.5
He-ionizing | 5-512 | 24 | 12
Ly$\alpha$ | 1214-1217 | 19 | 11.8 ±3.5
MUV | 1230-2588 | 105 ±7 | 105 ±7
Total∗ | 5-50,000 | 4.80 ±0.19 ×10^5 | 4.80 ±0.19 ×10^5

∗From Addison et al. (2019)

Since the EUV flux is uncertain to an order of magnitude, we analyze the
behavior of the escaping atmosphere over a wide range of EUV fluxes. For the
first stellar spectrum we constructed, which has low EUV flux compared to
literature values, we set up four 1D hydrodynamic simulations: at 1x, 1.3x,
2x, and 3x the nominal stellar XUV flux as seen by the planetary dayside
(which is half the flux seen by the substellar point). For the second
spectrum, which has high EUV flux compared to other literature values, we set
up five 1D hydrodynamic simulations: at 0.5x, 1x, 1.3x, 2x, and 3x the nominal
stellar XUV flux. Following Zhang et al. (2022b, a), we ran the simulation for
100 time units with advection off, where the time unit is calculated as the
planetary radius divided by 10 km/s (roughly the sound speed). Unlike in this
previous study, we did not run the models for another 100 time units with
advection on, due to numerical problems. Also unlike in the previous study, we
irradiate the atmosphere with the average flux experienced by the planetary
dayside, and not the flux experienced by the substellar point. The former is
half of the latter, and is more representative of the outflow as a whole.

Figure 7: Relationship between maximum excess absorption depth in our 1D
models and the stellar EUV flux at 1 AU. Blue points represent XUV stellar
spectra scaled from the spectrum of Lampón et al. (2021b), while red points
represent XUV spectra scaled from our own spectrum.

Figure 7 shows the results of the TPCI simulations as a function of EUV flux.
The simulation results are similar to those of other 1D models. In agreement
with Lampón et al. (2021a), we find that the the outflow is hot ($>\sim$10,000
K) and in the recombination regime, which is characterized by a sharp
transition from neutral to ionized hydrogen. For the 1/2 nominal XUV model
using the first stellar spectrum, which is the model with median EUV flux (see
Figure 7), the outflow is 95% neutral at 1.12 $R_{p}$ and drops to 5% neutral
by 1.32 $R_{p}$. In the recombination regime, most of the incident XUV is
radiated away via recombination, leading to low mass loss efficiency. In fact,
the mass loss rate we obtain for the aforementioned model, $5\times 10^{9}$
g/s, is substantially lower than even the $10^{11}$ g/s found by Lampón et al.
(2021a, b), and similar to the $10^{10}$ g/s found by Caldiroli et al. (2022)
for similar XUV flux. This corresponds to a mass loss efficiency of 1.3%,
similar to the 2–3% found by Caldiroli et al. (2022).

The absorption depth is higher than observed even in the simulation with the
lowest EUV. Increasing the EUV flux 34% above this lowest point increased the
helium absorption by 60%, but the depth increases roughly with the square root
of EUV flux thereafter, so that further 34% increases in flux only produce 15%
increases in depth.

HD 189733 has been extensively observed at XUV wavelengths, providing us with
good empirical constraints on the magnitude of the star’s variability.
Pillitteri et al. (2022) analyzed 25 XMM-Newton observations spanning a total
of 8 years, the last of which was taken before the first helium observation.
They found that the 25–75% range of quiescent flux varies by 23%, and that 44%
of the total observing time was occupied by flares, which typically lasted
several kiloseconds and increased the XUV flux by a few tens of percent (see
their Fig. B.2). Using the data behind their Figure B.3 (Pillitteri, private
comm.), we calculate that the typical variability in the X-ray flux, as
defined by the gap between the 84th and 16th percentiles, is 60%. EUV flux
increases with X-ray flux as $F_{\rm EUV}\propto F_{\rm X}^{0.86}$ (Sanz-
Forcada et al., 2011), so a 60% increase in X-rays corresponds to a 50%
increase in EUV, and a 23% increase in X-rays corresponds to a 19% increase in
EUV. These numbers suggest that the observational variability plotted in
Figure 4 can plausibly be due to variability in stellar EUV output if the EUV
flux is low ($\sim$6 erg s-1 cm-2 at 1 AU), but is less likely if the EUV flux
is high ($>\sim$10).

It is worthwhile to briefly consider other sources of variability that we
could not explore with our 1D model. The effect of differing stellar wind
conditions on the outflow has been extensively studied with 3D simulations
(e.g. McCann et al. 2019; Vidotto & Cleary 2020), but generally without
modelling the helium line. MacLeod & Oklopčić (2022) does model the helium
line, and find that given the same planetary mass loss rate, the equivalent
width of helium absorption changes by 10% when the stellar wind is increased
10x from weak to moderate, but changes far more drastically when the wind is
strong enough to “break through” the sonic surface and confine the outflow.
Zhang et al. (2022a) study the effect of different stellar wind conditions on
the helium absorption from TOI 560b, a young mini-Neptune, by alternately
halving the stellar wind density and velocity. They find changes in peak
absorption of up to 25%, but smaller changes in equivalent width of up to 12%.
These studies show that changes in stellar wind conditions can plausibly
explain a portion of the variability we see in HD 189733b.

Another source of variability is shear instability, which exists even with
constant irradiation and wind from the star. For the inflated sub-Saturn
WASP-107b, Wang & Dai (2021b) uses a 3D hydrodynamic model to predict
$\sim$10% fluctuations in helium absorption over hour-long timescales. These
fluctuations, however, change the planetary mass loss rate by less than 1%.
Shear instability, like the impact of the stellar wind, is not possible to
model in 1D.

One final source of variability that has been modelled in the literature is
flares. Flares increase the star’s XUV output, and a sufficiently long-lasting
flare would increase helium absorption, as predicted by our models. However,
depending on the spectrum of the flare, the immediate effect of a flare might
be to ionize metastable helium and reduce helium absorption. After a dynamical
timescale, when the surge of photoevaporative mass loss generated by the flare
has reached higher altitudes, helium absorption should increases above
baseline (Wang & Dai, 2021a). For WASP-69b, the planet simulated by Wang & Dai
(2021a), this takes about an hour; helium absorption then stays elevated for 2
days after the end of the flare. Thus, even though flares can both decrease
and increase helium absorption, the latter is far more likely to be observed
than the former. To explore the effect of flares that raise the XUV flux by
tens of percent on a high-gravity planet like HD 189733b–rather than a 10x
flare on a low gravity planet like WASP-69b, as Wang & Dai (2021a) did–new 3D
simulations are desirable.

##  6 Conclusion

Metastable helium observations provide us with a sensitive probe of planetary
outflow characteristics. In this paper, we present a new helium transit
observation of HD 189733b, a hot Jupiter that has been observed over a wide
range of wavelengths in the fifteen years since its discovery. This planet
will be also observed by 5 _James Webb Space Telescope_ programs in Cycle 1
alone; these programs will doubtlessly further expand our knowledge of its
atmospheric properties.

When combined with the set of eight published helium transit observations of
HD 189733b, our new observations add to the growing body of evidence
suggesting that this planet’s outflow properties may vary in time. HD 189733
is a relatively active K dwarf and it has been extensively observed at high
energies by XMM-Newton and Chandra. Our simulations suggest that stellar XUV
variability, by itself, can plausibly explain the observed variations in the
planet’s helium transit depth. However, this does not mean that other sources
of variability we do not model, such as stellar wind variations, shear
instability, and short flares, cannot also explain the variability.

Understanding the high-energy environment and hydrostatic atmosphere of the
planet is crucial to understanding mass loss, and the quantity of atmospheric
observations available for HD 189733b make it one of the most promising
planets for understanding mass loss physics in detail. In order to better
understand the amplitude, timescale, and likely cause of this planet’s
observed variability, it is important to continue monitoring of the planet in
the helium line, preferably with simultaneous H$\alpha$ and/or X-ray
measurements. This will in turn inform how much we should rely on single-epoch
observations of helium absorption to quantify the mass loss rates and
population-level properties of other close-in planets.

Software: numpy (van der Walt et al., 2011), scipy (Virtanen et al., 2020),
matplotlib (Hunter, 2007)

The data presented herein were obtained at the W. M. Keck Observatory, which
is operated as a scientific partnership among the California Institute of
Technology, the University of California and the National Aeronautics and
Space Administration. The Observatory was made possible by the generous
financial support of the W. M. Keck Foundation. AO gratefully acknowledges
support from the Dutch Research Council NWO Veni grant.

## References

  * Adams (2011) Adams, F. C. 2011, ApJ, 730, 27
  * Addison et al. (2019) Addison, B., Wright, D. J., Wittenmyer, R. A., et al. 2019, PASP, 131, 115003
  * Allart et al. (2019) Allart, R., Bourrier, V., Lovis, C., et al. 2019, A&A, 623, A58
  * Andretta et al. (2017) Andretta, V., Giampapa, M. S., Covino, E., Reiners, A., & Beeck, B. 2017, ApJ, 839, 97
  * Bouchy et al. (2005) Bouchy, F., Udry, S., Mayor, M., et al. 2005, A&A, 444, L15
  * Bourrier et al. (2020) Bourrier, V., Wheatley, P. J., Lecavelier des Etangs, A., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 559
  * Caldiroli et al. (2022) Caldiroli, A., Haardt, F., Gallo, E., et al. 2022, A&A, 663, A122
  * Cauley et al. (2018) Cauley, P. W., Kuckein, C., Redfield, S., et al. 2018, AJ, 156, 189
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137
  * France et al. (2013) France, K., Froning, C. S., Linsky, J. L., et al. 2013, ApJ, 763, 149
  * France et al. (2022) France, K., Fleming, B., Youngblood, A., et al. 2022, Journal of Astronomical Telescopes, Instruments, and Systems, 8, 014006
  * Fulton & Petigura (2018) Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109
  * Guilluy et al. (2020) Guilluy, G., Andretta, V., Borsa, F., et al. 2020, A&A, 639, A49
  * Harrington et al. (2022) Harrington, J., Himes, M. D., Cubillos, P. E., et al. 2022, The Planetary Science Journal, 3, 80
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Husser et al. (2013) Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Kasper et al. (2020) Kasper, D., Bean, J. L., Oklopčić, A., et al. 2020, AJ, 160, 258
  * Knutson et al. (2007) Knutson, H. A., Charbonneau, D., Allen, L. E., et al. 2007, Nature, 447, 183
  * Lampón et al. (2021a) Lampón, M., López-Puertas, M., Czesla, S., et al. 2021a, A&A, 648, L7
  * Lampón et al. (2021b) Lampón, M., López-Puertas, M., Sanz-Forcada, J., et al. 2021b, A&A, 647, A129
  * Linsky et al. (2014) Linsky, J. L., Fontenla, J., & France, K. 2014, ApJ, 780, 61
  * MacLeod & Oklopčić (2022) MacLeod, M., & Oklopčić, A. 2022, ApJ, 926, 226
  * McCann et al. (2019) McCann, J., Murray-Clay, R. A., Kratter, K., & Krumholz, M. R. 2019, The Astrophysical Journal, 873, 89
  * McCullough et al. (2014) McCullough, P. R., Crouzet, N., Deming, D., & Madhusudhan, N. 2014, ApJ, 791, 55
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, The Astrophysical Journal Supplement Series, 170, 228
  * Morello et al. (2014) Morello, G., Waldmann, I. P., Tinetti, G., et al. 2014, ApJ, 786, 22
  * Nortmann et al. (2018) Nortmann, L., Pallé, E., Salz, M., et al. 2018, Science, 362, 1388
  * Oklopčić (2019) Oklopčić, A. 2019, ApJ, 881, 133
  * Oklopčić & Hirata (2018) Oklopčić, A., & Hirata, C. M. 2018, ApJ, 855, L11
  * Pillitteri et al. (2022) Pillitteri, I., Micela, G., Maggio, A., Sciortino, S., & Lopez-Santiago, J. 2022, A&A, 660, A75
  * Poppenhaeger (2022) Poppenhaeger, K. 2022, MNRAS, 512, 1751
  * Poppenhaeger et al. (2013) Poppenhaeger, K., Schmitt, J. H. M. M., & Wolk, S. J. 2013, ApJ, 773, 62
  * Rosenthal et al. (2021) Rosenthal, L. J., Fulton, B. J., Hirsch, L. A., et al. 2021, ApJS, 255, 8
  * Salz et al. (2015a) Salz, M., Banerjee, R., Mignone, A., et al. 2015a, A&A, 576, A21
  * Salz et al. (2016) Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75
  * Salz et al. (2015b) Salz, M., Schneider, P. C., Czesla, S., & Schmitt, J. H. M. M. 2015b, A&A, 576, A42
  * Salz et al. (2018) Salz, M., Czesla, S., Schneider, P. C., et al. 2018, A&A, 620, A97
  * Sanz-Forcada et al. (2011) Sanz-Forcada, J., Micela, G., Ribas, I., et al. 2011, A&A, 532, A6
  * Seager & Sasselov (2000) Seager, S., & Sasselov, D. D. 2000, ApJ, 537, 916
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, The Astrophysical Journal, 699, 564
  * Sing et al. (2011) Sing, D. K., Pont, F., Aigrain, S., et al. 2011, MNRAS, 416, 1443
  * Spake et al. (2018) Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nature, 557, 68
  * Steinrueck et al. (2021) Steinrueck, M. E., Showman, A. P., Lavvas, P., et al. 2021, MNRAS, 504, 2783
  * Szabó & Kiss (2011) Szabó, G. M., & Kiss, L. L. 2011, The Astrophysical Journal, 727, L44
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22
  * Vidotto & Cleary (2020) Vidotto, A. A., & Cleary, A. 2020, Monthly Notices of the Royal Astronomical Society, 494, 2417
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wang & Dai (2021a) Wang, L., & Dai, F. 2021a, ApJ, 914, 98
  * Wang & Dai (2021b) —. 2021b, ApJ, 914, 99
  * Zhang et al. (2020) Zhang, M., Chachan, Y., Kempton, E. M. R., Knutson, H. A., & Chang, W. H. 2020, ApJ, 899, 27
  * Zhang et al. (2022a) Zhang, M., Knutson, H. A., Wang, L., Dai, F., & Barragán, O. 2022a, AJ, 163, 67
  * Zhang et al. (2021) Zhang, M., Knutson, H. A., Wang, L., et al. 2021, AJ, 161, 181
  * Zhang et al. (2022b) —. 2022b, AJ, 163, 68