# Learned Thresholds Token Merging and Pruning for Vision Transformers

Maxim Bonnaerens<EMAIL_ADDRESS>
IDLab-AIRO, Ghent University - imec Joni Dambre<EMAIL_ADDRESS>
IDLab-AIRO, Ghent University - imec

###### Abstract

Vision transformers have demonstrated remarkable success in a wide range of
computer vision tasks over the last years. However, their high computational
costs remain a significant barrier to their practical deployment. In
particular, the complexity of transformer models is quadratic with respect to
the number of input tokens. Therefore techniques that reduce the number of
input tokens that need to be processed have been proposed. This paper
introduces Learned Thresholds token Merging and Pruning (LTMP), a novel
approach that leverages the strengths of both token merging and token pruning.
LTMP uses learned threshold masking modules that dynamically determine which
tokens to merge and which to prune. We demonstrate our approach with extensive
experiments on vision transformers on the ImageNet classification task. Our
results demonstrate that LTMP achieves state-of-the-art accuracy across
reduction rates while requiring only a single fine-tuning epoch, which is an
order of magnitude faster than previous methods. Code is available at
https://github.com/Mxbonn/ltmp.

input image | patchified image | layer 1 | layer 2 | layer 3 | layer 4 | layer 5
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
layer 6 | layer 7 | layer 8 | layer 9 | layer 10 | layer 11 | layer 12
Figure 1: Visualization of the merging and pruning as applied to image
patches. In every layer the most similar tokens are merged and any unimportant
tokens are pruned. The visualizations show the remaining tokens after every
layer in DeiT-S.

##  1 Introduction

The adoption of transformers (Vaswani et al., 2017), originally developed for
natural language processing, in the field of computer vision with Vision
Transformers (ViT) (Dosovitskiy et al., 2021) has led to significant progress
in the field. But despite the impressive results of vision transformers, their
success has come with a cost; ViT models are computationally expensive and
require larger datasets (e.g. ImageNet-21k instead of ImageNet-1k (Deng et
al., 2009)) and prolonged training times (Dosovitskiy et al., 2021). In order
to benefit from their accuracy in downstream tasks and applications, the use
of pretrained models has become essential. However, their adoption on
resource-constrained devices such as mobile and embedded platforms remains
limited due to the high computational cost.

To reduce the computational cost of transformers, previous work has focused on
techniques such as distillation (Wu et al., 2022), quantization (Liu et al.,
2021b) and pruning. Pruning techniques have explored pruning model weights
(Gordon et al., 2020), attention heads (Voita et al., 2019), and input tokens
(Rao et al., 2021). This last approach is very effective, as the model
complexity of a transformer is quadratic to the number of tokens. In vision
transformers, the tokens are non-overlapping patches of an image, e.g. a token
may represent patches of $16\times 16$ pixels. Token pruning has attracted
research interest as it matches our intuition that not all parts of an image
are equally important. The self-attention mechanism in transformers, due to
its ability to process variable length inputs and its order-agnostic
characteristic, enables unstructured reduction of the number of tokens between
layers. This was previously non-trivial with convolutional networks.

Most token pruning approaches calculate an importance score for every token in
each layer and remove the least important tokens. While token pruning has been
shown to be an effective compression technique, removing tokens results in
loss of information which limits the amount of tokens that can be pruned. In
order to recover from this loss of information, most pruning approaches
require substantial retraining to be effective. Additionally, some recent
pruning techniques have incorporated token combining techniques where the
pruned tokens are combined into a single token that aggregates the information
that would otherwise be lost (Kong et al., 2022).

Token merging (ToMe) (Bolya et al., 2023) takes this combining technique one
step further. ToMe exclusively combines pairs of tokens into new tokens rather
than pruning them. This has as advantage that it does not discard but
summarizes information, leading to better accuracy while being equally
effective in reducing computational complexity.

In this work, we introduce Learned Thresholds token Merging and Pruning
(LTMP). Our approach combines the benefits of token merging, which allows us
to combine rather than discard token information, with pruning, which allows
us to remove uninformative tokens. To the best of our knowledge, this is the
first work to extensively combine these two reduction techniques, which leads
to improved accuracy compared to previous work. Our approach uses learned
threshold masking modules, which allow the model to learn thresholds that
determine which tokens to prune and which ones to merge in each layer. This
enables adaptive token reduction, while only requiring two learnable
parameters per transformer block. As a result, our approach converges within a
single epoch, reducing the fine-tuning cost by an order of magnitude compared
to other learnable pruning approaches.

Our contributions can be summarized as follows:

  * •

We propose to combine token merging with token pruning, enabling us to achieve
high token reduction rates with minimal loss of accuracy.

  * •

Our method introduces learned threshold masking modules, which require only
two learnable parameters per transformer block, allowing our approach to
converge within a single epoch.

  * •

We optimize the thresholds using a novel budget-aware training loss for which
we introduce a reduction target $r_{target}$ and an actual FLOPs reduction
factor $r_{FLOPs}$. This allows us to create models of any size and allows the
model to freely distribute the reduction operations across layers.

##  2 Related work

###  2.1 Efficient Vision Transformers

Initially, transformers (Vaswani et al., 2017) were adopted from NLP to
computer vision (Dosovitskiy et al., 2021) for their impressive accuracy. But
despite their success in many vision tasks, ViT-based models could not compete
with lightweight CNNs for deployment on resource-constrained platforms (Wang
et al., 2022). To create efficient ViT models, several architecture changes
have been proposed which modify the expensive attention mechanisms (Kitaev et
al., 2019; Chen et al., 2021; Liu et al., 2021a; Li et al., 2021). In this
paper, we look at token pruning, but other pruning approaches that have been
successfully applied to transformers are, among others, weight pruning (Gordon
et al., 2020) and attention heads pruning (Voita et al., 2019).

###  2.2 Token Pruning

The flexibility of transformers with respect to the sequence length and order
of the inputs allows token pruning, something which was previously non-trivial
to do in convolutional-based models. Token pruning methods can differ in
various ways, such as the score used to determine the importance of each
token. Pruning methods can also differ in the way that token reduction is
applied. In fixed rate pruning (Goyal et al., 2020; Rao et al., 2021; Bolya et
al., 2023; Liang et al., 2022; Xu et al., 2022) a predefined number of tokens
is removed per layer, while in adaptive approaches (Kim et al., 2022; Yin et
al., 2021; Liu et al., 2022) the tokens are pruned dynamically based on the
input.

The most recent pruning approaches do not only prune tokens but they also
create a single additional token after each pruning step. This token
aggregates the information of the pruned tokens and limits the loss of
accuracy while pruning. EViT (Liang et al., 2022), Evo-ViT (Xu et al., 2022)
and SPViT (Kong et al., 2022) use a weighted average based on the importance
score to create the new fused token.

###  2.3 Token merging

Token Merging (ToMe), as introduced by Bolya et al. (2023), introduces a
lightweight token matching algorithm to merge similar tokens. It is as fast as
pruning while being more accurate. ToMe partitions all tokens into two sets
$\mathbb{A}$ and $\mathbb{B}$ of roughly equal size by alternating and
calculates similarity scores for each token in $\mathbb{A}$ with every token
in $\mathbb{B}$. The similarity score is defined as the cosine similarity
between the key vectors ($\mathbf{K}$) used in the self-attention layer. The
final similarity score of a token in $\mathbb{A}$ is the highest similarity
score with any token in $\mathbb{B}$. Based on this score, ToMe merges the $k$
most similar tokens through averaging and concatenates the two sets back
together.

##  3 Learned thresholds token merging and pruning

###  3.1 Overview

An overview of our framework is shown in Figure 2. Given any vision
transformer, our approach adds merging (LTM) and pruning (LTP) components with
learned threshold masking modules in each transformer block between the Multi-
head Self-Attention (MSA) and MLP components. Based on the attention in the
MSA, importance scores for each token and similarity scores between tokens are
computed. Learned threshold masking modules then learn the thresholds that
decide which tokens to prune and which ones to merge.

Figure 2: Overview of our approach. LTMP contains a merging and pruning
component, each with a learned threshold masking module. The components are
added between the Multi-head Self-Attention and MLP components of each
transformer block.

###  3.2 Motivation

Figure 3: Visualization of token pruning (middle) compared to token merging
(right). The visualization shows the remaining tokens after the 10-th layer in
ViT-B when pruning or merging 16 tokens per layer.

Although token merging is generally more accurate than pruning as it combines
tokens instead of discarding them, it is not always better to merge tokens
instead of discarding them. In some cases, it may be more beneficial to prune
an unimportant token rather than merging the most similar tokens, as the
similarity between them may not be very high.

In this section, we explore whether token merging and token pruning are
techniques that can be combined. Figure 3 visualizes tokens kept by pruning
and by merging on one specific image, we observe that the kept tokens are
noticeably different between both approaches. To quantify the relation between
merging and pruning, we calculated the Kendall tau rank correlation (Kendall,
1938) between the _importance scores_ used in token pruning and the
_similarity scores_ used in token merging. However, as merging only calculates
similarity scores for tokens in set $\mathbb{A}$ (see Section 2.3), we compute
the Kendall tau correlation between the similarity scores and the importance
scores of only the tokens in set $\mathbb{A}$. We calculated the correlations
over 1000 images using a ViT-B where each layer pruned and merged 8 tokens in
a fixed top-$k$ manner and report the results in Table 1. We find that the
$\tau$ correlations between both scores are low, especially in the early
layers where the most important merging and pruning is done. We, therefore,
propose combining token merging and token pruning.

Table 1: Kendall $\tau$ correlation between _importance scores_ in token pruning and _similarity scores_ in token merging. The correlations are calculated over 1000 images using a ViT-B where each layer pruned and merged 8 tokens. layer | 1 | 3 | 5 | 7 | 9 | 11
---|---|---|---|---|---|---
$\tau$ | 0.01 | 0.14 | 0.21 | 0.19 | 0.15 | 0.22

###  3.3 Learned thresholds

