# REVISITING HETEROPHILY IN GRAPH CONVOLUTIONAL NETWORKS BY LEARNING
REPRESENTATIONS ACROSS TOPOLOGICAL AND FEATURE SPACES

###### Abstract

Graph convolution networks (GCNs) have been enormously successful in learning
representations over several graph-based machine learning tasks. Specific to
learning rich node representations, most of the methods have solely relied on
the homophily assumption and have shown limited performance on the
heterophilous graphs. While several methods have been developed with new
architectures to address heterophily, we argue that by learning graph
representations across two spaces i.e., topology and feature space GCNs can
address heterophily. In this work, we experimentally demonstrate the
performance of the proposed GCN framework over semi-supervised node
classification task on both homophilous and heterophilous graph benchmarks by
learning and combining representations across the topological and the feature
spaces. 11footnotetext: Code and additional implementation details can be
found here.

Index Terms—  Graph Convolutional Networks, Homophily, Heterophily.

##  1 Introduction

22footnotetext: This work is supported by SERB MATRICS grant.

Owing to their great capacity in node representation learning, graph
convolutional networks (GCNs) have widely been adopted across numerous graph-
based machine learning tasks ranging from graph classification [1], link
prediction [2], recommendation systems [3], fault diagnosis [4], network
embedding [5], and computer vision [6]. However, most of the existing methods
are centered around the homophily assumption for semi-supervised node
classification (SSNC) task, i.e, a pair of nodes tend to be connected if they
are similar or alike. The homophilous assumption is mainly attributed to the
recursive neighbourhood aggregation mechanism in GCNs wherein a node’s
representation is obtained from its neighbours which are likely to have the
same label. Several works [7, 8] advocate that GNNs are best suited to address
homophily and are not suitable for heterophilous graphs and have designed new
architectures to explicitly address heterophily in graphs. However, Ma _et
al._ [9] have shown that there exist certain specific types of heterophily
over which GCNs in their original form can show strong performance. As per
[9], GCNs learn equivalent embeddings for nodes that share same neighbourhood
distribution. Further, Wang _et al._ [10] studied the performance of multi-
channel graph convolutional networks in obtaining the information from both
topological and features space. In this work, we take a step further and
demonstrate that by learning and fusing features from topological and feature
spaces, GCNs can address heterophily in a more general sense, i.e., not
restricted to any specific kind of neighbourhood distributions.

Homophily and heterophily in graphs. When the connected nodes in a graph have
similar labels or properties, the graph is said to have homophily. Homophily
is found in several real-world graphs such as friendship networks [11],
political networks [12, 13], citation networks [14] and others. When the
connected nodes have different properties or labels, the graph is said to
possess heterophily. The examples of heterophilous graphs include dating
networks or molecular networks [7]. We observe that implicitly there is some
sort of commonality within a heterophilous connections. Consider a dating
network, where potential dates are more likely to be connected if they share
common interest or hobbies. Similarly, in molecular or protein networks,
different types of amino acids make up a protein based on their bio-chemical
compatibility. In both these examples, the common attributes between nodes of
different labels shows inherent homophily (commonality) in a heterophilous
connection.

GCNs for heterophily. Following these observations, we argue that GCNs can
show improved performance over heterophilous graphs if they can be modeled to
identify these hidden commonalities within a heterophilous connection. In
traditional GCNs, the node representation learning is guided solely by
topology graphs. However, both network topology and node features offer
holistic understanding for the graph representation learning. We cannot obtain
accurate neighborhood information by solely considering the topological view
as the node features are often smoothened out when the feature information is
propagated over network topological structure. Also, GCNs in their original
form fail to adaptively learn correlation information between topological and
feature views [10]. In this work, we propose a multi-channel GCN framework
that first extracts information individually from the topology and feature
graph by using two separate GCN encoders and then uses another encoder to
extract information common to both the spaces for the node classification task
and apply this to address heterophily. Remember that the role of the common
encoder is to fill in the missing information obtained from each of the
individual spaces.

##  2 Related Work

In this section, we briefly review the development of GCNs over the years and
understand how they have been applied to address heterophily. Bruna _et al._
[15] designed the graph convolution operation in Fourier domain using graph
Laplacian and Defferrard _et al._ [16] employed its Chebyshev expansion to
improve the efficiency. Later, Kipf and Welling [17] simplified the
convolution operation and proposed feature aggregation from the one-hop
neighbors. After this, a plenty of variants like GraphSAGE [18], GAT [19], SGC
[20], and GMNN [21]. were introduced. While they have have shown performance
improvement over node classification, most of them have suffered due to
inherent low-pass filtering of node features due to Laplacian smoothing. For a
more detailed review on different variants of GCNs, we refer the readers to
[22, 23]. Surprisingly, all these methods only use a single topology graph for
node aggregation and fail to completely utilize the information from the
feature space. A few methods like AM-GCN [10] have shown better performance by
extracting embeddings from both topology graph and feature graph. However,
these have mostly performed well on homophilous graphs. Several recent works
have methods to address heterophily in graphs by re-designing or modifying
model architectures such as Geom-GCN [24], H2GCN [7], GPR-GNN [25], and CPGNN
[8]. Further, [9] has established that GCN can potentially achieve good
performance on some specific types of heterophily in graphs. In this work, our
focus is on understanding whether the information learned and fused across
topological and feature graphs can be used to address heterophily in general.

##  3 Method

###  3.1 Mathematical Preliminaries.

Let us consider an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$
characterized by an adjacency matrix $\mathbf{A}$, a diagonal degree matrix
$\mathbf{D}$ ($\mathbf{\widetilde{A}}$ and $\mathbf{\widetilde{D}}$ when self-
loops are considered), $N$ nodes, and a node feature matrix $\mathbf{X}\in\rm
I\\!R^{N\times d}$ containing a $d$-dimensional feature vector for each node
in the graph.

GCN Layer. A single GCN layer takes the form of
$\mathbf{H}^{(l+1)}=\sigma(\widetilde{\mathbf{P}}\mathbf{H}^{(l)}\mathbf{W}^{(l)})$,
where $\mathbf{H}^{(l)}$ and $\mathbf{H}^{(l+1)}$ represent the input and
output features of the $l^{th}$ layer, $\mathbf{W}^{(l)}$ represents the
parameter (weight) matrix that transforms features, and $\sigma$ represents
ReLU activation. Further,
$\widetilde{\mathbf{P}}=\widetilde{\mathbf{D}}^{-1/2}\widetilde{\mathbf{A}}\widetilde{\mathbf{D}}^{-1/2}$
is the degree normalized graph Laplacian matrix.

Homophily. The edge homophily ratio is defined as the fraction of edges that
connect nodes with the same labels [9]. Mathematically,

$\centering\mathcal{H}(\mathcal{G},\mathcal{Y})=\frac{1}{|\mathcal{E}|}\sum_{(i,j)\in\mathcal{E}}\mathbbm{1}(y_{i}=y_{j})\@add@centering$
(1)

Here, $\mathcal{Y}$ is the set of all the node labels $y_{i}$ and
$\mathbbm{1}(\cdot)$ is an indicator function. $0.5\leq\mathcal{H}(\cdot)\leq
1$ corresponds to high homophily and $0\leq\mathcal{H}(\cdot)\leq 0.5$
corresponds to high heterophily.

###  3.2 Proposed Framework

The overall information flow in the proposed framework is shown in Figure 1.
The main idea is to propagate the node features through both the topological
and the feature space. We start with two 2-layer GCN encoders
$\mathcal{T}_{enc}$ and $\mathcal{F}_{enc}$ that take
$\mathcal{G}=(\mathbf{H^{(0)}},\mathbf{A})$ and
$\mathcal{G}_{f}=(\mathbf{H^{(0)}},\mathbf{A}_{f})$ as input, respectively
(with $\mathbf{H^{(0)}}=f_{MLP}(\mathbf{X})$), to obtain the feature maps
$\mathbf{Z}_{T}$ and $\mathbf{Z}_{F}$ corresponding to two views. Please note
that $\mathcal{G}$ is the regular topological graph and $\mathcal{G}_{f}$ is
the k-nearest neighbour feature graph based on node feature matrix
$\mathbf{X}$ and adjacency matrix $\mathbf{A}_{f}$ of the k-NN graph.
33footnotetext: We use cosine similarity to compute the similarity $s_{ij}$
between two node features $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$.
$s_{ij}=\frac{\mathbf{x}_{i}.\mathbf{x}_{j}}{|\mathbf{x}_{i}||\mathbf{x}_{j}|}$

Further, each GCN layer is designed to have a residual connection to initial
features $\mathbf{H}^{(0)}$ (or $\mathbf{X}$), as described in Equation 2.

$\centering\mathbf{H}^{(l+1)}=\alpha\mathbf{H}^{(l)}+(1-\alpha)\mathbf{H}^{(0)}\@add@centering$
(2)

Here, $\alpha=0.8$. We deploy a residual connection from the input features
$\mathbf{H}^{(0)}$ such that the final representation of each node retains at
least a fraction of information $(\alpha)$ from the input features. This helps
to slow down the rate of convergence of features to a point or smaller
subspace - a common cause of oversmoothing in GCNs [26, 27].

As discussed in Section 1, we learn the information with common
characteristics between these two spaces. We deploy a common GCN encoder
$\mathcal{C}_{enc}$ with shared parameter to learn the embeddings
$\mathbf{Z}_{CT}$ and $\mathbf{Z}_{CF}$ which in turn are combined together to
get the common embedding $\mathbf{Z}_{C}$ using a multi-layer perceptron
(MLP). However, unlike methods like AM-GCN [10], we use a weighted combination
of the learned features ($\mathbf{Z}_{T}$ and $\mathbf{Z}_{F}$) and the
initial features $\mathbf{X}$ as input to $\mathcal{C}_{enc}$. Specifically,
$\mathbf{Z}_{{C}_{in}}=\zeta\mathbf{H_{\star}}+(1-\zeta)\mathbf{Z}_{\star}$
Here, $\mathbf{H_{\star}}=f_{MLP}(\mathbf{H^{0}})$ and $\mathbf{Z}_{\star}$
represents either $\mathbf{Z}_{T}$ or $\mathbf{Z}_{F}$ and we set $\zeta=0.85$
This formulation was found to provide performance improvement over just using
the learned $\mathbf{Z}_{T}$ or $\mathbf{Z}_{F}$.

With $\mathbf{Z}_{T}$, $\mathbf{Z}_{F}$, and $\mathbf{Z}_{C}$ at hand, we now
wish to find an optimal combination of these features that decides the final
predicted labels, as governed by Equation 3.

$\centering\widetilde{\mathbf{Z}_{T}}=\bm{\alpha}_{T}\odot\mathbf{Z}_{T}+\bm{\alpha}_{C}\odot\mathbf{Z}_{C}\@add@centering$
$\centering\widetilde{\mathbf{Z}_{F}}=\bm{\alpha}_{F}\odot\mathbf{Z}_{F}+\bm{\alpha}_{C}\odot\mathbf{Z}_{C}\@add@centering$
(3)

