parameters of $f_{2}$ are updated based on an exponential moving average of
$f_{1}$ parameters. This design is independent of negative samples and avoids
the computational overheads of employing graph and node features across the
training batch.

Figure 11: BGRL design with online, target encoders $f_{1},f_{2}$, their node
projectors $z_{n1},z_{n2}$, g-l contrastive mode and BL objective.

###  D.6 Evaluation protocols

Generally, GCL frameworks differ in the choice of augmentors, encoder designs,
contrastive modes and objective functions. Thus, it is unclear how the
performance of an augmentor can be fairly tested. To address this issue, we
fix as many configurations as possible and focus solely on the impact of the
augmentor. Our reasoning is based on the well-established benchmark study by
Zhu et al. (2021b). Adhering to this approach, we leverage the designs of
widely adopted frameworks such as GRACE (Zhu et al., 2020), MVGRL (Hassani &
Khasahmadi, 2020), GraphCL (You et al., 2020) and BGRL (Thakoor et al., 2021),
where topological augmentation is necessary and evaluate on unsupervised node
and graph classification tasks. These frameworks were chosen so as to
extensively experiment on encoder designs, contrastive modes and objectives.
We leverage graph convolution network (GCN) (Kipf & Welling, 2016a) and graph
isomorphism network (GIN) (Xu et al., 2018) based encoders for node and graph-
based tasks respectively. Also, based on the empirical benchmark observations
in Zhu et al. (2021b), we perform feature masking with a constant masking
ratio of $0.3$ in all the experiments and perform an extensive grid-search
over the hyper-parameters to find their ideal values. The perturbation ratios
$\gamma_{1},\gamma_{2}$ are selected from $\\{0.0,0.1,0.2,0.3,0.4,0.5\\}$ to
prevent excessive graph corruption. The number of layers in the encoder are
chosen from $\\{2,4,8\\}$, hidden feature dimensions are chosen from
$\\{128,256,512\\}$, learning rate is chosen from
$\\{10^{-2},10^{-3},10^{-4}\\}$, weight decay is set to $10^{-5}$ and the
maximum epoch count is set to $2000$. An early stopping strategy with a
patience interval of 50 epochs with respect to contrastive loss is also
employed. Finally, we choose the Adam optimizer for learning the embeddings.
For testing, we follow a linear evaluation protocol (Velickovic et al., 2019;
Zhu et al., 2021b) where the embeddings learned by the encoders in
unsupervised settings are classified using a logistic regression classifier
with a train, validation and test split of $10\%,10\%,80\%$ respectively. The
process is repeated with $10$ random splits and the classification accuracies
are reported. The selected hyper-parameters based on a grid search are
reported in Table 12, 13, 14, 15 and can be reproduced with the provided code
and instructions.

Table 12: Hyper-parameters for GRACE-based design for node-classification
Dataset | mode | loss | $\tau$ | $\gamma_{1}$ | $\gamma_{2}$ | L | lr | wd | hidden dim
---|---|---|---|---|---|---|---|---|---
CORA | l-l | InfoNCE | $0.2$ | $0.5$ | $0.4$ | $2$ | $10^{-4}$ | $10^{-5}$ | $256$
AMAZON-PHOTO | l-l | InfoNCE | $0.2$ | $0.3$ | $0.4$ | $2$ | $10^{-2}$ | $10^{-5}$ | $256$
PUBMED | l-l | InfoNCE | $0.2$ | $0.4$ | $0.1$ | $2$ | $10^{-3}$ | $10^{-5}$ | $256$
COAUTHOR-CS | l-l | InfoNCE | $0.4$ | $0.5$ | $0.4$ | $2$ | $10^{-4}$ | $10^{-5}$ | $256$
COAUTHOR-PHY | l-l | InfoNCE | $0.4$ | $0.5$ | $0.5$ | $2$ | $10^{-2}$ | $10^{-5}$ | $128$
Table 13: Hyper-parameters for MVGRL-based design for node-classification
Dataset | mode | loss | $\gamma_{1}$ | $\gamma_{2}$ | L | lr | wd | hidden dim
---|---|---|---|---|---|---|---|---
CORA | g-l | JSD | $0.0$ | $0.4$ | $2$ | $10^{-3}$ | $10^{-5}$ | $512$
AMAZON-PHOTO | g-l | JSD | $0.0$ | $0.3$ | $2$ | $10^{-4}$ | $10^{-5}$ | $512$
PUBMED | g-l | JSD | $0.1$ | $0.3$ | $2$ | $10^{-3}$ | $10^{-5}$ | $512$
COAUTHOR-CS | g-l | JSD | $0.0$ | $0.2$ | $2$ | $10^{-4}$ | $10^{-5}$ | $512$
COAUTHOR-PHY | g-l | JSD | $0.1$ | $0.4$ | $2$ | $10^{-3}$ | $10^{-5}$ | $512$
Table 14: Hyper-parameters for GraphCL-based design for graph-classification
Dataset | mode | loss | $\tau$ | $\gamma_{1}$ | $\gamma_{2}$ | L | lr | wd | hidden dim
---|---|---|---|---|---|---|---|---|---
PROTEINS | g-g | InfoNCE | $0.5$ | $0.2$ | $0.2$ | $2$ | $10^{-2}$ | $10^{-5}$ | $128$
IMDB-BINARY | g-g | InfoNCE | $0.5$ | $0.3$ | $0.5$ | $2$ | $10^{-2}$ | $10^{-5}$ | $128$
MUTAG | g-g | InfoNCE | $0.5$ | $0.2$ | $0.2$ | $2$ | $10^{-3}$ | $10^{-5}$ | $128$
IMDB-MULTI | g-g | InfoNCE | $0.5$ | $0.2$ | $0.2$ | $2$ | $10^{-2}$ | $10^{-5}$ | $128$
NCI1 | g-g | InfoNCE | $0.2$ | $0.2$ | $0.3$ | $2$ | $10^{-3}$ | $10^{-5}$ | $128$
Table 15: Hyper-parameters for BGRL-based design for graph-classification
Dataset | mode | loss | $\gamma_{1}$ | $\gamma_{2}$ | L | lr | wd | hidden dim
---|---|---|---|---|---|---|---|---
PROTEINS | g-l | BL | $0.1$ | $0.1$ | $2$ | $10^{-2}$ | $10^{-5}$ | $256$
IMDB-BINARY | g-l | BL | $0.5$ | $0.3$ | $2$ | $10^{-3}$ | $10^{-5}$ | $128$
MUTAG | g-l | BL | $0.2$ | $0.2$ | $2$ | $10^{-3}$ | $10^{-5}$ | $128$
IMDB-MULTI | g-l | BL | $0.2$ | $0.2$ | $2$ | $10^{-3}$ | $10^{-5}$ | $128$
NCI1 | g-l | BL | $0.3$ | $0.1$ | $2$ | $10^{-3}$ | $10^{-5}$ | $256$

###  D.7 Baselines

We report the baseline node and graph classification performance from
previously published reports (Hassani & Khasahmadi, 2020; You et al., 2020;
Zhu et al., 2021c; Xu et al., 2021). For node classification, we leverage raw
features (Velickovic et al., 2019), DeepWalk (Perozzi et al., 2014), Graph
Auto Encoder (GAE) (Kipf & Welling, 2016b), Deep Graph Infomax (DGI)
(Velickovic et al., 2019), GCN networks and Graph Attention Networks (GAT)
(Veličković et al., 2017) (see Table 16). For graph classification, we
leverage Graphlet Kernel (GK) (Shervashidze et al., 2009), Weisfeiler-Lehman
Kernel (WL) (Shervashidze et al., 2011), Deep Graph Kernel (DGK) (Yanardag &
Vishwanathan, 2015), Multi-scale Laplacian Graph Kernel (MLG) (Kondor & Pan,
2016), node2vec (Grover & Leskovec, 2016), sub2vec (Adhikari et al., 2018),
graph2vec (Narayanan et al., 2017), InfoGraph (Sun et al., 2019), GIN, GCN and
GAT networks (see Table 17).