####  3.3.1 Learned thresholds pruning

Our learned thresholds approach is conceptually similar to learned token
pruning as introduced in Kim et al. (2022). In each transformer block an
_importance score_ is calculated for every token
$\mathbf{x}_{i},i\in\\{1,...,n\\}$, where $n=hw$ is the number of tokens111We
omit the [CLS] class token for simplicity, during pruning and/or merging we
always keep the [CLS] token.. A threshold
$\theta^{l}\in\mathbb{R},l\in\\{1,...,L\\}$, where $L$ is the number of
transformer blocks, determines which tokens to keep and which to prune in each
layer; only tokens with an importance score above the threshold are kept.

In order to prune tokens adaptively, we introduce a threshold masking module
that, given the importance scores $\mathbf{s}^{l}\in\mathbb{R}^{n}$, learns a
pruning threshold $\theta^{l}$ and outputs which tokens to keep.

$M(\mathbf{s}^{l}_{i},\theta^{l})=\begin{cases}1,&\text{if
}\mathbf{s}^{l}_{i}>\theta^{l}\\\ 0,&\text{otherwise}\end{cases}$ (1)

However, in order to make $\theta^{l}$ learnable during training, the
threshold masking module needs to be differentiable. We achieve this by
implementing the threshold masking module as a straight-through estimator
(Bengio et al., 2013), where we estimate the masking function during
backpropagation as

$M(\mathbf{s}^{l}_{i},\theta^{l})=\sigma(\frac{\mathbf{s}^{l}_{i}-\theta^{l}}{\tau})$
(2)

where $\sigma(x)$ is the sigmoid function and $\tau$ is the temperature
hyperparameter.

During inference we only keep the tokens in the $l$-th block where
$M(\mathbf{s}^{l}_{i},\theta^{l})=1$. However, during training, we can not
simply drop tokens as that does not allow the model to backpropagate the
influence of the threshold on the model performance. We, therefore, create a
mask indicating which tokens are kept and which ones are pruned. Every
threshold masking module only updates the entries of the mask for the tokens
that have not yet been removed prior to that layer, as tokens that are pruned
in an earlier layer have to remain pruned. We construct the pruning mask
$\mathbf{m}^{l}\in[0,1]^{n}$ as follows:

$\mathbf{m}^{l}_{i}=\begin{cases}M(\mathbf{s}^{l}_{i},\theta^{l}),&\text{if
}\mathbf{m}^{l-1}_{i}=1\\\ \mathbf{m}^{l-1},&\text{otherwise}\end{cases}$ (3)

The learned token pruning implementation in Kim et al. (2022) multiplies its
mask with the tokens in order to create zero valued tokens. However, these
tokens do not remain zero due to bias terms in MLP layers; furthermore adding
zero valued tokens changes attention calculations compared to removing those
tokens. Instead, our approach makes changes to the only place where tokens
influence each other: the attention mechanism. 222Technically, tokens also
influence each other during layer normalization, however as pruning is done on
pretrained models, we simply use the global statistics from pretraining during
normalization.

Recall the original formula for attention (Vaswani et al., 2017):

$\operatorname{Attention}(\mathbf{Q},\mathbf{K},\mathbf{V})=\operatorname{softmax}(\frac{\mathbf{QK}^{T}}{\sqrt{d_{k}}})\mathbf{V}$
(4)

In order to avoid that the masked tokens influence the attention mechanism, we
propose a modified function:

$\operatorname{Attention\\_with\\_mask}(\mathbf{Q},\mathbf{K},\mathbf{V},\mathbf{m})=\mathbf{S}\mathbf{V}$
(5)

where,

$\mathbf{S}_{ij}=\frac{\exp(\mathbf{A}_{ij})\mathbf{m}_{j}}{\sum_{k=1}^{N}\exp(\mathbf{A}_{ik})\mathbf{m}_{k}},1\leq
i,j,k\leq n$ (6)

and,

$\mathbf{A}=\mathbf{Q}\mathbf{K}^{T}/\sqrt{d_{k}}\in\mathbb{R}^{n\times n}$
(7)

Equation 6 computes a masked softmax, which is equivalent to a softmax
calculated with the pruned tokens removed.
$\operatorname{Attention\\_with\\_mask}$ is conceptually similar to the masked
attention as found in the transformer decoder of language models. However,
where the masking in transformer decoders is done by setting masked tokens to
$-\infty$, our approach requires the influence of the straight-through
estimator mask to propagate to the thresholds during backpropagation.

####  3.3.2 Learned thresholds merging

Token merging is originally a top-$k$ approach, meaning that it merges based
on a fixed rate and has no learnable parameters. We modify ToMe to use
thresholds instead of top-$k$ by applying the same techniques as introduced in
Section 3.3.1; this is by adding our learned threshold masking module, in
which similarity scores above these thresholds are selected for merging, and
by changing the attention function to Equation 5.

####  3.3.3 Learned thresholds merging and pruning

With learnable thresholds, it is trivial to combine merging and pruning, as we
can simply add a learned threshold masking module that learns thresholds for
importance scores and another module that learns thresholds for similarity
scores.

###  3.4 Training Strategy

####  3.4.1 Training objective

To effectively reduce the number of tokens in the transformer blocks, it is
necessary to include a regularization loss term in the training process.
Without this loss, the model has no incentive to prune any tokens and the
pruning thresholds will simply be set to $0$ as the most accurate model uses
all inputs. We propose a budget-aware training loss which introduces a
reduction target $r_{\text{target}}$ for the FLOPs of the vision transformer.

Let us denote $\phi_{\text{module}}(n,d)$ as a function that calculates the
FLOPs of a module based on the number of tokens $n$ and the embedding
dimension $d$. The actual FLOPs reduction factor $r_{\text{FLOPs}}$ of a ViT
can then be computed as:

$r_{\text{FLOPs}}=\frac{\phi_{\text{PE}}(n,d)}{\phi_{\text{ViT}}(n,d)}+\sum_{l=1}^{L}\frac{\phi_{\text{BLK}}(n,d)}{\phi_{\text{ViT}}(n,d)}\left(\frac{\phi_{\text{MSA}}(\bar{\mathbf{m}}^{l-1}n,d)}{\phi_{\text{BLK}}(n,d)}+\frac{\phi_{\text{MLP}}(\bar{\mathbf{m}}^{l}n,d)}{\phi_{\text{BLK}}(n,d)}\right)+\frac{\phi_{\text{HEAD}}(\bar{\mathbf{m}}^{l}n,d)}{\phi_{\text{ViT}}(n,d)}$
(8)

where $\bar{\mathbf{m}}^{l}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{m}^{l}_{i}$ is
the percentage of input tokens that are kept after the $l$-th threshold
masking operation and $\bar{\mathbf{m}}^{0}=1$. PE, BLK and HEAD denote the
different components of a vision transformer: the patch embedding module, the
transformer blocks and the classification head.

As the vast majority of the FLOPs in a vision transformer occurs in the
transformer blocks ($\approx 99\%$ in ViT-S), we ignore the FLOPs in the patch
embedding and classification head:
$\frac{\phi_{\text{PE}}(n,d)}{\phi_{\text{ViT}}(n,d)}=\frac{\phi_{\text{HEAD}}(n,d)}{\phi_{\text{ViT}}(n,d)}\approx
0$. That means that we can simplify

$\frac{\phi_{\text{BLK}}(n,d)}{\phi_{\text{ViT}}(n,d)}\approx\frac{1}{L},$ (9)

where $L$ is the number of transformer blocks.

The FLOPs of a transformer block and its two components, the MSA and MLP can
be computed as:

$\phi_{\text{MSA}}(n,d)=4nd^{2}+2n^{2}d$ (10) $\phi_{\text{MLP}}(n,d)=8nd^{2}$
(11)
$\phi_{\text{BLK}}(n,d)=\phi_{\text{MSA}}(n,d)+\phi_{\text{MLP}}(n,d)=12nd^{2}+2n^{2}d$
(12)

Substituting Equations 10, 11 and 12 into Equation 8 gives:

$r_{\text{FLOPs}}\approx{}\sum_{l=1}^{L}\frac{1}{L}\left(\frac{2\bar{\mathbf{m}}^{l-1}nd^{2}+(\bar{\mathbf{m}}^{l-1}n)^{2}d+4\bar{\mathbf{m}}^{l}nd^{2}}{6nd^{2}+n^{2}d}\right)$
(13)

Given this FLOPs reduction factor $r_{\text{FLOPs}}$ as a function of the
threshold masks, we define our regularization loss as the squared error
between the reduction target and the actual FLOPs reduction factor:

$\mathcal{L}_{\text{reg}}=(r_{\text{target}}-r_{\text{FLOPs}})^{2}$ (14)

This regularization loss is then combined with the classification loss, for
which we adopt the standard cross entropy loss.

$\mathcal{L}=\mathcal{L}_{\text{CE}}+\lambda\mathcal{L}_{\text{reg}}$ (15)

The overall training objective is to learn thresholds that optimize the model
while reducing the model complexity to a certain reduction target. The
combination of learned thresholds and our budget-aware loss enables the model
to optimally distribute merging and pruning across layers.

####  3.4.2 Training schedule

LTMP only adds two learnable parameters per transformer block (one for pruning
and one for merging). As is common in pruning it is applied to pretrained
models. We therefore only update the thresholds during training and keep all
other trainable parameters fixed, allowing LTMP to converge within a single
epoch.

Figure 4: Accuracy/FLOPs trade-offs for LTMP variants of DeiT-Tiny, -Small and
-Base. FLOPs are plotted in logarithmic scale. ToMe accuracy/FLOPs trade-offs
are shown as a comparison.

##  4 Experiments

In this section, we demonstrate our approach through extensive experiments on
the ImageNet-1k classification task (Deng et al., 2009) using various ViT
variants (Steiner et al., 2022; Touvron et al., 2021). All pretrained models
are taken from the timm PyTorch library (Wightman, 2022).