Here, $\bm{\alpha}_{T}$, $\bm{\alpha}_{F}$, and $\bm{\alpha}_{C}$ are the
feature-level attention weights for $\mathbf{Z}_{T}$, $\mathbf{Z}_{F}$, and
$\mathbf{Z}_{C}$, respectively, such that $\bm{\alpha}_{\star}=f_{\star
MLP}(Z_{\star})$ where $\star\in\\{T,F,C\\}$, $f_{\star MLP}$ is the
corresponding MLP, and $\odot$ is the element-wise Hadamard product.
Consequently, we obtain the final set of labels $\widehat{\mathbf{Y}}$ as
shown in Equation 4.

$\centering\widehat{\mathbf{Y}}=softmax(\mathbf{W}.\widehat{\mathbf{Z}}+\mathbf{b})\@add@centering$
(4)

where,
$\widehat{\mathbf{Z}}=(\widetilde{\mathbf{Z}_{T}}||\widetilde{\mathbf{Z}_{F}})$.

Fig. 1: The architecture of the proposed framework.

###  3.3 Training Details

We train and evaluate the proposed framework over both homophilous and
heterophilous graph datasets as listed in Table 1.

Closeness Constraint. We want the two output embeddings $\mathbf{Z}_{CT}$ and
$\mathbf{Z}_{CF}$ arising out of $\mathcal{C}_{enc}$ to be close as possible.
Therefore, we apply a closeness constraint enabling to search for commonality.
Following [10], we first obtain the similarity of $N$ nodes as,

$\centering\mathbf{S}_{T}=\mathbf{Z}_{CTnor}.{\mathbf{Z}^{T}}_{CTnor}\@add@centering$
$\centering\mathbf{S}_{F}=\mathbf{Z}_{CFnor}.{\mathbf{Z}^{T}}_{CFnor}\@add@centering$
(5)

Here, $\mathbf{Z}_{CTnor}$ and $\mathbf{Z}_{CFnor}$ are the $L_{2}$-normalized
version of $\mathbf{Z}_{CT}$ and $\mathbf{Z}_{CF}$, respectively. Using these,
we define the closeness constraint using the Frobenious norm, as shown in
Equation 6.

$\centering\mathcal{L}_{c}=||\mathbf{S}_{T}-\mathbf{S}_{F}||^{2}_{F}\@add@centering$
(6)

Disparity Constraint. An input graph from the same view generates two
representations, i.e.,
$\mathcal{G}\rightarrow(\mathbf{Z}_{T},\mathbf{Z}_{CT})$ and
$\mathcal{G}_{f}\rightarrow(\mathbf{Z}_{F},\mathbf{Z}_{CF})$. In order to
ensure that they capture different information, we introduce a disparity
regularization $\mathcal{L}_{d}$ (similar to [5]) to enlarge distances between
them using the negative cosine similarity, as described in Equation 7.

$\centering\mathcal{L}_{d}=-\frac{1}{N}\sum_{i=1}^{N}\Biggl{(}\frac{\mathbf{z}_{t}^{(i)}.\mathbf{z}_{ct}^{(i)}}{||\mathbf{z}_{t}^{(i)}||.||\mathbf{z}_{ct}^{(i)}||}+\frac{\mathbf{z}_{f}^{(i)}.\mathbf{z}_{cf}^{(i)}}{||\mathbf{z}_{f}^{(i)}||.||\mathbf{z}_{cf}^{(i)}||}\Biggl{)}\@add@centering$
(7)

Classification Loss. Let $\mathbf{Y}$ be the one-hot encoded label matrix,
$\mathcal{V}_{train}$ be the set of training nodes, and $\widehat{\mathbf{Y}}$
be the prediction, as defined in Equation 4. We minimize the standard cross-
entropy loss to train the network.

$\centering\mathcal{L}_{cl}=-\sum_{v\in\mathcal{V}_{train}}\mathbf{Y}ln\widehat{\mathbf{Y}}\@add@centering$
(8)

The overall optimization objective function is as described in Equation 9.

$\mathcal{L}_{T}=\lambda\mathcal{L}_{cl}+\beta\mathcal{L}_{c}+\gamma\mathcal{L}_{d}$
(9)

Here, $\lambda,\beta$, and $\gamma$ are the weights for the classification
loss, closeness, and disparity constraints, respectively. 111Due to space
constraints, we provide additional implementation details including the hyper-
parameter values in the associated config files here.

##  4 Experimental Evaluation

In this section, we evaluate the proposed approach on homophilous and
heterophilous graph datasets (see Table 1) and reason out several design
choices. Please note that we have chosen high performing standard GCN variants
and AM-GCN [10] for a fair comparison keeping in mind that our aim is to
establish the ability of GCNs in addressing heterophily with any major change
in its original formulation. All experiments have been performed with $40$
labeled nodes per class. We find that this is a reasonable approximation of
the real world scenario where we are less likely to provided with lots of
labeled nodes.

Dataset | # Classes | # Nodes | # Edges | # Features | Heterophily
---|---|---|---|---|---
ACM | 3 | 3,025 | 13,128 | 1,870 | 0.179
Citeseer | 6 | 3,327 | 4,732 | 3,703 | 0.253
CoraFull | 70 | 19,793 | 65,311 | 8710 | 0.412
BlogCatalog | 6 | 5,196 | 1,71,743 | 8,189 | 0.599
UAI2010 | 19 | 3,067 | 28,311 | 4,973 | 0.636
Flickr | 9 | 7,575 | 2,39,738 | 12,047 | 0.761

Table 1: Dataset Statistics.

###  4.1 Semi-supervised node classification

Table 2 shows the performance of the proposed approach with the most recent
and relevant baselines over six different graph datasets in terms of accuracy
and the F1-score. We observe that the proposed framework either performs best
or the second best among different baselines. It achieves relative
improvements of around $3.78\%$ (accuracy) and $4.61\%$ (F1-score) averaged
over the heterophilous datasets. Further, the proposed framework is observed
to outperform GCN and kNN-GCN on nearly all the datasets. This indicates the
effectiveness of extracting useful information from both the spaces rather
than only from the topological space (as in GCN) or the feature space (as in
kNN-GCN). Moreover, we find that GCN operating solely on a topology graph
fails to provide better result than on a feature space graph. This is due to
the inherent useful differences in information that each space carries, thus
establishing the need of introducing feature space graph in GCN. Overall, the
margin of improvement over heterophilous datasets is observed to be higher
than that over the homophilous datasets. This establishes our claim that
learning representation across topological and feature spaces does help in
better addressing the heterophily.

Network design choices We discuss the three important design choices in the
proposed framework. (a) Initial Residual Connection: We observed an
improvement of $2.14\%$ (accuracy) and $3.09\%$ (F1-score) on an average by
using initial residual connection (Equation 2) in GCNs versus vanilla GCN.
(b). Feature-wise attention values: Owing to a better feature selectivity for
node classification [27], we observed improved performance by learning
feature-wise attention values (weighing each feature of every node) instead of
just one attention value per node (as AM-GCN[10]). Finally, (c) Negative
cosine similarity as disparity constraint: It speeds up (almost $\times 4$)
the training (and offers better performance) when compared to Hilbert-Schmidt
Independence Crtierion (HSIC), as used in [10].

Datasets | Metrics | GCN | kNN-GCN | GAT | MixHop | AM-GCN | Ours
---|---|---|---|---|---|---|---
ACM | Acc | 89.06 | 81.66 | 88.60 | 82.34 | 90.76 | 90.71
F1 | 89.00 | 81.53 | 88.55 | 81.13 | 90.66 | 90.75
Citeseer | Acc | 73.10 | 61.54 | 73.04 | 71.48 | 74.70 | 74.70
F1 | 69.70 | 59.33 | 69.58 | 67.40 | 69.81 | 70.39
CoraFull | Acc | 58.12 | 44.80 | 60.61 | 57.20 | 58.60 | 59.90
F1 | 54.87 | 40.42 | 56.43 | 53.55 | 54.54 | 55.31
BlogCatalog | Acc | 71.28 | 80.84 | 67.40 | 71.66 | 84.94 | 88.83
F1 | 70.71 | 80.16 | 66.39 | 70.84 | 84.32 | 88.39
UAI2020 | Acc | 51.80 | 68.74 | 63.74 | 65.05 | 73.40 | 76.12
F1 | 33.80 | 54.45 | 45.08 | 53.86 | 62.19 | 64.46
Flickr | Acc | 45.48 | 75.08 | 38.44 | 55.19 | 79.81 | 82.90
F1 | 43.27 | 75.40 | 36.94 | 56.25 | 78.60 | 82.82

Table 2: Node classification results in terms of Accuracy and F1-score across
different real-world datasets. Best performance is reported in BOLD and second
best performance is (”underlined”).

1

Input: $\mathcal{G}=\\{\mathcal{V},\mathcal{E}\\}$, $\widehat{\mathcal{H}}$,
$N$ and $\mathcal{Y}=\\{y_{1},y_{2},...,y_{\mathcal{C}}\\}$

Output: $\mathcal{G^{\prime}}=\\{\mathcal{V},\mathcal{E^{\prime}}\\}$

2 Initialise $\mathcal{G^{\prime}}=\\{\mathcal{V},\mathcal{E}\\}$ and
$\widehat{\mathcal{H}}=h_{init}$; the original heterophily of the dataset.

3 Obtain $L$; the list containing $10$ heterophily ratios such that
$h_{init}\leq\widehat{\mathcal{H}}\leq 0.95$

4

5for  _$\widehat{\mathcal{H}}\in L$_ do

6 Compute $K$, the number of heterophilous edges to added

7 Sample $K$ nodes such that each node $i\sim Uniform(\mathcal{V})$

8

9 for  _$i=1$ to $K$_ do

10 Obtain the label $y_{i}$ for node $i$

11 Sample a label $y_{j}\in\mathcal{Y}-y_{i}$ with probability $N_{j}/N$

12 Sample a node $j\sim Uniform(\mathcal{V}_{j})$

13 Update the edge set
$\mathcal{E^{\prime}}=\mathcal{E^{\prime}}\cup\\{(i,j)\\}$

14

15return $\mathcal{G^{\prime}}=\\{\mathcal{V},\mathcal{E^{\prime}}\\}$

Algorithm 1 Adding Heterophilous Edges

###  4.2 Effect of extent of heterophily

In this section, we attempt to understand the effect of variation in the
extent of heterophily $\widehat{\mathcal{H}}$ over the model performance. We
synthetically alter the topological graph
$\mathcal{G}=(\mathbf{X},\mathbf{A})$ starting from the original heterophily
$h_{init}$ up-to $95\%$ heterophily in all the six graph datasets as per
Algorithm 1. Here, $N=|\mathcal{V}|$ and $\mathcal{V}_{j}$ is the set of nodes
with label $y_{j}$.

