# COEM: Cross-Modal Embedding for MetaCell Identification

Haiyi Mao  Minxue Jia  Jason Xiaotian Dou  Haotian Zhang  Panayiotis V. Benos

###### Abstract

Metacells are disjoint and homogeneous groups of single-cell profiles,
representing discrete and highly granular cell states. Existing metacell
algorithms tend to use only one modality to infer metacells, even though
single-cell multi-omics datasets profile multiple molecular modalities within
the same cell. Here, we present Cross-MOdal Embedding for MetaCell
Identification (COEM), which utilizes an embedded space leveraging the
information of both scATAC-seq and scRNA-seq to perform aggregation, balancing
the trade-off between fine resolution and sufficient sequencing coverage. COEM
outperforms the state-of-the-art method SEACells by efficiently identifying
accurate and well-separated metacells across datasets with continuous and
discrete cell types. Furthermore, COEM significantly improves peak-to-gene
association analyses, and facilitates complex gene regulatory inference tasks.

Machine Learning, ICML

##  1 Introduction

Recent advances in single-cell technologies have enabled measuring chromatin
accessibility and gene expression levels simultaneously at single-cell
resolution. Various sequencing technologies, such as sci-CAR (Cao et al.,
2018), SNARE-seq (Chen et al., ), SHARE-seq (Ma et al., 2020), and 10X
Genomics Multiome, have been developed and applied to comprehensive studies of
cell heterogeneity, developmental dynamics, and cis-regulatory elements (CREs)
(Chu et al., 2022). However, extreme data sparsity is observed in Single-cell
RNA sequencing (scRNA-seq) and Single-cell sequencing assay for transposase-
accessible chromatin (scATAC-seq) data. This impedes cis-regulation inference
at the single-cell level.

Figure 1: Violin plot of Pearson correlation between metacell-aggregated gene
expression and accessibility on (A) CD34+ bone marrow dataset and (B) T-cell
depleted bone marrow dataset. There are a considerable number of negative
peak-to-gene associations in metacells identified by SEACell with scATAC-seq,
while use of COEM-identified metacells improves peak-to-gene correlation.

High-throughput single-cell profiling of biological samples typically leads to
repetitive sampling of highly similar and statistically equivalent cells. The
concept of metacells (Baran et al., 2019) has been proposed to maintain
statistical confidence. Metacells are sets of scRNA-seq cell profiles
representing distinct, highly granular cell states. Metacell aggregates, with
their sufficient sequencing coverage, help reduce sparsity-derived problems in
downstream analyses.

The MetaCells algorithm (Baran et al., 2019) infers metacells based on the
partition of k-nearest neighbor (KNN) similarity graphs on scRNA-seq data
only, but it fails on scATAC-seq data. The SEACells approach (Single-cell
aggregation of cell-states) (Persad et al., 2022) is designed to identify
metacells in RNA or ATAC modality. It significantly improves peak-to-gene
associations, i.e. between chromatin accessibility peak and gene expression.
However, SEACells identifies metacells based on only one modality: holding
onto a strict biological assumption that accessible chromatin is associated
with active transcription consistently, despite different data modalities
possibly exhibiting a time lag between chromatin remodeling of a gene and its
transcription. The occurrence of this ”time lag” manifests when:

  1. 1.

Chromatin opens before transcription initiates

  2. 2.

mRNA degrades follows chromatin closing (Li et al., 2021).

So, it is not surprising that we observe negative peak-to-gene associations
across metacells, as Fig. 1 shows. We aim to investigate that integrated
analysis of scATAC-seq and scRNA-seq data to improve our ability to discover
and characterize cell states, while avoiding biased associations.

Our new algorithm, COEM, is designed to integratively analyze multi-omic
sequencing single-cell data (e.g., gene expression and chromatin accessibility
measured in the same cell). COEM addresses the fusion challenge (Baltrušaitis
et al., 2018) by first learning joint low-dimensional latent representations
for scRNA-seq and scATAC-seq data through a multi-view, multi-modal
variational auto-encoder (VAE) model; further utilizing a spectral clustering
approach, based on graph connectivity, to identify metacells representing
highly granular, distinct cell states on both RNA and chromatin levels. We
found COEM to be more accurate than SEACells in metacell identification across
datasets with discrete cell types and continuous developmental trajectories.
Also, COEM-identified metacells substantially improve the peak-to-gene
association analysis and empower CREs prediction.

Figure 2: The main steps of COEM. 1. The joint profiled scRNA-seq and scATAC-
seq are fed to the modified VAE to extract a common latent representations
$\mathbf{z}$. 2. Then a similarity matrix are constructed based on
$\mathbf{z}$ and KNN graphs. 3. Spectral clustering is applied to identify the
higher granularity clusters where Metacells belongs to. 4. The Metacells
features are aggregated for each clusters. 5. Metacells are availabe for
downstreaming tasks.

##  2 Method

###  2.1 Multi-Modal Joint Embedding

COEM learns a low dimensional representation from the joint scRNA-seq and
scATAC-seq data, measured in the same cell. We adopt a modified multi-view VAE
to extract a common latent space of the two data modalities. This is similar
to previous work in VAEs representing the single-cell sequence data (Lopez et
al., 2018; Li et al., 2022; Ding et al., 2018; Mao et al., 2020; Xu et al.,
2020) in less noisy low dimensions. The generative model of the VAEs is
derived as follows. Given a cell, the scRNA gene expression
$\mathbf{x}\in\mathcal{P}_{\mathbf{x}}$, scATAC chromatin accessibility data
$\mathbf{y}\in\mathcal{P}_{\mathbf{y}}$, cell type $c$, and the common latent
representation $\mathbf{z}$, we have the following probability densities,

$\displaystyle
p_{\theta}(\mathbf{x},\mathbf{z}|c)=p_{\theta}(\mathbf{x}|\mathbf{z},\mathbf{c})p_{\theta}(\mathbf{z}|c)]$
(1) $\displaystyle
p_{\theta}(\mathbf{y},\mathbf{z}|c)=p_{\theta}(\mathbf{y}|\mathbf{z},\mathbf{c})p_{\theta}(\mathbf{z}|c),$

where $\theta$ represents the parameters of the generative models. The
likelihood $p_{\theta}$ of $\mathbf{x}$ and $\mathbf{y}$ depends on the common
latent space $\mathbf{z}$ and the cell type $c$. With the assumption that the
distribution of observed scRNA-seq $\mathbf{x}$, scATAC-seq $\mathbf{y}$ are
independent given the common latent representations $\mathbf{z}$ and cell type
$c$, $\mathbf{x}\perp\\!\\!\\!\perp\mathbf{y}|\mathbf{z},c$. So, we have
following joint generative model,