All our experiments are trained for a single epoch, using SGD without momentum
and a batch size of 128. The remaining training settings such as augmentations
are set to the default values of timm. The hyperparameters that are introduced
in LTMP are set to $\tau=0.1$ and $\lambda=10$. As the importance scores and
similarity scores have values in different ranges we use separate learning
rates for the thresholds in the pruning modules and the merging modules:
$5\cdot 10^{-6}$ for the pruning thresholds and $5\cdot 10^{-3}$ for the
merging thresholds. To select these hyperparameters, we used a separate held-
out validation set containing 2% of the original ImageNet training set. The
final hyperparameters are determined based on their performance on DeiT-Small
with $r_{target}$ set to $0.65$. To show the robustness of these
hyperparameters, we did not change them between model variations and
reductions targets (i.e. ViT-Tiny with $r_{target}$ set to $0.45$, is fine-
tuned with the same hyperparameters).

###  4.1 Main results

As most other pruning approaches require extensive fine-tuning, the most
commonly used baseline vision transformer in other works is the data-efficient
vision transformer DeiT (Touvron et al., 2021), for this reason, we also
report on DeiT models.

In Figure 4, we show our approach applied to DeiT-Tiny, -Small and -Base. For
each model we vary $r_{target}$ such that we obtain a set of reduced models,
each with a different model size. Table 2 lists the detailed results for the
DeiT models. Models used in one of the other tables in the paper are
highlighted in gray. We included $r_{\text{target}}$, such that our results
can be reproduced. In Appendix A, we additionally included results on the
standard ViT models. All accuracy numbers are reported with 3 significant
figures; a brief analysis of the error ranges on the accuracy numbers can be
found in Appendix B.

Table 2: Detailed results for DeiT models.

Model | $r_{\text{target}}$ | FLOPs | Accuracy
---|---|---|---
DeiT-T (Baseline) | 1.258G | 72.2
DeiT-T | 0.85 | 1.091G | 72.1
DeiT-T | 0.80 | 1.029G | 72.0
DeiT-T | 0.75 | 0.974G | 72.0
DeiT-T | 0.70 | 0.918G | 71.8
DeiT-T | 0.65 | 0.866G | 71.7
DeiT-T | 0.60 | 0.813G | 71.5
DeiT-T | 0.55 | 0.752G | 71.0
DeiT-T | 0.50 | 0.703G | 70.5
DeiT-T | 0.45 | 0.652G | 69.8

Model | $r_{\text{target}}$ | FLOPs | Accuracy
---|---|---|---
DeiT-S (Baseline) | 4.608G | 79.9
DeiT-S | 0.85 | 3.940G | 79.8
DeiT-S | 0.80 | 3.753G | 79.8
DeiT-S | 0.75 | 3.518G | 79.8
DeiT-S | 0.70 | 3.326G | 79.7
DeiT-S | 0.65 | 3.155G | 79.7
DeiT-S | 0.62 | 3.016G | 79.6
DeiT-S | 0.60 | 2.955G | 79.5
DeiT-S | 0.55 | 2.722G | 79.3
DeiT-S | 0.50 | 2.523G | 79.0
DeiT-S | 0.45 | 2.314G | 78.6

Model | $r_{\text{target}}$ | FLOPs | Accuracy
---|---|---|---
DeiT-B (Baseline) | 17.583G | 81.99
DeiT-B | 0.85 | 15.007G | 81.9
DeiT-B | 0.80 | 14.152G | 81.9
DeiT-B | 0.75 | 13.403G | 81.7
DeiT-B | 0.70 | 12.571G | 81.6
DeiT-B | 0.65 | 11.754G | 81.3
DeiT-B | 0.60 | 11.022G | 81.0
DeiT-B | 0.55 | 10.276G | 80.3
DeiT-B | 0.50 | 9.608G | 79.4
DeiT-B | 0.45 | 9.003G | 78.8

###  4.2 Comparison to other work

In Table 3, we compare the reported top-1 accuracy, FLOPs333We follow the
convention of reporting FLOPs as multiply-adds., and fine-tune epochs of other
pruning approaches to our work. The other approaches we compare with are SPViT
(Kong et al., 2022), DynamicVit (Rao et al., 2021), EViT (Liang et al., 2022),
EvoViT (Xu et al., 2022) and ToMe (Bolya et al., 2023). Most works report on a
pruned model with around $3.0$ GFLOPs, which for DeiT-S corresponds to
$r_{target}\approx 0.65$ in our approach. Only SPViT reports on a different
model size, which is why it is compared separately. The results show that LTMP
reaches state-of-the-art accuracy at a fraction of the fine-tuning epochs
required by other learnable methods. The accuracy of LTMP matches or exceeds
the accuracy of other token pruning approaches which require a minimum of 30
fine-tune epochs. Only EViT is able to reach a higher accuracy than LTMP, but
only when drastically increasing the fine-tune epochs to 100, which is two
orders of magnitude more than our approach. Because ToMe requires no fine-
tuned checkpoints for comparisons, we are able to compare LTMP to ToMe more
extensively over a wide range of model sizes. Figure 4 shows the
accuracy/FLOPs trade-offs for LTMP and ToMe. Our experiments show that LTMP
consistently outperforms ToMe across model sizes.

Table 3: Comparison to other token reduction approaches on DeiT-S. Our method reaches state-of-the-art accuracy with significantly fewer fine-tune epochs than other learnable approaches. Method | FLOPs | Accuracy |  fine-tune epochs
---|---|---|---
DeiT-S (Baseline) | 4.6G | 79.8 | -
SPViT | 3.8G | 79.8 | 75
LTMP (Ours) | 3.8G | 79.8 | 1
DynamicViT | 2.9G | 79.3 | 30
EViT | 3.0G | 79.5 | 30
EViT | 3.0G | 79.8 | 100
Evo-ViT | 3.0G | 79.4 | 300
ToMe | 3.0G | 79.3 | 0
LTMP (Ours) | 3.0G | 79.6 | 1

###  4.3 Synergy between merging and pruning

In order to analyze the synergy between token merging and token pruning, we
have examined the distribution of merged and pruned tokens across each layer
of the vision transformer. The results, shown in Figure 5 for DeiT-S, reveal
that token merging is the dominant reduction operation in the early layers of
the transformer, while token pruning is more prevalent in the final layers.
This aligns with the findings of Bolya et al. (2023) that merging is more
effective than pruning, as it can summarize information. However, once all
similar tokens are merged, it is more beneficial to prune the least
informative tokens rather than merging tokens that are not as similar. In
other words, the first layers are mainly used to combine similar patches and
once this is mostly done, token pruning removes the unimportant parts of the
input.

Figure 5: Distribution of the number of tokens ($k$) removed by the merging
and pruning parts in each layer of LTMP DeiT-S $r_{FLOPs}\approx 0.5$.

###  4.4 Design Choices

####  4.4.1 Importance score

A key component of pruning approaches is the _importance score_ used to
determine which tokens to remove. The two most common choices for the
importance score are:

$\text{class attention score }s_{i}=\sum_{j=1}^{h}\mathbf{S}_{j0i}$ (16)

where $S\in\mathbb{R}^{h\times n\times n}$ is the multi-headed extension of
the attention softmax matrix (see Equation 6) where values at index $0$
correspond to the [CLS] token, and

$\text{mean column attention score }s_{i}=\frac{1}{h\cdot
n}\sum_{j=1}^{h}\sum_{k=1}^{n}\mathbf{S}_{jki}$ (17)

which can be interpreted as the normalized amount that all tokens
$\mathbf{x}_{k}$ attend to token $\mathbf{x}_{i}$ (Kim et al., 2022).

The results in Table 4 are obtained from DeiT-S LTP models and show that the
mean column attention score performs slightly better than the class attention
score, but not by a significant margin. For the remainder of this work, we use
the mean column attention score (Equation 17) as the importance score.

####  4.4.2 Merging and pruning order

Table 4: Ablation of the design choices regarding the pruning _importance score_ and the order in which to apply merging and pruning. All experiments are performed on DeiT-S variants. FLOPs | 2.3G | 2.7G | 3.1G
---|---|---|---
Accuracy |  |  |
class attention | 76.3 | 78.2 | 79.2
column mean attention | 76.5 | 78.2 | 79.2
LTPM | 78.4 | 79.2 | 79.5
LTMP | 78.6 | 79.3 | 79.7

As shown earlier in Table 1, the correlation between the pruning importance
score and the merging similarity score is low but not zero. This means that
for small reduction target values $r$, the same token might hit the thresholds
for both merging and pruning. In Table 4, we compare pruning followed by
merging (LTPM) with merging followed by pruning (LTMP). The results confirm
that the order has no noticeable influence when the reduction rate is small,
but once more tokens need to be removed LTMP is superior to LTPM. This is not
surprising as merging has been shown to be more accurate than pruning (Bolya
et al., 2023) and that in token merging and pruning, more tokens get merged
than pruned (see Figure 5).

###  4.5 Ablation

Our approach has two important components: the learned thresholds and the
combination of merging and pruning. In Table 5, we ablate the main components
of our approach on DeiT-S for two different model sizes. We compare top-$k$
pruning and merging, where $k$ tokens are pruned in each transformer block, to
learned thresholds variants. For merging, the top-$k$ approach is equal to
what is used in ToMe (Bolya et al., 2023). Additionally, we compare merging
and pruning individually with the combined approach.

Figure 6: Distribution of the number of tokens ($k$) removed in LTP and LTM
variants of DeiT-S $r_{FLOPs}\approx 0.5$.

The results show that learned thresholds improve the accuracy of pruning
significantly, while for merging the improvements are only marginal. This
difference in accuracy between learned thresholds and top-$k$ can be explained
by examining the distribution of removed tokens as shown in Figure 6. In this
box plot, we compare the distribution of the number of removed tokens with LTM
and LTP with uniform top-$k$ pruning (i.e. the dashed line in the figure). As
can be seen in the figure, the distribution of removed tokens with LTM is
close to the uniform top-$k$ distribution, which results in many of the same
tokens being merged, while for LTP the distribution is notably dissimilar.