Fig. 2: Effect of increasing heterophily on homophilous (first column) and
heterophilous (second column) datasets.

Figure 2 shows the performance variation with increasing heterophily. Each
point on plot in Figure 2 represents the performance of the proposed framework
on the generated graph corresponding to the heterophily value in x-axis with
the first point denoting the original graph, i.e, $K=0$. As depicted in Figure
2, as $K$ increases (thus, the heterophily) the classification performance
first decreases and then eventually begins to increase (except on the CORA
dataset). This shows that the the proposed GCN framework can perform under
higher heterophily as well. The shape of the curve is attributed to the “phase
transition” where at the beginning the initial topology is disturbed by added
edges making the topological space less informative. However, as we add more
edges, the network adapts itself to look for heterophilous properties in the
feature space and starts showing better performance. While this is similar to
that observed in [9] for specific kind of heterophily, our framework applies
to a general case by adapting itself to increasing heterophily. However, the
performance on CORA dataset reduces as $K$ increases. This can be attributed
to the fact that randomly adding heterophilous edges to a highly homophilous
graph can sometimes not be very meaningful and can turn out to be noisy,
similar to adding an edge between nodes that are actually never related in
real-world.

###  4.3 Analysis of attention values

In this section, we analyze the variation in attention values during the
training process. The main idea is to get an understanding of how does the
framework decide to choose between the information from topological and
feature space under different extent of heterophily. Here, we consider ACM and
Citeseer (homophilous) and Flickr and BlogCatalog (heterophilous) datasets for
visualization. Since the attention maps $(\bm{\alpha}_{T})$,
$(\bm{\alpha}_{F})$, and $(\bm{\alpha}_{C})$ are matrices, we first
L2-normalize them across rows and then compute the average of these row-wise
norms. Figure 3 shows that at the beginning, the attention values (in the
L2-norm sense) of topology $(\bm{\alpha}_{T})$, feature $(\bm{\alpha}_{F})$,
and common $(\bm{\alpha}_{C})$ embeddings are almost the same. However, as the
number of epochs increase, the attention values begin to differ. Specifically,
for homophilous datasets (Figure 3, first column), $(\bm{\alpha}_{T})$
increases more than $(\bm{\alpha}_{F})$ and vice-versa for the heterophilous
datasets (Figure 3, second column). Following these observations, we see that
for homophilous datasets, the framework focuses more on information from the
topological space and that on the feature space over heterophilous datasets.
The information from both the spaces augment network learning by adaptively
assigning appropriate attention values to information from each space and
perform well on both homophilous and heterophilous datasets.

Fig. 3: Variation of L2-normalised attention values learned over topological
and feature spaces.

##  5 Conclusion

We explore the ability of GCN in fusing information from network topological
and feature space, and additionally address heterophily in graphs. While there
is a method that shows GCNs in their original form addressing certain specific
type of heterophily, others have advocated that GCNs are implicitly designed
with homophily in mind. Further, they have not been suitable for heterophilous
graphs unless their architectures are modified and redesigned specifically to
address heterophily. We have shown that the proposed framework is able to
handle heterophily across different real-world datasets with varied amount of
heterophily not restricted solely to certain specific types. In the near
future, we would like to explore whether such a multi-channel framework can be
used to predict the missing node features in a graph and extend to dynamic
graphs.

## References

  * [1] Muhan Zhang, Zhicheng Cui, Marion Neumann, and Yixin Chen,  “An end-to-end deep learning architecture for graph classification,”  in Proceedings of the AAAI conference on artificial intelligence, 2018, vol. 32.
  * [2] Jiaxuan You, Rex Ying, and Jure Leskovec,  “Position-aware graph neural networks,”  in International conference on machine learning. PMLR, 2019, pp. 7134–7143.
  * [3] Rex Ying, Ruining He, Kaifeng Chen, Pong Eksombatchai, William L Hamilton, and Jure Leskovec,  “Graph convolutional neural networks for web-scale recommender systems,”  in Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2018, pp. 974–983.
  * [4] Zhiwen Chen, Jiamin Xu, Tao Peng, and Chunhua Yang,  “Graph convolutional network-based method for fault diagnosis using a hybrid of measurement and prior knowledge,”  IEEE transactions on cybernetics, 2021.
  * [5] Xiaolong Fan, Maoguo Gong, Yue Wu, AK Qin, and Yu Xie,  “Propagation enhanced neural message passing for graph representation learning,”  IEEE Transactions on Knowledge and Data Engineering, 2021.
  * [6] Long Zhao, Xi Peng, Yu Tian, Mubbasir Kapadia, and Dimitris N Metaxas,  “Semantic graph convolutional networks for 3d human pose regression,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2019, pp. 3425–3435.
  * [7] Jiong Zhu, Yujun Yan, Lingxiao Zhao, Mark Heimann, Leman Akoglu, and Danai Koutra,  “Beyond homophily in graph neural networks: Current limitations and effective designs,”  Advances in Neural Information Processing Systems, vol. 33, pp. 7793–7804, 2020.
  * [8] Jiong Zhu, Ryan A Rossi, Anup Rao, Tung Mai, Nedim Lipka, Nesreen K Ahmed, and Danai Koutra,  “Graph neural networks with heterophily,”  in Proceedings of the AAAI Conference on Artificial Intelligence, 2021, vol. 35, pp. 11168–11176.
  * [9] Yao Ma, Xiaorui Liu, Neil Shah, and Jiliang Tang,  “Is homophily a necessity for graph neural networks?,”  arXiv preprint arXiv:2106.06134, 2021.
  * [10] Xiao Wang, Meiqi Zhu, Deyu Bo, Peng Cui, Chuan Shi, and Jian Pei,  “Am-gcn: Adaptive multi-channel graph convolutional networks,”  in Proceedings of the 26th ACM SIGKDD International conference on knowledge discovery & data mining, 2020, pp. 1243–1253.
  * [11] Miller McPherson, Lynn Smith-Lovin, and James M Cook,  “Birds of a feather: Homophily in social networks,”  Annual review of sociology, pp. 415–444, 2001.
  * [12] Elisabeth R Gerber, Adam Douglas Henry, and Mark Lubell,  “Political homophily and collaboration in regional planning networks,”  American Journal of Political Science, vol. 57, no. 3, pp. 598–610, 2013.
  * [13] Mark Newman,  Networks,  Oxford university press, 2018.
  * [14] Valerio Ciotti, Moreno Bonaventura, Vincenzo Nicosia, Pietro Panzarasa, and Vito Latora,  “Homophily and missing links in citation networks,”  EPJ Data Science, vol. 5, pp. 1–14, 2016.
  * [15] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun,  “Spectral networks and locally connected networks on graphs,”  arXiv preprint arXiv:1312.6203, 2013.
  * [16] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst,  “Convolutional neural networks on graphs with fast localized spectral filtering,”  Advances in neural information processing systems, vol. 29, 2016\.
  * [17] Thomas N Kipf and Max Welling,  “Semi-supervised classification with graph convolutional networks,”  arXiv preprint arXiv:1609.02907, 2016.
  * [18] William L Hamilton, Rex Ying, and Jure Leskovec,  “Inductive representation learning on large graphs,”  arXiv preprint arXiv:1706.02216, 2017.
  * [19] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio,  “Graph attention networks,”  arXiv preprint arXiv:1710.10903, 2017.
  * [20] Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger,  “Simplifying graph convolutional networks,”  in International conference on machine learning. PMLR, 2019, pp. 6861–6871.
  * [21] Meng Qu, Yoshua Bengio, and Jian Tang,  “Gmnn: Graph markov neural networks,”  in International conference on machine learning. PMLR, 2019, pp. 5241–5250.
  * [22] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip,  “A comprehensive survey on graph neural networks,”  IEEE transactions on neural networks and learning systems, vol. 32, no. 1, pp. 4–24, 2020.
  * [23] Ziwei Zhang, Peng Cui, and Wenwu Zhu,  “Deep learning on graphs: A survey,”  IEEE Transactions on Knowledge and Data Engineering, 2020.
  * [24] Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang,  “Geom-gcn: Geometric graph convolutional networks,”  arXiv preprint arXiv:2002.05287, 2020.
  * [25] Benjamin Paul Chamberlain, James Rowbottom, Maria Gorinova, Stefan Webb, Emanuele Rossi, and Michael M Bronstein,  “Grand: Graph neural diffusion,”  arXiv preprint arXiv:2106.10934, 2021.
  * [26] Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li,  “Simple and deep graph convolutional networks,”  in International Conference on Machine Learning. PMLR, 2020, pp. 1725–1735.
  * [27] Ashish Tiwari, Richeek Das, and Shanmuganathan Raman,  “Exploring deeper graph convolutions for semi-supervised node classification,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 5463–5467.


# The Propagation of Strong Shocks into Planetary and Stellar Atmospheres

Almog Yalinewich  1 Andrey Remorov  1 1Canadian Institute for Theoretical
Astrophysics  60 St. George St  Toronto  ON M5S 3H8  Canada

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

###### Abstract

In this work we present a mathematical model for the propagation of the shock
waves that occur in graded density profiles. These waves can occur in a wide
range of astrophysical events, such as collisions in planetary and stellar
atmospheres, common envelope explosions and peculiar type Ia supernovae. The
behaviour of the shock wave and its evolution can be modelled using type II
self similar solutions. In such solutions the evolution of the shock wave is
determined by boundary conditions at the shock front and a singular point in
the shocked region. We show how the evolution can be determined for different
equations of state and density profiles, and compare these results to
numerical simulations. These findings are also applied to a variety of
astrophysical phenomena to further test their validity.

###### keywords:

shock waves – hydrodynamics – stars: atmospheres

††pubyear: 2021††pagerange: The Propagation of Strong Shocks into Planetary
and Stellar Atmospheres–8

##  1 Introduction

The problem of a surface explosion has been considered by many authors in the
past, using numerical simulations (Potter et al., 2012; Monteux & Arkani-
Hamed, 2019) and lab experiments (Holsapple & Housen, 2012; Mazzariol & Alves,
2019; Holsapple, 1987). The primarily incentive for these works was to study
cratering events due to planetary impacts. Among the terrestrial impact events
considered are the Chicxulub Cretaceous/Tertiary impact (Pope et al., 1997;
Newman et al., 1999) and the Tunguska event (Artemieva & Shuvalov, 2016).
Related phenomena occurring on other astronomical bodies are being explored as
well. Some of these include asteroid impacts on the atmosphere of Venus
(Korycansky et al., 2000), meteoroid collisions causing particle ejections
from the Bennu asteroid (Bottke et al., 2020), and various impacting objects
colliding with Jupiter causing flashes of certain light curves (Hueso et al.,
2013), the most famous of which is the light curve generated by the Shoemaker
Levy 9 collision (Zahnle & Mac Low, 1995). That being said, subsurface
explosions are not restricted to planets, and can also occur in star. For
example, it has been suggested that the 1954 precursor to the peculiar
supernova iptf2014hls was caused by an explosion in a common envelope
(Yalinewich & Matzner, 2019).