$\displaystyle p_{\theta}(\mathbf{x},\mathbf{y},\mathbf{z}|c)$
$\displaystyle=p_{\theta}(\mathbf{x},\mathbf{y}|\mathbf{z},c)p_{\theta}(\mathbf{z}|c)$
(2)
$\displaystyle=p_{\theta}(\mathbf{x},|\mathbf{z},c)p_{\theta}(\mathbf{y},|\mathbf{z},c)p_{\theta}(\mathbf{z}|c).$

Moreover, we have the following Evidence Lower Bound (ELBO) for two modalities
with conditional variational approximation $q_{\phi}$,

$\displaystyle\log p_{\theta}(\mathbf{x}|c)\geq$ $\displaystyle-
KL(q_{\phi}(\mathbf{z}|\mathbf{x},c))||p_{\theta}(\mathbf{z}|c))$ (3)
$\displaystyle+\mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x},c))}(\log
p_{\theta}(\mathbf{x}|\mathbf{z},c))$ $\displaystyle\log
p_{\theta}(\mathbf{y}|c)\geq$ $\displaystyle-
KL(q_{\phi}(\mathbf{z}|\mathbf{y},c))||p_{\phi}(\mathbf{z}|c))$
$\displaystyle+\mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{y},c)}(\log
p_{\theta}(\mathbf{y}|\mathbf{z},c)).$

Normally, we use an inference network (encoder) to learn the variational
approximation $q_{\phi}$ and a generative network (decoder) to learn
$p_{\theta}$. Particularly, we assume $p_{\theta}(\mathbf{x}|\mathbf{z},c)$ is
a zero-inflated Poisson (ZIP) distribution,
$p_{\theta}(\mathbf{y}|\mathbf{z},c)$ is a negative binomial (NB)
distribution, and $p_{\theta}(\mathbf{z}|c)$ is a Gaussian Mixture Model
$\mathbf{z}\sim\mathcal{N}(\mathbf{u}_{c},\Sigma_{c})$ (Li et al., 2022; Lopez
et al., 2018). Given the proposed joint-embedding generative model, we adopt
the scMVP (Li et al., 2022) architecture to model the generative network and
inference network. Specifically, there is a two-channel attention-based
encoder network for the scRNA-seq and scATAC-seq inputs, then concatenating
them to derive the posterior distribution
$p_{\phi}(\mathbf{z}|\mathbf{x},\mathbf{y},c)$. Next, the imputed scRNA and
scATAC profiles are reconstructed by an attention based two-channel decoder
network to compute the likelihood
$p_{\theta}(\mathbf{x},\mathbf{y}|\mathbf{z},c)$.

###  2.2 Metacells Identification

Similar to previous work (Baran et al., 2019; Persad et al., 2022), COEM is a
graph-based approach that employs spectral clustering to compute metacells. To
be more specific, COEM consists of the following four steps:

  1. 1.

Apply cross-modal joint embedding VAE to learn low common representations
$\mathbf{z}$ to encode the common information of scRNA-seq and scATAC-seq.

  2. 2.

Build a k-nearest neighbors (KNN) graph on the low dimensional space
$\mathbf{z}$ based on Euclidean distances.

  3. 3.

Then similarity matrix is constructed from the previous KNN graph by
connectivity. The similarity matrix is based on radial basis function (RBF)
kernel to encode the non-linear relations in the $\mathbf{z}$.

  4. 4.

Utilize spectral clustering to identify the cluster and to subsequently
aggregate metacells as in (Persad et al., 2022).

This is one of the differences between COEM and SEACells (Persad et al.,
2022), which uses archetypal analysis for metacell identification, instead of
spectral clustering. While they both are graph based methods, instead of the
convexity assumption of archetypal analysis of SEACells, spectral clustering
of COEM purely utilizes the graphs’ connectivity information and enjoys a
faster running speed. We empirically validate that spectral clustering can
improve metacell identification and running time.

methods | CD34 | sciCAR | 10xPBMC | Snare cellline
---|---|---|---|---
Spectral | 23.7 | 6.3 | 45.2 | 1.7
Archetypal | 191.5 | 103.6 | 277.4 | 28.0
Table 1: Running Time (seconds). We use the same low dimensional embeddings $\mathbf{z}$ to compare archetypal analysis (used in SEACells) and spectral clustering (used in COEM). | RNA_pca |
---|---|---
methods | separation | compactness
COEM | 0.235 $\tiny\pm 0.015$ |  0.193 $\tiny\pm 0.017$
SeaC-ATAC |  0.196 $\tiny\pm 0.007$ |  0.686 $\tiny\pm 0.035$
SeaC-RNA |  0.177 $\tiny\pm 0.051$ | 0.150 $\tiny\pm 0.062$
| ATAC_svd |
methods | separation | compactness
COEM | 0.091 $\tiny\pm 0.006$ | 0.028 $\tiny\pm 0.002$
SeaC-ATAC |  0.086 $\tiny\pm 0.002$ |  0.004 $\tiny\pm 0.001$
SeaC-RNA |  0.041 $\tiny\pm 0.015$ |  0.032 $\tiny\pm 0.014$
methods | Cell-type Purity |
COEM | 0.907 $\tiny\pm 0.006$ |
SeaC-ATAC |  0.666 $\tiny\pm 0.01$ |
SeaC-RNA |  0.734 $\tiny\pm 0.002$ |
Table 2: Metacacell Evaluation on sci-CAR celline dataset. Greater separation,
lower compactness, and higher cell-type purity indicate better performance.
COEM identified metacells with higher purity and separation
| RNA_pca |
---|---|---
methods | separation | compactness
COEM | 0.801 $\tiny\pm 0.022$ |  0.419 $\tiny\pm 0.001$
SeaC-ATAC |  0.723 $\tiny\pm 0.024$ |  1.061 $\tiny\pm 0.140$
SeaC-RNA |  0.704 $\tiny\pm 0.041$ | 0.292 $\tiny\pm 0.217$
| ATAC_svd |
methods | separation | compactness
COEM | 0.710 $\tiny\pm 0.013$ |  0.662 $\tiny\pm 0.019$
SeaC-ATAC |  0.665 $\tiny\pm 0.056$ |  0.640 $\tiny\pm 0.153$
SeaC-RNA |  0.566 $\tiny\pm 0.070$ | 0.547 $\tiny\pm 0.291$
Table 3: Metacacell Evaluation on T-cell depleted bone marrow dataset. COEM
metacells have higher separation than SEACell’s

##  3 Experiments

###  3.1 Evaluation Metrics

As in the SEACells paper (Persad et al., 2022), we use compactness,
separation, and cell-type purity for the metacell benchmark. Compactness
measures how homogeneous the cells within a metacell are. Separation assesses
how metacells are distinct from each other. Cell-type purity, instead,
measures the consistency of cell-types among cells that constitute a metacell.
We computed diffusion components using principal components from RNA assay
(RNA_PCA) or singular values singular values from ATAC assay (ATAC_SVD), and
then quantified the compactness and separation of metacells. Since SEACells
identifies metacells from one modality, we denote them as SeaC-ATAC and SeaC-
RNA.