We also observe that naively combining merging and pruning by applying both
techniques with an equal fixed rate is worse than only token merging.

This ablation shows that both components are essential in LTMP. Combining
merging and pruning outperforms the individual techniques but only when using
the learned thresholds to balance the merging and pruning.

Table 5: Ablation of the two main components of LTMP on DeiT-S: learned thresholds and combining merging with pruning. Method | setting | FLOPs | Accuracy
---|---|---|---
pruning |  |  |
Top-$k$ | $k=16$ | 2.3G | 74.8
LTP | $r_{FLOPs}=0.5$ | 2.3G | 76.5
Top-$k$ | $k=11$ | 3.0G | 77.9
LTP | $r_{FLOPs}=0.65$ | 3.0G | 79.2
Merging |  |  |
Top-$k$ (ToMe) | $k=16$ | 2.3G | 78.1
LTM | $r_{FLOPs}=0.5$ | 2.3G | 78.2
Top-$k$ (ToMe) | $k=11$ | 3.0G | 79.3
LTM | $r_{FLOPs}=0.65$ | 3.0G | 79.5
Merging & pruning |  |  |
Top-$k$ | $k=8+8$ | 2.3G | 77.6
LTMP | $r_{FLOPs}=0.5$ | 2.3G | 78.6
Top-$k$ | $k=6+6$ | 3.1G | 79.0
LTMP | $r_{FLOPs}=0.65$ | 3.0G | 79.6

###  4.6 Inference speed

Throughout this paper, we have reported FLOPs as complexity metric. While
FLOPs are often regarded as a poor proxy for latency, it has also been shown
that latency improvements on one type of hardware often do not translate to
similar improvements on other hardware, especially in mobile and embedded
devices (Bonnaerens et al., 2022). As our complexity improvements come from
reducing the input tokens and both our masking modules and pruning and merging
implementations are parallelized, we believe FLOPs to be the best available
metric to report complexity improvements.

Nevertheless, to demonstrate our approach, we have benchmarked it on a mobile
device using PyTorch’s optimize_for_mobile function and speed_benchmark
Android binary. Table 6 shows the latency of the baseline DeiT-S and our LTMP
(with $r_{FLOPs}\approx 0.5$) reduced variant. The benchmark is conducted on a
Google Pixel 7 and averaged over 200 runs (with 50 warm-up runs prior). The
results show that the latency improvements, which achieve a reduction of
$49.52\%$, are nearly identical to the theoretical FLOPs improvements, which
have a reduction of $50.12\%$.

LTMP is also faster than ToMe while not only merging but also pruning. This
likely comes from the argsort operator that is used in top-$k$ approaches such
as ToMe and which is not well supported in many frameworks (Prillo &
Eisenschlos, 2020). Unfortunately, despite Evit, Evo-Vit and DynamicVit having
an open-source PyTorch implementation, they use operations that are not
supported by TorchScript which is required for the mobile speed_benchmark
tool.

Table 6: Latency benchmark on a Google Pixel 7. Method | FLOPs | Latency | Accuracy
---|---|---|---
DeiT-S (Baseline) | 4.6G | 212 ms | 79.8
LTMP (Ours) | 2.3G | 107 ms | 78.6
ToMe | 2.3G | 118 ms | 76.9
input | patchified | layer 1 | layer 2 | layer 3 | layer 4 | layer 5 | layer 6 | layer 7 | layer 8 | layer 9 | layer 10 | layer 11 | layer 12
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |  |  |  |  |
Figure 7: More visualization of the merging and pruning as applied to image
patches. In every layer, the most similar tokens are merged and any
unimportant tokens are pruned. The visualizations show the remaining tokens
after every layer in DeiT-S.

###  4.7 Limitations

Our learned thresholds approach requires a batch size of 1 during inference as
each image is reduced differently. This is not a limitation for most resource-
constrained applications as they typically operate inference with a batch size
of 1. If desired, our method could be extended to accommodate larger batch
sizes by either incorporating masking to a common reduction size or by
converting thresholds to the average number of tokens removed per operation
and layer, and applying these values in a top-$k$ adaptation. The former
approach will result in higher computational complexity while the latter will
result in lower accuracy.

###  4.8 Visualizations

In Figure 7, we illustrate the merging and pruning of the tokens as they are
processed through the vision transformer. It can be observed how similar parts
of the image get merged and how unimportant parts of the image are pruned.

##  5 Conclusion

In this work, we introduced Learned Thresholds token Merging and Pruning
(LTMP) for vision transformers. LTMP makes it possible to reduce the
computational cost of a vision transformer to any reduction target value with
minimal loss in accuracy. LTMP adaptively reduces the number of input tokens
processed by merging similar tokens and pruning the unimportant ones. Our
implementation uses learned thresholds which enable different merging and
pruning rates between different images and allows the model to learn the
optimal trade-off between merging and pruning across layers. As LTMP only
introduces two learnable parameters per transformer block, our method is able
to converge within a single epoch, which is an order of magnitude quicker than
other learnable approaches.

## Acknowledgments

This research received funding through the Research Foundation Flanders (FWO-
Vlaanderen) under grant 1S47820N.

## References

  * Bengio et al. (2013) Yoshua Bengio, Nicholas Léonard, and Aaron Courville.  Estimating or propagating gradients through stochastic neurons for conditional computation.  _arXiv preprint arXiv:1308.3432_ , 2013.
  * Bolya et al. (2023) Daniel Bolya, Cheng-Yang Fu, Xiaoliang Dai, Peizhao Zhang, Christoph Feichtenhofer, and Judy Hoffman.  Token merging: Your vit but faster.  In _International Conference on Learning Representations_ , 2023.
  * Bonnaerens et al. (2022) Maxim Bonnaerens, Matthias Freiberger, Marian Verhelst, and Joni Dambre.  Hardware-aware mobile building block evaluation for computer vision.  _Applied Sciences_ , 12(24):12615, 2022.
  * Chen et al. (2021) Chun-Fu Richard Chen, Quanfu Fan, and Rameswar Panda.  Crossvit: Cross-attention multi-scale vision transformer for image classification.  In _Proceedings of the IEEE/CVF international conference on computer vision_ , pp. 357–366, 2021.
  * Deng et al. (2009) J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  ImageNet: A Large-Scale Hierarchical Image Database.  In _CVPR09_ , 2009.
  * Dosovitskiy et al. (2021) Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In _International Conference on Learning Representations_ , 2021.
  * Gordon et al. (2020) Mitchell Gordon, Kevin Duh, and Nicholas Andrews.  Compressing bert: Studying the effects of weight pruning on transfer learning.  In _Proceedings of the 5th Workshop on Representation Learning for NLP_ , pp. 143–155, 2020.
  * Goyal et al. (2020) Saurabh Goyal, Anamitra Roy Choudhury, Saurabh Raje, Venkatesan Chakaravarthy, Yogish Sabharwal, and Ashish Verma.  Power-bert: Accelerating bert inference via progressive word-vector elimination.  In _ICML_ , 2020.
  * Kendall (1938) Maurice G Kendall.  A new measure of rank correlation.  _Biometrika_ , 30(1/2):81–93, 1938.
  * Kim et al. (2022) Sehoon Kim, Sheng Shen, David Thorsley, Amir Gholami, Woosuk Kwon, Joseph Hassoun, and Kurt Keutzer.  Learned token pruning for transformers.  In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , KDD ’22, pp. 784–794. Association for Computing Machinery, 2022.
  * Kitaev et al. (2019) Nikita Kitaev, Lukasz Kaiser, and Anselm Levskaya.  Reformer: The efficient transformer.  In _International Conference on Learning Representations_ , 2019.
  * Kong et al. (2022) Zhenglun Kong, Peiyan Dong, Xiaolong Ma, Xin Meng, Wei Niu, Mengshu Sun, Bin Ren, Minghai Qin, Hao Tang, and Yanzhi Wang.  Spvit: Enabling faster vision transformers via soft token pruning.  In _Computer Vision – ECCV 2022_. Springer Nature Switzerland, 2022\.
  * Li et al. (2021) Yawei Li, Kai Zhang, Jiezhang Cao, Radu Timofte, and Luc Van Gool.  Localvit: Bringing locality to vision transformers.  _CoRR_ , abs/2104.05707, 2021.
  * Liang et al. (2022) Youwei Liang, Chongjian GE, Zhan Tong, Yibing Song, Jue Wang, and Pengtao Xie.  EVit: Expediting vision transformers via token reorganizations.  In _International Conference on Learning Representations_ , 2022.
  * Liu et al. (2022) Xiangcheng Liu, Tianyi Wu, and Guodong Guo.  Adaptive sparse vit: Towards learnable adaptive token pruning by fully exploiting self-attention.  _arXiv preprint arXiv:2209.13802_ , 2022.
  * Liu et al. (2021a) Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin transformer: Hierarchical vision transformer using shifted windows.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 10012–10022, 2021a.
  * Liu et al. (2021b) Zhenhua Liu, Yunhe Wang, Kai Han, Wei Zhang, Siwei Ma, and Wen Gao.  Post-training quantization for vision transformer.  _Advances in Neural Information Processing Systems_ , 34:28092–28103, 2021b.
  * Prillo & Eisenschlos (2020) Sebastian Prillo and Julian Eisenschlos.  Softsort: A continuous relaxation for the argsort operator.  In _International Conference on Machine Learning_ , pp. 7793–7802. PMLR, 2020.
  * Rao et al. (2021) Yongming Rao, Wenliang Zhao, Benlin Liu, Jiwen Lu, Jie Zhou, and Cho-Jui Hsieh.  Dynamicvit: Efficient vision transformers with dynamic token sparsification.  _NeurIPS_ , 2021.
  * Steiner et al. (2022) Andreas Peter Steiner, Alexander Kolesnikov, Xiaohua Zhai, Ross Wightman, Jakob Uszkoreit, and Lucas Beyer.  How to train your vit? data, augmentation, and regularization in vision transformers.  _Transactions on Machine Learning Research_ , 2022.
  * Touvron et al. (2021) Hugo Touvron, Matthieu Cord, Matthijs Douze, Francisco Massa, Alexandre Sablayrolles, and Herve Jegou.  Training data-efficient image transformers & distillation through attention.  In _International Conference on Machine Learning_ , volume 139, pp. 10347–10357, July 2021.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Voita et al. (2019) Elena Voita, David Talbot, Fedor Moiseev, Rico Sennrich, and Ivan Titov.  Analyzing multi-head self-attention: Specialized heads do the heavy lifting, the rest can be pruned.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pp. 5797–5808, 2019.
  * Wang et al. (2022) Xudong Wang, Li Lyna Zhang, Yang Wang, and Mao Yang.  Towards efficient vision transformer inference: a first study of transformers on mobile devices.  In _Proceedings of the 23rd Annual International Workshop on Mobile Computing Systems and Applications_ , pp. 1–7, 2022.
  * Wightman (2022) Ross Wightman.  timm: Pytorch image models, scripts, pretrained weights, and more.  https://github.com/rwightman/pytorch-image-models, 2022.
  * Wu et al. (2022) Kan Wu, Jinnian Zhang, Houwen Peng, Mengchen Liu, Bin Xiao, Jianlong Fu, and Lu Yuan.  Tinyvit: Fast pretraining distillation for small vision transformers.  In _European Conference on Computer Vision_ , pp. 68–85. Springer, 2022.
  * Xu et al. (2022) Yifan Xu, Zhijie Zhang, Mengdan Zhang, Kekai Sheng, Ke Li, Weiming Dong, Liqing Zhang, Changsheng Xu, and Xing Sun.  Evo-vit: Slow-fast token evolution for dynamic vision transformer.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, pp. 2964–2972, 2022.
  * Yin et al. (2021) Hongxu Yin, Arash Vahdat, Jose Alvarez, Arun Mallya, Jan Kautz, and Pavlo Molchanov.  Adavit: Adaptive tokens for efficient vision transformer.  _arXiv preprint arXiv:2112.07658_ , 2021.