In this work, we develop an analytic formalism to model the propagation of
strong shocks into an increasing density profiles, such as planetary and
stellar atmospheres. For simplicity, we will assume that the explosion is so
strong that we can neglect the ambient pressure, and that the radius is much
larger than the size of the initial hot spot. We will also neglect gravity and
assume an ideal gas equation of state for all materials. Under these
assumptions, we can use the method of self similarity to describe the
evolution of the shock wave (Barenblatt, 2014). The propagation of such a
shock wave will be markedly different from the propagation of a shock wave
into a medium with flat or shallow density profile. The latter can be
described by the celebrated Sedov Taylor solution (Sedov, 1946; Taylor, 1950),
where behaviour of the shock is governed by energy conservation. In contrast,
in the problem we discuss the energy available to drive the shock decreases
with time. This is because shocked fluid elements expand and accelerate away
from the shock front, and at some point cross a critical position called the
sonic point where information cannot travel back to the shock front. For this
reason, the behaviour of this shock is determined by the behaviour near this
sonic point. Such solutions to the hydrodynamic equations are termed “type
II”, to distinguish them from the “type I” solutions, which are governed by
energy conservation. These type II solutions have previously been used to
model accelerating shocks propagating into a steeply declining density
profiles, (Waxman & Shvarts, 1992; Schlichting, 2014). In contrast, in this
work we consider a decelerating type II shocks propagating into an increasing
density profile.

Our model is a novel generalisation to the classical impulsive piston problem
(Adamskii, 1956; Yalinewich & Remorov, 2020). In the classical problem, a thin
wafer collides with a semi infinite space filled with a cold gas with a
uniform density. One can consider this problem as the one dimensional analogue
of an impact event. Since this problem is one dimensional and slab symmetric,
it can be solved using the self similar method. The solution predicts a power
law relation between the shock velocity and the swept up mass $v\propto
m^{-\beta}$. Interestingly, the same power law relation also holds well in the
three dimensional case (i.e. a sphere hitting a flat surface head on). In this
study we consider a variation on the classical problem where the density of
target can vary as some power law of distance from the impact site
$\rho\propto x^{\omega}$. Not only will this description allow us to consider
an atmospheric graded density profile, but it will also help us understand why
the one dimensional analogue describes the three dimensional case well. This
is because one way to think about the difference between the slab symmetric
and three dimensional cases is how the mass increases with the distance. For
example, if the density is uniform, then in slab geometry the mass increases
linearly with the distance while in three dimensions the mass increases as a
cube of the distance. These differences can be accounted for in the slab
symmetric analogue by changing the density profile in the target.

The plan of the paper is as follows. In section 2 we show how self similarity
can be used to solve the generalised impulsive piston problem. In section 3 we
demonstrate how the solutions we found can be used to model a variety of
astrophysical problems. In section 4 we discuss the results and conclude.

##  2 Self Similar Model

###  2.1 Impulsive Piston Problem

Let us consider an explosion close to the surface of a large star or planet,
so that the radius of curvature is unimportant and the atmospheric layers can
be considered planar. This explosion will generate a shock wave that travels
into the atmosphere. While the distance between the centre of the shock and
the shock front is much smaller than the atmospheric scale height, the shock
front interacts with an ambient medium with the same density, and hence the
behaviour can be described by the classical Sedov Taylor explosion (Sedov,
1946; Taylor, 1950). When the distance between the shock front and explosion
centre becomes comparable to the scale height, the top of the shock encounters
a declining density profile and accelerates. A mathematical model describing
this part of the shock was developed in (Sakurai, 1960). In this work we focus
on the part of the shock that propagates deeper into the atmosphere, and since
the density the shock encounters increases the shock decelerates. When the
distance between the centre of the explosion and the deepest point on the
shock front is much larger than the atmospheric scale height the shock becomes
self similar. Therefore, the effective radius of the shock grows as some power
law of time $R\propto t^{\alpha}$, where We define the effective radius of the
shock as the distance between the deepest point of the shock front to the
centre of the explosion.

As was done in a previous paper (Yalinewich & Remorov, 2020), instead of
considering a three dimensional problem, we consider the one-dimensional
analogue called the impulsive piston problem (Adamskii, 1956; Zel’dovich &
Raizer, 1967). In this scenario, a thin wafer hits a much thicker slab of
material (with a certain density distribution) and both are perfectly cold
prior to the collision. As a result of this collision a shock wave emerges
from the contact surface and travels into the target. Once the shock wave has
travelled a distance much larger than the width of the wafer, there is no
other relevant length scale in this process other than the position of the
shock front, and hence we expect the shock wave to evolve in a self similar
way. In the original one dimensional impulsive piston problem, the target is
considered to have uniform density. However, in this work we consider a graded
density profile defined by a power law dependence on the distance from the
edge

$\rho_{a}=kx^{\omega}$ (1)

where $k$ and $\omega>0$ are constants and $x$ is the distance from the edge.

In the next sections we describe how self similarity allows us to solve the
impulsive piston problem. It does so by reduce the hydrodynamic equations,
which are partial differential equations, to ordinary differential equations.
We also describe how these equation can be numerically integrated, and how the
value of $\alpha$ can be obtained, as was done for the case $\omega=0$ in
(Yalinewich & Remorov, 2020).

###  2.2 Self Similar Equations

The mathematical machinery presented in the following subsections is very
similar to the methods used to analyse cratering events in planetary
collisions in a previous paper (Yalinewich & Remorov, 2020). The main
difference is that the previous work considered a uniform density profile,
compared to the graded density profile examined here.

The slab symmetric hydrodynamic equations in one dimension are given by
(Landau & Lifshitz, 1987), which are respectively the continuity, momentum,
and entropy conservation, where $\rho$ is density, $v$ is velocity, $p$ is
pressure, and $s$ is entropy:

$\frac{\partial\rho}{\partial t}+\rho\frac{\partial v}{\partial
x}+v\frac{\partial\rho}{\partial x}=0$ (2)

$\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial
x}+\frac{1}{\rho}\frac{\partial p}{\partial x}=0$ (3)

$\frac{\partial s}{\partial t}+v\frac{\partial s}{\partial x}=0\,.$ (4)

We assume an ideal gas equation of state

$s=\ln p-\gamma\ln\rho$ (5)

where $\gamma$ is the adiabatic index.

We can replace the pressure by the sound speed using the relation $p=\rho
c^{2}/\gamma$. After this substitution we have a system of partial
differential equations in $v,c$ and $\rho$. We make the assumption that the
one dimensional shock has a self similar behaviour, so the position of the
shock front $X(t)$ evolves as a power law in time

$X(t)\propto t^{\alpha}\,.$ (6)

where $\alpha$ is a constant. One of the main goals of the analysis is
determining $\alpha$. The partial differential equations can be reduced to
ordinary differential equations in the dimensionless position $\chi=x/X$.
Following (Landau & Lifshitz, 1987), we define dimensionless hydrodynamical
variables $V,C,D$ in the following way:

$v(x,t)=\frac{dX(t)}{dt}\chi V\left(\chi\right)$ (7)
$c(x,t)=\frac{dX(t)}{dt}\chi C\left(\chi\right)$ (8)
$\rho(x,t)=kX^{\omega}(t)D\left(\chi\right)\,.$ (9)

After this substitution, the equations contain terms involving different time
derivatives of $X$. These time derivatives can be factored out using equation
6. We introduce the following parameter

$\delta=\frac{\ddot{X}X}{\dot{X}^{2}}=1-\frac{1}{\alpha}\,.$ (10)

This parameter is the power law index of the shock velocity - position
relation $dX/dt\propto X^{\delta}$.

After making these adjustments and simplifying, the hydrodynamic equations
equations 2, 3 and 4 become the matrix equation
$\overleftrightarrow{M}d\vec{A}/\chi=\vec{B}$ where:

$M=\begin{bmatrix}\chi\left(V-1\right)&\chi D&0\\\ \chi
C^{2}&\gamma\chi(V-1)D&2\chi CD\\\
C\chi(\gamma-1)\left(1-V\right)&0&2D\chi\left(V-1\right)\end{bmatrix}$ (11)

$\vec{A}=\begin{bmatrix}D\\\ V\\\ C\end{bmatrix}$ (12)

$\vec{B}=\begin{bmatrix}-\omega D-DV\\\ -\delta\gamma
DV-\gamma(V^{2}-V)D-2C^{2}D\\\ -(-\gamma\omega+\omega)CD-2CDV\end{bmatrix}$
(13)

These equations are linear in the derivatives $dV/d\chi,dC/d\chi,dD/d\chi$ so
these can be isolated. Our choice of parameters allows us to go a step further
and obtain a single ordinary differential equation involving just $C$ and $V$.
We note that

$\frac{dV}{d\chi}=\frac{\tilde{N}_{1}}{\tilde{D}_{1}}$ (14)

and

$\frac{dC}{d\chi}=\frac{N_{1}}{D_{1}}$ (15)

where

$\tilde{N}_{1}=-C^{2}V\gamma-2C^{2}\delta+V^{3}\gamma+V^{2}\delta\gamma-2V^{2}\gamma-V\delta\gamma+V\gamma-\omega
C^{2}$ (16) $\tilde{D}_{1}=\chi\gamma\left(C^{2}-V^{2}+2V-1\right)$ (17)
$N_{1}=C\left((V+\omega)\gamma\left(V-1\right)\left(V\gamma-V-\gamma+1\right)-\left(C^{2}-\gamma\left(V-1\right)^{2}\right)\times\right.$
$\left.\times\left(2V+2\delta-\omega\gamma+\omega-2\right)-\left(2C^{2}+V^{2}\gamma+V\delta\gamma-V\gamma\right)\left(V\gamma-V-\gamma+1\right)\right)$
(18)
$D_{1}=2\chi\left(C^{2}\left(V-1\right)+C^{2}\left(V\gamma-V-\gamma+1\right)-\gamma\left(V-1\right)^{3}\right)$
(19)

Dividing one equation by the other, we get a single ODE

$\frac{dC}{dV}=\frac{N_{2}}{D_{2}}$ (20)

where:

$N_{2}=C\gamma\left(-(V+\omega)\gamma\left(V-1\right)\left(V\gamma-V-\gamma+1\right)+\right.$
$\left.\left(C^{2}-\gamma\left(V-1\right)^{2}\right)\left(2V+2\delta-\gamma\omega+\omega-2\right)+\right.$
$\left.+\left(2C^{2}+V^{2}\gamma+V\delta\gamma-V\gamma\right)\left(V\gamma-V-\gamma+1\right)\right)\times$
(21) $\left(C^{2}-V^{2}+2V-1\right)$