Table 16: Baseline node classification accuracies from published reports.
Method | CORA | Amazon-Photo | PubMed | Coauthor-CS | Coauthor-Phy
---|---|---|---|---|---
Raw Feat | $47.90\pm 0.40$ | $78.53\pm 0.00$ | $69.10\pm 0.30$ | $90.37\pm 0.00$ | $93.58\pm 0.00$
DeepWalk | $67.20\pm 0.00$ | $89.44\pm 0.11$ | $65.30\pm 0.00$ | $84.61\pm 0.22$ | $91.77\pm 0.15$
DeepWalk+feat | $70.70\pm 0.60$ | $90.05\pm 0.08$ | $74.30\pm 0.90$ | $87.70\pm 0.04$ | $94.90\pm 0.09$
GAE | $71.50\pm 0.40$ | $91.62\pm 0.13$ | $72.10\pm 0.50$ | $90.01\pm 0.71$ | $94.92\pm 0.07$
DGI | $82.30\pm 0.60$ | $91.61\pm 0.22$ | $76.80\pm 0.60$ | $92.15\pm 0.63$ | $94.51\pm 0.52$
GCN | $81.50\pm 0.00$ | $92.42\pm 0.22$ | $79.0\pm 0.00$ | $93.03\pm 0.31$ | $95.65\pm 0.16$
GAT | $83.0\pm 0.70$ | $92.56\pm 0.35$ | $79.0\pm 0.30$ | $92.31\pm 0.24$ | $95.47\pm 0.15$
Table 17: Baseline graph classification accuracies from published reports.
Method | MUTAG | PROTEINS | IMDB-BINARY | IMDB-MULTI | NCI1
---|---|---|---|---|---
GK | $81.70\pm 2.10$ | $-$ | $65.87\pm 0.98$ | $43.90\pm 0.40$ | $-$
WL | $80.72\pm 3.00$ | $72.92\pm 0.56$ | $72.30\pm 3.44$ | $47.0\pm 0.50$ | $80.01\pm 0.50$
DGK | $87.44\pm 2.72$ | $73.30\pm 0.82$ | $66.96\pm 0.56$ | $44.60\pm 0.50$ | $80.31\pm 0.46$
MLG | $87.9\pm 1.60$ | $-$ | $66.6\pm 0.30$ | $41.2\pm 0.00$ | $80.8\pm 1.30$
node2vec | $72.63\pm 10.20$ | $57.49\pm 3.57$ | $-$ | $-$ | $54.89\pm 1.61$
sub2vec | $61.05\pm 15.80$ | $53.03\pm 5.55$ | $55.26\pm 1.54$ | $36.70\pm 0.80$ | $52.84\pm 1.47$
graph2vec | $83.15\pm 9.25$ | $73.30\pm 2.05$ | $71.10\pm 0.54$ | $50.40\pm 0.90$ | $73.22\pm 1.81$
InfoGraph | $89.01\pm 1.13$ | $74.44\pm 0.31$ | $73.03\pm 0.87$ | $49.70\pm 0.50$ | $76.20\pm 1.06$
GIN | $89.4\pm 5.60$ | $76.2\pm 2.80$ | $75.1\pm 5.10$ | $52.3\pm 2.80$ | $82.7\pm 1.70$
GCN | $85.6\pm 5.80$ | $75.2\pm 3.60$ | $74.0\pm 3.4$ | $51.9\pm 3.8$ | $80.2\pm 2.0$
GAT | $89.4\pm 6.10$ | $74.7\pm 4.00$ | $70.5\pm 2.30$ | $47.8\pm 3.10$ | $66.6\pm 2.20$

##  Appendix E Augmentor Overheads

We establish an augmentor overhead benchmark in terms of memory usage and
latency for an extensive range of techniques on the benchmark datasets.
Specifically, we measure:

MEMORY: The memory (in MB) that is required solely by the augmentation phase.

LATENCY(CPU): The wall-clock time (in seconds) for performing the
augmentation.

LATENCY(w/ GPU): The wall-clock time (in seconds) for performing the
augmentation when GPU is available.

We run $10$ experiments for each augmentor-dataset combination with a fixed
perturbation ratio of $0.5$ (and use a batch size of $128$ for graph datasets)
to report the mean and standard deviation of the metrics in Table 18, 19. From
Table 18, we can observe that simple NodeDropping and EdgeDropping schemes are
computationally the most efficient, with EdgeDroppingEVC and diffusion-based
approaches being the least. The code for EdgeDroppingEVC has been adopted from
the original paper by Zhu et al. (2021c) and leverages the NetworkX library to
compute centrality scores. Furthermore, since this approach lacks GPU support,
the hardware acceleration is minimal. In fact, the overheads of data transfers
from GPU to CPU and the construction of NetworkX graphs significantly increase
the memory consumption and latencies for this technique. Similar patterns can
be observed for other augmentors as well where the GPU-CPU overheads dominate
the parallelization benefits. On the contrary, diffusion-based approaches
leverage GPU acceleration for heavy matrix operations and gain significant
speedups to achieve reasonable latencies. A similar pattern can be observed in
Table 19 for graph datasets with a batch size of $128$. Note that some of the
results show $0$ variance due to numerical rounding.

The overheads of $rLap$ in Table 18, 19 lie within a very narrow margin of the
most efficient techniques such as EdgeDropping, NodeDropping and
RandomWalkSubgraph. These low latencies are achieved by designing the
augmentor in C++ and minimizing the interaction with Python APIs.
Additionally, we compile Eigen without MKL optimizations and also do not
support a GPU kernel for $rLap$ yet, which explains the lack of performance
improvement when hardware acceleration is available. However, we tend to
achieve minimal resource consumption by just leveraging the C++ APIs to the
fullest. We leave the development of GPU support for $rLap$ as future work.