##  Appendix A ViT results

The main paper reports on DeiT models as they are most commonly used in
related pruning works. However the standard pretrained ViT models as found in
the timm library (Wightman, 2022; Steiner et al., 2022), are more accurate
than the DeiT models while having the same number of FLOPs. While DeiT models
are often chosen because of their more efficient training, LTMP requires only
a single training epoch making the more accurate ViT models the preferred
vision transformer to prune. We therefore also include results on ViT models
here. Figure 8 shows the Accuracy/FLOPs trade-off curve for LTMP reduced ViT
models. We also included the DeiT models as comparison. It can be observed
that ViT is indeed often the better choice, except for the heavily reduced
variants of the ‘small’ sized models, where the performance of ViT degrades
faster than DeiT. More detailed listings of the results can be found in Table
7.

Figure 8: Accuracy/FLOPs trade-offs for LTMP variants of ViT-Tiny, -Small and
-Base. DeiT variants are also plotted as comparison. FLOPs are plotted in
logarithmic scale. Table 7: Detailed results for ViT models.

Model | $r_{\text{target}}$ | FLOPs | Accuracy
---|---|---|---
ViT-T (Baseline) | 1.258G | 75.45
ViT-T | 0.85 | 1.094G | 75.4
ViT-T | 0.80 | 1.029G | 75.4
ViT-T | 0.75 | 0.970G | 75.3
ViT-T | 0.70 | 0.911G | 75.2
ViT-T | 0.65 | 0.855G | 74.9
ViT-T | 0.60 | 0.803G | 74.6
ViT-T | 0.55 | 0.737G | 74.0
ViT-T | 0.50 | 0.688G | 73.0
ViT-T | 0.45 | 0.639G | 71.7
ViT-T | 0.40 | 0.600G | 70.4

Model | $r_{\text{target}}$ | FLOPs | Accuracy
---|---|---|---
ViT-S (Baseline) | 4.608G | 81.37
ViT-S | 0.85 | 3.895G | 81.4
ViT-S | 0.80 | 3.689G | 81.1
ViT-S | 0.75 | 3.495G | 81.1
ViT-S | 0.70 | 3.309G | 80.9
ViT-S | 0.65 | 3.134G | 80.5
ViT-S | 0.60 | 2.915G | 80.2
ViT-S | 0.55 | 2.692G | 79.8
ViT-S | 0.50 | 2.501G | 79.0
ViT-S | 0.45 | 2.321G | 78.0
ViT-S | 0.40 | 2.248G | 76.4

Model | $r_{\text{target}}$ | FLOPs | Accuracy
---|---|---|---
ViT-B (Baseline) | 17.583G | 84.54
ViT-B | 0.85 | 15.024 | 84.5
ViT-B | 0.80 | 14.258G | 84.4
ViT-B | 0.75 | 13.461G | 84.3
ViT-B | 0.70 | 12.685G | 84.1
ViT-B | 0.65 | 11.975G | 83.8
ViT-B | 0.60 | 11.152G | 83.2
ViT-B | 0.55 | 10.280G | 82.8
ViT-B | 0.50 | 9.559G | 82.3
ViT-B | 0.45 | 8.841G | 81.3
ViT-B | 0.40 | 8.465G | 80.3

##  Appendix B Accuracy error range

In Table 8, we report the accuracy of 5 runs with different random seeds for
the main results in the paper (i.e. the data points highlighted in gray in
Table 2). It can be observed that the error range on the accuracy is about
$\pm 0.3$ for these data points. We, therefore, report the results in this
paper with 3 significant figures.

Model | $r_{\text{target}}$ | FLOPs | Acc. #1 | Acc. #2 | Acc. #3 | Acc. #4 | Acc. #5 | Acc. range
---|---|---|---|---|---|---|---|---
DeiT-S | 0.80 | 3.8G | 79.816 | 79.802 | 79.824 | 79.794 | 79.814 | 0.030
DeiT-S | 0.62 | 3.0G | 79.598 | 79.596 | 79.584 | 79.578 | 79.564 | 0.034
DeiT-S | 0.55 | 2.7G | 79.262 | 79.266 | 79.242 | 79.274 | 79.254 | 0.032
DeiT-S | 0.45 | 2.3G | 78.608 | 78.604 | 78.588 | 78.574 | 78.558 | 0.050
Table 8: Accuracy ranges for 5 different runs of the main data points in the
paper.


11institutetext: Max Planck Institute for Solar System Research, Justus-von-
Liebig-Weg 3, 37077 Göttingen, Germany
11email<EMAIL_ADDRESS>
# Scattering, absorption, and thermal emission by large cometary dust
particles

Synoptic numerical solution

Johannes Markkanen 11 Jessica Agarwal 11

###### Abstract

Context. Remote light scattering and thermal infrared observations provide
clues about the physical properties of cometary and interplanetary dust
particles. Identifying these properties will lead to a better understanding of
the formation and evolution of the Solar System.

Aims. We present a numerical solution for the radiative and conductive heat
transport in a random particulate medium enclosed by an arbitrarily shaped
surface. The method will be applied to study thermal properties of cometary
dust particles.

Methods. The recently introduced incoherent Monte Carlo radiative transfer
method developed for scattering, absorption, and propagation of
electromagnetic waves in dense discrete random media is extended for radiative
heat transfer and thermal emission. The solution is coupled with the
conductive Fourier transport equation that is solved with the finite-element
method.

Results. The proposed method allows the synoptic analysis of light scattering
and thermal emission by large cometary dust particles consisting of
submicrometer-sized grains. In particular, we show that these particles can
sustain significant temperature gradients resulting in the superheating factor
phase function observed for the coma of comet 67P/Churyumov-Gerasimenko.

###### Key Words.:

Radiative transfer, methods: numerical, comets: individual: 67P

##  1 Introduction

The physical properties of cosmic dust particles, which include size, shape,
porosity, and composition, can be constrained from remote light scattering and
thermal emission observations. Observations of scattered solar and thermally
emitted light have mostly been interpreted separately using different particle
models and numerical methods (Kolokolova et al., 2004; Kimura et al., 2016).
This may lead to contradictions in the retrieved particle properties as the
interpretation procedure is not self-consistent. Using the same numerical
method and the same particle model for scattered solar and thermally emitted
light (i.e., _synoptic modeling)_ results in more reliable estimates for the
physical properties of dust particles. Such modeling, however, introduces
various numerical challenges due to the multiscale nature of the problem. The
solution is needed from ultraviolet to thermal infrared. In addition, cometary
dust particles contain multiscale structures. They are thought to be large
aggregates consisting of submicrometer-sized grains (Güttler et al., 2019).

Numerically exact electromagnetic techniques allow for analyses of scattering
and radiative heat transfer properties of small aggregated particles.
Combining the standard electromagnetic methods with the statistical
fluctuation driven electrodynamic solvers, the radiative heat transfer problem
can be solved directly from Maxwell’s equations (Jin et al., 2017). This
approach is computationally heavy, and consequently it is not applicable to
particles that are much larger than the optical wavelengths or that contain a
large number of individual grains. For sufficiently large objects such as
asteroids and comet nuclei, a typical approach is to localize absorption and
emission sources on the surfaces of the object using a proper boundary
condition. Then, by introducing the effective conductive and radiative heat
transfer coefficient, heat transport inside the object can be modeled (Huebner
et al., 2006; Delbo et al., 2015). This simple model is a valid approximation
only for one-dimensional problems, and it introduces significant errors when
the object has locally three-dimensional features (on length scales on the
order of the mean free path length). Thus, large objects with small details
cannot be treated rigorously with the existing computational techniques, due
to the large computational time and memory requirement or inadequate
approximations.