###  3.2 Results

We use five datasets for evaluation: sci-CAR cell line (Cao et al., 2018),
SNARE-seq cell line (Chen et al., ), 10X Genomics Multiome PBMC, 10X Multiome
CD34+ bone marrow (hematopoietic stem and progenitor cells) and 10X Multiome
T-cell depleted bone marrow (Persad et al., 2022). As showed in Table 2 and 3,
and Table 4 ,5, and 6 in appendix, COEM outperforms SEACells in most cases.

The sci-CAR dataset profiles well-labeled cells, including cell lines 293T,
3T3, a 293T/3T3 cell mixture, and A549 cells after 0, 1, or 3 hours of
dexamethasone treatment. We evaluate cell-type purity in sci-CAR dataset, and
metacells from COEM algorithm has greater purity (0.91) than SeaC-ATAC (0.67)
and SeaC-RNA (0.74), which suggests that leveraging two modalities improves
our ability to define cellular states, whereas achetypal analysis on joint
embedding representations cannot guarantee high purity of metacells shown in
Table 2.

Then, we measure compactness and separation across datasets with discrete cell
types. Generally, COEM-based methods separate metacells better than SeaC-ATAC
and SeaC-RNA shown in Tables 3, 2, 4, 5 and 6. But they may have higher
compactness based on only one modality, since COEM balances cell variability
from two modalities at the same time. The metacells identified by SEACells
have lower compactness because metacell identification and compactness
evaluation are based on the diffusion components from a single modality. But
those metacells may not be well-separated and compacted based on the diffusion
components from the other modality. Therefore, the association analysis
between gene expression and peak accessibility may be biased using metacells
inferred from only one modality. Further, spectral clustering outperforms
archetypal analysis for identifying metacells with high cell-type purity and
separation based on low dimensional embedding $\mathbf{z}$. We also measure
running time with the same low dimensional embedding $\mathbf{z}$ for
archetypal analysis and spectral clustering shown in Table 1. We can see
spectral clustering has generally one order of magnitude shorter running time
in the four datasets.

Analyzing chromatin accessibility and gene expression jointly allows us to
infer the relationship between open-chromatin peaks and active transcription,
aiding COEM in CRE discovery. Metacells can be used as a wise strategy to have
high resolution and sufficient sequencing coverage for gene regulation
inference from single-cell data. We evaluate the co-variation in chromatin
accessibility and gene expression across metacells by computing Pearson
correlations for each peak within +/- 100kb of the gene. We observe a large
portion of peaks is negatively correlated with gene expression in SeaC-ATATC
metacells (CD34+ bone marrow: 25.1%; T-cell depleted bone marrow: 23.0%)
(Figure 1), while a small portion of negative peak-to-gene associations is
identified in COEM metacells (CD34+ bone marrow: 1.3%; T-cell depleted bone
marrow: 3.0%), which is consistent with reported ratio of negative peak-to-
gene associations (1.2 $\sim$ 11 %) at the single-cell level. This shows that
the bias from the “time lag” phenomenon during differentiation is exaggerated
in SEACells metacells but not in COEM. Besides, we could observe stronger
peak-to-gene correlations for core genes based on COEM metacells. For example,
GATA2 gene acts as a master regulator of proliferation and maintenance of
hematopoietic stem and progenitor cells. The peak-to-gene correlation of GATA2
in COEM metacells from CD34+ bone marrow dataset is 0.73, while its
correlation is 0.62 in SeaC-ATAC metacells. (The peak-to gene correlation of
GATA2 is 0.1 at single-cell level)

##  4 Discussion

Our experiments show that integrating information from both scRNA-seq and
scATAC-seq is beneficial for identifying accurate and robust metacells from
single-cell data. Especially joint-embedding, low-dimensional representations
can enable the identification of cellular states with high cell-type purity
and well-separated metacells. On top of that, the proposed COEM algorithm
enjoys a faster running time. Moreover, it avoids high numbers of negative
peak-to-gene links that are false. Hence, the joint profiled sequence data are
able to depict a comprehensive landscape of cells’ stages. Furthermore,
(Bilous et al., 2021) has demonstrated that metacells inferred from scRNA-seq
are compatible with RNA velocity model. Therefore, COEM is also able to be
applied to mutli-omic velocity model for estimating epigenomic and
transcriptomic dynamics.

In future work, we would explore more machine learning techniques for multi-
modal sequence data integration and metacell identification. For example,
scATAC-seq and scRNA-seq are naturally causally related. It is interesting to
utilize this unique relation to learn more robust, causally sufficient, and
efficient representations (Blei, 2020; Wang & Jordan, 2021; Veitch et al.,
2021; Dou et al., 2022; Tabib et al., 2021). Besides, the results are
promising for use of an optimal transport approach to model cellular stage
changes across different modalities on the metacell level (Schiebinger et al.,
2019).