Table 18: Statistics of resource consumption by augmentors on node
classification datasets.
augmentor | dataset | memory | latency(cpu) | latency(w/ gpu)
---|---|---|---|---
EdgeAddition | CORA | $3.23\pm 0.09$ | $0.0024\pm 0.0$ | $0.0033\pm 0.0001$
EdgeDropping | CORA | $0.51\pm 0.0539$ | $0.0004\pm 0.0$ | $0.0019\pm 0.0001$
EdgeDroppingDegree | CORA | $3.85\pm 1.1307$ | $0.003\pm 0.001$ | $0.0046\pm 0.0032$
EdgeDroppingEVC | CORA | $39.58\pm 0.14$ | $0.2214\pm 0.0038$ | $0.3542\pm 0.0021$
EdgeDroppingPR | CORA | $3.18\pm 0.0872$ | $0.0023\pm 0.0001$ | $0.0037\pm 0.0001$
MarkovDiffusion | CORA | $228.73\pm 38.1103$ | $0.7966\pm 0.1115$ | $0.2539\pm 0.0519$
NodeDropping | CORA | $1.31\pm 1.3729$ | $0.0006\pm 0.0$ | $0.0026\pm 0.0008$
PPRDiffusion | CORA | $140.0\pm 0.1342$ | $0.2531\pm 0.0017$ | $0.5933\pm 0.0589$
RandomWalkSubgraph | CORA | $3.26\pm 0.102$ | $0.001\pm 0.0$ | $0.0027\pm 0.0001$
rLap | CORA | $0.38\pm 0.04$ | $0.0037\pm 0.0005$ | $0.0048\pm 0.0002$
EdgeAddition | AMAZON-PHOTO | $20.94\pm 5.4471$ | $0.0368\pm 0.0017$ | $0.0071\pm 0.0003$
EdgeDropping | AMAZON-PHOTO | $5.64\pm 0.5953$ | $0.0032\pm 0.0015$ | $0.0022\pm 0.0001$
EdgeDroppingDegree | AMAZON-PHOTO | $26.13\pm 4.2223$ | $0.0391\pm 0.0007$ | $0.0041\pm 0.0001$
EdgeDroppingEVC | AMAZON-PHOTO | $2264.1\pm 0.004$ | $1.554\pm 0.012$ | $1.4576\pm 0.024$
EdgeDroppingPR | AMAZON-PHOTO | $11.32\pm 1.4999$ | $0.0094\pm 0.0001$ | $0.0043\pm 0.0003$
MarkovDiffusion | AMAZON-PHOTO | $158.42\pm 3.0019$ | $30.6898\pm 0.127$ | $1.7687\pm 0.0044$
NodeDropping | AMAZON-PHOTO | $1.53\pm 0.19$ | $0.0017\pm 0.0012$ | $0.0028\pm 0.0001$
PPRDiffusion | AMAZON-PHOTO | $389.33\pm 2.1808$ | $1.6942\pm 0.0247$ | $0.8229\pm 0.0112$
RandomWalkSubgraph | AMAZON-PHOTO | $7.26\pm 0.1356$ | $0.0024\pm 0.0001$ | $0.0029\pm 0.0007$
rLap | AMAZON-PHOTO | $32.6\pm 0.5119$ | $0.0848\pm 0.0023$ | $0.0891\pm 0.0014$
EdgeAddition | PUBMED | $10.93\pm 0.7785$ | $0.0126\pm 0.0002$ | $0.0045\pm 0.0001$
EdgeDropping | PUBMED | $1.89\pm 0.1136$ | $0.002\pm 0.0003$ | $0.0022\pm 0.0001$
EdgeDroppingDegree | PUBMED | $15.45\pm 2.1238$ | $0.015\pm 0.001$ | $0.0032\pm 0.0003$
EdgeDroppingEVC | PUBMED | $2236.23\pm 0.7537$ | $1.2787\pm 0.0152$ | $1.3101\pm 0.0156$
EdgeDroppingPR | PUBMED | $4.79\pm 0.2071$ | $0.0059\pm 0.0005$ | $0.0043\pm 0.0003$
MarkovDiffusion | PUBMED | $152.87\pm 2.0995$ | $44.4244\pm 0.1306$ | $2.5587\pm 0.0054$
NodeDropping | PUBMED | $1.78\pm 0.172$ | $0.0021\pm 0.0002$ | $0.0036\pm 0.0002$
PPRDiffusion | PUBMED | $438.42\pm 1.8595$ | $19.0975\pm 0.3605$ | $2.5488\pm 0.0135$
RandomWalkSubgraph | PUBMED | $1.74\pm 0.1625$ | $0.0044\pm 0.0003$ | $0.003\pm 0.0001$
rLap | PUBMED | $10.59\pm 3.1908$ | $0.0269\pm 0.0008$ | $0.0315\pm 0.0039$
EdgeAddition | COAUTHOR-CS | $16.71\pm 3.0409$ | $0.0286\pm 0.0015$ | $0.0057\pm 0.0003$
EdgeDropping | COAUTHOR-CS | $3.02\pm 0.5758$ | $0.007\pm 0.0018$ | $0.0025\pm 0.0003$
EdgeDroppingDegree | COAUTHOR-CS | $15.1\pm 2.3078$ | $0.0363\pm 0.0014$ | $0.0037\pm 0.0001$
EdgeDroppingEVC | COAUTHOR-CS | $1044.6\pm 0.5745$ | $10.2603\pm 0.0266$ | $10.6018\pm 0.1096$
EdgeDroppingPR | COAUTHOR-CS | $5.09\pm 0.359$ | $0.0114\pm 0.0003$ | $0.0052\pm 0.001$
MarkovDiffusion | COAUTHOR-CS | $177.9667\pm 7.0002$ | $66.5972\pm 0.1857$ | $3.3583\pm 0.1758$
NodeDropping | COAUTHOR-CS | $2.33\pm 0.1269$ | $0.0035\pm 0.0004$ | $0.0033\pm 0.0001$
PPRDiffusion | COAUTHOR-CS | $416.58\pm 3.9992$ | $46.9031\pm 1.4907$ | $2.2614\pm 0.1293$
RandomWalkSubgraph | COAUTHOR-CS | $5.63\pm 0.1792$ | $0.0066\pm 0.0002$ | $0.003\pm 0.0001$
rLap | COAUTHOR-CS | $13.6\pm 0.9798$ | $0.0508\pm 0.0007$ | $0.0602\pm 0.002$
EdgeAddition | COAUTHOR-PHY | $44.59\pm 8.3973$ | $0.0975\pm 0.0048$ | $0.0121\pm 0.0019$
EdgeDropping | COAUTHOR-PHY | $11.57\pm 1.5685$ | $0.0098\pm 0.0009$ | $0.002\pm 0.0001$
EdgeDroppingDegree | COAUTHOR-PHY | $20.7\pm 0.1483$ | $0.1219\pm 0.0025$ | $0.0046\pm 0.0003$
EdgeDroppingEVC | COAUTHOR-PHY | $2474.86\pm 0.4294$ | $28.2138\pm 0.0304$ | $28.9077\pm 0.2304$
EdgeDroppingPR | COAUTHOR-PHY | $16.96\pm 1.8205$ | $0.0311\pm 0.0002$ | $0.0041\pm 0.002$
MarkovDiffusion | COAUTHOR-PHY | $316.0\pm 0.4$ | $316.1449\pm 0.1584$ | $15.1789\pm 0.2203$
NodeDropping | COAUTHOR-PHY | $4.08\pm 1.7831$ | $0.0067\pm 0.0004$ | $0.0046\pm 0.0001$
PPRDiffusion | COAUTHOR-PHY | $695.0\pm 2.397$ | $124.368\pm 2.9031$ | $10.394\pm 0.2053$
RandomWalkSubgraph | COAUTHOR-PHY | $7.85\pm 0.9563$ | $0.014\pm 0.0003$ | $0.0057\pm 0.0047$
rLap | COAUTHOR-PHY | $43.9\pm 11.3004$ | $0.191\pm 0.0258$ | $0.2119\pm 0.0064$
Table 19: Statistics of resource consumption by augmentors on graph
classification datasets.
augmentor | dataset | memory | latency(cpu) | latency(w/ gpu)
---|---|---|---|---
EdgeAddition | IMDB-BINARY | $4.78\pm 0.3628$ | $0.0307\pm 0.0003$ | $0.0191\pm 0.0121$
EdgeDropping | IMDB-BINARY | $2.15\pm 0.75$ | $0.0079\pm 0.0038$ | $0.0038\pm 0.0004$
EdgeDroppingDegree | IMDB-BINARY | $7.51\pm 0.6978$ | $0.0346\pm 0.0002$ | $0.0121\pm 0.001$
EdgeDroppingEVC | IMDB-BINARY | $2101.35\pm 2.3127$ | $0.7977\pm 0.0068$ | $0.9294\pm 0.0138$
EdgeDroppingPR | IMDB-BINARY | $2.27\pm 0.0458$ | $0.0184\pm 0.0002$ | $0.0166\pm 0.0001$
MarkovDiffusion | IMDB-BINARY | $195.43\pm 60.3246$ | $8.2739\pm 0.0466$ | $0.4589\pm 0.006$
NodeDropping | IMDB-BINARY | $2.03\pm 0.064$ | $0.0046\pm 0.0001$ | $0.0032\pm 0.0001$
PPRDiffusion | IMDB-BINARY | $170.13\pm 22.767$ | $0.658\pm 0.0068$ | $0.5295\pm 0.0105$
RandomWalkSubgraph | IMDB-BINARY | $2.09\pm 0.1136$ | $0.0063\pm 0.0001$ | $0.0096\pm 0.0016$
rLap | IMDB-BINARY | $6.96\pm 0.3583$ | $0.0423\pm 0.0027$ | $0.0412\pm 0.0009$
EdgeAddition | IMDB-MULTI | $5.25\pm 1.0249$ | $0.0359\pm 0.0017$ | $0.0193\pm 0.0001$
EdgeDropping | IMDB-MULTI | $2.18\pm 0.1327$ | $0.0064\pm 0.0$ | $0.0046\pm 0.0001$
EdgeDroppingDegree | IMDB-MULTI | $8.72\pm 1.1232$ | $0.0365\pm 0.0001$ | $0.0157\pm 0.0001$
EdgeDroppingEVC | IMDB-MULTI | $2119.91\pm 2.3364$ | $0.806\pm 0.0104$ | $0.8579\pm 0.0062$
EdgeDroppingPR | IMDB-MULTI | $4.45\pm 0.7159$ | $0.0221\pm 0.0003$ | $0.0236\pm 0.0004$
MarkovDiffusion | IMDB-MULTI | $28.02\pm 97.7836$ | $5.5597\pm 0.0245$ | $0.386\pm 0.0046$
NodeDropping | IMDB-MULTI | $2.17\pm 0.1005$ | $0.0054\pm 0.0001$ | $0.0045\pm 0.0012$
PPRDiffusion | IMDB-MULTI | $98.67\pm 41.6684$ | $0.3907\pm 0.0183$ | $0.7974\pm 0.0134$
RandomWalkSubgraph | IMDB-MULTI | $4.07\pm 0.9144$ | $0.0075\pm 0.0001$ | $0.0107\pm 0.0006$
rLap | IMDB-MULTI | $7.93\pm 0.6165$ | $0.0421\pm 0.0023$ | $0.0415\pm 0.0009$
EdgeAddition | MUTAG | $3.22\pm 0.14$ | $0.0018\pm 0.0001$ | $0.0046\pm 0.0001$
EdgeDropping | MUTAG | $1.61\pm 0.0943$ | $0.0004\pm 0.0$ | $0.0025\pm 0.0001$
EdgeDroppingDegree | MUTAG | $4.01\pm 0.7713$ | $0.0021\pm 0.0001$ | $0.0047\pm 0.0004$
EdgeDroppingEVC | MUTAG | $2118.92\pm 5.3619$ | $0.0445\pm 0.0005$ | $0.0683\pm 0.0061$
EdgeDroppingPR | MUTAG | $3.35\pm 0.1118$ | $0.0025\pm 0.0001$ | $0.006\pm 0.0004$
MarkovDiffusion | MUTAG | $164.0\pm 24.8769$ | $0.3848\pm 0.0067$ | $0.1774\pm 0.0075$
NodeDropping | MUTAG | $1.61\pm 0.0943$ | $0.0008\pm 0.0$ | $0.0033\pm 0.0001$
PPRDiffusion | MUTAG | $117.04\pm 0.2973$ | $0.0829\pm 0.0019$ | $0.602\pm 0.0101$
RandomWalkSubgraph | MUTAG | $2.7\pm 0.3821$ | $0.0009\pm 0.0$ | $0.0044\pm 0.0008$
rLap | MUTAG | $1.62\pm 0.3682$ | $0.0021\pm 0.0$ | $0.0041\pm 0.0001$
EdgeAddition | PROTEINS | $6.68\pm 0.3124$ | $0.0288\pm 0.0014$ | $0.0153\pm 0.0013$
EdgeDropping | PROTEINS | $2.58\pm 0.06$ | $0.0058\pm 0.0012$ | $0.004\pm 0.0001$
EdgeDroppingDegree | PROTEINS | $10.48\pm 1.6024$ | $0.0319\pm 0.0003$ | $0.0127\pm 0.0002$
EdgeDroppingEVC | PROTEINS | $2095.08\pm 0.6554$ | $1.1132\pm 0.0209$ | $1.3092\pm 0.2126$
EdgeDroppingPR | PROTEINS | $5.27\pm 0.064$ | $0.0195\pm 0.0018$ | $0.0182\pm 0.0001$
MarkovDiffusion | PROTEINS | $203.03\pm 84.7664$ | $25.7943\pm 0.1671$ | $1.1789\pm 0.0087$
NodeDropping | PROTEINS | $2.0\pm 0.6678$ | $0.0079\pm 0.0019$ | $0.0062\pm 0.0001$
PPRDiffusion | PROTEINS | $295.79\pm 26.6581$ | $4.8792\pm 0.0261$ | $1.6547\pm 0.549$
RandomWalkSubgraph | PROTEINS | $2.23\pm 0.064$ | $0.0074\pm 0.0016$ | $0.0099\pm 0.0001$
rLap | PROTEINS | $6.45\pm 0.3442$ | $0.0405\pm 0.0024$ | $0.0407\pm 0.0009$
EdgeAddition | NCI1 | $4.26\pm 0.5122$ | $0.0517\pm 0.0006$ | $0.046\pm 0.0015$
EdgeDropping | NCI1 | $1.72\pm 0.6447$ | $0.0108\pm 0.0001$ | $0.0085\pm 0.0002$
EdgeDroppingDegree | NCI1 | $7.04\pm 0.6296$ | $0.0574\pm 0.0013$ | $0.0389\pm 0.0009$
EdgeDroppingEVC | NCI1 | $2071.43\pm 5.9791$ | $1.9465\pm 0.1448$ | $2.3067\pm 0.0682$
EdgeDroppingPR | NCI1 | $4.11\pm 0.1136$ | $0.0548\pm 0.0064$ | $0.0603\pm 0.0017$
MarkovDiffusion | NCI1 | $141.825\pm 6.632$ | $37.7575\pm 0.3818$ | $2.0278\pm 0.0061$
NodeDropping | NCI1 | $2.22\pm 0.1249$ | $0.019\pm 0.0027$ | $0.0219\pm 0.0008$
PPRDiffusion | NCI1 | $157.16\pm 8.8786$ | $8.8779\pm 0.0595$ | $2.4582\pm 0.011$
RandomWalkSubgraph | NCI1 | $2.68\pm 0.5862$ | $0.0181\pm 0.0008$ | $0.0229\pm 0.0014$
rLap | NCI1 | $2.82\pm 0.5474$ | $0.074\pm 0.0003$ | $0.0801\pm 0.0011$