and

$D_{2}=2\left(C^{2}\left(V-1\right)+C^{2}\left(V\gamma-V-\gamma+1\right)-\gamma\left(V-1\right)^{3}\right)\times$
(22)
$\times\left(C^{2}V\gamma+2C^{2}\delta+C^{2}\omega-V^{3}\gamma-V^{2}\delta\gamma+2V^{2}\gamma+V\delta\gamma-V\gamma\right)\,.$

###  2.3 Boundary Conditions

In order to integrate equation 20, boundary conditions and the value of
$\delta$ are needed. The boundary conditions at the shock front are given by
the Rankine Hugoniot conditions for a strong shock (Zel’dovich & Raizer, 1967;
Landau & Lifshitz, 1987).

$V_{f}=\frac{2}{\gamma+1}$ (23)

and

$C_{f}=\frac{\sqrt{2\gamma\left(\gamma-1\right)}}{\gamma+1}\,.$ (24)

In some self similar problems, like the Sedov Taylor explosion, the parameter
$\delta$ can be be inferred directly from conservation laws. Such problems are
known as Type I solutions. In our case, however, $\delta$ is determined by the
condition that the hydrodynamic trajectory passes smoothly through some
singularity. These are called Type II solutions (Waxman & Shvarts, 1992). The
singularity occurs at a point where information cannot propagate back to the
shock front, and is therefore also referred to as the sonic point. Energy that
flows through the sonic point cannot travel back to the shock front, and so
the energy available to sustain the shock is not conserved.

Equation 14 has a singularity when

$C=1-V$ (25)

This curve is known as the sonic line. On the sonic line, the denominator
(equation 22) vanishes. To prevent a divergence, the numerator also has to
vanish. This happens on a specific point on the sonic line, which we call the
sonic point

$V_{s}=\frac{2\delta+\omega}{\delta(2-\gamma)+\omega}$ (26)
$C_{s}=-\frac{\delta\gamma}{\delta(2-\gamma)+\omega}$ (27)

Knowing the boundary conditions and $\frac{dC}{dV}$, we can integrate from the
sonic point to the shock front. To avoid a numerical divergence, the starting
point is slightly shifted from the sonic point in the positive $V$ direction
by $dV_{i}$, and in $C$ direction by $dV_{i}\frac{dC}{dV}\Bigr{|}_{sonic}$ for
some small $dV_{i}\ll 1$.

The correct value of $\delta$ is such that the curve $C\left(V\right)$
satisfies both boundary conditions at the sonic point and the shock front. We
use the shooting method to find the value of $\delta$. We guess a value for
$\delta$, numerically integrate equation 20 w.r.t to $V$ from the sonic point
to the shock front, and note the value of $C$ at the shock front. We then use
the bisection method to refine the value of $\delta$ to minimise the distance
between the value of $C$ obtained from numeric integration and its theoretical
value (equation 24). An example for some hydrodynamic $C\left(V\right)$
trajectories for the same value of $\gamma,\omega$ but different values of
$\delta$ can be seen in figure 1.

In the case $C_{s}>0$, numerical integration is straightforward since both the
sonic point and the shock front are on the same side of the impact site. An
example for several hydrodynamic trajectories for fixed $\gamma,\omega$ such
that $C_{s}>0$ can be seen in figure 1. The case $C_{s}<0$ is more complicated
because the integration goes through $x=0$. Because of the way we defined the
self similar variables (equation 9), they diverge at $x=0$, even though the
physical quantities $v$ and $c$ remain finite. We can circumvent this
difficulty by noting that the Mach number, i.e. the ratio between the velocity
and the speed of sound, remains finite and changes continuously across $x=0$.
The sonic point in the case $C_{s}<0$ lies in the fourth quadrant of the of
the $C-V$ plane (i.e. positive $V$ and negative $C$). Numerical integration
proceeds to even higher values of $V$, moving away from the shock front. Far
away from the sonic point, the curve attains some asymptotic slope, and this
slope is the Mach number at $x=0$, which we’ll denote by $\mathcal{M}_{0}$.
Once we have this piece of information, we can restart the numerical
integration in the second quadrant ($C$ is positive and $V$ is negative),
beginning from some arbitrarily highly negative $V=V_{r}$ (subscript r for
restart) and $C_{r}=V_{r}/\mathcal{M}_{0}$. From this point we can continue
the numerical integration all the way to the shock front. An example for
several hydrodynamic trajectories for fixed $\gamma,\omega$ such that
$C_{s}<0$ can be seen in figure 2.

Figure 1: First type of integration, for $C_{s}>0$. Here $\gamma=3.5$ and
$\omega=3$. The numbers in the legend show the values of $\delta$ used in each
integration. This integration is straightforward and there are no coordinate
singularities.

Figure 2: Second type of integration, for $C_{s}<0$, and
$\gamma=1.47,\omega=0.3$. In this situation, the shock and sonic points are on
opposite sides of the boundary separating the fluid and the vacuum, so a
mathematical singularity occurs when integrating from the sonic point to the
shock front.

In this manner we calculated $\delta$ for different values of $\gamma$ and
$\omega$. The results are displayed in the top panel of figure 3. This figure
shows that $\delta$ is mostly dependent on $\omega$, and weakly dependent on
$\gamma$. For completeness, the asymptotic cases of $\gamma\rightarrow\infty$
; $\omega\rightarrow\infty$ ; and
$\omega\rightarrow\infty,\gamma\rightarrow\infty$ are covered in the appendix.

This result prompted us to consider another power law index, $\beta=d\ln
v/d\ln m$, which relates the velocity of the shock to the swept up mass. This
parameter can be related to $\delta$ via

$\beta=\frac{\delta}{1+\omega}$ (28)

A contour plot of $\beta(\gamma,\omega)$ is shown in figure 3, in which one
can see that $\beta$ does not vary much as a function of $\gamma$ and
$\omega$.

Conservation of energy and momentum laws allows us to set upper and lower
limits on the power law parameters (Zel’dovich & Raizer, 1967)

$-\frac{1}{2}>\beta>-1$ (29) $-\frac{\omega+1}{2}>\delta>-\omega-1$ (30)

and

$\frac{2}{3+\omega}>\alpha>\frac{1}{2+\omega}\,.$ (31)

###  2.4 Comparison with Simulations

In order to verify our results, we compare them to the numerical simulations
performed in (Yalinewich & Matzner, 2019). The authors used hydrodynamic
simulations to obtain the original power law parameter $\alpha$, defined in 6.
They obtained $\alpha\left(\gamma=5/3,\omega=3\right)=0.19$ and
$\alpha\left(\gamma=5/3,\omega=3/2\right)=0.25$. We note that these
calculations were performed in spherical geometry, so in order to calculate
$\beta$ we use

$\beta=\frac{\delta}{3+\omega}=\frac{1-1/\alpha}{3+\omega}\,.$ (32)

Using equation 32 we calculate the numerical values
$\beta\left(\gamma=5/3,\omega=3\right)=-0.71$ and
$\beta\left(\gamma=5/3,\omega=3/2\right)=-0.67$. Using the numerical
integration methods, we obtained
$\beta\left(\gamma=5/3,\omega=3\right)=-0.671$ and
$\beta\left(\gamma=5/3,\omega=3/2\right)=-0.663$. Despite using different
geometric setups and approaches to determining the power law scaling
parameters, the values for $\beta$ obtained here are reasonably close to those
obtained previously with hydrodynamic simulations.

Figure 3: Contour plots of $\delta(\gamma,\omega)$ and $\beta(\gamma,\omega)$.

##  3 Applications

A number of previous studies have already considered the role of shock waves
from a surface explosions in different astrophysical phenomena. (Yalinewich &
Schlichting, 2019) studied atmospheric mass loss from planetary collisions.
(Yalinewich & Matzner, 2019) studied the optical transient that results from a
rapid release of energy by a neutron star while inside the envelope of a giant
companion star. Finally, (Yalinewich & Caplan, 2021) showed that the shapes of
craters due to collisions with primordial black holes would differ from
typical craters, and hence the moon can serve as a dark matter detector.

In addition to these cases, in this section we present more possible
astrophysical applications to the formalism developed here. In particular, we
discuss how the analytical formalism presented above can reproduce results
from numerical hydrodynamic simulations and observations.

###  3.1 Bolides

Even when objects burn up completely in the atmosphere, they produce a
pressure wave that can cause damage on the ground. The most recent example is
the Chelyabinsk event (Emel’yanenko et al., 2013). The damage caused by the
explosion (measured by the reports on damage per capita) declined with
distance according to $r^{-2.4}$, similar to the measured decline in over-
pressure $r^{-2.6}$. We can interpret this decline in pressure using the model
developed above. After the shock wave reaches the ground, it expands further
only in the horizontal direction. For this reason the swept up mass scales
with the area $m\propto r^{2}$. The density at ground level is roughly
constant, and so the pressure scales as $p\propto\rho v^{2}\propto
m^{-2\beta}\propto r^{-4\beta}$. Choosing $\beta=2/3$, we find $p\propto
r^{-2.7}$, which is relatively close to the empirical result.

The damage caused by the shock depends on whether the shock is strong by the
time it reaches the ground, i.e. whether the shock pressure is much larger
than the atmospheric pressure when the shock reaches the ground.

As the bolide moves downward, it interacts with the atmosphere, decelerates
and heats up. This heat causes it to expand, which increases the deceleration
and heating rate and so forth. The air density at the the altitude at which
the bolide decelerates considerably is (Chyba et al., 1993)

$\rho_{d}\approx\rho_{i}\left(\frac{R_{i}}{h}\right)^{2}$ (33)

where $\rho_{i}$ is the density of the impactor, $R_{i}$ its radius and $h$ is
the scale height. The bolide then deposits the energy in the gas and creates a
shock wave that expands outward. At early times all parts of the shock
interact with roughly the same density, so the shock evolves as a classical
Sedov Taylor solution (Taylor, 1950; Sedov, 1946). The difference between
density at different points on the shock becomes significant when the shock
radius is comparable to the scale height. At that point the swept up mass is

$M_{ST}\approx\rho_{d}h^{3}\approx\rho_{i}R_{i}^{2}h\,.$ (34)

Since energy is conserved, the velocity at the end of the Sedov Taylor phase
is

$v_{s}\approx v_{i}\sqrt{\frac{m_{i}}{M_{ST}}}\approx
v_{i}\sqrt{\frac{R_{i}}{h}}$ (35)

where $m_{i}\approx\rho_{i}R_{i}^{3}$ is the mass of the impactor. As the
shock moves to lower altitudes, it follows the solution found in the previous
section, so the shock velocity on the ground is

$v_{g}\approx v_{ST}\left(\frac{M_{ST}}{\rho_{0}h^{3}}\right)^{\beta}\approx
v_{i}\sqrt{\frac{R_{i}}{h}}\left(\frac{R_{i}^{2}\rho_{i}}{h^{2}\rho_{0}}\right)^{\beta}\,.$
(36)

