11institutetext: Department of Biomedical Engineering, School of Biomedical
Engineering and Imaging Sciences, King’s College London, London, SE1 7EH, UK
22institutetext: Centre for the Developing Brain, Department of Perinatal
Imaging and Health, School of Biomedical Engineering and Imaging Sciences,
King’s College London, London, SE1 7EH, UK 33institutetext: Biomedical Image
Analysis Group, Department of Computing, Imperial College London, London, SW7
2AZ, UK 44institutetext: Klinikum rechts der Isar, Technical University of
Munich, Munich, 81675, Germany

# Improving Phenotype Prediction using Long-Range Spatio-Temporal Dynamics of
Functional Connectivity

Simon Dahan 11 Logan Z. J. Williams 1122 Daniel Rueckert 3344 Emma C. Robinson
1122

###### Abstract

The study of functional brain connectivity (FC) is important for understanding
the underlying mechanisms of many psychiatric disorders. Many recent analyses
adopt graph convolutional networks, to study non-linear interactions between
functionally-correlated states. However, although patterns of brain activation
are known to be hierarchically organised in both space and time, many methods
have failed to extract powerful spatio-temporal features. To overcome those
challenges, and improve understanding of long-range functional dynamics, we
translate an approach, from the domain of skeleton-based action recognition,
designed to model interactions across space and time. We evaluate this
approach using the Human Connectome Project (HCP) dataset on sex
classification and fluid intelligence prediction. To account for subject
topographic variability of functional organisation, we modelled functional
connectomes using multi-resolution dual-regressed (subject-specific) ICA
nodes. Results show a prediction accuracy of 94.4% for sex classification (an
increase of 6.2% compared to other methods), and an improvement of correlation
with fluid intelligence of 0.325 vs 0.144, relative to a baseline model that
encodes space and time separately. Results suggest that explicit encoding of
spatio-temporal dynamics of brain functional activity may improve the
precision with which behavioural and cognitive phenotypes may be predicted in
the future.

###### Keywords:

Human Connectome Project Graph Convolution Networks Temporal Dynamics
Functional MRI Phenotyping.

## Introduction

Functional MRI (fMRI) is widely recognised as a cornerstone technique for
relating brain processes to behaviour [24, 15]. In particular, the correlation
between regional time-series activations, referred to as functional
connectivity (FC), is known to be subject-specific [4, 11] and discriminate
between populations [23, 13]. This has led to it being popularly investigated
as a potential clinical biomarker for many psychiatric and neurological
disorders [13, 9]. Despite several promising studies on large cohorts, for
example, the Autism Spectrum Disorder study ABIDE [13], several recent works
have emphasised the need to better account for morphological, or topographic,
heterogeneity of cortical architecture when building functional networks [2,
12, 6]. This is considered especially important for disease modelling [17],
where phenotypes are often heterogeneous and may be classified into sub-types
[28]. Under such circumstances fitting all data to a single global-average
model of healthy functional organisation, stands to mask subtle features of
disease [17, 28], as well as obscure understanding of mechanisms of disorders
that usually co-occur [28]. For these reasons, increasing efforts are being
made to look towards individual subject-level analyses and move away from the
case-control approach [18].

To this end, recent approaches have looked towards better capture of inter-
subject variation when deriving functional networks. For example,
probabilistic functional modes (or PROFUMO) [8] uses a variational form of
independent component analysis (ICA) to explicitly account for subject-
specific variability during group-wise factorisation of fMRI data into
temporal and spatial modes; [11] implements a multi-session hierarchical
Bayesian model (MS-HBM) in order to cluster functional connectivity data into
group-wise comparable, but subject-specific, functional parcellations; and [6]
projects a hand-annotated, population-average multimodal parcellation of the
cortex onto individuals by training a multi-layer perceptron (MLP) classifier;
in the process revealing considerable topographic variation of regions when
observed across the group.

At the same time, many studies are moving away from modelling static
functional connectivity (averaged over all time) towards explicit modelling of
temporal dynamics, as a means of better capturing temporal markers of
behaviour [27, 5, 10, 19]. Of particular note, for this paper, are methods
which translate models used for video understanding [29], especially graph
convolution networks (GCN) used in skeleton-based action recognition [5]. In
this paper we specifically investigate the potential for translation of new
techniques which explicitly encode interactions, across nodes, in both space
and time [16]. The specific contributions of this paper are as follows:

  1. 1.

Explicit modelling of FC spatio-temporal dynamics improves prediction of
behavioural phenotypes: We demonstrate improvements over existing GCN models
of dynamic FC [5] by adopting a new model [16], which aggregates functional
signals in both space and time.

  2. 2.

Accounting for inter-subject cortical heterogeneity improves performance: We
compare models trained using a 22 region group-average parcellation and multi-
resolution dual-regressed (subject-specific) ICA nodes, and demonstrate gains
in performance by accounting for topographic variability of functional
organisation.

## Related Works

The use of GCNs to study FC supports learning of complex, non-linear, long-
range interactions between correlated functional network states. While early
approaches focused on the analysis of static FC [13, 14, 3], demonstrating
improved performance in metric learning for behavioural prediction of
neuropsychiatric disorders such as Autism Spectrum Disorder [13, 14]; more
recent papers have started to investigate the potential of GCNs to model FC
temporal dynamics [5, 10]. Of these, [5] takes inspiration from a skeleton-
based GCN model for video understanding [29], translating this to the domain
of dynamic FC modelling, by building a spatial graph from static FC (estimated
across all node time series) but self-connecting nodes temporally along a
regular grid. By applying spatial-only and temporal-only convolutions, the
method successfully captured long-range (spatial and temporal) dependencies
but failed to extract powerful space-time features by decoupling space and
time graphs. By contrast, the more recent MS-G3D of [16] goes further to
capture spatio-temporal semantics between nodes by encoding them through novel
graph convolutional units which aggregate across space and time; in this way
encoding complex actions, for example, the action of sitting down or standing
up, in which the top half of a skeleton moves before the bottom. In this
paper, we hypothesise that such long-range spatio-temporal behaviour is also
observed in brain functional dynamics. We, therefore, benchmark this method
against the method proposed by [5] as a means for quantifying this
effect.111The code for this experiment will be made available at:
http://www.github.com/metrics-lab/ST-fMRI/

## Methods

### Data

This study used functional data from the first session (15min - 1200 frames -
0.72s/frame) of the Human Connectome Project (HCP) S1200 release [25], where
individual subject node timeseries for 1003 participants were derived
following group-wise independent component analysis (Group-ICA) [1, 24]. In
brief, individual functional MRI data was processed through the HCP pipelines
[7, 24] and aligned using MSMAll [22, 21]. Then, spatial Group-ICA (with 468
subjects) was applied at different dimensionalities (15 to 300), factorising
the data into a set of spatial maps (functional nodes) and associated time
courses. These were propagated into individuals using dual regression, in this
way tailoring the brain parcellations and time courses at the individual level
(illustrated in figure 1). We used all available data for the sex
classification task (534 males and 469 females), but had to exclude 4 subjects
for the fluid intelligence prediction task (due to missing data). We used Penn
Progressive Matrices as scores of fluid intelligence, referred as
_PMAT24_A_CR_ in the HCP behavioral data. Similarly to [5], we also benchmark
against use of group-average functional parcellation, where the group average
HCP multimodal cortical parcellation (180 regions per hemisphere) [6] was
aggregated into 22 regions of interest (ROIs) per hemisphere by combining
nodes with common functions and large contiguous areas. In this latter
configuration, nodes time series correspond to the average BOLD signal in each
ROI. Subject-specific ICA nodes and group-averaged 22 ROIs time series were
normalised by z-scores.

Figure 1:  Illustration of the pipeline used in this study. Group-ICA
parcellation regressed into individuals is combined with FC to build a spatio-
temporal graph sequence representing time courses of brain activation.

### Model

In this paper, long-range, spatio-temporal dynamics of FC are modelled by
extending the ST-GCN model [29] (utilised in [5]) to incorporate information
over multiple spatial and temporal scales, by adapting the multi-scale,
spatio-temporal graph operator (MS-G3D) approach of [16]. Following the ST-GCN
model of FC used by [5], functional connectomes are modelled as graphs
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, with
$\mathcal{V}=\\{v_{1},v_{2}....v_{N}\\}$ representing the set of functional
regions (or nodes); and $\mathcal{E}$ representing spatial and temporal
connections (or edges). In this framework, the spatial graph
$\mathbf{A}\in\mathbb{R}^{\mathbb{N}\times\mathbb{N}}$:

$\mathbf{A}_{i,j}=\left\\{\begin{array}[]{c}corr(i,j)\text{ if }i\neq j\\\
0\text{ otherwise }\end{array}\right.$ (1)

represents static FC estimated over all time, with $corr(i,j)$ measuring the
magnitude of correlation over the full timeseries (1200 frames). Spatio-
temporal interactions are then modelled in the form of a feature tensor
$\mathbf{X}\in\mathbb{R}^{T\times N\times C}$:
$\mathcal{X}=\allowbreak\left\\{x_{n,t}\in\mathbb{R}^{C}|t,n\in\mathbb{Z},\break
1\leq t\leq T,1\leq n\leq N\right\\}$, where the channel dimension of the
input feature tensor represents the node time courses ($C=1$). The MS-G3D
framework [16] models spatio-temporal graph convolutions (ST-GC) as blocks
composed of two streams: a _factorised pathway_ which decouples space and
time, and a _G3D pathway_ which aggregates space-time information. In _the
factorised pathway_ , space-only convolutions are represented by:

$\mathbf{X}_{t}^{(l+1)}=\sigma\left(\mathbf{\widetilde{D}}^{-1/2}(\mathbf{A}+\mathbf{I})\mathbf{\widetilde{D}}^{-1/2}\mathbf{X}_{t}^{(l)}\Theta^{(l)}\right)$
(2)

with $\Theta^{(l)}$ the learnable weights at layer l and
$\mathbf{\widetilde{D}}$ the degree matrix of $(\mathbf{A}+\mathbf{I})$. Here,
$\sigma$ is a ReLU activation function, except for the last ST-GC block of the
model. Still in _the factorised pathway_ , temporal-only convolutions (TCN)
are implemented with a bottleneck design, including four convolutional
filters, all with kernels of size (3,1) along the temporal dimension but with
various dilation factors ($d\in[1,4]$), and residual connections. In parallel,
the _G3D pathway_ splits the full frame sequence into sub time-windows of
$\tau\in\left\\{3,5\right\\}$, and combines regional spacetime information
using a bespoke spatio-temporal graph convolution operation (G3D), applied on
every single sub time-window:

$\left[\mathbf{X}_{(\tau)}^{(l+1)}\right]_{t}=\sigma\left(\mathbf{\widetilde{D}}^{-1/2}_{(\tau)}\mathbf{\tilde{A}}_{(\tau)}\mathbf{\widetilde{D}}^{-1/2}_{(\tau)}\left[\mathbf{X}_{(\tau)}^{(l)}\right]_{t}\Theta^{(l)}\right)$
(3)

with

$\mathbf{\tilde{A}}_{(\tau)}=\left[\begin{array}[]{ccc}\mathbf{(A+I)}&\ldots&\mathbf{(A+I)}\\\
\vdots&\ddots&\vdots\\\
\mathbf{(A+I)}&\cdots&\mathbf{(A+I)}\end{array}\right]\in\mathbb{R}^{\mathbb{\tau
N}\times\mathbb{\tau N}}$ (4)

The G3D operation (Eq (3)), therefore extends the spatial convolution (Eq (2))
into a space-time operation by stacking the adjacency matrix $\mathbf{(A+I)}$
into $\mathbf{\tilde{A}}_{(\tau)}$ connecting nodes across $\tau$ frames, in
space and time: edges being extended to connect nodes to their neighbours
across the entire sub time-window; in this way constructing a spatio-temporal
subgraph $\mathcal{G}_{(t)}=(\mathcal{V}_{(t)},\mathcal{E}_{(t)})$.
$\mathbf{\widetilde{D}}_{(\tau)}$ and $\mathbf{X}_{(\tau)}$ are natural
extensions of the previously introduced notations.

In addition, a multi-scale (MS) operator can be embedded into any of the
previously described convolutions. This aims at aggregating more structural
information, by relating a node to its K-closest neighbours, through
substituting the convolution operation in (Eq 2) by
$\sum_{k=0}^{K}\mathbf{\widetilde{D}^{-1/2}_{(k)}}\mathbf{\tilde{A}_{(k)}}\mathbf{\widetilde{D}^{-1/2}_{(k)}}\mathbf{X_{t}^{(l)}}\Theta^{(l)}_{(k)}$
where:

$\mathbf{\tilde{A}}_{(k)}=\mathbf{I}+\mathbbm{1}\left(\mathbf{\tilde{A}}^{k}\geq
1\right)-\mathbbm{1}\left(\mathbf{\tilde{A}}^{k-1}\geq
1\right)|\mathbf{\tilde{A}}_{(1)}=\mathbf{A}+\mathbf{I};\mathbf{\tilde{A}}_{(0)}=\mathbf{I}$
(5)

In summary, where the ST-GCN model [29] adopted by [5] factorises spatio-
temporal convolutions into spatial-only and temporal-only blocks, useful for
modelling information flow along space and time (but limited in their capacity
to model space-time interactions), our method adds in parallel a space-time
convolution module connecting all nodes in a short temporal window.
Importantly, the TCN block adopted here from MS-G3D [16] also improves
capacity for more long-range temporal modelling by aggregating multi-scale
temporal convolutions, over several temporal scales.

### Architecture

Network implementation (illustrated in figure 2) follows the design principles
of MS-G3D [16] but adapts it to the case of FC. At each training iteration, a
time window of T consecutive frames is sampled at random from the entire
sequence and is passed as input to the model. In the present architecture
(figure 2), three ST-GC blocks were stacked with 96, 192, 384 outputs
channels, respectively. In each of these blocks, two _G3D pathways_
convolutions (Eq (3)) were implemented; these were passed spatio-temporal
graphs formed from time-windows of consecutive frames at two temporal lengths
($\tau\in\left\\{3,5\right\\}$). In parallel, the factorised pathway applies
spatial-only GCN (Eq (2)), followed by temporal-only convolutions TCN, on the
entire time-window sequence. Then, a multi-scale temporal module aggregates
the two pathways, increasing long-range temporal modelling by connecting
regional sub-windows. Finally, predictions are made either as classification
or regression, by simply replacing a sigmoid activation layer (for
classification) with a regression head for the fluid intelligence task. We
refer to this architecture as _Brain-MS-G3D_.

Figure 2: Architecture of our Brain-MS-G3D. A sequence of T frames is used as
input of the model. Each G3D pathway processes T sub time-windows of lenght
$\tau\in\left\\{3,5\right\\}$, while the factorised pathway process the entire
sequence T with spatial-only (GCN) and temporal-only (TCN) blocks. A final TCN
block aggregates all pathways.

### Experimental Setup

We first trained models, the ST-GCN [5] and our Brain-MS-G3D, on the 22 ROIs
group-average cortical parcellation used by [5]. Then both models were also
compared on the subject-specific node time series derived from dual-regressed
group-ICA at different dimensionalities (15 to 200 nodes). Following the
approach used in [5], we benchmark the proposed MS-G3D [16] model against the
ST-GCN [29] used in [5], by training models using 5-fold cross-validation. In
both cases models were trained on different length time-windows T, from the
entire 1200 time-points sequence. At each training iteration, a temporal
window was randomly extracted from the entire sequence to be processed by the
model. At inference time, predictions from $V=64$ voters were averaged. In all
cases, an Adam optimiser was used with binary cross-entropy loss (for
classification) and mean squared error (MSE) for regression; learning rate was
set to $1e^{-3}$ by default and with a weight decay of $1e^{-3}$. All
experiments were run on a single NVIDIA TITAN RTX 24GB GPU. The batch size was
maximised to fit into memory depending on the node resolution: for instance,
32 elements per batch for 200 nodes but 512 for 25 nodes.

## Results

Results obtained are summarised in figure 3 and Table 1. In Table 1, we
compared our method to the ST-GCN model [5], on the 22 ROIs and the ICA
subject-specific node timeseries at different dimensionalities for sex
classification and fluid intelligence prediction, and to results reported in a
recent method which models FC dynamics using spatio-temporal attention on
graphs [10]. In figure 3, performance on sex classification is evaluated
across different ICA (and group parcellation) dimensionalities, on models
trained and optimised for a maximum number of 2000 iterations. Results show
that our Brain-MS-G3D model, inspired by [16], outperforms the ST-GCN model
[29], used in [5], at every dimensionality, with an average increase in
performance of $+6.8\%$, and a maximum of $+11.1\%$ at dimensionality 100. In
Table 1, we report results for optimised trainings with more iterations, using
the best time-window lenght obtained from figure 3. Our Brain-MS-G3D notably
reached a 94.4% sex classification accuracy at dimensionality 200. We also
note in Table 1 an increase in classification for both models when moving from
the group-average parcellation (22 ROIs) to a similar ICA dimensionality (25
nodes). For sex classification, we experimented the impact of changing the
length of time-windows $T$ for three node dimensionalities (15,25,50): with
$T\in\left\\{50,75,100\right\\}$, representing respectively 4.2%, 6.3% and
8.3% of the entire time sequence. Range of performance across time windows is
reported through the vertical bars in figure 3, where plotted values
correspond to the highest accuracy for $T\in\left\\{50,75,100\right\\}$. In
the case of Brain-MS-G3D, all experiments suggest an increase in performance
with larger time-windows (best results with $T=100$), while the gap reduces
when spatial resolution increases. By comparison, the ST-GCN model
demonstrates a larger variance in performance, and the 100-time window network
always under-performing the others. As widely reported in the literature [11,
12, 19], fluid intelligence was a much more difficult task to predict. As
detailed in [20], performances for fluid intelligence performance are reported
as correlation values in Table 1. Our results quantitatively outperform
established methodology such as MAGE and HMM on fluid intelligence prediction
(see [19, 26]).

Figure 3: Comparison of sex classification results between ST-GCN and our Brain-MS-G3D on ICA subject-specific data at different resolution (15,25,50,100,200).  Methods | Data |  Sex (% acc) |  Fluid Intelligence (corr)
---|---|---|---
ST-GCN [5] | 22 Group ROIs | 81.8 [83.7] | 0.144
STAGIN-SERO [10] | 400 ROIs | [88.20] | N/A
STAGIN-GARO [10] | 400 ROIs | [87.01] | N/A
Brain-MS-G3D | 22 Group ROIs | 84.7 | 0.269
Brain-MS-G3D | ICA-15 | 81.5 | 0.286
Brain-MS-G3D | ICA-25 | 86.1 | 0.313
ST-GCN [5] | ICA-25 | 82.1 | N/A
Brain-MS-G3D | ICA-50 | 90.9 | 0.325
Brain-MS-G3D | ICA-100 | 93.9 | 0.317
Brain-MS-G3D | ICA-200 | 94.4 | 0.324
Table 1: Results of sex classification and fluid intelligence prediction on
the HCP dataset. Results presented were optimised for a maximum of 10k
iterations. In square brackets, result reported in original publications.

## Discussion

While our approach introduced a spatio-temporal modelling of functional
connectivity that has proved to increase phenotype predictions, additional
experiments revealed that the multi-scale operator (Eq (5)), used in [16], has
only a small impact on Brain-MS-G3D performances. In the case of brain
connectomes, the spatial graph is by construction a weighted-undirected graph
that connects all nodes together, while in skeleton analysis the spatial graph
is much sparser. Thresholding correlation values in the adjacency matrix could
increase the impact of spatial scaling (Eq (5)). Additionally, in figure 3, we
illustrated that increasing the spatial resolution of the ICA tends to improve
performances. However, the length of time sequence has also been pointed as
crucial to understand the temporal hierarchy of FC [27]. We had to limit the
study to a relatively small number of time frames per sequence (maximum of
128), mostly due to hardware limitations, but experiments on larger time
windows with smaller ICA resolution showed dramatic increases in prediction
accuracy, which should be further investigated. Finally, figure 4 illustrates
the impact of parcellation for modelling subject-variability in functional
connectivity: atypical subject-specific activation might be averaged out by
parcellation.

In conclusion, this study improved behavioural and demographic predictions on
the HCP dataset by introducing a state-of-the-art method for skeleton-based
action recognition to explicitly modelling spatio-temporal dynamics of FC, and
accounting for inter-subject variability in functional organisation through
the use of dual-regressed ICA maps. Future work will explore the potential for
improving performance further by increasing the spatial resolution of the
model.

## Acknowledgements

Data were provided by the Human Connectome Project, WU-Minn Consortium
(Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657)
funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for
Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at
Washington University [25].

(a) 22 ROIs group-average parcellation derived from [6]

(b) Group-average activation map

(c) Subject-specific activation map

Figure 4: Pink arrows highlight Area 55b, which shows topographic variability
between subject. In (a), the functional connectivity of this region is
averaged across the larger ROI (light green). (b) group average functional
connectivity of the HCP language task contrast ”Story vs. Baseline”, and (c)
shows atypical functional connectivity of the same language task contrast in a
single subject.

## References

  * [1] Beckmann, C., Smith, S.: Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Transactions on Medical Imaging 23(2), 137–152 (2004)
  * [2] Bijsterbosch, J.D., Woolrich, M.W., Glasser, M.F., Robinson, E.C., Beckmann, C.F., et al.: The relationship between spatial configuration and functional connectivity of brain regions. Elife (2018)
  * [3] Dsouza, N.S., Nebel, M.B., Crocetti, D., Robinson, J., Mostofsky, S., et al.: M-gcn: A multimodal graph convolutional network to integrate functional and structural connectomics data to predict multidimensional phenotypic characterizations (2021)
  * [4] Finn, E.S., Shen, X., Scheinost, D., Rosenberg, M.D., Huang, J., et al.: Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity. Nature neuroscience 18(11), 1664–1671 (2015)
  * [5] Gadgil, S., Zhao, Q., Pfefferbaum, A., Sullivan, E.V., Adeli, E., et al.: Spatio-temporal graph convolution for resting-state fmri analysis (2021)
  * [6] Glasser, M.F., Coalson, T.S., Robinson, E.C., Hacker, C.D., Harwell, J., et al.: A multi-modal parcellation of human cerebral cortex. Nature (7615), 171–178 (2016)
  * [7] Glasser, M.F., Sotiropoulos, S.N., Wilson, J.A., Coalson, T.S., Fischl, B., et al.: The minimal preprocessing pipelines for the human connectome project. NeuroImage 80, 105–124 (2013)
  * [8] Harrison, S.J., Bijsterbosch, J.D., Segerdahl, A.R., Fitzgibbon, S.P., Farahibozorg, S.R., et al.: Modelling subject variability in the spatial and temporal characteristics of functional modes. NeuroImage 222, 117226 (2020)
  * [9] Huang, Z.A., Zhu, Z., Yau, C.H., Tan, K.C.: Identifying Autism Spectrum Disorder From Resting-State fMRI Using Deep Belief Network. IEEE Transactions on Neural Networks and Learning Systems (2020)
  * [10] Kim, B.H., Ye, J.C., Kim, J.J.: Learning dynamic graph representation of brain connectome with spatio-temporal attention (2021)
  * [11] Kong, R., Li, J., Orban, C., Sabuncu, M.R., Liu, H., et al.: Spatial topography of individual-specific cortical networks predicts human cognition, personality, and emotion. Cerebral cortex 29(6), 2533–2551 (2019)
  * [12] Kong, R., Yang, Q., Gordon, E., Xue, A., Yan, X., et al.: Individual-specific areal-level parcellations improve functional connectivity prediction of behavior. bioRxiv (2021)
  * [13] Ktena, S.I., Parisot, S., Ferrante, E., Rajchl, M., Lee, M., et al.: Metric learning with spectral graph convolutions on brain connectivity networks. NeuroImage 169, 431–442 (2018)
  * [14] Li, X., Zhou, Y., Gao, S., Dvornek, N., Zhang, M., et al.: Braingnn: Interpretable brain graph neural network for fmri analysis. bioRxiv (2020)
  * [15] Liégeois, R., Li, J., Kong, R., Orban, C., Van De Ville, D., et al.: Resting brain dynamics at different timescales capture distinct aspects of human behavior. Nature Communications 10(1), 1–9 (2019)
  * [16] Liu, Z., Zhang, H., Chen, Z., Wang, Z., Ouyang, W.: Disentangling and unifying graph convolutions for skeleton-based action recognition (2020)
  * [17] Marquand, A.F., Kia, S.M., Zabihi, M., Wolfers, T., Buitelaar, J.K., et al.: Conceptualizing mental disorders as deviations from normative functioning. Molecular Psychiatry (10), 1415–1424 (2019)
  * [18] Marquand, A.F., Rezek, I., Buitelaar, J., Beckmann, C.F.: Understanding Heterogeneity in Clinical Cohorts Using Normative Models: Beyond Case-Control Studies. Biological Psychiatry 80(7), 552–561 (2016)
  * [19] Pervaiz, U., Vidaurre, D., Gohil, C., Smith, S.M., Woolrich, M.W.: Multi-dynamic modelling reveals strongly time-varying resting fmri correlations. bioRxiv (2021)
  * [20] Pervaiz, U., Vidaurre, D., Woolrich, M.W., Smith, S.M.: Optimising network modelling methods for fMRI. NeuroImage 211, 116604 (2020)
  * [21] Robinson, E.C., Garcia, K., Glasser, M.F., Chen, Z., Coalson, T.S., et al.: Multimodal surface matching with higher-order smoothness constraints. NeuroImage 167, 453–465 (2018)
  * [22] Robinson, E.C., Jbabdi, S., Glasser, M.F., Andersson, J., Burgess, G.C., et al.: Msm: A new flexible framework for multimodal surface matching. NeuroImage 100, 414–426 (2014)
  * [23] Smith, S.M., Nichols, T.E., Vidaurre, D., Winkler, A.M., Behrens, T.E., et al.: A positive-negative mode of population covariation links brain connectivity, demographics and behavior (2015)
  * [24] Smith, S.M., Vidaurre, D., Beckmann, C.F., Glasser, M.F., Jenkinson, M., et al.: Functional connectomics from resting-state fmri. Trends in Cognitive Sciences 17(12), 666–682 (2013)
  * [25] Van Essen, D.C., Smith, S.M., Barch, D.M., Behrens, T.E., Yacoub, E., et al.: The wu-minn human connectome project: An overview. NeuroImage 80, 62–79 (2013)
  * [26] Vidaurre, D., Abeysuriya, R., Becker, R., Quinn, A.J., Alfaro-Almagro, F., et al.: Discovering dynamic brain networks from big data in rest and task. NeuroImage 180, 646–656 (2018)
  * [27] Vidaurre, D., Smith, S.M., Woolrich, M.W.: Brain network dynamics are hierarchically organized in time. Proceedings of the National Academy of Sciences 114(48), 12827–12832 (2017)
  * [28] Wolfers, T., Rokicki, J., Alnæs, D., Berthet, P., Agartz, I., et al.: Replicating extensive brain structural heterogeneity in individuals with schizophrenia and bipolar disorder. Human Brain Mapping 42(8), 2546–2555 (2021)
  * [29] Yan, S., Xiong, Y., Lin, D.: Spatial temporal graph convolutional networks for skeleton-based action recognition (2018)


# Self-Similar solution of rotating eruptive outflows on its equatorial plane

Takatoshi Ko  Research Center for the Early Universe (RESCEU), Graduate School
of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033,
Japan  Department of Astronomy, Graduate School of Science, The University of
Tokyo, Tokyo, Japan  Kotaro Fujisawa  Department of Physics, Graduate School
of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033,
Japan  Toshikazu Shigeyama  Research Center for the Early Universe (RESCEU),
Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan  Department of Astronomy, Graduate School of Science,
The University of Tokyo, Tokyo, Japan

(Received Sep 01, 2022; Revised Sep 27, 2022; Accepted Nov 05, 2022)

###### Abstract

We construct axisymmetric self-similar solutions of transonic outflows
emanating from a point source including the effect of the rotation. The
solutions are constructed exclusively on the equatorial plane. The features of
solutions are determined by three parameters; the adiabatic index $\gamma$,
the dimensionless coordinate of the transonic point, and the dimensionless
azimuthal velocity at the transonic point. We classify the solutions into five
groups according to the asymptotic behaviors. We find that the behaviors of
the self-similar solutions change at $\gamma=11/9$. In addition, some
solutions show double-power-law density profiles, which are usually seen in
ejecta from a binary merger or nova-like explosion. Thus, our self-similar
solutions can be applied not only to the outflow blowing from the central
spinning objects, but also to the ejecta erupted from the binary merger or
nova-like explosion.

XXX

††journal: ApJ

##  1 Introduction

All objects in the universe rotate and have angular momentum. Among them,
there are many objects formed via binary stellar mergers, for which the
effects of rotation cannot be ignored. Some of such objects emit outflows
containing high angular momentum, and the outflows are very important in the
evolution of the system.

In addition to many objects that are spinning around and emitting rotational
outflow at high velocities on their own, there are also objects that emit
rotating ejecta by binary merger. Both observational researches and numerical
simulations showed that binary mergers erupt a lot of such rotation-rich
ejecta (e.g., Tanikawa et al., 2015; Schneider et al., 2019; Shibata &
Hotokezaka, 2019), which bring essential information on the original objects.
Schneider et al. (2019) reported three-dimensional simulations of the
coalescence of two massive stars and showed there were some ejecta from the
merger. Furthermore, binary white dwarf (WD) merger is known to cause a wide
variety of phenomena such as type Ia supernovae (e.g., Webbink, 1984; Iben &
Tutukov, 1984; Bildsten et al., 2007; Fink et al., 2007; Guillochon et al.,
2010; Dan et al., 2011; Pakmor et al., 2013), the formation of a massive WD
(e.g., Schwab et al., 2016; Gvaramadze et al., 2019; Oskinova et al., 2020),
and collapse into a neutron star through accretion-induced collapse (AIC)
(e.g., Saio & Nomoto, 1985; Taam & van den Heuvel, 1986; Piro & Kulkarni,
2013; Moriya, 2016). The ejecta released from the binary WD merger is one of
the important parameters that determines the physics of these phenomena. Some
numerical results of binary WD mergers involving the ejecta have been reported
(e.g., Tanikawa et al., 2015; Sato et al., 2015, 2016). The ejecta from binary
neutron star (NS) merger are very important because the ejecta become the
electromagnetic counterpart known as kilonova (e.g., Li & Paczyński, 1998;
Metzger et al., 2010; Abbott et al., 2017) and a plausible origin of elements
heavier than iron (Tsujimoto & Shigeyama, 2014; Lattimer & Schramm, 1974).
There are a lot of reports about the numerical simulations of binary NS merger
and the properties of the ejecta (for review, e.g., Shibata & Hotokezaka,
2019). Some simulations showed the double-power-law density profile of the
ejecta (e.g., Hotokezaka et al., 2018).

Many stationary wind models have been constructed to understand the nature of
these winds or ejecta blowing from the central objects without rotation (e.g.,
Parker, 1965; Kato & Hachisu, 1994) and with rotation (e.g., Weber & Davis,
1967; Lyutikov, 2018; Kashiyama et al., 2019).

Parker (1965) gave a review of the one-dimensional steady wind model that
includes the gravity of the central object. This model explains the solar wind
in the vicinity of the sun very successfully. Kato & Hachisu (1994)
constructed a steady wind model including the effects of radiation transfer
and explicitly showed that the wind is accelerated by the radiation pressure.
Weber & Davis (1967) constructed an axisymmetric steady wind model including
the central gravity, rotation, and magnetic field. The model is formulated as
one-dimensional problems by focusing on the equatorial plane. Kashiyama et al.
(2019), by combining the models of Weber & Davis (1967) and Kato & Hachisu
(1994), constructed an axisymmetric steady wind model including the central
gravity, rotation, radiation transfer, and magnetic fields on the equatorial
plane. This reported the nature of the wind from a fast rotating massive WD
produced via binary merger, which is known to be blowing fast wind at a speed
of up to $\sim 16,000~{}\mathrm{km\,s^{-1}}$ from optical observations
(Gvaramadze et al., 2019). The model attempts to reproduce the observed
features of this object by the behavior of a rotating magnetic wind on the
equatorial plane emitted from the central WD.

Recent short cadence observations will detect emissions from these events in
the near future at a very early stage when the outflow has not yet reached a
steady state. To prepare for this situation, we aim to construct a solution
which can describe the early stage of the evolution of the outflow including
the effect of rotation in this study.

In order to describe the evolution up to the steady state, a useful way is to
assume the self-similarity of the evolving blast waves (so-called self-similar
solution). Taylor (1950), Sedov (1959) and Parker (1961) calculated self-
similar outflow models to study a sudden ejection of matter such as a coronal
mass ejection. Their models assumed the existence of a shock formed by a
collision of the ejected material with the ambient medium. Thus, their self-
similar models are strongly affected by the surroundings, and it is difficult
to understand the pure nature of the ejection. In addition, the self-similar
solutions assuming the existence of one shock are effective only after the
reverse shock disappeared. In order to describe the evolution of the outflow
until a shock is generated, a transonic outflow model is an important tool.
Cheng (1977) and Fukue (1984) reported such one-dimensional self-similar
transonic outflow models. Cheng (1977) constructed a self-similar transonic
outflow model including the central gravity by assuming a polytropic equation
of state. That is, the flow is assumed to be isentropic. Fukue (1984) took an
approach to the self-similar transonic outflow different from Cheng (1977).
Fukue (1984) did not adopt the polytropic equation of state but assumed the
adiabatic evolution. None of these self-similar models took into account
rotation.

Kashiyama et al. (2019) reported that a rapid rotation plays an important role
in the outflow properties. Thus, in order to apply to objects such as the
fast-spinning WD, the transonic outflow model should be constructed by taking
into account of rotation. Therefore, in this work we construct self-similar
transonic outflow solutions including the central gravity and the high
azimuthal velocity on the equatorial plane. It should be noted that the
solutions constructed in this work include the solution of Cheng (1977)
because of the same setup except for rotation. Thus, we report the solutions
with a topological examination including Cheng (1977) solutions.

This paper is organized as follows. In Sect. 2, we introduce the self-similar
setup and the integration of our work. In Sect. 3.1, we present the results of
our calculations and discuss them. In Sect. 4, we conclude the property of our
self-similar calculations.

##  2 Methods

###  2.1 Model

We consider rotating stellar outflow emanating from a central point source
that exerts the gravity. We assume that the system is axisymmetric around the
rotational axis and consider the outflow exclusively on the equatorial plane.
Under these assumptions, the continuity equation and the Euler equation become

$\frac{\partial\rho}{\partial t}=-\frac{1}{r^{2}}\frac{\partial}{\partial
r}(r^{2}\rho v_{r}),$ (1) $\frac{\partial v_{r}}{\partial
t}=\frac{v_{\phi}^{2}}{r}-v_{r}\frac{\partial v_{r}}{\partial
r}-\frac{1}{\rho}\frac{\partial p}{\partial
r}-\frac{G{\color[rgb]{0,0,0}{\color[rgb]{0,0,0}M_{*}}}}{r^{2}},$ (2)
$\frac{\partial v_{\phi}}{\partial
t}=-\frac{1}{r}v_{r}\frac{\partial}{\partial r}(rv_{\phi}),$ (3)

where, $r$ is the radius, $t$ is the time, $G$ is the gravitational constant,
${\color[rgb]{0,0,0}M_{*}}$ is the mass of the central object and
$\rho,v_{r},v_{\phi},p$ are the density, the radial velocity, the azimuthal
velocity, and the pressure, respectively. They are functions of $r$ and $t$.
Here we assume that the flow is isentropic because the cooling timescale is
much longer than the dynamical timescale in the eruptive outflow context:

$p=K\rho^{\gamma},$ (4)

where $K$ is a constant, and $\gamma$ is the adiabatic index.

###  2.2 Self-Similar Formalism

We introduce the similarity variables as follows:

$x=\frac{r^{3}}{At^{2}},$ (5)
$\rho=\left(\frac{r^{2}}{Kt^{2}}\right)^{\frac{1}{\gamma-1}}\Omega(x),$ (6)
$v_{r}=\frac{r}{t}V_{r}(x),$ (7) $v_{\phi}=\frac{r}{t}V_{\phi}(x),$ (8)

where $A=G{\color[rgb]{0,0,0}M_{*}}$.

Next, we rewrite equations (1)-(3) in terms of similarity variables:

$W\dotproduct df=B,$ (9)

where $W$, $df$ and $B$ can be written as follows:

$W=\matrixquantity(x\left(\gamma-1\right)\Omega&0&x(\gamma-1)v\\\
x(\gamma-1)v&0&x\gamma(\gamma-1)\Omega^{\gamma-2}\\\ 0&x(\gamma-1)v&0),$ (10)
$df=\matrixquantity(\frac{dV_{r}}{dx}\\\ \frac{dV_{\phi}}{dx}\\\
\frac{d\Omega}{dx}),$ (11) $B=\matrixquantity(\frac{\left(-3\gamma
V_{r}+V_{r}+2\right)\Omega{}}{3}\\\
\frac{(1-\gamma)\left(V_{r}^{2}-V_{\phi}^{2}+V_{r}+1/x\right)}{3}-\frac{2\gamma\Omega^{\gamma-1}}{3}\\\
-\frac{(\gamma-1)(2V_{r}-1)V_{\phi}}{3}).$ (12)

Here, $v(x)$ has been introduced as

$v(x)=V_{r}\left(x\right)-2/3.$ (13)

To solve these equations numerically, we rewrite the equation as follows:

$df=W^{-1}\dotproduct B.$ (14)

To obtain a transonic solution, we require that $df$ should not diverge even
at a critical point where the determinant of $W$ becomes $0$, i.e.,

$\mathrm{det}\ W=(\gamma-1)^{3}x^{3}\Omega^{3}v\left(v^{2}-C^{2}\right)=0,$
(15)

where $C=\sqrt{\gamma\Omega^{\gamma-1}}$ is the dimensionless sound speed.
Otherwise, equation (14) becomes singular at the critical point where $v=C$
and the flow truncates at this point. It should be noted from this equation
that the density distribution can also be singular at a point where $v=0$,
which describes a contact surface.

###  2.3 Singularity Analysis and Boundary Conditions

In order to obtain a transonic solution from equation (9) without numerical
divergence at the critical point, both of the denominator and numerator of
equation (14) must be equal to $0$, i.e. the following relation should be
satisfied at the critical point:

$18C^{2}x+\frac{15\gamma-21}{\gamma-1}Cx+2x-9+9V_{\phi}^{2}x=0.$ (16)

Here, we write $V_{\phi}{\color[rgb]{0,0,0}(x=x_{\rm crit})}=\alpha
V_{K}{\color[rgb]{0,0,0}(x=x_{\rm crit})}$, where $V_{K}$ is the dimensionless
azimuthal velocity in Kepler rotation and $\alpha$ is defined at the critical
point. Since the azimuthal velocity in Kepler rotation is
$v_{K}=\sqrt{G{\color[rgb]{0,0,0}M_{*}}/r}$, $V_{K}$ can be expressed as
follows:

$V_{K}=\frac{t}{r}\sqrt{\frac{G{\color[rgb]{0,0,0}M_{*}}}{r}}=\frac{1}{\sqrt{x}}.$
(17)

Substituting this to equation (16), we obtain the following relation at the
critical point $x=x_{\rm{crit}}$:

$18C^{2}x+\frac{15\gamma-21}{\gamma-1}Cx+2x-9(1-\alpha^{2})=0,$ (18)

or

$x_{\rm
crit}(\alpha,C)=\frac{9(1-\alpha^{2})}{18C^{2}+\frac{15\gamma-21}{\gamma-1}C+2}.$
(19)

Figure 1 shows this relation for several $\gamma$. It can be easily seen from
the equation and the figure that the position of the critical point decreases
when $|\alpha|$ (the rotation) is large for a constant $C$. Here, let $f(C)$
be the denominator of equation (19), the discriminant $D$ of the quadratic
equation $f(C)=0$ :

$D=\frac{9\left(9\gamma-11\right)\left(\gamma-3\right)}{(\gamma-1)^{2}},$ (20)

indicates that $f(C)$ and equation (19) change the behavior depending on
$\gamma$ greater or less than $11/9$. When $\gamma>11/9$, $D$ becomes negative
and $f(C)$ does not vanish for any real $C$. Thus, $x_{\rm{crit}}$ takes the
maximum value $x_{\rm{crit}}=9(1-\alpha^{2})/2$ at $C=0$, decreases
monotonically with increasing $C$, and approaches to $0$ as $C$ approaches
infinity (left panel of Fig. 1). When $\gamma\leq 11/9$, $x_{\rm{crit}}$ can
take all positive values because $f(C)=0$ has a real solution of $C$ and then
$x_{\rm{crit}}$ can be infinite (center and right panels in Fig. 1).

From equations (15) and (18), all the physical quantities at the singular
point are obtained by giving $x_{\rm{crit}}$ and $\alpha$. Therefore, we set
the values of $x_{\rm{crit}}$ and $\alpha$ and initial values at the critical
point and integrate equations (9) inwardly and outwardly avoiding numerical
divergence.

###  2.4 Integration

We use the Runge–Kutta method to numerically integrate equation (14) from the
critical point, using the `scipy.integrate` module (Virtanen et al. 2020). As
mentioned in Sect.2.3, the integration starts from points near the critical
point for some parameter sets of $(\gamma,x_{\mathrm{crit}},\alpha)$.

##  3 Results and Discussion

In this section, we present self-similar solutions for several parameter sets
and the analysis of the self-similar solutions. In addition, we present the
solutions with dimensions and discuss the features of the solutions.

###  3.1 Dimensionless solutions

Figure 2 shows the distribution of the dimensionless Mach number defined as

$M(x)=\frac{v(x)}{C(x)}=\frac{V_{r}(x)-2/3}{\sqrt{\gamma\Omega^{\gamma-1}(x)}},$
(21)

for each parameter set $(\gamma,x_{\mathrm{crit}},\alpha)$. Where $M(x)=1$
(the sonic point) at a place other than $x=x_{\rm crit}$, the derivatives in
equation (14) diverge and the flow truncates as shown in Figure 2. This means
that the flow is in contact with the vacuum at the sonic singular point with
some exceptions (see section 3.1.6).

####  3.1.1 Classification

We found that solutions with rotation ($\alpha\neq 0$) can be classified into
five types based on the number of sonic points: solutions with $\gamma\leq
11/9$ being classified into four types (I, II, III, and IV) and solutions with
$\gamma>11/9$ into one (V). The flow with $\gamma\leq 11/9$ has a sonic point
or a contact surface inside the transonic point depending on whether
$|\alpha|$ is larger than a threshold value or not. For $\gamma>11/9$ (type
V), on the other hand, an inner sonic point always appears unless the flow
does not rotate at all (i.e., $\alpha\neq 0$).

####  3.1.2 Solutions with $\gamma\leq 11/9$

For $\gamma\leq 11/9$, the solutions are classified based on whether a
singular point exists inside ($x<x_{\rm crit}$) or outside ($x>x_{\rm crit}$)
of the transonic point. Type I solution has sonic points on both sides of the
transonic point (see the line with a parameter set of
$(\gamma,x_{\mathrm{crit}},\alpha)=(7/6,0.05,0.2)$ in Fig. 2, for example) and
the flow truncates on both sides of the sonic points. Type II solution
truncates at the sonic point outside the transonic point and also truncates at
the contact surface inside of the transonic point where $V_{r}=2/3$ ($v=0$)
instead of $M=1$ (e.g., $(\gamma,x_{\mathrm{crit}},\alpha)=(7/6,0.05,0.1)$).
On the other hand, type III solution has an inner sonic point but no outer
sonic point and extends to infinity (e.g.,
$(\gamma,x_{\mathrm{crit}},\alpha)=(7/6,0.0005,0.2)$). Type IV solution (e.g.,
$(\gamma,x_{\mathrm{crit}},\alpha)=(7/6,0.0005,0.1)$) has a contact surface
inside the transonic point and has no sonic point on either side of the
transonic point. At the inner sonic point in type I and III solutions, $V_{r}$
and $\Omega$ diverge to infinity as shown in Figure 3. Type II and type IV
solutions appear when the value of $|\alpha|$ is small or the rotation at the
transonic point is slow.

Solutions of types III and IV extend to infinity and approach homologous
expansion. This is indicated from Figure 3 in which $V_{r}$ of these solutions
approach unity at large $x$. We recognize such flows approaching homologous
expansion ($V_{r}=1$) as outflow type solutions. We can obtain the asymptotic
behaviors of the density and the azimuthal velocity listed in table 1 by
substituting $V_{r}=1$ ($v=1/3$) into the equation (9). Equations of the
conservation of the mass and the angular momentum yield

$\frac{x}{3}\frac{d\Omega}{dx}=-\Omega,$ (22)

and

$\frac{x}{3}\frac{dV_{\phi}}{dx}=-\frac{V_{\phi}}{3},$ (23)

respectively. Then $\Omega\propto x^{-3}$ and $V_{\phi}\propto x^{-1}$ satisfy
these equations. In the inner region where $x_{\rm{min}}\ll x\ll
x_{\rm{crit}}$, $V_{r}\gg 1$ from Figure 3. Then we can show $V_{\phi}\propto
x^{-2/3}$ because

$xV_{r}\frac{dV_{\phi}}{dx}=\frac{2V_{\phi}}{3},$ (24)

from equation (14).

All solutions with $\gamma\leq 11/9$ and a finite $\alpha$ have rotational
velocities greater than the Kepler rotation in the innermost region. Note that
the boundary surfaces between the four types (I, I, III and IV) in the
parameter space of $(\gamma,x_{\mathrm{crit}},\alpha)$ can only be determined
numerically.

####  3.1.3 Solutions with $\gamma>11/9$

Type V solution with $\gamma>11/9$ always has an inner sonic point but no
outer sonic point. As shown in Table 1, the asymptotic power-law behaviors of
physical quantities have the same exponents as types III and IV, but there is
a difference from type III solutions. Type V without rotation extends to the
center, while type III never reaches the center.

The fact that the behavior of the solution changes across $\gamma=11/9$ is
consistent with the change in behavior across $n=1/(\gamma-1)=9/2$ in equation
(A11) in Cheng (1977).

####  3.1.4 Outer critical points

Figure 4 exemplifies transitions from type II to type IV solutions with
decreasing $x_{\rm crit}$ by showing the relation between the critical point
$x_{\rm{crit}}$ and the ratio of the coordinates of the outer sonic point
$x_{\rm{max}}$ to $x_{\rm{crit}}$. The divergence of the values of
$x_{\rm{max}}/x_{\rm{crit}}$ to infinity in this plot indicates that the type
of the solution shifts from type II to type IV. Another noticeable feature
observed from Figures 2 and 4 is that the ratio $x_{\rm{max}}/x_{\rm crit}$
decreases with increasing $x_{\rm{crit}}$ following a scaling law:
$x_{\rm{max}}/x_{\rm{crit}}\propto x_{\rm{crit}}^{18(1-\gamma)}$ as indicated
in Figure 4.

####  3.1.5 Inner critical points

Figures 2 and 5 show that the faster the rotation (the larger $|\alpha|$), the
closer the distance between the transonic point $x_{\rm{crit}}$ and the inner
sonic point for the flow with $\gamma>11/9$. Thus the ratio
$x_{\rm{min}}/x_{\rm{crit}}$ increases with increasing $|\alpha|$ and type V
solution with a finite value of $\alpha$ has an inner sonic point (at
$x=x_{\rm{min}}$) and the flow truncates there. On the other hand, the flow
with $\gamma\leq 11/9$ has a different trend because the flow with a small
$|\alpha|$ can truncate at the contact surface where $v$ becomes 0 rather than
at a sonic point. When the inner boundary is a contact surface (i.e., $v=0$),
the ratio of $x_{\rm{min}}/x_{\rm{crit}}$ decreases with increasing $\alpha$
(dashed line in Fig. 5). In contrast, the flow with a larger $|\alpha|$
truncates at the point where $v=C$, and the ratio $x_{\rm{min}}/x_{\rm{crit}}$
increases with increasing $\alpha$. This trend can be approximated as
$x_{\rm{min}}/x_{\rm{crit}}\propto\exp(b|\alpha|)$, where $b$ is a
proportional coefficient. For solutions with an inner sonic point, the values
of $x_{\rm{min}}/x_{\rm{crit}}$ monotonically decrease with increasing
$x_{\rm{crit}}(\alpha)$. Since $x_{\rm{crit}}(\alpha)$ is a monotonically
decreasing function of $\alpha$, the subsonic region shrinks with increasing
$|\alpha|$. This feature can be seen in the steady outflow model as well
(e.g., see figure 3 in Lyutikov, 2018). In addition, the value of
$x_{\rm{min}}/x_{\rm{crit}}$ increases with decreasing $x_{\rm{crit}}$, which
means that the flow truncates immediately close to the transonic point for
small $x_{\rm{crit}}$. These two facts reflect the relative importance of the
gravity and the rotation. A small $x_{\rm{crit}}$ indicates a significant
centrifugal force relative to the gravity since
$x={r^{3}}\big{/}{G{\color[rgb]{0,0,0}M_{*}}t^{2}}$.

####  3.1.6 The effect of rotation

Figure 3 compares distributions of the dependent variables for different types
of solutions. In each type of solutions, the flow truncates at a point inside
the transonic point with $M=0$ or $M=1$, and the azimuthal velocity there
sometimes exceeds that of the Kepler rotation where the flows truncate.
Furthermore, since $V_{\phi}$ becomes much smaller than $V_{r}$ and the Kepler
rotation outside the transonic point ($x>x_{\rm{crit}}$), there is little
effect of rotation in any type of solution.

###  3.2 Self-Similar Solutions with Dimension

In this section, we attach the solutions obtained in Sect. 3.1 with dimensions
by using equations (5) - (8). For simplicity, we show solutions with
dimensions setting $A=1$ and $K=1$. We have classified solutions into five
types. Both of type I and type II solutions terminate at a certain point
outside the critical point. Type III and type IV solutions have density
distributions increasing with radius. This feature would not be feasible in
reality. Furthermore, some of these I-V types of solutions rotate faster than
the Kepler rotation inside the transonic point. Such a high azimuthal velocity
may not be feasible without taking into account the effects of strong magnetic
fields (e.g., Kashiyama et al., 2019).

Thus, we focus on type V solutions in which the azimuthal velocity never
exceeds that of the Kepler rotation and show their features in the rest of
this section. In the limit of large $r$ or small $t$, the velocities and the
density of type V solution have the following asymptotic temporal and radial
dependence;

$v_{r}(r,t)\propto\frac{r}{t},$ (25) $v_{\phi}(r,t)\propto\frac{t}{r^{2}},$
(26) $\rho(r,t)\propto
t^{\frac{6\gamma-8}{\gamma-1}}r^{-\frac{9\gamma-11}{\gamma-1}}\propto
t^{-3}\left(\frac{r}{t}\right)^{-\frac{9\gamma-11}{\gamma-1}}.$ (27)

Equations (25) and (27) show that matter in the outer region freely expands.

In the inner region ($x_{\rm min}\ll x\ll x_{\rm crit}$), the velocities and
the density have the following asymptotic temporal and radial dependences;

$v_{r}(r,t)\propto
t^{\frac{3\gamma-5}{3\gamma-3}}r^{-\frac{2\gamma-3}{\gamma-1}},$ (28)
$v_{\phi}(r,t)\propto t^{1/3}r^{-1},$ (29) $\rho(r,t)\propto
t^{0}r^{-\frac{1}{\gamma-1}}.$ (30)

Equations (28) and (30) indicate that the solution in the inner region is in a
nearly hydrostatic equilibrium state as long as $\gamma<5/3$ and $v_{\phi}\ll
v_{K}$, because the radial velocity becomes much smaller than the sound speed
there. Thus,

$\gamma
K\rho^{\gamma-2}\frac{d\rho}{dr}+\frac{G{\color[rgb]{0,0,0}M_{*}}}{r^{2}}=0,$
(31)

holds there. However, in Type V solution, it is numerically found that if the
$|\alpha|$ is smaller than a certain threshold, the azimuthal velocity exceeds
the Kepler velocity in the inner region, and this $|\alpha|$ threshold becomes
smaller for larger $\gamma$. For example, $\gamma=1.5$ yields the threshold of
$|\alpha|\sim 0.1$ and $|\alpha|\sim 0.007$ for $\gamma=1.6$. Thus, type V
solutions with lower $|\alpha|$ rotate faster than the Kepler rotation at the
innermost points while with higher $|\alpha|$ do not exhibit such peculiar
features. Thus, type V solution has the inner region that can be regarded as a
rotating steady state as long as $v_{\phi}\ll v_{K}$ holds. The steady state
can be attained with large $t$’s as well as with small $x$’s because a large
$t$ also indicates a small $x$. Thus the inner region of type V solution
attains a rotating steady state and this region expands over time to a larger
radius. In this type V solutions, the total angular momentum increases over
time and the angular momentum continuously enters from the inner boundary.

The above calculations show that the flow with $\gamma>11/9$ converges to
homologous expansion at large $x$, and its density profile has a double-power-
law distribution as follows:

$\rho(r,t)\propto r^{-\frac{1}{\gamma-1}}\ \rm{(inner)},$ (32)
$\rho(r,t)\propto
t^{\frac{6\gamma-8}{\gamma-1}}r^{-\frac{9\gamma-11}{\gamma-1}}\ \rm{(outer)}.$
(33)

Since the double-power-law of the density profile is realized when
$\gamma>11/9$, the density profile is always double-power-law when the
pressure is dominated by radiation ($\gamma\sim 4/3$). Since such double-
power-law of the density profile is commonly observed in ejecta erupted from
binary mergers (e.g., Shibata & Hotokezaka, 2019), the model constructed in
this work may explain the ejecta from the binary-merger system when the ejecta
mass can be ignored compared with the central object mass. Note that, the
density profile becomes a ”single-power law” when $\gamma=4/3$. In this work,
we assume that the ejecta mass can be neglected, and this assumption holds
when

${\color[rgb]{0,0,0}M_{*}\gg 4\pi\int r^{2}\rho dr}$ (34)

holds if we assume the spherical symmetry. Using our solutions, we can
calculate the right hand side:

${\color[rgb]{0,0,0}4\pi\int r^{2}\rho
dr=4\pi\int\left(At^{2}x\right)^{2/3}\left(\frac{r^{2}}{Kt^{2}}\right)^{\frac{1}{\gamma-1}}\Omega(x)dr.}$
(35)

With fixed time $t$, this equation can be expressed as

${\color[rgb]{0,0,0}4\pi\int r^{2}\rho
dr=A^{\frac{3\gamma-1}{3\gamma-3}}K^{-\frac{1}{\gamma-1}}t^{\frac{6\gamma-8}{3\gamma-3}}\int
x^{\frac{2}{3\gamma-3}}\Omega(x)dx.}$ (36)

Thus, for the solutions with $\gamma>4/3$, this assumption holds at small $t$,
and the solutions with $\gamma<4/3$, this holds at large $t$. However, for the
solutions with $\gamma<4/3$, the ejecta mass is monotonically decreasing and
it is not an outflow solution. In type V solution with
$(\gamma,x_{\mathrm{crit}},\alpha)=(5/3,0.05,0.1)$, for example, this relation
can be calculated as

${\color[rgb]{0,0,0}4\pi\int r^{2}\rho dr\sim 1.4A^{2}K^{-3/2}t.}$ (37)

In the context of massive white dwarf observed by Gvaramadze et al. (2019),
$M_{*}\sim 1.5M_{\odot}$ and the mass loss rate $\dot{M}\sim 3.5\times
10^{-6}M_{\odot}$/yr are estimated by Gvaramadze et al. (2019) and Kashiyama
et al. (2019). If the adiabatic index is $\gamma=5/3$, the adiabatic constant
can be estimated as $K\sim 4\times 10^{31}$ by differentiating equation (37)
with $t$. Substituting these values to equation (34), we can obtain a
criterion that the self-gravity of ejecta mass can be neglected as

${\color[rgb]{0,0,0}t\ll 4\times 10^{5}\mathrm{yr}.}$ (38)

This is much larger than the evolution timescale of the massive white dwarf
system which is estimated to be a few $10^{3}$ year (Kashiyama et al., 2019).
Thus this work may explain the outflow evolution.

Moreover, Figure 6 shows that the radial velocity is below the escape velocity
in the inner region. By contrast, matter in the outer most region is freely
expanding. The double power law in the density distribution is a result of a
”weak” explosion in which the explosion energy is a fraction of the binding
energy of the original matter. Ejecta with such a feature have been observed
from numerical simulations for the eruption of the envelope of a massive star
prior to a supernova Tsuna et al. (2021). Thus there is potential for
applications in the formation of the circumstellar matter of an interaction-
powered supernova even if there is rotation.

In this work, we only consider the outflows on the equatorial plane, however,
we can use the same basic equations (1)-(3) to study the angular dependency of
the outflows (e.g., Müller & Vink, 2014). Müller & Vink (2014) provided two-
dimensional outflow models with rotation by assuming zero angular flow
velocity. Suppose there is no interaction between radial flows with different
polar angles. In that case, it is possible to investigate the properties of
rotating outflows in a two-dimensional plane with each polar angle. Those are
future works.

##  4 Conclusions

In this work, we have constructed a rotating transonic self-similar outflow in
the equatorial plane. We find that outflow solutions exist in the case of
$\gamma>11/9$, and these solutions have double-power-law density profiles as
follows:

$\rho(r,t)\propto r^{-\frac{1}{\gamma-1}}\rm{(inner)},$ (39) $\rho(r,t)\propto
t^{\frac{6\gamma-8}{\gamma-1}}r^{-\frac{9\gamma-11}{\gamma-1}}\rm{(outer)}.$
(40)

Such double-power-law density profiles are also common in the ejecta erupted
from binary merger systems.

The effect of rotation is noticeable in the subsonic region. In the subsonic
region, the inner singular point appears and the location of the point
approaches the transonic point as the rotation increases. This behavior is
consistent with the rotating steady outflow. On the other hand, the rotation
does not affect the solution of the supersonic region. Though our solutions
have some common features to steady rotating outflow models, our solutions
never reach any steady states and continue to expand because of the existence
of the inner sonic point.

The self-similar solutions in this paper do not include the effect of the
magnetic fields. Magnetic fields are keys to understanding an outflow from
merger remnant (Kashiyama et al., 2019; Schneider et al., 2019). In order to
study an outflow from a fast rotating white dwarf with strong magnetic fields,
it is necessary to construct a self-similar outflow model with magnetic
fields, which will be our future works.

This work is also supported by JSPS KAKENHI grant Nos. 20K14512 (KF),
20H05639, 22K03688, and 22K03671 (TS), MEXT, Japan.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Bildsten et al. (2007) Bildsten, L., Shen, K. J., Weinberg, N. N., & Nelemans, G. 2007, ApJ, 662, L95, doi: 10.1086/519489
  * Cheng (1977) Cheng, A. F. 1977, ApJ, 213, 537, doi: 10.1086/155185
  * Dan et al. (2011) Dan, M., Rosswog, S., Guillochon, J., & Ramirez-Ruiz, E. 2011, ApJ, 737, 89, doi: 10.1088/0004-637X/737/2/89
  * Fink et al. (2007) Fink, M., Hillebrandt, W., & Röpke, F. K. 2007, A&A, 476, 1133, doi: 10.1051/0004-6361:20078438
  * Fukue (1984) Fukue, J. 1984, PASJ, 36, 87
  * Guillochon et al. (2010) Guillochon, J., Dan, M., Ramirez-Ruiz, E., & Rosswog, S. 2010, ApJ, 709, L64, doi: 10.1088/2041-8205/709/1/L64
  * Gvaramadze et al. (2019) Gvaramadze, V. V., Gräfener, G., Langer, N., et al. 2019, Nature, 569, 684, doi: 10.1038/s41586-019-1216-1
  * Hotokezaka et al. (2018) Hotokezaka, K., Kiuchi, K., Shibata, M., Nakar, E., & Piran, T. 2018, ApJ, 867, 95, doi: 10.3847/1538-4357/aadf92
  * Iben & Tutukov (1984) Iben, I., J., & Tutukov, A. V. 1984, ApJS, 54, 335, doi: 10.1086/190932
  * Kashiyama et al. (2019) Kashiyama, K., Fujisawa, K., & Shigeyama, T. 2019, ApJ, 887, 39, doi: 10.3847/1538-4357/ab4e97
  * Kato & Hachisu (1994) Kato, M., & Hachisu, I. 1994, ApJ, 437, 802, doi: 10.1086/175041
  * Lattimer & Schramm (1974) Lattimer, J. M., & Schramm, D. N. 1974, ApJ, 192, L145, doi: 10.1086/181612
  * Li & Paczyński (1998) Li, L.-X., & Paczyński, B. 1998, ApJ, 507, L59, doi: 10.1086/311680
  * Lyutikov (2018) Lyutikov, M. 2018, arXiv e-prints, arXiv:1811.10777.  https://arxiv.org/abs/1811.10777
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, MNRAS, 406, 2650, doi: 10.1111/j.1365-2966.2010.16864.x
  * Moriya (2016) Moriya, T. J. 2016, ApJ, 830, L38, doi: 10.3847/2041-8205/830/2/L38
  * Müller & Vink (2014) Müller, P. E., & Vink, J. S. 2014, A&A, 564, A57, doi: 10.1051/0004-6361/201323031
  * Oskinova et al. (2020) Oskinova, L. M., Gvaramadze, V. V., Gräfener, G., Langer, N., & Todt, H. 2020, A&A, 644, L8, doi: 10.1051/0004-6361/202039232
  * Pakmor et al. (2013) Pakmor, R., Kromer, M., Taubenberger, S., & Springel, V. 2013, ApJ, 770, L8, doi: 10.1088/2041-8205/770/1/L8
  * Parker (1961) Parker, E. N. 1961, ApJ, 133, 1014, doi: 10.1086/147105
  * Parker (1965) —. 1965, Space Sci. Rev., 4, 666, doi: 10.1007/BF00216273
  * Piro & Kulkarni (2013) Piro, A. L., & Kulkarni, S. R. 2013, ApJ, 762, L17, doi: 10.1088/2041-8205/762/2/L17
  * Saio & Nomoto (1985) Saio, H., & Nomoto, K. 1985, A&A, 150, L21
  * Sato et al. (2015) Sato, Y., Nakasato, N., Tanikawa, A., et al. 2015, ApJ, 807, 105, doi: 10.1088/0004-637X/807/1/105
  * Sato et al. (2016) —. 2016, ApJ, 821, 67, doi: 10.3847/0004-637X/821/1/67
  * Schneider et al. (2019) Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211, doi: 10.1038/s41586-019-1621-5
  * Schwab et al. (2016) Schwab, J., Quataert, E., & Kasen, D. 2016, MNRAS, 463, 3461, doi: 10.1093/mnras/stw2249
  * Sedov (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics
  * Shibata & Hotokezaka (2019) Shibata, M., & Hotokezaka, K. 2019, Annual Review of Nuclear and Particle Science, 69, 41, doi: 10.1146/annurev-nucl-101918-023625
  * Taam & van den Heuvel (1986) Taam, R. E., & van den Heuvel, E. P. J. 1986, ApJ, 305, 235, doi: 10.1086/164243
  * Tanikawa et al. (2015) Tanikawa, A., Nakasato, N., Sato, Y., et al. 2015, ApJ, 807, 40, doi: 10.1088/0004-637X/807/1/40
  * Taylor (1950) Taylor, G. 1950, Proceedings of the Royal Society of London Series A, 201, 159, doi: 10.1098/rspa.1950.0049
  * Tsujimoto & Shigeyama (2014) Tsujimoto, T., & Shigeyama, T. 2014, A&A, 565, L5, doi: 10.1051/0004-6361/201423751
  * Tsuna et al. (2021) Tsuna, D., Takei, Y., Kuriyama, N., & Shigeyama, T. 2021, PASJ, 73, 1128, doi: 10.1093/pasj/psab063
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355, doi: 10.1086/161701
  * Weber & Davis (1967) Weber, E. J., & Davis, Leverett, J. 1967, ApJ, 148, 217, doi: 10.1086/149138

Figure 1: The relation between the dimensionless sound velocity $C$ and $x_{\rm crit}$ at the transonic point for a few $\gamma$’s. From equation (18), this relation has distinct features depending on whether the value of $\gamma$ is greater than $11/9$ or not. The three panel shows the three cases: $\gamma>11/9$ (left panel), $\gamma=11/9$ (middle panel) and $\gamma<11/9$ (right panel). For $\gamma>11/9$, $x_{\rm{crit}}$ has a maximum value, but for $\gamma\leq 11/9$, $x_{\rm{crit}}$ can take all positive values. Table 1: The asymptotic forms of the dimensionless variables $\Omega$ and $V_{\phi}$ for each type. In types I and II, the exponent is for $\gamma=7/6$. In types III, IV and V, the exponent is for general $\gamma$. While type III, IV, and V solutions exhibit the same asymptotic behavior of the dimensionless variables, the density distributions ($\rho(r,t)$) of these types are different because of different $\gamma$. It should be noted that when the flow does not truncate at $x>x_{\rm{crit}}$, we set $x_{\rm{max}}=\infty$. Type | Asymptotic forms for $x_{\rm{min}}\ll x\ll x_{\rm{crit}}$ | Asymptotic forms for $x_{\rm{crit}}\ll x\ll x_{\rm{max}}$ | Examples of the parameter sets
---|---|---|---
I | $V_{\phi}\propto x^{-0.68}$ | $V_{\phi}\propto x^{-0.72}$ | $(\gamma,x_{\mathrm{crit}},\alpha)$
| $\Omega\propto x^{-4.9}$ | $\Omega\propto x^{-4.2}$ | $=(7/6,0.05,0.2)$
II | $V_{\phi}\propto x^{-0.72}$ | $V_{\phi}\propto x^{-0.72}$ | $(\gamma,x_{\mathrm{crit}},\alpha)$
| $\Omega\propto x^{-4.9}$ | $\Omega\propto x^{-4.2}$ | $=(7/6,0.05,0.1)$
III, IV | $V_{\phi}\propto x^{-2/3}$ | $V_{\phi}\propto x^{-1}$ | $(\gamma,x_{\mathrm{crit}},\alpha)$
| $\Omega\propto x^{-\frac{1}{\gamma-1}}$ | $\Omega\propto x^{-3}$ | $=(7/6,0.0005,0.2)$
V | $V_{\phi}\propto x^{-2/3}$ | $V_{\phi}\propto x^{-1}$ | $(\gamma,x_{\mathrm{crit}},\alpha)$
| $\Omega\propto x^{-\frac{1}{\gamma-1}}$ | $\Omega\propto x^{-3}$ | $=(5/3,0.0005,0.1)$
Figure 2: The dimensionless Mach numbers as functions of $x$ for several
parameter sets. Some solutions truncate at the point where $v=C$ (i.e., $M=1$)
again inside and/or outside of the transonic point $x=x_{\rm{crit}}$. The
other solutions truncate at the point where $v=0$ or $M=0$ inside of the
transonic point. The left two panels show the distribution of $M(x)$ for
$\gamma=7/6$. The middle two show that for $\gamma=11/9$ and the right two
show that for $\gamma=5/3$. The upper and lower panels show solutions with
different coordinates $x_{\rm{crit}}$ of transonic points. Figure 3: A
comparison of dimensionless variables as functions of $x$ between different
types of solutions. In this figure, the following parameter sets
$(\gamma,x_{\mathrm{crit}},\alpha)$ are used: type I (7/6,0.05,0.2), type II
(7/6,0.05,0.1), type III (7/6,0.0005,0.2), type IV (7/6,0.0005,0.1) and type V
(5/3,0.0005,0.1). Figure 4: The relation between $x_{\rm{crit}}$ and
$x_{\rm{max}}/x_{\rm{crit}}$ for two different adiabatic indices denoted in
the legend. Here, $x_{\rm{max}}$ denotes the coordinate of the outer sonic
point. Both of the lines show the type of the flow changes from II to IV with
decreasing $x_{\rm{crit}}$. A large $x_{\rm{crit}}$ inhibits the integration
toward the outer region and thus results in $x_{\rm{max}}$ very close to
$x_{\rm{crit}}$. Note that the dashed horizontal line indicates $x_{\rm
max}/x_{\rm crit}=1$. Figure 5: The relation between the rotation strength
$\alpha$ and $x_{\rm{min}}/x_{\rm{crit}}$. Here, $x_{\rm{min}}$ denotes the
coordinate of the inner truncation point. The flow truncates because $v=C$
(solid lines) or $v=0$ (dashed lines) at $x=x_{\rm{min}}$. For the truncation
with $v=C$, $x_{\rm{min}}$ is found to monotonically increase with increasing
rotation ($|\alpha|$). For large $\alpha$,
$x_{\rm{min}}/x_{\rm{crit}}\propto\exp(b\alpha)$ is satisfied as far as
$\gamma\leq 11/9$. When the parameter set is
$(\gamma,x_{\rm{crit}})=(7/6,0.05)$, the value of $b$ becomes about 28. Figure
6: Time evolution of physical quantities of type V solutions. The upper two
panels and the bottom left panel show the time evolution of the radial
velocity (the top left), the azimuthal velocity (the top right), and the
density (the bottom left) of type V solution
$(\gamma,x_{\mathrm{crit}},\alpha)=(5/3,0.05,0.1)$ where the azimuthal
velocity does not exceed the Kepler rotation anywhere. The lower right panel
shows the time evolution of the density of type V solution
$(\gamma,x_{\mathrm{crit}},\alpha)=(3/2,0.05,0.005)$ where the azimuthal
velocity exceeds the Kepler rotation (the bottom right panel) in the innermost
region. Each dot indicates a transonic point where $x=x_{\rm{crit}}$. For the
former solution, quasi-hydrostatic equilibrium is established because the flow
velocity is subsonic in the inner region. For the latter solution, the
rotational effect is not negligible and the radial velocity approaches the
sound speed in the inner region.