## References

  * Baltrušaitis et al. (2018) Baltrušaitis, T., Ahuja, C., and Morency, L.-P.  Multimodal machine learning: A survey and taxonomy.  _IEEE transactions on pattern analysis and machine intelligence_ , 2018.
  * Baran et al. (2019) Baran, Y., Bercovich, A., Sebe-Pedros, A., Lubling, Y., Giladi, A., Chomsky, E., Meir, Z., Hoichman, M., Lifshitz, A., and Tanay, A.  Metacell: analysis of single-cell rna-seq data using k-nn graph partitions.  _Genome biology_ , 20(1):1–19, 2019.
  * Bilous et al. (2021) Bilous, M., Tran, L., Cianciaruso, C., Gabriel, A., Michel, H., Carmona, S., Pittet, M., and Gfeller, D.  Metacells untangle large and complex single-cell transcriptome networks (preprint).  2021\.
  * Blei (2020) Blei, D. M.  Representation learning: A probabilistic perspective.  2020\.
  * Cao et al. (2018) Cao, J., Cusanovich, D. A., Ramani, V., Aghamirzaie, D., Pliner, H. A., Hill, A. J., Daza, R. M., McFaline-Figueroa, J. L., Packer, J. S., Christiansen, L., Steemers, F. J., Adey, A. C., Trapnell, C., and Shendure, J.  Joint profiling of chromatin accessibility and gene expression in thousands of single cells.  _Science_ , 361(6409):1380–1385, 2018.
  * (6) Chen, S., Lake, B. B., and Zhang, K.  High-throughput sequencing of the transcriptome and chromatin accessibility in the same cell.  _Nature Biotechnology_ , 37(12):1452–1457.  ISSN 1087-0156, 1546-1696.
  * Chu et al. (2022) Chu, T., Wang, Z., Pe’er, D., and Danko, C. G.  Cell type and gene expression deconvolution with bayesprism enables bayesian integrative analysis across bulk and single-cell rna sequencing in oncology.  _Nature Cancer_ , 2022.
  * Ding et al. (2018) Ding, J., Condon, A., and Shah, S. P.  Interpretable dimensionality reduction of single cell transcriptome data with deep generative models.  _Nature communications_ , 9(1):1–13, 2018.
  * Dou et al. (2022) Dou, J. X., Luo, L., and Yang, R. M.  An optimal transport approach to deep metric learning (student abstract).  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 36(11):12935–12936, Jun. 2022.
  * Li et al. (2021) Li, C., Virgilio, M., Collins, K. L., and Welch, J. D.  Single-cell multi-omic velocity infers dynamic and decoupled gene regulation.  _bioRxiv_ , 2021.
  * Li et al. (2022) Li, G., Fu, S., Wang, S., Zhu, C., Duan, B., Tang, C., Chen, X., Chuai, G., Wang, P., and Liu, Q.  A deep generative model for multi-view profiling of single-cell rna-seq and atac-seq data.  _Genome biology_ , 23(1):1–23, 2022.
  * Lopez et al. (2018) Lopez, R., Regier, J., Cole, M. B., Jordan, M. I., and Yosef, N.  Deep generative modeling for single-cell transcriptomics.  _Nature methods_ , 15(12):1053–1058, 2018.
  * Ma et al. (2020) Ma, S., Zhang, B., LaFave, L. M., Earl, A. S., Chiang, Z., Hu, Y., Ding, J., Brack, A., Kartha, V. K., Tay, T., Law, T., Lareau, C., Hsu, Y.-C., Regev, A., and Buenrostro, J. D.  Chromatin potential identified by shared single-cell profiling of rna and chromatin.  _Cell_ , 183(4):1103–1116.e20, 2020.  ISSN 0092-8674.
  * Mao et al. (2020) Mao, H., Broerman, M. J., and Benos, P. V.  Interpretable factors in scrna-seq data with disentangled generative models.  In _2020 IEEE 20th International Conference on Bioinformatics and Bioengineering_ , pp. 85–88. IEEE, 2020.
  * Persad et al. (2022) Persad, S., Choo, Z.-N., Dien, C., Masilionis, I., Chaligné, R., Nawy, T., Brown, C. C., Pe’er, I., Setty, M., and Pe’er, D.  Seacells: Inference of transcriptional and epigenomic cellular states from single-cell genomics data.  _bioRxiv_ , 2022.
  * Schiebinger et al. (2019) Schiebinger, G., Shu, J., Tabaka, M., Cleary, B., Subramanian, V., Solomon, A., Gould, J., Liu, S., Lin, S., Berube, P., et al.  Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming.  _Cell_ , 176(4):928–943, 2019.
  * Tabib et al. (2021) Tabib, T., Huang, M., Morse, N., Papazoglou, A., Behera, R., Jia, M., Bulik, M., Monier, D. E., Benos, P. V., Chen, W., et al.  Myofibroblast transcriptome indicates sfrp2hi fibroblast progenitors in systemic sclerosis skin.  _Nature communications_ , 12(1):1–13, 2021.
  * Veitch et al. (2021) Veitch, V., D' Amour, A., Yadlowsky, S., and Eisenstein, J.  Counterfactual invariance to spurious correlations in text classification.  In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), _Advances in Neural Information Processing Systems_ , volume 34, pp. 16196–16208, 2021.
  * Wang & Jordan (2021) Wang, Y. and Jordan, M. I.  Desiderata for representation learning: A causal perspective.  _arXiv preprint arXiv:2109.03795_ , 2021.
  * Xu et al. (2020) Xu, H., Luo, D., Henao, R., Shah, S., and Carin, L.  Learning autoencoders with relational regularization.  In _International Conference on Machine Learning_ , pp. 10576–10586, 2020.

##  Appendix A Appendix

In this section, we show the results of additional experiments on datasets
including Snare Cellline, 10x_PBMC, and 10X Multiome CD34+ bone marrow.

Table 4: Metacacell Evaluation on SNARE-seq celline dataset | RNA_pca |
---|---|---
methods | separation | compactness
COEM | 0.280 $\tiny{\pm 0.015}$ |  0.856 $\tiny{\pm 0.185}$
SeaC-ATAC |  0.191 $\tiny{\pm 0.013}$ |  0.617 $\tiny{\pm 0.112}$
SeaC-RNA |  0.067 $\tiny{\pm 0.011}$ | 0.067 $\tiny{\pm 0.004}$
| ATAC_svd |
methods | separation | compactness
COEM | 0.159 $\tiny{\pm 0.021}$ |  0.124 $\tiny{\pm 0.008}$
SeaC-ATAC |  0.097 $\tiny{\pm 0.010}$ | 0.004 $\tiny{\pm 0.001}$
SeaC-RNA |  0.095 $\tiny{\pm 0.018}$ |  0.037 $\tiny{\pm 0.012}$
Table 5: Metacacell Evaluation on CD34+ bone marrow dataset | RNA_pca |
---|---|---
methods | separation | compactness
COEM | 0.266 $\tiny{\pm 0.011}$ |  0.035 $\tiny{\pm 0.006}$
SeaC-ATAC |  0.203 $\tiny{\pm 0.013}$ |  0.037 $\tiny{\pm 0.007}$
SeaC-RNA |  0.256 $\tiny{\pm 0.037}$ | 0.026 $\tiny{\pm 0.003}$
| ATAC_svd |
methods | separation | compactness
COEM | 0.318 $\tiny{\pm 0.014}$ |  0.136 $\tiny{\pm 0.010}$
SeaC-ATAC |  0.284 $\tiny{\pm 0.015}$ | 0.025 $\tiny{\pm 0.003}$
SeaC-RNA |  0.244 $\tiny{\pm 0.123}$ |  0.115 $\tiny{\pm 0.021}$
Table 6: Metacacell Evaluation on PBMC dataset | RNA_pca |
---|---|---
methods | separation | compactness
COEM | 0.515 $\tiny{\pm 0.127}$ |  0.349 $\tiny{\pm 0.115}$
SeaC-ATAC |  0.406 $\tiny{\pm 0.129}$ |  3.315 $\tiny{\pm 1.21}$
SeaC-RNA |  0.428 $\tiny{\pm 0.133}$ | 0.294 $\tiny{\pm 0.046}$
| ATAC_svd |
methods | separation | compactness
COEM | 0.181 $\tiny{\pm 0.038}$ |  0.030 $\tiny{\pm 0.013}$
SeaC-ATAC |  0.169 $\tiny{\pm 0.024}$ | 0.022$\tiny{\pm 0.008}$
SeaC-RNA |  0.034 $\tiny{\pm 0.009}$ |  0.141 $\tiny{\pm 0.005}$


# Doppler Dimming and Brightening Effects in Solar Prominences