It is also possible to use the radiative transfer equation (RTE) in
combination with the Fourier heat transport equation to deal with heat
transfer in participating media (Howell et al., 2010; Modest, 2013). A large
number of different solution techniques to solve the RTE exist, for example
the discrete ordinates (Hardy et al., 2016), finite-volume (Coelho, 2014),
finite-element methods (Richling et al., 2001), and Monte Carlo ray tracing
methods (Howell, 1998). Nonetheless, the standard RTE is only applicable to
sparse discrete random media where the grains are in each other’s far zone.
Recently, a numerical framework for electromagnetic scattering, absorption,
and propagation in dense discrete random media was introduced which extends
the applicability range of the standard RTE (Muinonen et al., 2018; Markkanen
et al., 2018b; Väisänen et al., 2019). Furthermore, Markkanen et al. (2018a)
augmented the method with the coherent field contribution to account for
surface reflections and refractions allowing the scattering analysis of
millimeter-sized cometary dust particles consisting of submicrometer-sized
grains.

In this paper we extend the numerical method developed in Markkanen et al.
(2018a) for the radiative heat transfer problems. The radiative heat transfer
solution is coupled with the transient conductive heat transfer equation,
which is solved by the finite-element method (FEM). The presented approach
allows us to analyze the scattering and thermal properties of particles
consisting of aggregates of small grains over a wide wavelength range.
Finally, we show that large cometary dust particles can hold significant
temperature gradients, which explains the superheating phase function of the
coma of comet 67P/Churyumov-Gerasimenko measured by the Rosetta/VIRTIS
(Visible and Infrared Thermal Imaging Spectrometer) instrument.

##  2 Methods

In this section, we present our numerical framework for solving radiative and
conductive heat transfer in a particulate random medium enclosed by an
arbitrarily shaped boundary. We neglect heat convection and assume that
external radiation and thermal fluctuations are the only heat sources. We also
assume that the radiative heat transport is much faster than the conductive
heat transport. Thus, we consider the radiative part as instantaneous time
harmonic radiation.

###  2.1 Energy equation

Let us consider a finite body $\Omega\in\mathbb{R}^{3}$ characterized by the
position $\boldsymbol{r}$, dependent density $\rho$, specific heat capacity
$c_{p}$, and conductivity $\kappa$ in free space. Our goal is to solve the
time $t$ dependent temperature $T$ under the influence of the external
radiation source. The energy balance equation in $\Omega$ for the combined
conductive and radiative heat transfer reads as

$\rho c_{\rm p}\frac{\partial T}{\partial t}-\nabla\cdot\kappa\nabla
T=-\nabla\cdot\boldsymbol{q}_{\rm r}\,\,\,\text{in}\,\Omega$ (1)

with the boundary condition

$\boldsymbol{n}\cdot\kappa\nabla T=\boldsymbol{n}\cdot\boldsymbol{q}_{\rm
r}\,\,\,\text{on}\,\partial\Omega,$ (2)

where $\boldsymbol{q}_{\rm r}$ denotes the total radiative flux,
$\partial\Omega$ is the boundary of $\Omega$, and $\boldsymbol{n}$ is the
inward pointing normal vector on $\partial\Omega$.

###  2.2 Finite-element solution for the energy equation

To solve the energy balance equation (1) and (2), we apply the standard
finite-element method to discretize the spatial component of the energy
equation. Taking the inner product with the testing functions $w$ and
integrating by parts, the weak formulation reads as follows:

Find $T\in H^{1}(\Omega)$ such that

$\int_{\Omega}w\rho c_{\rm p}\frac{\partial T}{\partial
t}\,\text{d}V+\int_{\Omega}\nabla w\cdot\kappa\nabla
T\,\text{d}V=-\int_{\Omega}w\nabla\cdot\boldsymbol{q}_{\rm r}\,\text{d}V$ (3)

is valid for all $w\in H^{1}(\Omega)$. Here $H^{1}(\Omega)$ denotes the space
of square integrable functions whose derivatives are also square integrable.
We note that equation (3) automatically satisfies the Neumann boundary
condition (2).

Next we discretize the time derivative using the finite-difference formula

$\frac{\partial T}{\partial t}\approx\frac{1}{\tau}(T_{t+1}-T_{t}),$ (4)

where the subscript $t$ indicates the time step and $\tau$ is its size.
Substituting (4) into (3) and using the implicit backward Euler time
integration scheme, we obtain

$\displaystyle\rho c_{\rm
p}\tau^{-1}(\int_{\Omega}wT_{t+1}\,\text{d}V-\int_{\Omega}wT_{t}\,\text{d}V)+\int_{\Omega}\nabla
w\cdot\kappa\nabla T_{t+1}\,\text{d}V=$
$\displaystyle-\int_{\Omega}w(\nabla\cdot\boldsymbol{q}_{\rm r})_{t+1}.$ (5)

We discretize the domain $\Omega$ with linear tetrahedral elements and expand
the unknown temperature $T$ into the linear nodal basis functions $u^{m}$ as
$T\approx\sum_{m}x^{m}u^{m}$, where $x^{m}$ are the unknown coefficients.
Using Galerkin’s method, in which the testing functions and basis functions
are identical $w^{m}=u^{m}$ , we can write the mass and stiffness matrices as

$M=\rho c_{\rm p}\int_{\Omega}w^{n}u^{m}\,\text{d}V,$ (6)
$S=\int_{\Omega}\nabla w^{n}\cdot\kappa\nabla u^{m}\,\text{d}V,$ (7)

and the force vector as

$F=-\int_{\Omega}w^{n}\nabla\cdot\boldsymbol{q}_{\rm r}\,\text{d}V.$ (8)

Finally, we arrive at the expression for the temperature coefficient vector
$x$ given by

$x_{t+1}=(M+\tau S)^{-1}(Mx_{t}+\tau Fx_{t+1}).$ (9)

The force vector is strongly nonlinear with respect to temperature ($F\sim
T^{4}$). Hence, estimating $\tau Fx_{t+1}$ with $\tau Fx_{t}$ would require
using a very small time step $\tau$ to reach a converging numerical solution.
Therefore, we solve equation (9) for $x_{t+1}$ by using an iterative method.
In the next section we discuss how to compute the force vector $F$ in a densly
packed particulate random medium enclosed by an arbitrarily shaped boundary
$\partial\Omega$ with the Monte Carlo radiative transfer approach.

###  2.3 Radiative transfer solution

To calculate the absorbed and thermally emitted power (i.e., the force vector
$F$ in Eq. (9)), we employ an approximate version of the recently introduced
radiative transfer with the reciprocal transactions (R2T2) method for
electromagnetic scattering and propagation in dense discrete random media
(Muinonen et al., 2018; Markkanen et al., 2018b; Väisänen et al., 2019). R2T2
is based on the Monte Carlo ray tracing of the order-of-scattering
interactions in which the fundamental scatterers are volume elements
containing a large number of small grains. Incoherent electromagnetic
interactions within the volume elements and the order-of-scattering
interactions between the volume elements are computed via the exact
electromagnetic solver, namely the fast superposition T-matrix method (FaSTMM)
(Markkanen & Yuffa, 2017).

Here we use the far-field approximation to compute the order-of-scattering
interactions. Yet, the interactions within the volume elements are computed
exactly. We also neglect the reciprocal scattering diagrams (i.e., those
responsible for the coherent backscattering effect). This allows us to include
the coherent field contribution in terms of reflections and refractions on the
surfaces of a particle using geometric optics (Markkanen et al., 2018a;
Muinonen et al., 2019).

Computations proceed as follows. First, we compute scattering properties for
volume elements, as described in Section 2.3.1. This step is done separately
for each wavelength and for each material region, but it does not depend on
the shape and size of the entire particle. Second, the computed scattering
properties from the first step are used as input parameters in the combined
Monte Carlo geometric optics radiative transfer solver, as discussed in
Section 2.3.2. This step is computationally independent of the first step for
different particle shapes and sizes.

####  2.3.1 Scattering properties of volume elements

Input parameters for our radiative transfer algorithm are the incoherent
scattering phase function $M_{11}^{\rm ic}$, incoherent albedo $\omega$,
incoherent mean free path $\ell$, coherent effective refractive index $m_{\rm
eff}$, and absorption cross section of the ensemble averaged volume element
$C_{\rm abs}$.

We compute these parameters for spherical ensemble-averaged volume elements
that contain a large number of small spherical monomers. The radius of the
volume element $R_{0}$ should be large enough to capture the statistics of the
random medium. First, we generate a sample medium into a periodic box by
randomly depositing monomers into the box until the target volumetric filling
factor $v$ is reached. Second, we draw a sample spherical volume element from
the box, as explained in detail by Väisänen et al. (2019). Third, we apply
FaSTMM to compute the Jones matrix $S_{i}$ of such a volume element $i$.
Repeating the process $N$ times ($N=512$ is used in this work), we compute the
coherent Jones scattering matrix as

$S^{\rm c}=\frac{1}{N}\sum_{i=1}^{N}S_{i},$ (10)

and the incoherent Jones scattering matrix for the $i$th state as

$S^{\rm ic}_{i}=S_{i}-S^{\rm c}.$ (11)

Averaging the square of the incoherent Jones matrix over the ensemble of all
$N$ sample volume elements, we obtain the ensemble-averaged one-one element of
the Mueller matrix, composed of the incoherent scattering phase function

$M_{11}^{\rm ic}=\frac{1}{N}\sum_{i}^{N}|S^{\rm ic}_{i}|^{2}$ (12)

and the incoherent scattering cross section

$C_{\rm sca}^{\rm ic}=\int_{S^{2}}M_{11}^{\rm ic}\,\text{d}S,$ (13)

where $S^{2}$ is the surface of the unit sphere. The incoherent albedo is
defined as

$\omega=\frac{C_{\rm sca}^{\rm ic}}{C_{\rm sca}^{\rm ic}+C_{\rm abs}},$ (14)

where $C_{\rm abs}$ is the ensemble-average absorption cross section, and the
incoherent mean free path length as

$\ell=\frac{V}{C^{\rm ic}_{\rm sca}+C_{\rm abs}},$ (15)

in which $V$ is the volume of the volume element.