where $\rho_{0}$ is the air density on the ground so $\rho_{0}h^{3}$ is the
air mass enclosed within a radius $h$ from any point on the ground. The
pressure on the ground is

$p_{g}\approx\rho_{0}v_{g}^{2}\approx\rho_{0}v_{i}^{2}\frac{R_{i}}{h}\left(\frac{R_{i}^{2}\rho_{i}}{h^{2}\rho_{0}}\right)^{2\beta}\,.$
(37)

By comparing this pressure to the atmospheric pressure at ground level we can
obtain a critical value for the impactor radius

$R_{i}\approx
h\left(\frac{\rho_{0}}{\rho_{i}}\right)^{\frac{2\beta}{4\beta+1}}\left(\frac{p_{0}}{\rho_{0}v_{i}^{2}}\right)^{\frac{1}{4\beta+1}}\,.$
(38)

Substituting typical values and using $\beta=2/3$ yields

$R_{i}\approx
66\tilde{h}\tilde{\rho}_{0}^{1/11}\tilde{\rho}_{i}^{-4/11}\tilde{p}_{0}^{3/11}\tilde{v}_{i}^{-6/11}\,\rm
m$ (39)

where $\tilde{h}=h/8\rm\,km$,
$\tilde{\rho}_{0}=\rho_{0}/10^{-3}\rm\,\frac{g}{cm^{3}}$,
$\tilde{\rho}_{i}=\rho_{i}/3\rm\,\frac{g}{cm^{3}}$,
$\tilde{p}_{0}=p_{0}/10^{5}\,\rm Pa$ and
$\tilde{v}_{i}=v_{i}/10\,\rm\frac{km}{s}$. A impactor with a larger radius
creates a shock that remains strong when it reaches the ground. In numerical
simulations, the critical impactor size is around 50 metres (Shuvalov et al.,
2013). Incidentally, the over-pressure threshold for substantial damage to
brick and concrete surfaces is also comparable to the ambient atmospheric
pressure (Glasstone & Dolan, 1977). Therefore, if the shock wave is no longer
strong by the time it reaches the ground, it is unlikely to cause substantial
damage to buildings. Weak shocks can still cause damage to more brittle
materials such as glass and wood, as was the case with the Chelyabinsk bolide.

Finally, it is possible to estimate the atmospheric mass loss from bolides. An
impactor moving faster than the escape velocity will create a decelerating
shock wave, so that after it has swept up a certain amount of mass the
velocity will drop below the escape velocity. Since the velocity does not
change considerably after the shock, the amount of mass lost in this way will
be comparable to the mass accelerated to velocities exceeding the escape
velocity $v_{e}$. Hence, we can approximate the amount of atmospheric mass
loss $M$ by

$M\approx
m_{i}\left(\frac{h}{R_{i}}\right)^{\frac{2\beta-1}{\beta}}\left(\frac{v_{i}}{v_{e}}\right)^{1/\beta}\,.$
(40)

Verification of this relation and discussion of its implication for planet
formation is relegated to a future work.

###  3.2 Shoemaker Levy 9

The impact of the fragments of the comet Shoemaker Levy 9 produced bright
visible and infrared flashes that lasted, depending on the band, between tens
of seconds and tens of minutes. A direct comparison between the observational
data and theory is problematic because of a number of reasons. First, some of
the fragments landed on the far side of Jupiter, so that only the spacecraft
Galileo was able to observe the early phases of the impact (Neukum et al.,
1995), while ground based instruments could only begin to observe after the
plumes rose above the horizon (although some did detect a faint precursor due
to the passage of the fragments through the tenuous upper layers of Jupiter’s
atmosphere). Second, some of the ground based telescopes became saturated and
thus missed the peak (Graham et al., 1995). Nevertheless, such a comparison
has been performed (Zahnle & Mac Low, 1995), using a similar model for the
shock propagation, but without the analytic justification. We will not repeat
that analysis here, but instead focus on the analytic predictions that can be
obtained from our formalism. Also, for simplicity in this section we assume
that $\beta=2/3$.

Most of the light observed after the impact of Shoemaker Levy 9 on Jupiter did
not originate from the bolide itself, but from ejecta from the impact site
that collided with more distant parts on the surface of Jupiter (Zahnle & Mac
Low, 1995). The range of every fluid shell depends on the velocity with which
it emerges from the basin

$r\approx\frac{v^{2}}{g}\approx\frac{R_{i}^{\frac{10}{3}}h^{\frac{2}{3}}\rho_{i}^{\frac{4}{3}}v_{i}^{2}}{gm^{\frac{4}{3}}}\,.$
(41)

where $g$ is the Jovian surface gravity. The time it takes a fluid shell to
fall down to the surface is

$t_{f}\approx\frac{v}{g}\approx\frac{R_{i}^{\frac{5}{3}}\sqrt[3]{h}\rho_{i}^{\frac{2}{3}}v_{i}}{gm^{\frac{2}{3}}}\,.$
(42)

At any moment an observer can only see fluid elements that had time to fall
$t>t_{f}$. When the ejecta collides with the surface, it is shocked, and its
temperature is given by

$T\approx\mu v^{2}/k\approx\frac{g^{2}\mu t^{2}}{k}\,.$ (43)

where $\mu$ is the atomic mass and $k$ is the Boltzmann constant. The
luminosity is given by blackbody emission

$L\approx\sigma T^{4}r^{2}\approx\frac{g^{10}\mu^{4}\sigma t^{12}}{k^{4}}\,.$
(44)

We find that the luminosity rises very steeply, and so attains its peak value
almost instantaneously, in accordance with observations (Chapman, 1996). This
peak luminosity comes from the fastest shell, which is the Sedov mass. The
corresponding peak bolometeric luminosity is

$L_{p}\approx\frac{R_{i}^{12}\mu^{4}\sigma
v_{i}^{12}}{g^{2}h^{12}k^{4}}\approx
10^{22}\tilde{R}_{i}^{12}\tilde{\mu}^{4}\tilde{v}_{i}^{12}\tilde{g}^{-2}\tilde{h}^{-12}\,\rm
erg/s$ (45)

where $\tilde{R}_{i}=R_{i}/24\,\rm km$, $\tilde{v}_{i}=v_{i}/60\rm\,km/s$,
$\tilde{g}=g/60\rm m/s^{2}$, $\tilde{h}=h/24\,\rm km$, $\tilde{\mu}=\mu/m_{p}$
and $m_{p}$ is the proton mass. This result is a factor of a few larger than
the observed luminosity (Chapman, 1996), however the results are very
sensitive to the parameters. The duration of the transient will be comparable
to the fallback time

$t_{p}\approx\frac{v_{i}R_{i}}{gh}\approx
40\tilde{v}_{i}\tilde{R}_{i}\tilde{g}^{-1}\tilde{h}^{-1}\,\rm s.$ (46)

This estimate is also in agreement with the observations. The temperature at
the peak is given by

$T_{p}\approx\frac{R_{i}^{2}\mu v_{i}^{2}}{h^{2}k}\approx
750\tilde{R}_{i}^{2}\tilde{\mu}\tilde{v}_{i}^{2}\tilde{h}^{-2}\,\rm K\,.$ (47)

The observed temperature is larger than our result by a factor of three
(Nicholson et al., 1995).

###  3.3 Asteroid Impact Avoidance

The impact of rocky bodies on Earth poses a danger to the life and safety of
its inhabitants. For this reason, over the years different strategies for
asteroid impact avoidance have been studied (National Research Council, 2010,
and references therein). Asteroids whose diameter is less than a few hundred
meters can be deflected using kinetic impactors (sometimes called mass
drivers). The net change in the momentum of the target is larger than just the
momentum of the impactor, since upon impact debris from the target fly in the
opposite direction. We can use the formalism above to estimate the net change
in the momentum of the target.

Let us consider an impactor of mass $m$ that collides with a much larger
target at a velocity $v_{i}$. As a result, a shock wave emerges from the
impact site and excavates a crater. Crater excavation stops when the shock
velocity drops below some critical velocity $v_{e}$, which is comparable to
the velocity with which elastic waves propagate. The mass of the excavated
material is

$M\approx m\left(\frac{v_{i}}{v_{e}}\right)^{1/\beta}\,.$ (48)

Most of the ejected material moves with a velocity $v_{e}$, so the net change
in momentum is $Mv_{e}$. For convenience, we normalise this change in target
momentum with the momentum of the incident impactor

$\xi\approx\frac{Mv_{e}}{mv_{i}}\approx\left(v_{i}/v_{e}\right)^{\frac{1-\beta}{\beta}}\,.$
(49)

For $\beta=2/3$, we get $\xi\propto v_{i}^{0.5}$, which is similar to the
behaviour found in numerical simulations $\xi\propto v_{i}^{0.44}$ (Bruck Syal
et al., 2016).

Larger asteroid can only be deflected using nuclear surface explosions (Ahrens
& Harris, 1992). We can use our formalism to relate the energy and depth of
the explosion to the net momentum change. If the initial depth is $h$ and the
ambient density of the asteroid is $\rho_{0}$, then the shock wave conserves
energy until it sweeps up a mass $\rho_{0}h^{3}$, at which point the shock
velocity is $\sqrt{E/\rho_{0}h^{3}}$. At later times the shock evolves
according to our new formalism. Thus, the mass of the ejecta is

$M\approx\rho_{0}h^{3}\left(\sqrt{\frac{E}{\rho_{0}h^{3}v_{e}^{2}}}\right)^{1/\beta}$
(50)

where $v_{e}$ is the escape velocity from the asteroid. The net change
momentum is

$I\approx v_{e}M\approx
v_{e}\rho_{0}h^{3}\left(\sqrt{\frac{E}{\rho_{0}h^{3}v_{e}^{2}}}\right)^{1/\beta}\propto
h^{3-3/2\beta}\,.$ (51)

If $\beta=2/3$, we have $I\propto h^{0.75}$, which is similar to the
predictions based on empirical scaling laws (Meshcheryakov & Lipnitskii,
2015). If the explosion is too deep, then by the time the explosion makes it
to the surface the shock velocity drops below the escape velocity and the net
change in momentum is zero. The optimal depth is determined by the condition
that the shock is at the escape velocity when it reaches the surface

$h_{o}\approx\left(\frac{E}{\rho_{0}v_{e}^{4}}\right)^{1/3}\,.$ (52)

###  3.4 Hypervelocity White Dwarfs

Recent observations revealed an intriguing class of celestial objects. These
are white dwarfs with velocities that are markedly different than the other
stars in their neighbourhoods, and also feature unusual surface composition
(Kepler et al., 2016; Shen et al., 2018). It has been suggested that these are
white dwarfs that underwent some episode of runaway thermonuclear burning that
stopped before it consumed the entire white dwarf (Gänsicke et al., 2020). If
the reaction did not stop, then the white dwarf would have explode as a type
Ia supernova and leave no remnant behind. For this reason, a partially burned
white dwarf is sometimes referred to as a supernova survivor. It has also been
suggested that these events give rise to subluminous supernovae (Foley et al.,
2013).