Aaron W. Peat,1,2 Christopher M. J. Osborne,2 and Petr Heinzel1,3
1University of Wrocław, Centre of Scientific Excellence – Solar and Stellar
Activity, Kopernika 11, 51-622 Wrocław, Poland
2SUPA School of Physics and Astronomy, University of Glasgow, Glasgow, G12
8QQ, UK
3Astronomical Institute, Czech Academy of Sciences, 25165 Ondřejov, Czech
Republic E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We explored the impact that Doppler dimming and brightening effects from bulk
motions of solar prominences have on the formation of Ly$~{}\alpha$,
H$~{}\alpha$, and Mg ii h line profiles. We compared two schemes in which
these effects manifest; when the prominence is moving radially away from the
solar surface (radial case), and when the prominence is moving parallel to the
solar surface (horizontal case). To do this, we analysed 13,332 model profiles
generated through the use of the 1D NLTE (i.e. departures from Local
Thermodynamic equilibrium) radiative transfer (RT) code Promweaver, built on
the Lightweaver NLTE RT framework to mimic the behaviour and output of the 1D
NLTE RT code PROM. We found that horizontal velocities are just as, or more
important than radial velocities. This demonstrates that horizontal velocities
need to be accounted for when attempting to do any sort of forward modelling.

###### keywords:

Sun: filaments, prominence – radiative transfer – line: profiles

††pubyear: 2024††pagerange: Doppler Dimming and Brightening Effects in Solar
Prominences–Doppler Dimming and Brightening Effects in Solar Prominences

##  1 Introduction

Solar prominences are cool, dense structures suspended in the solar corona
which have been observed for centuries (e.g. Vassenius, 1733; Secchi, 1875;
Hale, 1929; Vyssotsky, 1949) and simulated for half a century (e.g. Ishizawa,
1971; Heasley et al., 1974; Heinzel et al., 1987; Gouttebroze, 2004). It has
been long shown that eruptive prominences are affected by Doppler dimming
(DDE) and brightening (DBE) effects – where the incident radiation from the
solar disc appears Doppler shifted in the frame-of-reference of the solar
prominence, causing a brightening or dimming effect in the scattered line
radiation. From here, we will use the acronym DXE to refer to both DDE and
DBE. The former occurs when the incident radiation is in emission, and the
latter when in absorption (Hyder & Lites, 1970; Heinzel & Rompolt, 1987;
Labrosse et al., 2010; Vial & Engvold, 2015).

More generally, when a prominence undergoes some bulk motion, DBE occurs when
the incident radiation increases in the wavelength range of the absorption
profile of the line in question, and DDE occurs when the incident radiation
decreases in the wavelength range of the absorption profile of the line
(Tandberg-Hanssen, 1995).

An underexplored avenue which produces DXE is when a prominence has zero
radial velocity, but non-zero velocity parallel to the solar surface (e.g.
Rompolt, 1967; Morimoto & Kurokawa, 2003). We believe that no large scale
quantitative study has been performed to understand the impact of such a
velocity field on the radiation observed from solar prominences. In the case
of prominences, this is commonly referred to as line-of-sight (LOS) velocity –
however, in observations, true LOS velocities commonly contain some radial
component. Here, we wish to demonstrate that, even with a zero radial velocity
component, prominences (and filaments) undergoing bulk motions experience
significant DXE.

To begin, we will briefly introduce the NLTE (i.e. departures from Local
Thermodynamic equilibrium) radiative transfer (RT) code which we use to
simulate solar prominence emission. Following this, we briefly discuss our
results and compare them to the effect when caused by non-zero radial
velocities. We then discuss the implications this has for NLTE inversions of
solar prominence atmospheres. Finally, we offer our conclusions.

##  2 Promweaver

Figure 1: How PROM and Promweaver model a prominence. Where $H$ is the
altitude of the prominence, and $L$ is the slab width. This figure is adapted
from Peat (2023).

Figure 2: The incident radiation of Ly $\alpha$, H $\alpha$, and Mg ii h&k as
seen by the prominence for $\mu=1$ (in the frame of the prominence, i.e. the
radiation that would be incident on the prominence from directly below). In
the Mg ii h&k panel, the velocities are calculated using Mg ii h (2803.53 Å)
as the rest wavelength.

To demonstrate how non-zero velocities parallel to the solar surface in solar
prominences produce DXE, we employ the use of the 1D NLTE RT code Promweaver
(Osborne, 2022), built on top of the NLTE RT framework, Lightweaver (Osborne &
Milić, 2021). Promweaver is written to efficiently and expediently replicate
the results of the 1D NLTE RT code PROM (Gouttebroze et al., 1993; Heinzel et
al., 1994; Labrosse & Gouttebroze, 2004; Levens & Labrosse, 2019). However,
unlike PROM, Promweaver has the ability to include any arbitrary bulk velocity
field, meanwhile PROM can only include radial velocities. Promweaver is not
the first prominence code that can include arbitrary velocities; for example,
MALI (Heinzel, 1995; Heinzel et al., 2014) includes an implementation of the
technique discussed in Heinzel & Rompolt (1987). However, in practice, it has
only been used for radial velocities in the context of erupting prominences.
Additionally, PROM produces half line profiles which are mirrored to produce
the full profiles, but Promweaver synthesises the entire line profile handling
arbitrary velocity fields and interactions between overlapping transitions.

Promweaver allows us to simulate any of the commonly modelled and observed
lines; here we focus on Ly $\alpha$, H $\alpha$, and Mg ii h using a 5-level
plus continuum H model and a 10-level plus continuum Mg model, following
Leenaarts et al. (2013). Partial frequency redistribution (PRD) is of critical
importance to the formation of the Ly $\alpha$ line (Heinzel et al., 1987).
Additionally, PRD plays a key role in the formation of the Mg ii resonance
lines (Paletou et al., 1993; Heinzel et al., 2014), but this only affects the
far wings of the emission, which may be below the detection limit (Heinzel et
al., 2014). Here we treat all resonance lines in PRD. These atoms are solved
simultaneously, allowing for mutual radiative interactions to occur.

Promweaver, like PROM or MALI, models a solar prominence as a monolithic slab
suspended in the solar corona (see Figure 1). It solves for statistical
equilibrium of the considered atoms in the given slab, additionally conserving
charge, and modifying the mass density of the structure to follow the
prescribed pressure stratification (which is affected by ionisation).
Promweaver can use one of two model atmospheres; isothermal and isobaric, and
one with a prominence-to-corona transition region (PCTR). The PCTR is the
region where the cool and dense prominence transitions into the hot and
tenuous corona. This region is modelled by the following equations (Kippenhahn
& Schlüter, 1957; Anzer & Heinzel, 1999; Labrosse & Gouttebroze, 2004),