##  Appendix F Additional Experiments

###  F.1 PPR Diffusion and rLap

In our experiments, we observed that PPRDiffusion on medium scale graphs and
MVGRL design leads to OOM and requires sub-graph sampling to proceed with the
GCL phases. To this end, we observed a decrease in unsupervised node
classification performance when the size of the sub-graph decreased (see Table
20). As the nodes for a sub-graph are randomly selected, one can’t guarantee
that the sub-graph will capture useful structural information for the GNN
encoders to exploit. In a worst case scenario, if all the nodes of a sub-graph
are disconnected, the GNN encoder essentially acts as an MLP. Thus, as the
size of the sub-graph reduces, we can expect a heavy corruption of topological
information in the augmented views.

Table 20: Evaluation (in accuracy) on benchmark node datasets with MVGRL based
design and PPR Diffusion with varying sub-graph sizes. PPRDiffusion indicates
the full graph view and PPRDiffusion-$x$ indicates a sampled sub-graph with
$x$ nodes as the graph view.
Augmentor | CORA | Amazon-Photo | PubMed | Coauthor-CS | Coauthor-Phy
---|---|---|---|---|---
PPRDiffusion | $84.12\pm 2.72$ | $90.65\pm 1.01$ | OOM | OOM | OOM
PPRDiffusion-$8192$ | $-$ | $-$ | $82.7\pm 0.85$ | $90.9\pm 1.06$ | $94.03\pm 0.5$
PPRDiffusion-$4096$ | $-$ | $87.38\pm 0.91$ | $80.79\pm 1.4$ | $89.0\pm 0.83$ | $92.52\pm 0.39$
PPRDiffusion-$2048$ | $80.93\pm 2.33$ | $85.65\pm 1.41$ | $76.4\pm 0.88$ | $88.43\pm 0.77$ | $89.24\pm 0.69$

####  F.1.1 Reducing computational overheads

For graphs such as COAUTHOR-PHY, the computational overheads of PPRDiffusion
are relatively significant when compared to other augmentors (see Table 18).
The standard practice to avoid such bottlenecks is to cache the output of the
diffusion operator ${\mathbf{S}}^{PPR}$ and reuse it for sampling sub-graphs.
To further reduce the overheads of computing ${\mathbf{S}}^{PPR}$ prior to
caching, we leverage Theorem 3.2.

Sketch: The COAUTHOR-PHY graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$
represents $|{\mathcal{V}}|=34,493$ nodes and $|{\mathcal{E}}|=495,924$ edges.
To randomly sample a sub-graph of size $8192$ from its diffused variant
${\mathbf{S}}^{PPR}$, one can compute a randomized Schur complement
${\mathbf{R}}_{\gamma|{\mathcal{V}}|}$ using $rLap$ with $\gamma=0.5$, apply
the PPR diffusion operator on ${\mathbf{R}}_{\gamma|{\mathcal{V}}|}$ to obtain
$\widetilde{{\mathbf{S}}}^{PPR}_{\gamma}$ and cache it.
$\widetilde{{\mathbf{S}}}^{PPR}_{\gamma}$ represents $\approx 17,000$ nodes
from which a sub-graph of size $8192$ can be sampled with relatively less
overheads.

We benchmark the computational overheads of:

  * •

PPRDiffusion on COAUTHOR-PHY graph followed by random sampling of a sub-graph
with $8192$ nodes.

  * •

$rLap$ on COAUTHOR-PHY graph with $\gamma=0.5,o_{v}=rand,o_{n}=asc$ followed
by PPRDiffusion and random sampling of a sub-graph with $8192$ nodes.

The results in Figure 12 demonstrate the effectiveness of the latter approach
in achieving $\approx 0.25\times$ reduction in memory consumption, $\approx
10\times$ reduction in CPU only latencies and $\approx 5\times$ reduction in
latencies when GPU is available. We also compared the unsupervised node
classification accuracy and observed a slight performance improvement of
$\approx 1\%$.

Limitations and Tradeoffs: From Algorithm 1, notice that after every outer
loop iteration, the number of edges in the graph representing the current
Schur complement state will be less than or equal to the number of edges at
the beginning of the iteration. To quantify the impact of diffusion on these
randomized Schur complements, we measured the edge counts of sub-graphs of
size $8192$ randomly sampled from the following graphs 777We consider the
$o_{v}=rand,o_{n}=asc$ variant of $rLap$, which is the default setting. :

  * •

COAUTHOR-CS:PPR: PPRDiffusion applied on COAUTHOR-CS.

  * •

COAUTHOR-CS:rLapPPR: $rLap$ on COAUTHOR-CS with $\gamma=0.5$ followed by
PPRDiffusion.

  * •

COAUTHOR-PHY:PPR: PPRDiffusion applied on COAUTHOR-PHY.

  * •

COAUTHOR-PHY:rLapPPR: $rLap$ on COAUTHOR-PHY with $\gamma=0.5$ followed by
PPRDiffusion.

The results in Figure 13 show that the randomly sampled sub-graph from
COAUTHOR-CS:rLapPPR has $\approx 2.5\times$ the number of edges in COAUTHOR-
CS:PPR. Similarly, the sub-graph from COAUTHOR-PHY:rLapPPR has $\approx
1.5\times$ the number of edges in COAUTHOR-PHY:PPR. Thus, although one can
save memory and time to address the bottlenecks of PPRDiffusion, the GNN
encoders would eventually consume additional memory due to the message passing
operations on relatively more edges.

In such situations, if training latencies take precedence over memory
constraints, then applying $rLap$ and proceeding with PPRDiffusion will lead
to better training throughput. Else, based on the edge count pattern in Figure
14, the practitioner can use a lower $\gamma$ value and consider the trade-
offs between augmentation latencies and available memory as per their setup.

|
---|---
(a) Memory (MB) overhead | (b) Latency (sec) overhead
|
(c) Latency w/ GPU (sec) overhead | (d) Node classification accuracy
Figure 12: Computational overheads and unsupervised node classification
accuracy of PPRDiffusion and $rLap$ \+ PPRDiffusion with $\gamma=0.5$ on
COAUTHOR-PHY dataset.

Figure 13: Plots of edge counts of randomly sampled sub-graphs of size $8192$
with $\gamma=0.5$ for $rLap$ based techniques.
|
---|---
(a) $\gamma=0.1$ | (b) $\gamma=0.2$
|
(c) $\gamma=0.3$ | (d) $\gamma=0.4$
Figure 14: Plots of edge counts of randomly sampled sub-graphs of size $8192$
with varying $\gamma$ for $rLap$ based techniques.


appears in the ‘after’ image between MM4 and MM5, though probably associated
with the latter. Bayandina et al. (2022) specluate that the H2O maser flare is
associated with this new source, noting that the vacuum light-travel time from
MM1 to the new source (in the sky plane) of 78 d is comparable to the 100 d
delay in the maser flares. If this association is correct, we offer two
possibilities for the generation of the new H2O maser cluster. The first is
that the radiation from the accretion event in MM1 has a spectral energy
distribution (SED) that is IR-dominated, and very similar to those in Stecklum
et al. (2021). This radiation impinges on gas containing very dense ($\sim
10^{10}$ H2 cm-3) and cold ($<$50 K) water molecules. Under these conditions,
a radiative pumping scheme can operate for the 22-GHz transition (Gray et al.,
2022). This possibility is attractive in that the IR burst is destructive for
the existing masers in hotter, more rarefied conditions, but can explain the
rapid rise of the flare (over 20-30 d, Bayandina et al. 2022) in a previously
maser-free site. The second possibility is that the dust shrouding of the
young stellar object in MM1 is patchy, allowing radiation approximating to its
raw SED to escape along certain lines of sight. Assuming most radiation with
$\lambda<111$ nm, is lost in dissociating or ionizing hydrogen, there is still
a band $145-111$ nm that can heat gas efficiently due to photoelectric
emission from grains. The 9 best burst models in Stecklum et al. (2021) (their
Table A6) have a mean stellar temperature, $T_{*}$ of 22 570 K and mean radius
of 10.31 R☉. From formulae in Section 30.2 of Draine (2011) such photoelectric
emission deposits energy at an approximate rate of $9.28\times
10^{-27}n_{\mathrm{H}}(n_{\gamma}/3\times 10^{-3}\mathrm{cm^{-3}})$ erg cm-3
s-1, where $n_{\gamma}$ is the number density of the UV photons and
$n_{\mathrm{H}}$ is the number density of H-nuclei. The radiation field, at
$T_{*}$, diluted to a distance of $13\,500$ au, has $n_{\gamma}\simeq 270$
cm-3 of photons in the 145-111 nm band, so for a number density
$n_{\mathrm{H}}=2\times 10^{8}$ cm-3, typical of the pre-shock gas in the
C-type models used in the present work, the photoelectric heating rate is
$1.64\times 10^{-14}$ J m-3 s-1. An increase, $\Delta T$, in the gas thermal
energy due to this heating rate, $\Gamma_{\mathrm{p}e}$, follows $\Delta
T=2\Gamma_{\mathrm{p}e}t/(3k_{\mathrm{B}}n_{\mathrm{H}_{2}})$ after $t$
seconds. From this formula, a temperature increase of 800 K requires only a
little over 100 s. There remains the question of how quickly this energy can
be transferred from the photoelectrons to the H2 gas. For electron knietic
energies of a few eV, the elastic collision cross-section for collisions with
H2 are typically $>$10-19 m2 (Yoon et al., 2008), so the mean-free time for
the electrons is $1.7\times 10^{5}/E_{\rm eV}^{1/2}$ s. Energy can therefore
be transferred from the radiation field to the H2 gas within a few days. The
rate of energy deposition will, of course, fall as the radiation progresses
into the gas. We use an extinction cross section due to dust of $8.16\times
10^{-26}$ m2 at 100 nm to calculate an optical depth at this wavelength of
$\tau_{1000}=2.441L_{\mathrm{AU}}$, so the energy is mostly deposited a layer
of astronomical-unit scale, linking ultraviolet extinction by dust to the
scale of model clouds studied in the present work. Since the time required for
heating is much less than the sound crossing time of the heated region, the
heated gas will become hot and overpressured with respect to its surroudings.
Under these conditions it is very likely that shock expansion will follow into
cooler gas further from the radiation source. In this way, an initial
radiative event caused by accretion can ultimately drive a collision-based
pumping scheme after a delay little longer than the light travel time to the
new maser zone.