Many models have been proposed to explain type Ia supernova (Hillebrandt &
Niemeyer, 2000), and in most of them, the explosion is initiated at a point
that is offset from the centre. Numerical simulations have shown that
sometimes the reaction will only run through a portion of the white dwarf
rather and not destroy it entirely (Jordan et al., 2012). We can use the
formalism developed here to calculate the aftermath of such a failed
supernova. For simplicity, let us consider the white dwarf to be a sphere of
uniform density $\rho$. Suppose further that the runaway reaction has burned
through a hemisphere of radius $h$ below the surface. The mass of the burnt
material is therefore $\rho h^{3}$. We also assume that energy released per
unit mass is $\varepsilon$, so the burnt material is accelerated to a velocity
$\sqrt{\varepsilon}$. From that hot spot a shock wave emerges and travels
through the white dwarf, but does not ignite more material. If the escape
velocity from the white dwarf is $v_{e}$, then the net mass loss is

$\Delta M\approx\rho
h^{3}\left(\frac{\varepsilon}{v_{e}^{2}}\right)^{1/2\beta}\,.$ (53)

The mass ejected from the explosion moves at the escape velocity, so the
velocity kick to the survivor is

$\frac{\Delta
v}{v_{e}}\approx\left(\frac{h}{R}\right)^{3}\left(\frac{\varepsilon}{v_{e}^{2}}\right)^{1/2\beta}$
(54)

where $R$ is the radius of the white dwarf. For $\beta=2/3$ we find

$\frac{\Delta
M}{M_{WD}}\approx\left(\frac{h}{R}\right)^{3}\left(\frac{\varepsilon}{v_{e}^{2}}\right)^{3/4}\approx
0.05\left(\frac{h/R}{0.3}\right)^{3}\left(\varepsilon/\frac{{\rm
MeV}}{m_{p}}\right)^{0.75}\left(\frac{v_{e}}{6\cdot
10^{3}\,\rm\frac{km}{s}}\right)^{-1.5}$ (55)

where $M_{WD}$ is the mass of the white dwarf and $m_{p}$ is the proton mass.
The kick velocity is

$\Delta v\approx
340\left(\frac{h/R}{0.3}\right)^{3}\left(\varepsilon/\frac{\rm
MeV}{m_{p}}\right)^{0.75}\left(\frac{v_{e}}{6\cdot
10^{3}\,\rm\frac{km}{s}}\right)^{-0.5}\,\rm\frac{km}{s}\,.$ (56)

These values are similar to those found in numerical simulations (Jordan et
al., 2012).

##  4 Conclusion

In this paper, we develop a self similar model for shocks propagating down
stellar or planetary atmospheres with graded density profiles as a result of a
rapid, localised deposition of energy. We achieve this by generalising the
impulsive piston problem to account for a graded density profile. In this
model, instead of considering the full three dimensional problem, we consider
a one dimensional slab symmetric analogue. In the analogue problem, a thin
wafer hits a half space target whose density is allowed to vary in the
direction normal to the surface. The evolution of the resulting shock is
determined by the conditions at the sonic point, a singularity beyond which
information cannot travel to the shock front.

Under these assumptions, we can apply the self similar method to solve the
hydrodynamic equations (Landau & Lifshitz, 1987). We do so by reducing the
partial differential flow equations to a single ordinary differential
equation. We determine the deceleration parameter of the shock
$\beta=d\ln\dot{X}/d\ln m$ (where $\dot{X}$ is the shock velocity and $m$ is
the swept up mass) by imposing boundary conditions at the shock front and at
the sonic point, and using the shooting method.

We find that the parameter $\beta$ is weakly dependent on the density profile.
For this reason, we can use this scaling relation even for different
geometries, since the main difference between slab, cylindrical and spherical
shock is the relation between the swept up mass and the distance from the
centre of the explosion.

We also demonstrate that the method we developed can be used to study a wide
range of astrophysical scenarios. These scenarios include the transient
resulting from the impact of meteors on Earth and other planets, mass loss due
to giant planetary collisions, crater morphology, explosions in a common
envelope and failed type Ia supernovae.

## Acknowledgements

We would like to thank Jonathan Fortney for suggesting the title for the
paper. AY is supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC), funding reference #CITA 490888-16. This work made
use of the sympy (Meurer et al., 2017), numpy (Oliphant, 2015) and matplotlib
(Hunter, 2007) python packages.

## Data Availability

The equations and results for Section 2 and the appendix can be found in: .

## References

  * Adamskii (1956) Adamskii V. B., 1956, Sov. Phys. Acoust, 2
  * Ahrens & Harris (1992) Ahrens T. J., Harris A. W., 1992, Deflection and fragmentation of near-Earth asteroids, doi:10.1038/360429a0, https://www.nature.com/articles/360429a0
  * Artemieva & Shuvalov (2016) Artemieva N. A., Shuvalov V. V., 2016, Annual Review of Earth and Planetary Sciences, 44, 37
  * Barenblatt (2014) Barenblatt G. I., 2014, Scaling, self-similarity, and intermediate asymptotics. Cambridge University Press, doi:10.1017/CBO9781107050242, /core/books/scaling-selfsimilarity-and-intermediate-asymptotics/3B56096C3B7E822794C81B51F7370B82
  * Bottke et al. (2020) Bottke W. F., et al., 2020, Journal of Geophysical Research: Planets, 125, e06282
  * Bruck Syal et al. (2016) Bruck Syal M., Michael Owen J., Miller P. L., 2016, Icarus, 269, 50
  * Chapman (1996) Chapman C. R., 1996, International Astronomical Union Colloquium, 156, 121
  * Chyba et al. (1993) Chyba C. F., Thomas P. J., Zahnle K. J., 1993, Nature, 361, 40
  * Emel’yanenko et al. (2013) Emel’yanenko V. V., et al., 2013, Solar System Research, 47, 240
  * Foley et al. (2013) Foley R. J., et al., 2013, Astrophysical Journal, 767, 57
  * Gänsicke et al. (2020) Gänsicke B. T., Koester D., Raddi R., Toloza O., Kepler S. O., 2020, Monthly Notices of the Royal Astronomical Society, 496, 4079
  * Glasstone & Dolan (1977) Glasstone S., Dolan P. J., 1977, Technical report, The Effects of Nuclear Weapons. Third edition. Technical Information Center, U.S. Department of Energy, doi:10.2172/6852629
  * Graham et al. (1995) Graham de Pater I., Jernigan J. G., Liu M. C., Brown M. E., 1995, Science, 267, 1320
  * Hillebrandt & Niemeyer (2000) Hillebrandt W., Niemeyer J. C., 2000, Annual Review of Astronomy and Astrophysics, 38, 191
  * Holsapple (1987) Holsapple K. A., 1987, International Journal of Impact Engineering, 5, 343
  * Holsapple & Housen (2012) Holsapple K. A., Housen K. R., 2012, Icarus, 221, 875
  * Hueso et al. (2013) Hueso R., et al., 2013, Astronomy and Astrophysics, 560, 55
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Jordan et al. (2012) Jordan G. C., Perets H. B., Fisher R. T., van Rossum D. R., 2012, Astrophysical Journal Letters, 761
  * Kepler et al. (2016) Kepler S. O., Koester D., Ourique G., 2016, Science, 352, 67
  * Korycansky et al. (2000) Korycansky D. G., Zahnle K. J., Law M. M. M., 2000, Icarus, 146, 387
  * Landau & Lifshitz (1987) Landau L. D., Lifshitz E. M., 1987, Fluid Mechanics, doi:10.1007/b138775
  * Mazzariol & Alves (2019) Mazzariol L. M., Alves M., 2019, International Journal of Impact Engineering, 133, 103364
  * Meshcheryakov & Lipnitskii (2015) Meshcheryakov S. A., Lipnitskii Y. M., 2015, Technical Physics, 60, 26
  * Meurer et al. (2017) Meurer A., et al., 2017, PeerJ Computer Science, 3, e103
  * Monteux & Arkani-Hamed (2019) Monteux J., Arkani-Hamed J., 2019, Icarus, 331, 238
  * National Research Council (2010) National Research Council 2010, Defending Planet Earth: Near-Earth-Object Surveys and Hazard Mitigation Strategies: Final Report
  * Neukum et al. (1995) Neukum G., et al., 1995, in European SL-9/Jupiter Workshop. pp 63–68
  * Newman et al. (1999) Newman W. I., Symbalisty E. M., Ahrens T. J., Jones E. M., 1999, Icarus, 138, 224
  * Nicholson et al. (1995) Nicholson P. D., et al., 1995, Geophysical Research Letters, 22, 1613
  * Oliphant (2015) Oliphant T. E., 2015, Guide to NumPy. Continuum Press, https://books.google.ca/books?id=g58ljgEACAAJ
  * Pope et al. (1997) Pope K. O., Baines K. H., Ocampo A. C., Ivanov B. A., 1997, Journal of Geophysical Research E: Planets, 102, 21645
  * Potter et al. (2012) Potter R. W., Collins G. S., Kiefer W. S., McGovern P. J., Kring D. A., 2012, Icarus, 220, 730
  * Sakurai (1960) Sakurai A., 1960, Communications on Pure and Applied Mathematics, 13, 353
  * Schlichting (2014) Schlichting H. E., 2014, The Astrophysical Journal Letters, Volume 795, Issue 1, article id. L15, 5 pp. (2014)., 795
  * Schlichting et al. (2015) Schlichting H., Sari R., Yalinewich A., 2015, Icarus, 247
  * Sedov (1946) Sedov L. I., 1946, Prikl. Mat. Mekh, 10, 241
  * Shen et al. (2018) Shen K. J., et al., 2018, The Astrophysical Journal, 865, 15
  * Shuvalov et al. (2013) Shuvalov V. V., Svettsov V. V., Trubetskaya I. A., 2013, Solar System Research, 47, 260
  * Taylor (1950) Taylor G. I., 1950, Proc. R. Soc. Lond. A, 201, 159
  * Waxman & Shvarts (1992) Waxman E., Shvarts D., 1992, Physics of Fluids A, 5, 1035
  * Yalinewich & Caplan (2021) Yalinewich A., Caplan M. E., 2021, Monthly Notices of the Royal Astronomical Society: Letters, 505, L115
  * Yalinewich & Matzner (2019) Yalinewich A., Matzner C. D., 2019, Monthly Notices of the Royal Astronomical Society, 490, 312
  * Yalinewich & Remorov (2020) Yalinewich A., Remorov A., 2020, Atmosphere, 11, 445
  * Yalinewich & Schlichting (2019) Yalinewich A., Schlichting H., 2019, MNRAS, 486, 2780
  * Zahnle & Mac Low (1995) Zahnle K., Mac Low M.-M., 1995, Journal of Geophysical Research, 100, 16885
  * Zel’dovich & Raizer (1967) Zel’dovich Y. B., Raizer Y. P., 1967, Physics of shock waves and high-temperature hydrodynamic phenomena. Dover