Finally, we solve the coherent effective medium parameter $m_{\rm eff}$. This
is done by matching the coherent scattering cross section $C^{\rm c}_{\rm
sca}$ to the $C_{\rm sca}$ of the equal-sized sphere using the Mie solution.
$C^{\rm c}_{\rm sca}$ is defined as in (13), but the incoherent Mueller matrix
is replaced with the coherent matrix.

####  2.3.2 Monte Carlo ray tracing

We trace the rays in the same tetrahedral mesh as is used in the FEM. This
makes it trivial to couple the radiative transfer solution (i.e.,
$\nabla\cdot\boldsymbol{q}_{\rm r}$) to the FEM solution. We decompose the
divergence of the total flux into two parts $\nabla\cdot\boldsymbol{q}_{\rm
r}=q_{\rm sol}+q_{\rm em}$ in which $q_{\rm sol}$ is the absorption rate due
to the solar radiation, and $q_{\rm em}$ is the absorption plus emission rate
due to the thermally fluctuating sources.

To compute $q_{\rm sol}$, we launch $N^{e}_{\rm ray}$ rays for a given
wavelength band $\Delta_{\lambda}$ from outside $\Omega$. The power of each
incident ray is given by

$E^{e}=\int_{\lambda_{1}}^{\lambda_{2}}\frac{P(\lambda)G}{N^{e}_{\rm
ray}}\,\text{d}\lambda,$ (16)

where $P(\lambda)$ is the flux density of the incident radiation and $G$ is
the geometric cross section of the particle. Each ray has a specific position
$\boldsymbol{r}$ and propagation direction $\hat{\boldsymbol{k}}$ depending on
the source.

For thermally emitted radiation, to compute $q_{\rm em}$, $N^{i}_{\rm ray}$
rays are launched from each tetrahedron $C_{t}$ with a random position and
direction. The emitted power per ray for a given wavelength is given by

$E^{i}_{C_{t}}=\frac{4N_{\rm d}}{N^{i}_{\rm
ray}}\int_{\lambda_{1}}^{\lambda_{2}}\pi B(\lambda,T_{c_{t}},m_{\rm
eff})C_{\rm abs}\,\text{d}\lambda,$ (17)

where $T_{c_{t}}$ is the average temperature of the tetrahedron $C_{t}$ and
$N_{\rm d}$ is the number density of the volume elements in the tetrahedron
$C_{t}$. The Planck function in a medium with the refractive index $m$ written
in terms of the free space wavelength $\lambda_{0}$ reads as

$B_{\lambda}(\lambda_{0},T,m)=\frac{2hc^{2}Re\\{m\\}^{2}}{\lambda_{0}^{5}}\frac{1}{\exp(\frac{hc}{\lambda_{0}k_{\rm
b}T})-1},$ (18)

where $h$ is the Planck constant, $c$ is the speed of light in vacuum, and
$k_{\rm b}$ is the Stefan-Boltzmann constant.

Once we have the initial position $\boldsymbol{r}$, direction
$\hat{\boldsymbol{k}}$, and power of the ray $E$, we start tracing. If the ray
hits the boundary $\partial\Omega$, it reflects and refracts according to
Snel’s law and the power is updated from Fresnel’s coefficients which are
calculated using the effective refractive index $m_{\rm eff}$. The position
and direction of the reflected and refracted rays are updated.

If the ray is in $\Omega$, the distance to the next scattering event is
generated as

$d=-\ell\ln\Phi,$ (19)