##  6 Conclusions

Shock-wave modification of gaseous clouds of scale a few au to a few tens of
au provides a mechanism that can explain maser flares with the largest
observed brightness temperatures ($\gtrsim 10^{17}$ K). The mechanism is
limited to transitions with a mainly collisional pumping mechanism, and the
output is highly directional, dominated by directions parallel to the shock
front, but this effect is less pronounced for C-type shocks. Magnetically
dominated (C-type) shocks can generate masers with similar flux densities to
the ideal hydrodynamic variety, and the magnetic types typically produce
flares with an extremely large variability index, since these models have a
minimal inversion in the 22-GHz transition of the maser species (ortho-H2O) in
the pre-shock fluid. C-type shocks also produce flares of shorter duration,
particular during the rise to peak, because faster shocks are required to
achieve the same (ultimate) compression. Density of the maser species in the
post-shock gas is important, but only in the sense that it needs to be in a
‘sensible’ range that is sufficient to achieve saturation, but low enough to
avoid quenching the pump. Maser beaming pattens show a strong asymmetry, being
considerably more extended parallel to the shock front. We find consistency
with earlier work that suggests that the observed image, perpendicular to the
shock front, in a shock-generated maser is much closer to the true size of the
supporting cloud than in the quiescent spherical case.

## Acknowledgments

MDG and SE acknowledge funding from the UK Science and Technology Facilities
Council (STFC) as part of the consolidated grant ST/P000649/1 to the Jodrell
Bank Centre for Astrophysics at the University of Manchester. MDG acknowledges
financial support from the National Astronomical Research Institute of
Thailand (NARIT) whilst on sabbatical at their HQ in Chiang Mai, Thailand.
This work was performed, in part, using the DiRAC Data Intensive service at
Leicester, operated by the University of Leicester IT Services, which forms
part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was
funded by BEIS capital funding via STFC capital grants ST/K000373/1 and
ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of
the National e-Infrastructure. Data used in this work was generated under
DiRAC award dp124. The authors would also like to thank the management of the
Chalawan supercomputer at NARIT.

## Data Availability

The data underlying this article will be shared on reasonable re- quest to the
corresponding author.