##  Appendix A Asymptotic Behaviour

###  A.1 Asymptotic Behaviour for an Infinite Adiabatic Index

In the limit $\gamma\rightarrow\infty$, the equations used in Section 2.2
break down. This is because the integration interval
$V\in\left[\frac{2\delta+\omega}{\delta(2-\gamma)+\omega},\frac{2}{\gamma+1}\right]$
shrinks to zero width in this limit. To overcome this difficulty, we define a
new variable $W=\gamma V$, as was done in (Yalinewich & Remorov, 2020). In
this new variable, the integration domain becomes
$W\in\left[-2-\frac{\omega}{\delta},2\right]$, and so remains finite in the
limit $\gamma\rightarrow\infty$. The boundary conditions for $C$ simplifies as
well, and become: $C=1$ at the sonic point and $C=\sqrt{2}$ at the shock
front. The expression for the derivative involving the new variable is given
by

$\frac{dC}{dW}=\frac{C(W\delta+C^{2}\omega+2C^{2}+2\delta-2)}{2(WC^{2}+W\delta-W+2C^{2}\delta+C^{2}\omega)}\,$
(57)

At the sonic point, the slope is given by

$\left.\frac{dC}{dW}\right|_{s}=-\frac{\delta(\delta-\omega+\sqrt{9\delta^{2}+2\delta\omega-12\delta+\omega^{2}+4}-2)}{8\delta^{2}+4\delta\omega-8\delta-4\omega}\,$
(58)

We can determine $\delta(\omega)$ and $\beta(\omega)$ using the shooting
integration method described in the previous section, since we know the
boundary conditions and $\frac{dC}{dW}$. A set of sample trajectories, for
$\omega=0.5$, is shown in figure 4. Applying the shooting integration method
with varying $\omega$ allows one to obtain the dependencies of $\delta$ and
$\beta$ on $\omega$, shown in figure 5. The values of $\beta$ and $\delta$ in
this asymptotic case fall below their respective energy conservation limits in
equations 29, 30, unlike in the opposite limit $\gamma=1$, where the values of
these parameters are equal to the momentum conservation limits.

Figure 4: Hydrodynamic trajectories for the case
$\gamma\rightarrow\infty,\omega=0.5$. Different coloured lined represent
trajectories with different values of $\delta$. The sonic line is represented
by the black line, and the shock front by the black cross.

Figure 5: Using the shooting integration method with a variety of $\omega$,
plots of $\delta(\omega)$ (at the top) and $\beta(\omega)$ (at the bottom) for
$\gamma\rightarrow\infty$ can be obtained. Although $\delta$ appears to be
linearly dependent on $\omega$, this is actually not the case, since
$\frac{d\delta}{d\omega}$ is not a constant function of $\omega$. It is noted
that $\beta$ has a much lower dependence on $\omega$ than $\delta$.

###  A.2 Asymptotic behaviour for an Infinite Density Slope

The limit $\omega\rightarrow\infty$ is equivalent to an exponential density
profile. In this scenario, different equations for the density profile,
dimensionless hydrodynamic quantities and shock front trajecory have to be
used. A similar problem has been considered (Schlichting et al., 2015), so
many of the expressions are similar. We note however that (Schlichting et al.,
2015) considered the portion of the shock wave that travels in the opposite
direction, towards lower densities, and so their power law index is different.

The density profile is defined by
$\rho(x)=\rho_{0}\exp\left(\frac{x}{h}\right)$ where $\rho_{0}$ is the density
at a reference point (where the explosion started), x is the distance from the
reference point, and h is the atmospheric scale height. In our situation, the
shock wave gets slower as time progresses, while in (Schlichting et al., 2015)
it gets faster and faster and traverses an infinite distance after a finite
time. The dimensionless position is $\chi=\frac{X(t)-x}{h}$ and the
dimensionless hydrodynamic variables are

$v(x,t)=\frac{dX(t)}{dt}V\quad\quad
c(x,t)=\frac{dX(t)}{dt}C\quad\quad\rho(x,t)=\rho_{0}e^{X(t)/h}D$ (59)

In the new framework, we start again with the same slab symmetric hydrodynamic
equations 2, 3 and 4 and ideal gas equation of state, and then substitute in
the new dimensionless parameters defined in 59.

Like in section 2.2, after these substitutions, the equations contain terms
involving different time derivatives of $X$. The time derivatives can be
factored out using equation, by eliminating the second derivative using

$\frac{d^{2}X}{dt^{2}}=-\frac{\delta}{h}\left(\frac{dX}{dt}\right)^{2}\,.$
(60)

This equation looks different from equation 10 because the trajectory of the
shock front is not a power law of time, but instead has a logarithmic
dependence on time given by $X(t)=\frac{h}{\delta}\ln(1+t/t_{0}$) where h is
the atmosphere scale height and $t_{0}$ is a scale time (equation 11 in
(Schlichting et al., 2015)).

In this limit, the hydrodynamic equations equations 2, 3 and 4 can be
represented as a matrix equation $\overleftrightarrow{M}d\vec{A}/\chi=\vec{B}$
where:

$M=\begin{bmatrix}1-V&-D&0\\\ -C^{2}&\gamma(1-V)D&-2CD\\\ -\gamma C-(-\gamma
C+C)V+C&0&-2DV+2D\end{bmatrix}$ (61)

$\vec{A}=\begin{bmatrix}D\\\ V\\\ C\end{bmatrix}$ (62)

$\vec{B}=\begin{bmatrix}-D\\\ -\gamma\delta DV\\\ \gamma CD-CD-2\delta
CD\end{bmatrix}$ (63)

Once again, these equations are linear in the derivatives
$\frac{dV}{d\chi},\frac{dC}{d\chi},\frac{dD}{d\chi}$, so we can isolate them,
and obtain a single differential equation of $\frac{dC}{dV}$ by going through
the same process as was done for equations 11 to 14:

$\frac{dC}{dV}=\frac{C(\frac{C^{2}\gamma}{\delta}-\frac{C^{2}}{\delta}-2C^{2}-V^{2}\gamma^{2}+3V^{2}\gamma+V\gamma^{2}-5V\gamma+2\gamma)}{-\frac{2C^{2}V}{\delta}-4C^{2}V+\frac{2C^{2}}{\delta}+4C^{2}+2V^{3}\gamma-4V^{2}\gamma+2V\gamma}$
(64)

Like in section 2.3, in order integrate equation 64, boundary conditions and
the value of $\alpha$ are needed. Once again, the boundary conditions at the
shock front are:

$V_{f}=\frac{2}{\gamma+1}$ (65)

and

$C_{f}=\frac{\sqrt{2\gamma\left(\gamma-1\right)}}{\gamma+1}\,.$ (66)

The denominator of equation 64 vanishes on the sonic line, given by equation
25. At the sonic point, the numerator also vanishes. This allows us to obtain
the sonic point boundary conditions for this scenario:

$V_{s}=\frac{\frac{1}{\delta}+2}{\frac{1}{\delta}+2-\gamma}$ (67)
$C_{s}=-\frac{\gamma}{\frac{1}{\delta}+2-\gamma}$ (68)

Knowing the boundary conditions and $\frac{dC}{dV}$, we can integrate from the
sonic point to the shock front. The integration process, with the exception of
the exact formulas for $\frac{dC}{dV}$ and boundary conditions, is equivalent
to section 2.3 above, so the details will not be repeated here. However, one
thing to note is that the $C_{s}<0$ never arises, because
$\frac{1}{\delta}+2<\gamma$ for all physically possible ranges of
$\delta,\gamma$. This means that only the first integration case, where the
sonic point and shock front are on the same side of the impact site, occurs.
An example of a hydrodynamic trajectory for this scenario is shown in figure
7. The bisection solving method can be used to obtain $\delta$ for each
possible value of $\gamma$, and then taking the reciprocal gets $\delta$. This
relation $\delta(\gamma)$ is plotted in figure 8. In the limit of $\gamma=1$,
we get $\delta=-1$, which corresponds to momentum conservation. In the limit
$\gamma\rightarrow\infty$, which corresponds to energy conservation, the
solution converges to approximately $\delta=-0.597$. Obtaining this value
requires the same substitution $W=\gamma V$ present in section A1.

###  A.3 Asymptotic behaviour for Infinite Adiabatic Index and Density Slope

Once again, when $\gamma\rightarrow\infty$, the integration along
$V\in\left[\frac{\frac{1}{\delta}+2}{\frac{1}{\delta}+2-\gamma},\frac{2}{\gamma+1}\right]$
shrinks to zero width, and using the new variable $W=\gamma V$ instead of $V$
remedies this. Now the integration range for $W$ is
$\left[-2-\frac{1}{\delta},2\right]$, and this is not empty, so the
integration process can occur. The boundary conditions in $C$ are $C=1$ at the
sonic point and $C=\sqrt{2}$ at the shock front. The expression for the
derivative in this situation is:

$\frac{dC}{dW}=\frac{C(W+\frac{C^{2}}{\delta}+2)}{2(W+\frac{C^{2}}{\delta}+2C^{2})}\,$
(69)

And at the sonic point this slope is:

$\left.\frac{dC}{dW}\right|_{s}=-\frac{\frac{3}{\delta}+\sqrt{\frac{9}{\delta^{2}}+\frac{8}{\delta}+16}}{4(\frac{1}{\delta}+2)}\,$
(70)

We can now determine $\delta(\gamma)$ using the shooting integration method,
since the boundary conditions and $\frac{dC}{dW}$ are known. Some hydrodynamic
trajectories for different values of $\delta$ close to the solution are shown
in figure 6. Using a bisection solver, we find that
$\lim_{\gamma\rightarrow\infty}\delta\approx-0.597$.

Figure 6: Hydrodynamic trajectories for the case
$\omega\rightarrow\infty,\gamma\rightarrow\infty$. Different coloured lined
represent trajectories with different values of $\frac{1}{\delta}$. The sonic
line is represented by the black line, and the shock front by the black cross.

Figure 7: Hydrodynamic trajectories for the exponential density profile and
$\gamma=3.5$. Fortunately, in this scenario, the integration is always
straightforward and there are no mathematical singularities. Once a
$\frac{1}{\delta}$ value is found whose integration curve ends close enough to
the shock front boundary condition, its reciprocal is taken to obtain the
corresponding $\delta$ value. Figure 8:  A plot of $\delta(\gamma)$ for the
exponential density profile scenario ($\omega\rightarrow\infty$) can be
obtained by applying the shooting integration method with a variety of
$\gamma$ and taking the reciprocal of the respective $\frac{1}{\delta}$
values.