$T(m)=T_{\text{cen}}+(T_{\text{tr}}-T_{\text{cen}})\left(1-4\frac{m}{M}\left(1-\frac{m}{M}\right)\right)^{\gamma},$
(1) $p(m)=4p_{c}\frac{m}{M}\left(1-\frac{m}{M}\right)+p_{\text{tr}},$ (2)

for $\gamma\geq 2$, and where $p_{c}=p_{\text{cen}}-p_{\text{tr}}$. Here $m$
is the column mass, $M$ is the total column mass, $T(m)$ is the temperature as
a function of column mass, $T_{\text{cen}}$ is the central (prominence)
temperature of the plasma, $T_{\text{tr}}$ is the (coronal) temperature at the
edge of the PCTR, and $\gamma$ is a dimensionless factor that defines the
extent of the PCTR. $p(m)$ is the pressure as a function of column mass,
$p_{\text{cen}}$ is the central (prominence) pressure, and $p_{\text{tr}}$ is
the (coronal) pressure at the edge of the PCTR. Eq. 1 is an empirical
expression described in Anzer & Heinzel (1999), while Eq. 2 is an analytical
formula from the Kippehahn-Schlüter model (KS; Kippenhahn & Schlüter, 1957;
Heasley & Mihalas, 1976) for magnetohydrostatic (MHS) equilibrium. These
equations are used as Promweaver is built to mimic PROM and this is how the
PCTR is constructed in PROM (Labrosse & Gouttebroze, 2004). As Eq. 2 comes
from MHS equilibrium, it should be stressed that the applicability of this
model for erupting prominences is questionable.

The incident radiation is generated using the FAL C atmosphere (Fontenla et
al., 1993) and ‘light cone’ method discussed in the appendix of Jenkins et al.
(2023). The incident radiation at $\mu=1$ (i.e. radially from the solar
surface) in the lines considered is shown in Fig. 2. In all situations we
employ fully angularly dependent boundary conditions including the effects of
limb darkening (self-consistently in a plane-parallel approximation as per
Jenkins et al. (2023)) and arbitrary prominence motions.

The emergent intensity is computed by a final pass through the model using the
computed level populations and ratio of the PRD emission to absorption
profiles. The cones in the boundary conditions are disabled to compute the
radiation along a specific ray (rather than the spherical integrals used to
solve the non-LTE problem). This method can handle arbitrary viewing angles
through the prominence slab, including when it is partially back-lit against
the solar disc.

Parameter | Unit | Iso Models | PCTR Models
---|---|---|---
$T_{\text{cen}}$ | kK | 8 | 6
$T_{\text{tr}}$ | kK | – | 100
$p_{\text{cen}}$ | dyn cm-2 | 0.05 | 0.1
$p_{\text{tr}}$ | dyn cm-2 | – | 0.01
L | km | 400 | 407.69a
M | g cm-2 | 2$\times 10^{-6}$ a | 5$\times 10^{-6}$
$v_{T}$ | km s-1 | 5 | 5
$\gamma$ |  | – | 4
Table 1: The base model parameters. For the isothermal and isobaric (iso)
models, $T_{\text{cen}}$ and $p_{\text{cen}}$ are taken as $T$ and $p$,
respectively. The iso models are built using the slab width, while the PCTR
models use maximum column mass to limit the extent of the structure.
a For this model type, this parameter is not strictly an input. It instead is
found when solving statistical equilibrium.

Using Promweaver, we generated a large grid of models. These models are all
based on one of two base models, which can be seen in Table 1. For isothermal
models, slab width (L) is an input and the mass density (M) is calculated when
solving statistical equilibrium. Meanwhile, for the PCTR models, mass density
is an input, and the slab width is calculated during the solving of
statistical equilibrium. Due to these differences, we attempted to select
variables such that the slab widths and mass densities of these model types
were of similar orders of magnitudes. To investigate the effect of bulk motion
on the line profiles, we create a modified set of these base models such that
they have horizontal velocities in the range 0–80 km s-1 with a step size of
0.8 km s-1. Another set was also created with radial velocities in the same
range in order to directly compare the two effects. As DXE is affected by
altitude (Rompolt, 1967, 1980a, 1980b, 1980c), we further modified these
models to include altitudes in the range 5–70 Mm with a step size of 2 Mm.
These combinations resulted in a total of 13,332 models.

Figure 3: How DD/DB effects from horizontal (vhor) and radial velocities
(vrad) are affected by a change in height. The top row is the results from the
isothermal and isobaric atmospheres, and the bottom the PCTR atmospheres.

Figure 4: How DD/DB effects from horizontal (vhor) and radial velocities
(vrad) are affected by a change in velocity. The top row is the results from
the isothermal and isobaric atmospheres, and the bottom the PCTR atmospheres.

These models were generated on a 36-core/72-thread dual socket motherboard
machine with two Intel Xeon Gold 6140s with a base clock speed of 2.30 GHz,
and a turbo speed of 3.70 GHz. The wall time111Real world time (i.e. measured
on a clock on the wall). taken to produce the 13,332 models was approximately
70 hours with a CPU time222The total time contributed by all CPU threads. of
approximately 5,000 hours.

##  3 Results

DXE is effectively produced when a prominence undergoes some sort of bulk
motion, causing the incident radiation to appear Doppler shifted in the frame
of the prominence. For incident radiation seen in absorption, this causes a
DBE as the shifted wings of the incident line profile increase the radiation
scattered by the prominence. For incident radiation seen in emission, the
wings also illuminate the prominence, but since the incident wings are dimmer
than the incident line core, the net incident radiation is decreased and a DDE
is observed. With double peaked emission profiles, at low velocities, we
expect to see a DBE as the horns of the incident radiation are now impacting
the outgoing radiation. At higher velocities, the wings of incident line
profiles take precedence resulting in a DDE.

All 13,332 of our statistical equilibrium calculations successfully converged.
The general relationships between velocity and height for the different model
atmospheres and velocity types can be seen in Figs. 4 and 4. A value of 40 km
s-1 was chosen as the fixed velocity for Fig. 4 simply because it was the
middle of our velocity grid. A value of 11 Mm was chosen as the fixed altitude
for Fig. 4 as altitudes of the order of 10 Mm are commonly used in papers
which use PROM (e.g. Zhang et al., 2019; Levens & Labrosse, 2019; Ruan et al.,
2019; Peat et al., 2021). Due to our step size, however, the closest altitude
to 10 Mm was 11 Mm.