## References

  * Abramowitz & Stegun (1972) Abramowitz M., Stegun I. A., 1972, Handbook of Mathematical Functions
  * Araya et al. (2010) Araya E. D., Hofner P., Goss W. M., Kurtz S., Richards A. M. S., Linz H., Olmi L., Sewiło M., 2010, ApJ, 717, L133
  * Asanok et al. (2010) Asanok K., Etoka S., Gray M. D., Thomasson P., Richards A. M. S., Kramer B. H., 2010, MNRAS, 404, 120
  * Ashimbaeva et al. (2020a) Ashimbaeva N. T., Colom P., Krasnov V. V., Lekht E. E., Pashchenko M. I., Rudnitskii G. M., Tolmachev A. M., 2020a, Astronomy Reports, 64, 586
  * Ashimbaeva et al. (2020b) Ashimbaeva N. T., Colom P., Krasnov V. V., Lekht E. E., Pashchenko M. I., Rudnitskii G. M., Tolmachev A. M., 2020b, Astronomy Reports, 64, 839
  * Ashimbaeva et al. (2016) Ashimbaeva N. T., Colom P., Lekht E. E., Pashchenko M. I., Rudnitskii G. M., Tolmachev A. M., 2016, Astronomy Letters, 42, 652
  * Ashimbaeva et al. (2018) Ashimbaeva N. T., Colom P., Lekht E. E., Pashchenko M. I., Rudnitskii G. M., Tolmachev A. M., 2018, Astronomy Reports, 62, 609
  * Ashimbaeva et al. (2020) Ashimbaeva N. T., Krasnov V. V., Lekht E. E., Pashchenko M. I., Rudnitskii G. M., Tolmachev A. M., 2020, Astronomy Reports, 64, 15
  * Barlow et al. (1996) Barlow M. J., Nguyen-Q-Rieu et al. T., 1996, A&A, 315, L241
  * Bayandina et al. (2022) Bayandina O. S., Brogan C. L., Burns R. A., Caratti o Garatti A., Chibueze J. O., van den Heever S. P., Kurtz S. E., MacLeod G. C., Moscadelli L., Sobolev A. M., et al. 2022, A&A, 664, A44
  * Bayandina et al. (2019) Bayandina O. S., Burns R. A., Kurtz S. E., Shakhvorostova N. N., Val’tts I. E., 2019, ApJ, 884, 140
  * Bergman & Humphreys (2020) Bergman P., Humphreys E. M. L., 2020, A&A, 638, A19
  * Burns et al. (2016) Burns R. A., Handa T., Nagayama T., Sunada K., Omodaka T., 2016, MNRAS, 460, 283
  * Burns et al. (2023) Burns R. A., Uno Y., Sakai N., Blanchard J., Rosli Z., Orosz G., Yonekura Y., Tanabe Y., Sugiyama K., Hirota T., et al. 2023, Nature Astronomy, 7, 557
  * Chen et al. (2020) Chen X., Sobolev A. M., Ren Z.-Y., Parfenov S., Breen S. L., Ellingsen S. P., Shen Z.-Q., Li B., MacLeod G. C., Baan W., et al. 2020, Nature Astronomy, 4, 1170
  * Chen & Cai (2001) Chen Y., Cai D., 2001, App. Maths. & Computation, 124, 351
  * Chibueze et al. (2021) Chibueze J. O., MacLeod G. C., Vorster J. M., Hirota T., Brogan C. L., Hunter T. R., van Rooyen R., 2021, ApJ, 908, 175
  * Colom et al. (2019) Colom P., Ashimbaeva N. T., Lekht E. E., Pashchenko M. I., Rudnitskii G. M., Krasnov V. V., Tolmachev A. M., 2019, Astronomy Reports, 63, 814
  * Colom et al. (2021) Colom P., Ashimbaeva N. T., Lekht E. E., Pashchenko M. I., Rudnitskij G. M., Krasnov V. V., Tolmachev A. M., 2021, MNRAS, 507, 3285
  * Colom et al. (2016) Colom P., Lekht E. E., Pashchenko M. I., Rudnitskii G. M., Tolmachev A. M., 2016, Astronomy Reports, 60, 730
  * Coppola et al. (2016) Coppola C. M., Mizzi G., Bruno D., Esposito F., Galli D., Palla F., Longo S., 2016, MNRAS, 457, 3732
  * Cragg et al. (1992) Cragg D. M., Johns K. P., Godfrey P. D., Brown R. D., 1992, MNRAS, 259, 203
  * De Ceuster et al. (2022) De Ceuster F., Ceulemans T., Srivastava A., Homan W., Bolte J., Yates J., Decin L., Boyle P., Hetherington J., 2022, The Journal of Open Source Software, 7, 3905
  * De Ceuster et al. (2020) De Ceuster F., Homan W., Yates J., Decin L., Boyle P., Hetherington J., 2020, MNRAS, 492, 1812
  * de Jong (1973) de Jong T., 1973, A&A, 26, 297
  * Deguchi & Watson (1989) Deguchi S., Watson W. D., 1989, ApJ, 340, L17
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  * Elitzur & Asensio Ramos (2006) Elitzur M., Asensio Ramos A., 2006, MNRAS, 365, 779
  * Elitzur et al. (1989) Elitzur M., Hollenbach D. J., McKee C. F., 1989, ApJ, 346, 983
  * Elitzur et al. (1992) Elitzur M., Hollenbach D. J., McKee C. F., 1992, ApJ, 394, 221
  * Elitzur et al. (1991) Elitzur M., McKee C. F., Hollenbach D. J., 1991, ApJ, 367, 333
  * Faure & Josselin (2008) Faure A., Josselin E., 2008, A&A, 492, 257
  * Field et al. (1994) Field D., Gray M. D., de St. Paer P., 1994, A&A, 282, 213
  * Frisch (1988) Frisch H., 1988, in Saas-Fee Advanced Course 18: Radiation in Moving Gaseous Media Radiative Transfer with Frequency Redistribution - Asymptotic methods scaling laws and approximate solutions. p. 339
  * Goedhart et al. (2004) Goedhart S., Gaylard M. J., van der Walt D. J., 2004, MNRAS, 355, 553
  * Goldreich & Kwan (1974) Goldreich P., Kwan J., 1974, ApJ, 190, 27
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
  * Gray (2012) Gray M. D., 2012, Maser Sources in Astrophysics. Cambridge University Press, Cambridge, UK
  * Gray et al. (2019) Gray M. D., Baggott J., Westlake J., Etoka S., 2019, MNRAS, 486, 4216 (Paper 2)
  * Gray et al. (2016) Gray M. D., Baudry A., Richards A. M. S., Humphreys E. M. L., Sobolev A. M., Yates J. A., 2016, MNRAS, 456, 374
  * Gray et al. (2020) Gray M. D., Etoka S., Pimpanuwat B., 2020, MNRAS, 498, L11
  * Gray et al. (2022) Gray M. D., Etoka S., Richards A. M. S., Pimpanuwat B., 2022, MNRAS, 513, 1354
  * Gray et al. (2020) Gray M. D., Etoka S., Travis A., Pimpanuwat B., 2020, MNRAS, 493, 2472 (Paper 3)
  * Gray et al. (2018) Gray M. D., Mason L., Etoka S., 2018, MNRAS, 477, 2628 (Paper 1)
  * Gwinn et al. (1992) Gwinn C. R., Moran J. M., Reid M. J., 1992, ApJ, 393, 149
  * Hirota et al. (2014) Hirota T., Tsuboi M., Kurono Y., Fujisawa K., Honma M., Kim M. K., Imai H., Yonekura Y., 2014, PASJ, 66, 106
  * Hollenbach et al. (2013) Hollenbach D., Elitzur M., McKee C. F., 2013, ApJ, 773, 70
  * Imai et al. (1999) Imai H., Iwata T., Miyoshi M., 1999, PASJ, 51, 473
  * Inayoshi et al. (2013) Inayoshi K., Sugiyama K., Hosokawa T., Motogi K., Tanaka K. E. I., 2013, ApJ, 769, L20
  * Kaufman & Neufeld (1996) Kaufman M. J., Neufeld D. A., 1996, ApJ, 456, 250
  * Kim et al. (2013) Kim J.-S., Kim S.-W., Kurayama T., Honma M., Sasao T., Surcis G., Cantó J., Torrelles J. M., Kim S. J., 2013, ApJ, 767, 86
  * King & Florance (1964) King J. I. F., Florance E. T., 1964, ApJ, 139, 397
  * Kobayashi et al. (2000) Kobayashi H., Shimoikura T., Omodaka T., Diamond P. J., 2000, in Hirabayashi H., Edwards P. G., Murphy D. W., eds, Astrophysical Phenomena Revealed by Space VLBI Monitoring of the Orion-KL Water Maser Outburst. pp 109–112
  * Krasnov et al. (2018) Krasnov V. V., Lekht E. E., Minnebaev V. M., Pashchenko M. I., Rudnitskii G. M., Tolmachev A. M., 2018, Astronomy Reports, 62, 200
  * Krasnov et al. (2015) Krasnov V. V., Lekht E. E., Rudnitskii G. M., Pashchenko M. I., Tolmachev A. M., 2015, Astronomy Letters, 41, 517
  * Lees (1973) Lees R. M., 1973, ApJ, 184, 763
  * MacLeod et al. (2018) MacLeod G. C., Smits D. P., Goedhart S., Hunter T. R., Brogan C. L., Chibueze J. O., van den Heever S. P., Thesner C. J., Banda P. J., Paulsen J. D., 2018, MNRAS, 478, 1077
  * MacLeod et al. (2021) MacLeod G. C., Smits D. P., Green J. A., van den Heever S. P., 2021, MNRAS, 502, 5658
  * McCarthy et al. (2023) McCarthy T. P., Breen S. L., Kaczmarek J. F., Chen X., Parfenov S., Sobolev A. M., Ellingsen S. P., Burns R. A., MacLeod G. C., Sugiyama K., et al. 2023, MNRAS, 522, 4728
  * Melnick et al. (2000) Melnick G. J., et al. 2000, ApJ, 539, L87
  * Moscadelli et al. (2022) Moscadelli L., Sanna A., Beuther H., Oliva A., Kuiper R., 2022, Nature Astronomy, 6, 1068
  * Moscadelli et al. (2019) Moscadelli L., Sanna A., Goddi C., Krishnan V., Massi F., Bacciotti F., 2019, A&A, 631, A74
  * Nakamura et al. (2021) Nakamura M., Motogi K., Fujisawa K., 2021, in 13th East Asian VLBI Workshop 2021 A kinematic study of the disk-outflow system around a high-mass protostar G59.783+0.065 probed by CH3OH and H2O masers. p. 3
  * Nedoluha & Watson (1991) Nedoluha G. E., Watson W. D., 1991, ApJ, 367, L63
  * Olech et al. (2020) Olech M., Szymczak M., Wolak P., Gérard E., Bartkiewicz A., 2020, A&A, 634, A41
  * Parfenov & Sobolev (2014) Parfenov S. Y., Sobolev A. M., 2014, MNRAS, 444, 620
  * Prozesky & Smits (2020) Prozesky A., Smits D. P., 2020, MNRAS, 491, 2536
  * Rajabi & Houde (2017) Rajabi F., Houde M., 2017, Science Advances, 3, e1601858
  * Rajabi et al. (2019) Rajabi F., Houde M., Bartkiewicz A., Olech M., Szymczak M., Wolak P., 2019, MNRAS, 484, 1590
  * Richards et al. (2005) Richards A. M. S., Cohen R. J., Crocker M., Lekht E. E., Mendoza E., Samodurov V. A., 2005, Ap&SS, 295, 19
  * Richards et al. (2011) Richards A. M. S., Elitzur M., Yates J. A., 2011, A&A, 525, A56
  * Richards et al. (2012) Richards A. M. S., Etoka S., Gray M. D., Lekht E. E., Mendoza-Torres J. E., Murakawa K., Rudnitskij G., Yates J. A., 2012, A&A, 546, A16
  * Royer et al. (2010) Royer P., Decin L., Wesson R. e. a., 2010, A&A, 518, L145
  * Salii et al. (2022) Salii S. V., Zinchenko I. I., Liu S.-Y., Sobolev A. M., Aberfelds A., Su Y.-N., 2022, MNRAS, 512, 3215
  * Sanna et al. (2010) Sanna A., Moscadelli L., Cesaroni R., Tarchi A., Furuya R. S., Goddi C., 2010, A&A, 517, A78
  * Schöier et al. (2005) Schöier F. L., van der Tak F. F. S., van Dishoeck E. F., Black J. H., 2005, A&A, 432, 369
  * Shakhvorostova et al. (2018) Shakhvorostova N. N., Vol’vach L. N., Vol’vach A. E., Dmitrotsa A. I., Bayandina O. S., Val’tts I. E., Alakoz A. V., Ashimbaeva N. T., Rudnitskii G. M., 2018, Astronomy Reports, 62, 584
  * Shimoikura et al. (2005) Shimoikura T., Kobayashi H., Omodaka T., Diamond P. J., Matveyenko L. I., Fujisawa K., 2005, ApJ, 634, 459
  * Slysh et al. (2010) Slysh V. I., Alakoz A. V., Migenes V., 2010, MNRAS, 404, 1121
  * Sobolev et al. (2018) Sobolev A. M., Moran J. M., Gray M. D., Alakoz A., Imai H., Baan W. A., Tolmachev A. M., Samodurov V. A., Ladeyshchikov D. A., 2018, ApJ, 856, 60
  * Stecklum et al. (2021) Stecklum B., Wolf V., Linz H., Caratti o Garatti A., Schmidl S., Klose S., Eislöffel J., Fischer C., Brogan C., Burns R. A., et al. 2021, A&A, 646, A161
  * Strelnitski (2007) Strelnitski V., 2007, in Haverkorn M., Goss W. M., eds, SINS - Small Ionized and Neutral Structures in the Diffuse Interstellar Medium Vol. 365 of Astronomical Society of the Pacific Conference Series, AU Structures in Turbulent Outflows from YSOs Revealed by H2O Masers. p. 196
  * Strelnitski et al. (2017) Strelnitski V. S., Holder B. P., Shishov V. I., Nezhdanova N. I., 2017, Astronomical and Astrophysical Transactions, 30, 173
  * Strelnitskii (1982) Strelnitskii V. S., 1982, Soviet Astronomy Letters, 8, 86
  * Surcis et al. (2011) Surcis G., Vlemmings W. H. T., Curiel S., Hutawarakorn Kramer B., Torrelles J. M., Sarma A. P., 2011, A&A, 527, A48
  * Surcis et al. (2014) Surcis G., Vlemmings W. H. T., van Langevelde H. J., Goddi C., Torrelles J. M., Cantó J., Curiel S., Kim S. W., Kim J. S., 2014, A&A, 565, L8
  * Tennyson et al. (2016) Tennyson J., Hulme K., Naim O. K., Yurchenko S. N., 2016, Journal of Physics B Atomic Molecular Physics, 49, 044002
  * Uscanga et al. (2005) Uscanga L., Cantó J., Curiel S., Anglada G., Torrelles J. M., Patel N. A., Gómez J. F., Raga A. C., 2005, ApJ, 634, 468
  * van der Walt (2011) van der Walt D. J., 2011, AJ, 141, 152
  * van Dishoeck et al. (2021) van Dishoeck E. F., Kristensen L. E., Mottram J. C., Benz A. O., Bergin E. A., Caselli P., Herpin F., Hogerheijde M. R., Johnstone D., Liseau e. a., 2021, A&A, 648, A24
  * Volvach et al. (2023) Volvach A. E., Volvach L. N., Larionov M. G., 2023, ApJ, 955, 10
  * Vol’vach et al. (2019) Vol’vach L. N., Vol’vach A. E., Larionov M. G., MacLeod G. C., van den Heever S. P., Wolak P., Olech M., Ipatov A. V., Ivanov D. V., Mikhailov A. G., 2019, Astronomy Reports, 63, 49
  * Volvach et al. (2019) Volvach L. N., Volvach A. E., Larionov M. G., MacLeod G. C., Wolak P., Olech M., Kramer B., Menten K., Kraus A., Brand J., Zanicelli A., Poppi S., Rigini S., 2019, Astronomy Letters, 45, 321
  * Volvach et al. (2019) Volvach L. N., Volvach A. E., Larionov M. G., Wolak P., Kramer B., Menten K., Kraus A., Brand J., Zanichelli A., Poppi S., Rigini S., Ipatov A. V., Ivanov D. V., Mikhailov A. G., Mel’nikov A., 2019, Astronomy Reports, 63, 652
  * Voronkov et al. (2005) Voronkov M. A., Sobolev A. M., Ellingsen S. P., Ostrovskii A. B., 2005, MNRAS, 362, 995
  * Xu et al. (2013) Xu Y., Li J. J., Reid M. J., Menten K. M., Zheng X. W., Brunthaler A., Moscadelli L., Dame T. M., Zhang B., 2013, ApJ, 769, 15
  * Yates et al. (1997) Yates J. A., Field D., Gray M. D., 1997, MNRAS, 285, 303
  * Yoon et al. (2008) Yoon J.-S., Song M.-Y., Han J.-M., Hwang S. H., Chang W.-S., Lee B., Itikawa Y., 2008, Journal of Physical and Chemical Reference Data, 37, 913