where $\Phi$ is a uniform random number within $[0,1[$. The scattering event
happens if the ray does not cross the boundary $\partial\Omega$. Absorbed
power $E_{\rm abs}=E(1-\omega)$ in the scattering event is added into the
total energy of the tetrahedron where the scattering event occurs, and the
ray’s power is updated ($E=E\omega$). A new propagation direction
$\hat{\boldsymbol{k}}$ is drawn from the cumulative distribution function
generated from the phase function of the ensemble averaged incoherent volume
element $M_{11}^{\rm ic}$. Finally, if the ray does not scatter and it does
not cross the boundary $\partial\Omega$ , intensity is collected and added to
the total scattering phase function. This process is repeated until the ray’s
power has decreased under the predefined threshold.

The radiative transfer computation gives us
$\int_{C_{t}}\nabla\cdot\boldsymbol{q}_{\rm r}\,\text{d}V$ for each
tetrahedron, which is then used to compute the force vector $F$ using (8). In
addition, since the escaped rays are recorded, the total scattering and
thermal emission phase functions and cross sections are obtained for each
wavelength.

##  3 Numerical results

In this section, we present some numerical results. To validate the method, we
compare our solution with the known approximate solution for small loosely
packed grains computed by employing the effective medium approximation (EMA)
and the Mie theory. Then, we apply the method to explain the superheating
factor phase function measured for the coma of comet 67P/Churyumov-
Gerasimenko.

###  3.1 Comparison to the effective medium approximation

Let us first consider a spherical medium consisting of small spherical
silicate grains of radius $r=$20\text{\,}\mathrm{nm}$$. The grains are
randomly positioned in a spherical domain $\Omega$ with radius $R$. The
volumetric filling factor is $\nu=0.05$. The wavelength dependent refractive
index for magnesium-iron silicate, from Dorschner et al. (1995), is plotted in
Fig. 1. Here we assume that the refractive index is independent of temperature
in the studied temperature range.

Figure 1: Real and imaginary parts of the complex refractive index for
magnesium-iron silicate.

The external radiation source is the Sun, which is assumed to radiate as a
perfect blackbody of temperature $T_{\odot}=5777$K. The solar flux $P$ in the
given wavelength band ($\lambda_{1}-\lambda_{2}$) at the distance
$\Delta=1.3$AU from the Sun is given by

$P=(\frac{r_{\odot}}{\Delta})^{2}\int_{\lambda_{1}}^{\lambda_{2}}\pi
B_{\lambda}(\lambda,T_{\odot})\,\text{d}\lambda,$ (20)

where $r_{\odot}$ is the radius of the Sun.

Since the grain size is small, we can calculate a reference solution at the
equilibrium temperature using the steady-state energy balance equation, EMA,
and the Mie theory. To compute the absorption cross section $C_{\rm abs}$, we
use the Maxwell-Garnett mixing rule (Sihvola, 2000) and the Mie theory. For an
isothermal spherical particle, the steady-state energy balance equation reads
as

$(\frac{r_{\odot}}{\Delta})^{2}\int_{0}^{\infty}\pi
B_{\lambda}(\lambda,T_{\odot})C_{\rm
abs}(\lambda)\,\text{d}\lambda=4\int_{0}^{\infty}\pi
B_{\lambda}(\lambda,T)C_{\rm abs}(\lambda)\,\text{d}\lambda.$ (21)

To compare the proposed method to the steady-state analysis with EMA (21), we
ran the transient simulation until the equilibrium was reached. Figure 2 shows
the total absorbed solar radiation and thermally emitted power as a function
of wavelength calculated by the proposed method and EMA for three different
particle sizes $R=0.01,0.1$, and $1.0$ mm. We observed excellent agreement
between the two methods; however, we expect a decrease in the accuracy of our
method as the particle size becomes smaller. This is evident as we used
geometric optics to approximate the coherent reflections and refractions on
the surface of a particle.

In the computations the volume element radius is
$r_{0}=$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and the time step is
$\tau=$0.1\text{\,}\mathrm{s}$$, $\tau=$1.0\text{\,}\mathrm{s}$$, and
$\tau=$10\text{\,}\mathrm{s}$$ for $0.01\text{\,}\mathrm{mm}$,
$0.1\text{\,}\mathrm{mm}$, and $1.0\text{\,}\mathrm{mm}$ particles,
respectively. We also assumed high thermal conductivity $\kappa=100$ W/m/K to
remove possible temperature gradients inside the body. Although it is not a
realistic assumption for real dust, it serves as a validation by simplifying
the problem. It is worth noting that the density $\rho=1000$ kg/m3 and the
specific heat capacity $c_{p}=500$ J/kg/K do not affect the equilibrium
temperature as $\frac{\partial T}{\partial t}=0$ in (1) for the steady-state
solution; they only affect the time it takes to reach equilibrium. Figure 3
shows the temperature evolution of the particles when the initial temperature
$T_{0}=200$K. We see that the transient temperature reaches the equilibrium
temperature computed by using the steady-state approach with EMA.

Figure 2: Total absorbed power of solar radiation (top) and thermally emitted
power (bottom) by particles with radii $10\text{\,}\mathrm{\SIUnitSymbolMicro
m}$, $100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{,}$ and
$1.0\text{\,}\mathrm{mm}$ as a function of wavelength computed by the proposed
method and the effective medium approximation (EMA) are shown. Grain size is
$20\text{\,}\mathrm{nm}$ and the volumetric filling factor $v=0.05$. Figure 3:
Temperature time evolution for porous ($v=0.05$) aggregates of 20 nm silicate
grains, with initial temperature $T_{0}=200$K. The unit on the time axis is
scaled by the radii of the particles $R$.

Next we studied larger grains with $r=100$ nm and increased the filling factor
to $v=0.3$. These values are commonly used to model cometary dust (Kolokolova
et al., 2004). It is evident from Fig. 4 and Fig. 5 that EMA leads to a
different absorption rate, which consequently affects the equilibrium
temperature and thermal emission; the EMA is no longer a valid approximation
at shorter wavelengths as the grain size is approximately the same as the
wavelength. Thus, care should be taken when applying EMAs to interpret
infrared observations of cometary dust particles consisting of 100 nm-sized
grains.

Figure 4: As in Fig. 2, but for grain size of $100\text{\,}\mathrm{nm}$ and
filling factor $v=0.3$.  Figure 5: As in Fig. 3, but for grain size of
$100\text{\,}\mathrm{nm}$ and filling factor $v=0.3$.

###  3.2 Application to the coma of 67P

We recently introduced a particle model that explains the scattering phase
functions of the coma of comet 67P/Chruyumov-Gerasimenko measured by the
Rosetta/OSIRIS (Optical, Spectrocopic and Infrared Remote Imaging System)
instrument (Markkanen et al., 2018a). The particle model consists of
aggregated submicrometer-sized organic grains and micrometer-sized silicate
grains. The shapes of the aggregates are Gaussian random spheres. It is thus
interesting to study whether the same model can also explain the thermal
infrared observations of the same comet. Since there is no refractive index
available for cometary organic material over a wide wavelength range, we used
the refractive index of carbonaceous dust analogues (Jager et al., 1998). For
silicates, we used the refractive index presented in Fig. 1 that corresponds
to magnesium-iron silicate mineral (Dorschner et al., 1995).

We concentrate on the superheating factor phase function reported by Bockelée-
Morvan et al. (2019) for the coma of the comet 67P. The superheating factor is
defined as the ratio of the color temperature $T_{\rm c}$ to the equilibrium
temperature of an ideal blackbody $T_{\rm b}=278/\Delta^{2}$[K] (i.e.,
$S=T_{\rm c}/T_{\rm b}$). To obtain the color temperature $T_{\rm c}$ for a
given phase angle, Bockelée-Morvan et al. (2019) fitted the Planck function to
the Rosetta/VIRTIS near-infrared spectra
($3-$5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$). They found that the
superheating factor phase function has a phase dependence indicating that the
particles can hold temperature gradients if the particles are assumed to be
randomly oriented.

Figure 6: Computed superheating factors as a function of phase angle for
different particle sizes. The heat conduction coefficient is $\kappa=0.00025$
W/m/K. Also shown are the measured values.

Computed superheating factors for different particle sizes together with the
measured values (Fig. 8 in Bockelée-Morvan et al. (2019)) as a function of
phase angle are presented in Fig. 6. We see that only large particles can
sustain temperature gradients resulting in the phase functions with
significant slopes. The difference between the maximum and minimum temperature
inside the particle with respect to the particle size is plotted in Fig. 7. It
is evident that in order to model the measured superheating factor phase
function slopes, the dust particles must sustain significant temperature
gradients. Large temperature gradients, in turn, may have interesting effects
on the dust dynamics in the inner coma via rocket forces caused by the
asymmetrical sublimation of ices. Thermal break-down of particles may also
play a role as the dust is transported from the inner to outer coma.
Nevertheless, these interesting topics are beyond of the scope of this paper
and are left for future research.

Figure 7: Difference between the maximum and minimum temperature inside the
particle as a function of particle radius for an example particle.

In our computations, the heat conduction coefficient is assumed to be
$\kappa=0.00025$ W/m/K, which corresponds to a realistic value for a porous
dust aggregate (Krause et al., 2011; Arakawa et al., 2017; Sakatani et al.,
2017). The heat conduction coefficient depends on porosity, grain size, and
composition. Increasing the heat conduction coefficient flattens the
superheating phase functions. Decreasing $\kappa$ increases the slope until
the radiative heat transport becomes more dominant. We also assumed that the
angular velocities of the particles in the direction perpendicular to the Sun
are slow ($\omega=\frac{2\pi}{50}\frac{0.1\text{mm}}{R}$ rad/s), hence the
contribution of rotation to the temperature gradients is negligible. The slow
rotation rate was introduced in order to speed up numerical convergence.
Higher angular velocities decrease temperature gradients and flatten the
superheating factor phase functions. However, it is important to note that if
the rotation axis points towards the Sun, it does not affect the temperature
gradients.

Finally, we averaged the modeled superheating phase functions over a
differential power-law size distribution of index -3. The averaging results in
a superheating factor phase function that is consistent with the measured
value if the conductive heat coefficient $\kappa=0.00025$ W/m/K, as
demonstrated in Fig 8. Such a low $\kappa$ value suggests that the particles
must be porous. The superheating factor phase function also indicates that the
dust particles in the inner coma must be large. The dominating particle size
range is around 10-$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Smaller
particles would show flatter superheating factor phase functions with higher
absolute values whereas larger particles would result in a steeper slope and
smaller absolute values. The size range is also consistent with other dust
models explaining the intensity phase functions measured by Rosetta/OSIRIS
(Moreno et al., 2018; Markkanen et al., 2018a).

Figure 8: Modeled superheating factor phase function for different conductive
heat transfer coefficients $\kappa$ averaged over the differential power-law
size distribution of index -3. The minimum cutoff is
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and the maximum
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

##  4 Conclusions

We presented a synoptic numerical solution for scattering, absorption, and
thermal emission by large aggregated dust particles consisting of
submicrometer-sized grains. The numerical solution is self-consistent and
allows for analysis of multi-instrument data sets. We solved the transient
energy equation using the finite-element method in which the radiative heat
transfer part is included as an additional forcing term. We developed a new
algorithm, based on the Monte Carlo radiative transfer with reciprocal
transactions framework, to compute the forcing term in the dense discrete
random medium enclosed by an arbitrarily shaped surface. The method explicitly
computes the radiative heat transfer part in contrast to the standard thermal
models in cometary science where it is included in the heat conduction
coefficient by assuming that heat transport is one dimensional. Thus, the
developed method provides a full three-dimensional solution for radiative heat
transfer problems.

We compared the solution of the proposed method to that obtained by employing
the effective medium approximation and Mie theory. The solutions match if the
grains are small and loosely packed, but deviate for grain sizes and packing
densities that are typically assumed to form cometary dust particles. This
implies that using EMAs to interpret thermal infrared observations of cometary
dust may lead to unreliable results. Finally, we showed that particles
$10-$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ in size can reproduce the
superheating factor phase function of the coma of comet 67P if the dust
particles hold significant temperature gradients. This means that the heat
conduction coefficient must be very low, hence particles must be porous. The
presented analysis also suggests that the visible phase functions measured by
Rosetta/OSIRIS and the superheating phase function from Rosetta/VIRTIS are
consistent with each other.

###### Acknowledgements.

This work has been funded by the ERC Starting Grant No. 757390 Comet and
Asteroid Re-Shaping through Activity (CAstRA). Computational resources have
been provided by Gesellschaft für Wissenschaftliche Datenverarbeitung mbH
Göttingen (GWDG).

## References

  * Arakawa et al. (2017) Arakawa, S., Tanaka, H., Kataoka, A., & Nakamoto, T. 2017, A&A, 608, L7
  * Bockelée-Morvan et al. (2019) Bockelée-Morvan, D., Leyrat, C., Erard, S., et al. 2019, A&A
  * Coelho (2014) Coelho, P. J. 2014, JQSRT, 145, 121
  * Delbo et al. (2015) Delbo, M., Mueller, M., Emery, J. P., Rozitis, B., & Capria, M. T. 2015, Asteroids IV, 107
  * Dorschner et al. (1995) Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503
  * Güttler et al. (2019) Güttler, C., Mannel, T., Rotundi, A., et al. 2019, A&A
  * Hardy et al. (2016) Hardy, D. L., Favennec, Y., Domingues, G., & Rousseau, B. 2016, J. Appl. Math. Phys., 4, 1488
  * Howell et al. (2010) Howell, J., Siegel, R., & Mengüc, M. P. 2010, Thermal Radiation Heat Transfer
  * Howell (1998) Howell, J. R. 1998, J. Heat Trans., 120, 547
  * Huebner et al. (2006) Huebner, W. F., Benkhoff, J., Capria, M.-T., et al., eds. 2006, Heat and Gas Diffusion in Comet Nuclei
  * Jager et al. (1998) Jager, C., Mutschke, H., & Henning, T. 1998, A&A, 332, 291
  * Jin et al. (2017) Jin, W., Messina, R., & Rodriguez, A. W. 2017, Phys. Rev. B, 95, 161409
  * Kimura et al. (2016) Kimura, H., Kolokolova, L., Li, A., & Lebreton, J. 2016, Light Scattering and Thermal Emission by Primitive Dust Particles in Planetary Systems, ed. A. Kokhanovsky (Berlin, Heidelberg: Springer Berlin Heidelberg), 363–418
  * Kolokolova et al. (2004) Kolokolova, L., Hanner, M. S., Levasseur-Regourd, A.-C., & Gustafson, B. Å. S. 2004, Physical properties of cometary dust from light scattering and thermal emission, ed. M. C. Festou, H. U. Keller, & H. A. Weaver, 577–604
  * Krause et al. (2011) Krause, M., Blum, J., Skorov, Y. V., & Trieloff, M. 2011, Icarus, 214, 286
  * Markkanen et al. (2018a) Markkanen, J., Agarwal, J., Väisänen, T., Penttilä, A., & Muinonen, K. 2018a, APJL, 868, L16
  * Markkanen et al. (2018b) Markkanen, J., Väisänen, T., Penttilä, A., & Muinonen, K. 2018b, Opt. Lett., 43, 2925
  * Markkanen & Yuffa (2017) Markkanen, J. & Yuffa, A. J. 2017, JQSRT, 189, 181
  * Modest (2013) Modest, M. F. 2013, in Radiative Heat Transfer (Third Edition), third edition edn., ed. M. F. Modest (Boston: Academic Press), 279 – 302
  * Moreno et al. (2018) Moreno, F., Guirado, D., Muñoz, O., et al. 2018, AJ, 156, 237
  * Muinonen et al. (2018) Muinonen, K., Markkanen, J., Väisänen, T., Peltoniemi, J., & Penttilä, A. 2018, Opt. Lett., 43, 683
  * Muinonen et al. (2019) Muinonen, K., Väisänen, T., Martikainen, et al. 2019, J. Vis. Exp, 149
  * Richling et al. (2001) Richling, S., Meinköhn, E., Kryzhevoi, N., & Kanschat, G. 2001, A&A, 380, 776
  * Sakatani et al. (2017) Sakatani, N., Ogawa, K., Iijima, Y., et al. 2017, AIP Advances, 7, 015310
  * Sihvola (2000) Sihvola, A. 2000, Subsurface sensing technologies and applications, 1, 393
  * Väisänen et al. (2019) Väisänen, T., Markkanen, J., Penttilä, A., & Muinonen, K. 2019, PLOS ONE, 14, 1