Fig. 4 shows how the wavelength integrated intensity of the emission changes
with altitude. The irregularities in the Ly $\alpha$ curves here are simply
due to numerical errors because of a slight misalignment between the
wavelength grids of the incident radiation and the emitted radiation from the
prominence. Unsurprisingly, as the height increases, the integrated intensity
falls – the Sun presents a smaller solid angle to the prominence. As the
height increases the radial velocity curves were expected to produce a greater
DXE effect than in the horizontal case. This is due to the lessening of the
horizontal component of the incident rays as height increases. However, Mg ii
h does not exhibit this behaviour until higher altitudes. H $\alpha$ and Ly
$\alpha$ both appear to diverge instead, with the two curves for Ly $\alpha$
in an isothermal and isobaric atmosphere appearing essentially parallel. The
integrated intensities of the lines in PCTR atmospheres are almost three times
higher than in isothermal and isobaric atmospheres for the Hydrogen lines, and
approximately only 10 per cent higher in Mg ii h.

Fig. 4 shows how the integrated intensity of the emission changes with
velocity for a fixed altitude. In H $\alpha$, for both atmosphere types,
radial and horizontal velocities produce a similar brightening effect in both
behaviour and value. For Ly $\alpha$ and Mg ii however, we see much more
complicated behaviour. Up to 30 km s-1, horizontal velocities produce a much
larger DBE in Mg ii h than the radial velocities. At higher velocities, the
horizontal velocities also produce a much greater DDE than radial. The
behaviour of Mg ii in the radial case is consistent with that found by Heinzel
et al. (2014). Ly $\alpha$ sees a similar effect, but only in the isothermal
and isobaric atmosphere with the horizontal velocity dimming at a much greater
rate than the radial velocity after its peak brightening. In the PCTR
atmosphere, there is a similar, but much less pronounced, relationship.

The asymmetry of the lines was briefly investigated. Gontikakis et al. (1997)
showed that when treated in PRD, the Ly $\alpha$ line profiles from moving
solar prominences are asymmetric. This is not possible under CRD, and another
reason why we treated all of the resonance lines in PRD. Here, the asymmetry
was calculated via the quantile method (Kerr et al., 2015). This works by
calculating the cumulative distribution function (CDF) of a line profile, and
normalising them such that the values in the function range from 0–1. The
wavelength at which the CDF is 0.5 ($\lambda_{50}$) is defined as the line
core. If the line is of a Gaussian shape, the difference between
$\lambda_{88}$ and $\lambda_{12}$ is the full width half maximum. From this,
we calculate the asymmetry by the following relationship,

$\text{Asymmetry}=\frac{\left(\lambda_{88}-\lambda_{50}\right)-\left(\lambda_{50}-\lambda_{12}\right)}{\lambda_{88}-\lambda_{12}}=\frac{\lambda_{88}-2\lambda_{50}+\lambda_{12}}{\lambda_{88}-\lambda_{12}},$
(3)

where a positive asymmetry means there is more emission in the blue, and a
negative asymmetry means there is more emission in the red. The asymmetry in
all of the line profiles is not appreciably affected in either the horizontal
or radial velocity modes, never exceeding $\pm$0.02. However, in observations,
asymmetry is usually measured on the order of $\pm$0.5 (Ruan et al., 2018;
Peat et al., 2021). Therefore, this effect would be negligible in observations
and cannot be a source of the asymmetry on its own.

##  4 Discussion

The way in which DXE occurs in the radial and horizontal velocity cases are
subtly different. In the much studied radial case, as the prominence moves
away from the solar surface, the incident radiation appears red-shifted in the
frame of the prominence, and if it is moving towards the the solar surface,
the radiation appears blue-shifted. However, for the horizontal case, the
leading edge of the prominence sees the incident radiation as blue-shifted,
while the trailing edge sees the incident radiation as red-shifted. The radial
case is only affected by one side of the wings of the incident radiation,
whereas the horizontal case is affected by both sides of the wings.

The shape and intensity of the incident radiation profiles is extremely
important to the total intensity of the line profiles produced by a prominence
whether it is undergoing some bulk motion or not. Here we used the Ly
$\alpha$, H $\alpha$, and Mg ii h emission from the aforementioned FALC
atmospheric model as the incident radiation. Whereas other studies modelling
the emission from prominences use observations of the lines in question (e.g.
Heinzel & Rompolt, 1987; Levens & Labrosse, 2019). However, if the incident
radiation were to change from FALC to that found in observations, the
qualitative relationships between the intensities and velocities would remain
the same.

Much of the focus of DXE has been contextualised in the form of eruptive
prominences moving predominantly radially (e.g. Heinzel & Rompolt, 1987; Vial
& Engvold, 2015; Zhang et al., 2019), but this work shows that horizontal
velocities also play a role. In observation, eruptive prominences are never
going to erupt perfectly perpendicular to the solar surface. This will
introduce a small horizontal component to its velocity, which is when the
horizontal effect is at its greatest for lines with incident radiation in
emission (see Fig. 4). Fig. 4 shows that radial and horizontal DXE appear to
dissipate with height at the same rate, illustrating that this horizontal
component will continue to be important as the prominence moves away from the
solar disc. Furthermore, it is likely that the solar spectra incident on an
erupting prominence are no longer symmetrical, increasing the complexity of
the DXE effect. If it is assumed that an erupting prominence only has radial
velocity, the radial velocity will be overestimated to account for the DXE
introduced by the horizontal effect.

##  5 Conclusions

In this study we demonstrated an application of the NLTE Lightweaver RT
framework, in the context of Promweaver, to explore DXE caused by different
orthogonal motions. To meet this goal, we generated 13,332 models of differing
parameters.

We showed the importance of all translational bulk motions of solar
prominences on the formation of observed line profiles. Past studies mainly
focused on radial velocities to better understand erupting prominences, but
neglected horizontal velocities entirely. This neglect may have caused derived
radial velocities of eruptive prominences from forward modelling to be
overestimated.

In both the isboaric/isothermal and PCTR atmospheres, Ly $\alpha$ and Mg ii h,
as velocity increases, the horizontal case produces a much stronger DXE effect
than its radial counterpart, but H $\alpha$ sees the greater effect caused by
the radial case. With respect to increasing altitude, the intgrated
intensities in the radial and horizontal cases both appear to diminish at
equivalent rates in the two hydrogen lines. However, with respect to Mg ii h,
the radial and horizontal cases converge. Perhaps at larger velocities, Ly
$\alpha$ and H $\alpha$ will also display this behaviour. In general, the
magnitude of the DXE produced by the horizontal case is comparable to that
produced in the radial case. In future studies, care should be taken when
forward modelling to account for effects introduced by horizontal velocities
and not to erroneously attribute it to only radial velocities.

## Acknowledgements

AWP and PH acknowledge financial support from the ‘Excellence
Initiative–Research University’ program for the years 2020–2026 at the
University of Wrocław, project No. BPIDUB.4610.96.2021.KG. AWP would also like
to thank Dr. Nicolas Labrosse for having the time for many impromptu
conversations.

CMJO is grateful for the support of the Royal Astronomical Society’s Norman
Lockyer fellowship, and the support of the University of Glasgow’s College of
Science and Engineering and Lord Kelvin/Adam Smith fellowship.