##  Appendix A Variable Velocity Development

With reference to equation (2), we define here the symbols used. The velocity
width parameter is $W(\tau)=\sqrt{2k_{\mathrm{B}}T_{\mathrm{K}}(\tau)/m_{X}}$,
and its domain average by node is $\bar{W}$, where $k_{\mathrm{B}}$ is
Boltzmann’s constant, $T_{\mathrm{K}}$ is the kinetic temperature and $m_{X}$
is the molecular mass of species $X$. Scaled distance along a ray, propagating
in direction $\hat{\mn@boldsymbol{n}}$, is represented by the optical depth
$\tau$. The parameter $\eta(\tau)$ is $\bar{W}/W(\tau)$, and dimensionless
frequency and velocity are respectively denoted by $\bar{\nu}$ and
$\mn@boldsymbol{u}$. The background dimensionless specific intensities,
$i_{\mathrm{BG}}(\Omega)$, are also multiples of the saturation intensity, and
here we allow for the possibility that the background level may vary with
direction ($\Omega$) on the sky surrounding the domain.

As in Papers 1-3, a nodal solution for the inversion is found through a
saturation term that involves a frequency-and-angle averaged intensity that is
itself eliminated in favour of an expression in terms of the inversions in
nodes that border the ray paths converging on the node in question. In the
earlier papers, it proved possible to perform the frequency averaging
analytically, but this cannot be applied with variable velocity: although the
analytical frequency integral can still be carried out, the resulting
expression cannot be factored to restore multiple integrals (from the power-
series expansion of an exponential) to the product of identical integrals,
each in a single spatial variable. It is therefore not possible to reduce the
expression to a single power series with one spatial integral as the argument.
We therefore adopt the rather different strategy described here.

A formal frequency average of equation (2) over a gaussian molecular response
generates the averaged intensity,

$\displaystyle i(\tau)$
$\displaystyle=\frac{i_{\mathrm{BG}}(\Omega)}{\mathrm{\pi}^{1/2}}\int_{-\infty}^{\infty}d\bar{\nu}e^{-(\bar{\nu}-\hat{\mn@boldsymbol{n}}\cdot\mn@boldsymbol{u}(\tau))^{2}}$
$\displaystyle\times\exp\left\\{\int_{\tau_{0}}^{\tau}d\tau^{\prime}\delta^{\prime}(\tau^{\prime})e^{-(\bar{\nu}-\hat{\mn@boldsymbol{n}}\cdot\mn@boldsymbol{u}(\tau^{\prime}))^{2}}\right\\},$
(20)

where we have now assumed that $\eta$, the width parameter, is 1 everywhere,
equivalent to assuming that the domain is isothermal. The frequency integral
in equation (20) is now modified by defining the new dimensionless frequency
$\varpi=\bar{\nu}-\hat{\mn@boldsymbol{n}}_{q}\cdot\mn@boldsymbol{u}(\mn@boldsymbol{r})$.
This change of variable is applied to make the gaussian in equation (20) zero-
centered, as convenient for numerical quadrature. With the further definition
$\Delta\mn@boldsymbol{u}(\tau^{\prime})=\mn@boldsymbol{u}(\mn@boldsymbol{r})-\mn@boldsymbol{u}(\tau^{\prime})$,
and specifying the intensity of a particular ray, $q$, equation (20) becomes

$\displaystyle i_{q}(\tau)$
$\displaystyle=\frac{i_{\mathrm{BG},q}}{\mathrm{\pi}^{1/2}}\\!\\!\int_{-\infty}^{\infty}\\!\\!\exp\left\\{\int_{\tau_{0}}^{\tau}\delta^{\prime}(\tau^{\prime})e^{-(\varpi+\hat{\mn@boldsymbol{n}}_{q}\cdot\Delta\mn@boldsymbol{u}(\tau^{\prime}))^{2}}d\tau^{\prime}\right\\}$
$\displaystyle\times e^{-\varpi^{2}}d\varpi.$ (21)

The exponential in braces in equation (21) contains a spatial integral that is
discretized according to a finite element scheme similar to that used in Paper
2 and Paper 3. We first define

$Y_{q}=\sum_{j=1}^{J(q)}\int_{\tau_{\mathrm{in},j}}^{\tau_{\mathrm{out},j}}\delta^{\prime}(\tau^{\prime})e^{-(\varpi+\hat{\mn@boldsymbol{n}}_{q}\cdot\mn@boldsymbol{u}(\tau^{\prime}))^{2}}d\tau^{\prime},$
(22)

where the expression in equation (22) shows the integral along ray $q$
separated into the $J(q)$ elemental contributions that appear between the
entry point to the domain and arrival at a target node. The integral in
equation (22) is now over the path through a single element, $j$, with entry
point $\tau_{\mathrm{in},j}$, and the corresponding exit point
$\tau_{\mathrm{out},j}$. The expression for $Y_{q}$ is equivalent to the
argument of the main exponential in equation (21), so this can now be written
in the much simpler form,

$i_{q}(\tau)=(i_{\mathrm{BG},q}/\mathrm{\pi}^{1/2})\int_{-\infty}^{\infty}e^{-\varpi^{2}}\exp\left\\{Y_{q}(\varpi)\right\\}d\varpi.$
(23)

If the $j$th element has a shape function for its $i$th node of the form,

$f_{i}(x,y,z)=a_{i}+b_{i}x+c_{i}y+d_{i}z,$ (24)

where $x,y,z$ are dimensionless coordinates within the element, and the
quantities with the subscript $i$ are shape-function coefficients, then
$f_{i}$ can be parameterized in terms of a single variable $0\leq t\leq 1$,
for example,

$f_{i}(t)=a_{i}+b_{i}(x_{\mathrm{in}}+t\Delta x)+c_{i}(y_{\mathrm{in}}+t\Delta
y)+d_{i}(z_{\mathrm{in}}+t\Delta z),$ (25)

where $(x_{\mathrm{in}},y_{\mathrm{in}},z_{\mathrm{in}})$ are the coordinates
of the entry point $\tau_{\mathrm{in},j}$, and $(\Delta x,\Delta y,\Delta x)$
are the coordinate differences between the exit and entry points. Summing
terms that do not multiply $t$ as
$A_{i}=a_{i}+b_{i}x_{\mathrm{in}}+c_{i}y_{\mathrm{in}}+d_{i}z_{\mathrm{in}}$,
and those that do, as $B_{i}=b_{i}\Delta x+c_{i}\Delta y+d_{i}\Delta z$, the
$i$th shape function becomes simply,

$f_{i}(t)=A_{i}+B_{i}t.$ (26)

Any spatially variable quantity can then be found, in the interior of element
$j$, by summing equation (26), multiplied by nodal values of the quantity,
over the nodes. For example, for some ray distance $\tau_{\mathrm{in},j}\leq
t\leq\tau_{\mathrm{out},j}$, the dimensionless bulk velocity is approximated
by

$\mn@boldsymbol{u}(t)=\sum_{m=1}^{M}\mn@boldsymbol{u}_{m}(A_{m}+B_{m}t),$ (27)

where $\mn@boldsymbol{u}_{m}$ is the velocity at the $m$th mode, and each
element has $M$ nodes. A similar expansion to equation (27) is made for the
unknown, $\delta^{\prime}(t)$.

We now use equation (27) and its analogue for $\delta^{\prime}(t)$ to write
$Y_{q}$ in terms of nodal values of quantities that vary across the domain.
With the aid of the further definitions,

$\alpha_{qj}^{(\varpi)}=\varpi+\hat{\mn@boldsymbol{n}}_{q}\cdot\mn@boldsymbol{u}_{i}-\sum_{m=1}^{M}A_{m}\hat{\mn@boldsymbol{n}}_{q}\cdot\mn@boldsymbol{u}_{m},$
(28)

where $\mn@boldsymbol{u}_{i}$ is the velocity at the target node, and

$\beta_{qj}=\sum_{m=1}^{M}B_{m}\hat{\mn@boldsymbol{n}}_{q}\cdot\mn@boldsymbol{u}_{m},$
(29)

equation (22) reduces to the form, with $\tau_{j}$ as the distance covered in
passing through element $j$,

$Y_{q}=\sum_{j=1}^{J(q)}\tau_{j}\sum_{\mu=1}^{M}\delta^{\prime}_{\mu}\int_{0}^{1}(A_{j\mu}+B_{j\mu}t)e^{-(\alpha_{qj}^{(\varpi)}-\beta_{qj}t)^{2}}dt.$
(30)