PH was also supported by the grant No. 22-34841S of the Czech Funding Agency.
He would like to express his deep gratitude to Prof. Bogdan Rompolt from the
University of Wrocław for introducing him to problems of DXE.

This research used version 3.8.2 of Matplotlib (Hunter, 2007), version 1.24.4
of NumPy (Harris et al., 2020), and version 5.3.3 of Astropy
(http://www.astropy.org) a community-developed core Python package for
Astronomy (The Astropy Collaboration et al., 2013, 2018).

## Data Availability

This research presented an application of version 0.1.1 of Promweaver
available at https://doi.org/10.5281/zenodo.6546897. The most recent version
of Promweaver is available to download via pip or at
https://github.com/Goobley/Promweaver.
All of the models used in this study can be found in the form of a Python
pickle at https://doi.org/10.5281/zenodo.12083238.

## References

  * Anzer & Heinzel (1999) Anzer U., Heinzel P., 1999, Astronomy & Astrophysics, 349, 974
  * Fontenla et al. (1993) Fontenla J. M., Avrett E. H., Loeser R., 1993, ApJ, 406, 319
  * Gontikakis et al. (1997) Gontikakis C., Vial J. C., Gouttebroze P., 1997, A&A, 325, 803
  * Gouttebroze (2004) Gouttebroze P., 2004, Astronomy & Astrophysics, 413, 733
  * Gouttebroze et al. (1993) Gouttebroze P., Heinzel P., Vial J.-C., 1993, Astronomy & Astrophysics, 99
  * Hale (1929) Hale G. E., 1929, ApJ, 70, 265
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Heasley & Mihalas (1976) Heasley J. N., Mihalas D., 1976, ApJ, 205, 273
  * Heasley et al. (1974) Heasley J. N., Mihalas D., Poland A. I., 1974, The Astrophysical Journal, 192, 181
  * Heinzel (1995) Heinzel P., 1995, A&A, 299, 563
  * Heinzel & Rompolt (1987) Heinzel P., Rompolt B., 1987, Sol. Phys., 110, 171
  * Heinzel et al. (1987) Heinzel P., Gouttebroze P., Vial J. C., 1987, Astronomy & Astrophysics, 183, 351
  * Heinzel et al. (1994) Heinzel P., Gouttebroze P., Vial J.-C., 1994, Astronomy & Astrophysics, 292, 656
  * Heinzel et al. (2014) Heinzel P., Vial J.-C., Anzer U., 2014, Astronomy & Astrophysics, 564, A132
  * Hunter (2007) Hunter J. D., 2007, Comp. Sci, Eng., 9, 90
  * Hyder & Lites (1970) Hyder C., Lites B., 1970, Solar Physics, 14
  * Ishizawa (1971) Ishizawa T., 1971, PASJ, 23, 75
  * Jenkins et al. (2023) Jenkins J. M., Osborne C. M. J., Keppens R., 2023, A&A, 670, A179
  * Kerr et al. (2015) Kerr G. S., Simões P. J. A., Qiu J., Fletcher L., 2015, Astronomy & Astrophysics, 582, A50
  * Kippenhahn & Schlüter (1957) Kippenhahn R., Schlüter A., 1957, Z. Astrophys., 43, 36
  * Labrosse & Gouttebroze (2004) Labrosse N., Gouttebroze P., 2004, The Astrophysical Journal, 617, 614
  * Labrosse et al. (2010) Labrosse N., Heinzel P., Vial J.-C., Kucera T., Parenti S., Gunár S., Schmieder B., Kilper G., 2010, Space Science Reviews, 151, 243
  * Leenaarts et al. (2013) Leenaarts J., Pereira T. M. D., Carlsson M., Uitenbroek H., De Pontieu B., 2013, The Astrophysical Journal, 772, 89
  * Levens & Labrosse (2019) Levens P. J., Labrosse N., 2019, Astronomy & Astrophysics, 625, A30
  * Morimoto & Kurokawa (2003) Morimoto T., Kurokawa H., 2003, PASJ, 55, 503
  * Osborne (2022) Osborne C. M. J., 2022, Goobley/Promweaver: v0.1.1, doi:10.5281/zenodo.6546897, https://doi.org/10.5281/zenodo.6546897
  * Osborne & Milić (2021) Osborne C. M. J., Milić I., 2021, The Astrophysical Journal, 917, 14
  * Paletou et al. (1993) Paletou F., Vial J. C., Auer L. H., 1993, Astronomy & Astrophysics, 274, 571
  * Peat (2023) Peat A. W., 2023, PhD thesis, University of Glasgow, UK
  * Peat et al. (2021) Peat A. W., Labrosse N., Schmieder B., Barczynski K., 2021, Astronomy & Astrophysics, 653, A5
  * Rompolt (1967) Rompolt B., 1967, Acta Astron., 17, 329
  * Rompolt (1980a) Rompolt B., 1980a, Hvar Observatory Bulletin, 4, 39
  * Rompolt (1980b) Rompolt B., 1980b, Hvar Observatory Bulletin, 4, 49
  * Rompolt (1980c) Rompolt B., 1980c, Hvar Observatory Bulletin, 4, 55
  * Ruan et al. (2018) Ruan G., Schmieder B., Mein P., Mein N., Labrosse N., Gunár S., Chen Y., 2018, The Astrophysical Journal, 865, 123
  * Ruan et al. (2019) Ruan G., Jejčič S., Schmieder B., Mein P., Mein N., Heinzel P., Gunár S., Chen Y., 2019, The Astrophysical Journal, 886, 134
  * Secchi (1875) Secchi A., 1875, Le Soleil. Gauthier-Villars, doi:10.3931/e-rara-14748
  * Tandberg-Hanssen (1995) Tandberg-Hanssen E., 1995, The Nature of Solar Prominences. Astrophysics and Space Science Library Vol. 199, Springer Netherlands, Dordrecht, doi:10.1007/978-94-017-3396-0, http://link.springer.com/10.1007/978-94-017-3396-0
  * The Astropy Collaboration et al. (2013) The Astropy Collaboration et al., 2013, Astronomy & Astrophysics, 558, A33
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration et al., 2018, The Astronomical Journal, 156, 123
  * Vassenius (1733) Vassenius B., 1733, Philosophical Transactions of the Royal Society of London, 38, 134
  * Vial & Engvold (2015) Vial J.-C., Engvold O., eds, 2015, Solar Prominences. Astrophysics and Space Science Library Vol. 415, Springer International Publishing, Cham, doi:10.1007/978-3-319-10416-4, http://link.springer.com/10.1007/978-3-319-10416-4
  * Vyssotsky (1949) Vyssotsky A. N., 1949, Meddelanden fran Lunds Astronomiska Observatorium Serie II, 126, 3
  * Zhang et al. (2019) Zhang P., Buchlin É., Vial J.-C., 2019, Astronomy & Astrophysics, 624, A72