There are various ways of carrying out the integral in equation (30): perhaps
the most satisfactory is to extract a factor of
$\exp-[\alpha_{qj}^{(\varpi)}]^{2}$ from the integral, while defining new
constants $a=\beta_{qj}^{2}$ and $b=\alpha_{qj}^{(\varpi)}\beta_{qj}$,
reducing equation (30) to the form

$Y_{q}=\\!\sum_{j=1}^{J(q)}\tau_{j}e^{-[\alpha_{qj}^{(\varpi)}]^{2}}\\!\\!\sum_{\mu=1}^{M}\delta^{\prime}_{\mu}\\!\\!\int_{0}^{1}\\!\\!(A_{j\mu}+B_{j\mu}t)e^{-(at^{2}-2bt)}dt,$
(31)

that is tractable in terms of standard integrals. Of particular interest is
the formula from Abramowitz & Stegun (1972), reproduced here as

$\int
e^{-(ax^{2}+2bx+c)}dx=\frac{\mathrm{\pi}^{1/2}e^{-c}}{2a^{1/2}}e^{b^{2}/a}\mathrm{erf}\left\\{\sqrt{a}x+\frac{b}{\sqrt{a}}\right\\}+\kappa,$
(32)

where $\kappa$ is a constant of integration. One integral is exactly of this
form, and the other is reducible to this same integral through integration by
parts. Results are in terms of exponentials and error functions, and the
integrated form of equation (31) is

$\displaystyle Y_{q}$
$\displaystyle=\sum_{j=1}^{J(q)}\sum_{\mu=1}^{M}\frac{\tau_{j}\delta^{\prime}_{\mu}B_{j\mu}}{2\beta_{qj}}\left\\{\frac{e^{-(\alpha_{qj}+\varpi)^{2}}-e^{-(\gamma_{qj}+\varpi)^{2}}}{\beta_{qj}}\right.$
$\displaystyle\left.+\frac{\mathrm{\pi}^{1/2}}{\beta_{qj}}(\varepsilon_{qj}+\varpi)\left[\mathrm{erf}(\alpha_{qj}+\varpi)-\mathrm{erf}(\gamma_{qj}+\varpi)\right]\right\\},$
(33)

where further new constants are $\alpha_{qj}=\alpha_{qj}^{(\varpi)}-\varpi$,
$\gamma_{qj}=\alpha_{qj}-\beta_{qj}$ and
$\varepsilon_{qj}=\alpha_{qj}+A_{j\mu}\beta_{qj}/B_{j\mu}$. If the double sum
over elements and local nodes in equation (33) is replaced by a single sum
over the global node numbers of those nodes that bound ray $q$, noting that
such nodes may be members of more than one element in $J(q)$, then we can pre-
compute the coefficients, $\Phi_{j,k}^{q,i}$, where $j$ is now the global node
index. These coefficients include everything in equation (33) except the nodal
values of the inversion itself, $\delta^{\prime}_{\mu}$. It is now possible to
write equation (33) in the simple form
$Y_{q}=\sum_{j=1}^{J(q)}\delta^{\prime}_{j}\Phi_{j,k}^{q,i}$. We substitute
this expression for $Y_{q}$ in equation (23), noting that $Y_{q}$ is a
function of $\varpi$ through the cofficients $\Phi_{j,k}^{q,i}$. We now
replace the frequency integral with a Gauss-Hermite numerical quadrature with
$K$ abcissae, $\varpi_{k}$, and corresponding weights $\zeta_{k}$, with $k$ in
the range 1 to $K$. We generally use the Gauss-Hermite quadrature scheme, as
this is designed to approximate integrals that are of a form where the
integrand is a product of a gaussian and another function, such that

$\int_{-\infty}^{\infty}e^{-\varpi^{2}}f(\varpi)d\varpi\simeq\sum_{k=1}^{K}\zeta_{k}f(\varpi_{k}).$
(34)

The resulting form of equation (23), discretized in frequency, is

$i_{q}(\tau)=(i_{\mathrm{BG},q}/\mathrm{\pi}^{1/2})\sum_{k=1}^{K}\zeta_{k}\exp\left\\{\sum_{j=1}^{J(q)}\delta^{\prime}_{j}\Phi_{j,k}^{q,i}\right\\}.$
(35)

A formal average of equation (35) over the $q$ rays with their associated
solid angles results in the mean intensity, $\bar{j}(\mn@boldsymbol{r})$, at
the target node, which is given by

$\bar{j}(\mn@boldsymbol{r})\\!=\\!\frac{1}{4\mathrm{\pi}^{3/2}}\\!\sum_{q=1}^{Q}\frac{i_{\mathrm{BG},q}{\cal
A}_{q}}{l_{q}^{2}}\\!\times\\!\sum_{k=1}^{K}\zeta_{k}\exp\\!\left[\tau_{\mathrm{M}}\\!\sum_{j=1}^{J(q)}\\!\delta^{\prime}_{j}\Phi_{j,k}^{q,i}\right],$
(36)

where a global depth multiplier, $\tau_{\mathrm{M}}$, has been extracted from
every saturation coefficient.

##  Appendix B Absolute Saturation Parameters

We begin with the definition of the saturation intensity from equation (16) of
the main text. As in Appendix A of Paper 3, we use a form of the definition
where we assume that the A-value of the maser and collisional rates across the
maser transition itself contribute negligibly compared to the overall loss
rate, represented by the parameter $\Gamma$.

The main difficulty in equation (16) is evaluating $\Gamma$ to a reasonable
accuracy, since the other parameters, comprising the line-centre frequency
(21.5549 GHz), the Einstein A-value of the maser transition
($A_{ij}=1.835\times 10^{-9}$ s-1) and the statistical weight ratio
$g_{i}/g_{j}=13/11$, are all easily available from the ortho-H2O radex
datafile provided by the LAMBDA database (Schöier et al., 2005). The necessary
numbers to calculate $\Gamma$ can be extracted from the same database, noting
that $\Gamma$ is essentially the non-maser loss rate from the upper level of
the maser transition, $J_{K_{a},K,c}=6_{2,6}$, in the vibrational ground
state.

There is one radiative transition from the upper maser level (to $5_{1,5}$)
with a significant A-value of 0.758 Hz. However, it is clear that with number
densities typical of the shocked gas in the models used in the present work,
significant contributions to $\Gamma$ will also come from collisional rate
coefficients. We ignore upward radiative transitions from $6_{2,6}$ on the
grounds that there is negligible external pumping radiation. The file from the
LAMBDA database contains two tables of downward second-order collisional rate
coefficients: the first is for collisions with H2, the second, with electrons.
These tables are derived from work by Faure & Josselin (2008). Although the
rate coefficients in the electron table are typically much stronger than those
from H2, by approximately 3 orders of magnitude, we do not consider electron
collisions further, since their effect is marginal, even with an electron
abundance as high as 10-4 with respect to H2. The sum of the downward rate
coefficients from $6_{2,6}$ is $3.168\times 10^{-10}$ cm3 s-1, using the
figures for 800 K. A similar total of $3.379\times 10^{-10}$ cm3 s-1 arises
from summing the first 20 upward rate coefficients. Transitions to higher
levels add little due to a general decay in the magnitude of the coefficients
and increasingly adverse Boltzmann exponentials. Summing these collisional
terms and the single strong A-value yields

$\Gamma=0.758+0.06547n_{8}\;\;\mathrm{s^{-1}}$ (37)

where $n_{8}$ is the number density of H2 in multiples of 108 cm-3.

Equation 37 has the slightly unfortunate effect of making the loss rate
density, and therefore model, dependent. As an example, we consider the model
described in Section 4.4. In this case $\tau_{\mathrm{M}}=3$ and the
corresponding number density of o-H2O in the unshocked gas from Fig. 2 is 440
cm-3. As the models in Gray et al. (2016) use a standard ratio of $3\times
10^{-5}$ for o-H2O to H2, we obtain a value of $n_{8}=0.147$. However, as
saturation is only important in the shocked gas we multiply this figure by 9,
the compression factor, noting that the post-shock inversion follows the
density in this model, but not in the others. Therefore, we arrive at
$n_{8}=1.32$. With the aid of equation (37) we then obtain $\Gamma=0.844$ s-1
and a saturation intensity from equation (16) of $3.11\times 10^{-11}$ W m-2
Hz-1.

###  B.1 Conversion of flux density to Jy

The starting point here is equation (A3) of Paper 3, which may be written
$f_{\mathrm{kpc}}=f_{\nu}(R_{\mathrm{AU}}/206265d_{\mathrm{kpc}})^{2}$, where
$f_{\mathrm{kpc}}$ is the dimensionless flux density of a source of size
$R_{\mathrm{AU}}$ in astronomical units at a distance $d_{\mathrm{kpc}}$
kiloparsec, and $f_{\nu}$ is the flux density at the standard distance of 1000
domain units and scaled by the saturation intensity of the model. To obtain a
dimensioned flux density, we need only multiply the above equation by the
saturation intensity of a given model and divide by $10^{-26}$ to put the
result in Jy. For the example given above where $I_{\mathrm{s}}=3.11\times
10^{-11}$ W m-2 Hz-1, the dimensioned flux density becomes

$F_{\mathrm{k}pc}=7.31\times
10^{4}f_{\nu}(R_{\mathrm{AU}}/d_{\mathrm{kpc}})^{2}\;\;\mathrm{Jy},$ (38)

which is equation (18) of the main text, noting that the leading constant is
significantly larger for the 22-GHz H2O transition in the present work than
for the 6.7-GHz methanol transition considered in Paper 3 because of the
higher saturation intensity in the water transition.