# Uncertainty Estimations by Softplus normalization in Bayesian Convolutional
Neural Networks with Variational Inference

Kumar Shridhar
Department of Computer Science
TU Kaiserslautern
MindGarage
<EMAIL_ADDRESS>
Authors contributed equally. Felix Laumann11footnotemark: 1
Department of Mathematics
Imperial College London
NeuralSpace
<EMAIL_ADDRESS>
Marcus Liwicki
Department of Computer Science
Luleå University of Technology
MindGarage
<EMAIL_ADDRESS>
###### Abstract

We introduce a novel uncertainty estimation for classification tasks for
Bayesian convolutional neural networks with variational inference. By
normalizing the output of a Softplus function in the final layer, we estimate
aleatoric and epistemic uncertainty in a coherent manner. The intractable
posterior probability distributions over weights are inferred by Bayes by
Backprop. Firstly, we demonstrate how this reliable variational inference
method can serve as a fundamental construct for various network architectures.
On multiple datasets in supervised learning settings (MNIST, CIFAR-10,
CIFAR-100), this variational inference method achieves performances equivalent
to frequentist inference in identical architectures, while the two desiderata,
a measure for uncertainty and regularization are incorporated naturally.
Secondly, we examine how our proposed measure for aleatoric and epistemic
uncertainties is derived and validate it on the aforementioned datasets.

##  1 Introduction

Convolutional neural networks excel at tasks in the realm of image
classification (e.g. [1, 2, 3]). However, from a probability theory
perspective, it is unjustifiable to use single point-estimates as weights to
base any classification on. Cnns with frequentist inference require
substantial amounts of data examples to train on and are prone to overfitting
on datasets with few examples per class.
In this work, we apply Bayesian methods to Cnns in order to add a measure for
uncertainty and regularization in their predictions, respectively their
training. This approach allows the network to express uncertainty via its
parameters in form of probability distributions (see Figure 1). At the same
time, by using a prior probability distribution to integrate out the
parameters, we compute the average across many models during training, which
gives a regularization effect to the network, thus preventing overfitting.
We build our Bayesian Cnn upon Bayes by Backprop [4, 5], and approximate the
intractable true posterior probability distributions $p(w|\mathcal{D})$ with
variational probability distributions $q_{\theta}(w|\mathcal{D})$, which
comprise the properties of Gaussian distributions $\mu\in\mathbb{R}^{d}$ and
$\sigma\in\mathbb{R}^{d}$, denoted $\mathcal{N}(\theta|\mu,\sigma^{2})$, where
$d$ is the total number of parameters defining a probability distribution. The
shape of these Gaussian variational posterior probability distributions,
determined by their variance $\sigma^{2}$, expresses an uncertainty estimation
of every model parameter. The main contributions of our work are as follows:

  1. 1.

We present how Bayes by Backprop can be efficiently applied to Cnns. We
therefore introduce the idea of applying two convolutional operations, one for
the mean and one for the variance.

  2. 2.

We empirically show how this generic and reliable variational inference method
for Bayesian Cnns can be applied to various Cnn architectures without any
limitations on their performances, but with intrinsic regularization effects.
We compare the performances of these Bayesian Cnns to Cnns which use single
point-estimates as weights, i.e. which are trained by frequentist inference.

  3. 3.

We explain and implement Softplus normalization by means of examining how to
estimate aleatoric and epistemic uncertainties without employing an additional
Softmax function in the output layer [6] which brings an inconsistency of
activation functions into the model.

This paper is structured as subsequently outlined: after we have introduced
our work here, we secondly review Bayesian neural networks with variational
inference, including previous works, an explanation of Bayes by Backprop and
its implementation in Cnn. Thirdly, we examine aleatoric and epistemic
uncertainty estimations with an outline of previous works and how our proposed
method directly connects to those. Fourthly, we present our results and
findings through experimental evaluation of the proposed method on various
architectures and datasets before we finally conclude our work.

Figure 1: Input image with exemplary pixel values, filters, and corresponding
output with point-estimates (left) and probability distributions (right) over
weights.

##  2 Bayesian convolutional neural networks with variational inference

Recently, the uncertainty afforded by Bayes by Backprop trained neural
networks has been used successfully to train feedforward neural networks in
both supervised and reinforcement learning environments [5, 7, 8], for
training recurrent neural networks [9], and for Cnns [10, 11]. Here, we review
this method for Cnns to construct a common foundation on which be build on in
section 3.

###  2.1 Related work

Applying Bayesian methods to neural networks has been studied in the past with
various approximation methods for the intractable true posterior probability
distribution $p(w|\mathcal{D})$. [12] started to propose various maximum-a-
posteriori (MAP) schemes for neural networks. They were also the first who
suggested second order derivatives in the prior probability distribution
$p(w)$ to encourage smoothness of the resulting approximate posterior
probability distribution. In subsequent work by [13], the first variational
methods were proposed which naturally served as a regularizer in neural
networks. [14] suggest taking an information theory perspective into account
and utilizing a minimum description length (MDL) loss. This penalises non-
robust weights by means of an approximate penalty based upon perturbations of
the weights on the outputs. [15] and [16] investigated the posterior
probability distributions of neural networks by using Laplace approximations.
As a response to the limitations of Laplace approximations, [17] investigated
the use of hybrid Monte Carlo for training neural networks, although it has so
far been difficult to apply these to the large sizes of neural networks built
in modern applications. More recently, [4] derived a variational inference
scheme for neural networks and [5] extended this with an update for the
variance that is unbiased and simpler to compute. [18] derives a similar
algorithm in the case of a mixture posterior probability distribution.
Several authors have derived how Dropout [19] and Gaussian Dropout [20] can be
viewed as approximate variational inference schemes [21, 22]. We compare our
results to [21]. Furthermore, structured variational approximations [23],
auxiliary variables [24], and Stochastic Gradient MCMC [25] have been proposed
to approximate the intractable posterior probability distribution. Recently,
the uncertainty afforded by Bayes by Backprop trained neural networks has been
used successfully to train feedforward neural networks in both supervised and
reinforcement learning environments [5, 7, 8], for training recurrent neural
networks [9], and convolutional neural networks [10, 11].

###  2.2 Bayes by Backprop

Bayes by Backprop [4, 5] is a variational inference method to learn the
posterior distribution on the weights $w\sim q_{\theta}(w|\mathcal{D})$ of a
neural network from which weights $w$ can be sampled in backpropagation. Since
the true posterior is typically intractable, an approximate distribution
$q_{\theta}(w|\mathcal{D})$ is defined that is aimed to be as similar as
possible to the true posterior $p(w|\mathcal{D})$, measured by the Kullback-
Leibler (KL) divergence [26]. Hence, we define the optimal parameters
$\theta^{opt}$ as

$\displaystyle\theta^{opt}$ $\displaystyle=\underset{\theta}{\text{arg min}}\
\text{KL}\ [q_{\theta}(w|\mathcal{D})\|p(w|\mathcal{D})]$ (1)
$\displaystyle=\underset{\theta}{\text{arg min}}\ \text{KL}\
[q_{\theta}(w|\mathcal{D})\|p(w)]-\mathbb{E}_{q(w|\theta)}[\log
p(\mathcal{D}|w)]+\log p(\mathcal{D})$

where

$\text{KL}\ [q_{\theta}(w|\mathcal{D})\|p(w)]=\int
q_{\theta}(w|\mathcal{D})\log\frac{q_{\theta}(w|\mathcal{D})}{p(w)}dw.$ (2)

This derivation forms an optimization problem with a resulting cost function
widely known as variational free energy [27, 28, 29] which is built upon two
terms: the former, $\text{KL}\ [q_{\theta}(w|\mathcal{D})\|p(w)]$, is
dependent on the definition of the prior $p(w)$, thus called complexity cost,
whereas the latter, $\mathbb{E}_{q(w|\theta)}[\log p(\mathcal{D}|w)]$, is
dependent on the data $p(\mathcal{D}|w)$, thus called likelihood cost. The
term $\log p(\mathcal{D})$ can be omitted in the optimization because it is
constant.
Since the KL-divergence is also intractable to compute exactly, we follow a
stochastic variational method [4, 5]. We sample the weights $w$ from the
variational distribution $q_{\theta}(w|\mathcal{D})$ since it is much more
probable to draw samples which are appropriate for numerical methods from the
variational posterior $q_{\theta}(w|\mathcal{D})$ than from the true posterior
$p(w|\mathcal{D})$. Consequently, we arrive at the tractable cost function (3)
which is aimed to be optimized, i.e. minimised with respect to $\theta$,
during training:

$\mathcal{F}(\mathcal{D},\theta)\approx\sum_{i=1}^{n}\log
q_{\theta}(w^{(i)}|\mathcal{D})-\log p(w^{(i)})-\log p(\mathcal{D}|w^{(i)})$
(3)

where $n$ is the number of draws. We sample $w^{(i)}$ from
$q_{\theta}(w|\mathcal{D})$.

###  2.3 Bayes by Backprop for convolutional neural networks

In this section, we explain our algorithm of building Cnns with probability
distributions over weights in each filter, as seen in Figure 1, and apply
Bayes by Backprop to compute the intractable posterior probability
distributions $p(w|\mathcal{D})$, as described in the previous section 2.2.
Notably, a fully Bayesian perspective on a Cnn is for most Cnn architectures
not accomplished by merely placing probability distributions over weights in
convolutional layers; it also requires probability distributions over weights
in fully-connected layers.

####  2.3.1 Local reparameterization trick for convolutional layers

We utilize the local reparameterization trick [22] and apply it to Cnns.
Following [22, 11], we do not sample the weights $w$, but we sample layer
activations $b$ instead due to its consequent computational acceleration. The
variational posterior probability distribution
$q_{\theta}(w_{ijhw}|\mathcal{D})=\mathcal{N}(\mu_{ijhw},\alpha_{ijhw}\mu^{2}_{ijhw})$
(where $i$ and $j$ are the input, respectively output layers, $h$ and $w$ the
height, respectively width of any given filter) allows to implement the local
reparamerization trick in convolutional layers. This results in the subsequent
equation for convolutional layer activations $b$:

$b_{j}=A_{i}\ast\mu_{i}+\epsilon_{j}\odot\sqrt{A^{2}_{i}\ast(\alpha_{i}\odot\mu^{2}_{i})}$
(4)

where $\epsilon_{j}\sim\mathcal{N}(0,1)$, $A_{i}$ is the receptive field,
$\ast$ signalises the convolutional operation, and $\odot$ the component-wise
multiplication.

####  2.3.2 Applying two sequential convolutional operations (mean and
variance)

The crux of equipping a Cnn with probability distributions over weights
instead of single point-estimates and being able to update the variational
posterior probability distribution $q_{\theta}(w|\mathcal{D})$ by
backpropagation lies in applying two convolutional operations whereas filters
with single point-estimates apply one. As explained in the previous section
2.3.1, we deploy the local reparametrization trick and sample from the
activations $b$. Since activations $b$ are functions of mean $\mu_{ijwh}$ and
variance $\alpha_{ijhw}\mu^{2}_{ijhw}$ among others, we are able to compute
the two variables determining a Gaussian probability distribution, mean
$\mu_{ijhw}$ and variance $\alpha_{ijhw}\mu^{2}_{ijhw}$, separately.
We pursue this in two convolutional operations: in the first, we treat the
output $b$ as an output of a Cnn updated by frequentist inference. We optimize
with Adam [30] towards a single point-estimate which makes the accuracy of
classifications in the validation dataset increasing. We interpret this single
point-estimate as the mean $\mu_{ijwh}$ of the variational posterior
probability distributions $q_{\theta}(w|\mathcal{D})$. In the second
convolutional operation, we learn the variance $\alpha_{ijhw}\mu^{2}_{ijhw}$.
As this formulation of the variance includes the mean $\mu_{ijwh}$, only
$\alpha_{ijhw}$ needs to be learned in the second convolutional operation
[31]. In this way, we ensure that only one parameter is updated per
convolutional operation, exactly how it would have been with a Cnn updated by
frequentist inference.
In other words, while we learn in the first convolutional operation the MAP of
the variational posterior probability distribution
$q_{\theta}(w|\mathcal{D})$, we observe in the second convolutional operation
how much values for weights $w$ deviate from this MAP. This procedure is
repeated in the fully-connected layers. In addition, to accelerate
computation, to ensure a positive non-zero variance
$\alpha_{ijhw}\mu^{2}_{ijhw}$, and to enhance accuracy, we learn
$\log\alpha_{ijhw}$ and use the Softplus activation function as further
described in the Experiments section 4.

##  3 Uncertainty estimation in Bayesian Cnns

In classification tasks, we are interested in the predictive distribution
$p_{\mathcal{D}}(y^{*}|x^{*})$, where $x^{*}$ is an unseen data example and
$y^{*}$ its predicted class. For a Bayesian neural network, this quantity is
given by:

$\displaystyle p_{\mathcal{D}}(y^{*}|x^{*})=\int p_{w}(y^{*}|x^{*})\
p_{\mathcal{D}}(w)\ dw$ (5)

In Bayes by Backprop, Gaussian distributions
$q_{\theta}(w|\mathcal{D})\sim\mathcal{N}(w|\mu,\sigma^{2})$, where
$\theta=\\{\mu,\sigma\\}$ are learned with some dataset
$\mathcal{D}=\\{x_{i},y_{i}\\}_{i=1}^{n}$ as we explained previously in 2.2.
Due to the discrete and finite nature of most classification tasks, the
predictive distribution is commonly assumed to be a categorical. Incorporating
this aspect into the predictive distribution gives us

$\displaystyle p_{\mathcal{D}}(y^{*}|x^{*})$
$\displaystyle=\int\text{Cat}(y^{*}|f_{w}(x^{*}))\mathcal{N}(w|\mu,\sigma^{2})\
dw$ (6)
$\displaystyle=\int\prod_{c=1}^{C}f(x_{c}^{*}|w)^{y_{c}^{*}}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(w-\mu)^{2}}{2\sigma^{2}}}\
dw$

where $C$ is the total number of classes and $\sum_{c}f(x_{c}^{*}|w)=1$.
As there is no closed-form solution due to the lack of conjugacy between
categorical and Gaussian distributions, we cannot recover this distribution.
However, we can construct an unbiased estimator of the expectation by sampling
from $q_{\theta}(w|\mathcal{D})$:

$\mathbb{E}_{q}[p_{\mathcal{D}}(y^{*}|x^{*})]=\int q_{\theta}(w|\mathcal{D})\
p_{w}(y|x)\ dw\approx\frac{1}{T}\sum_{t=1}^{T}p_{w_{t}}(y^{*}|x^{*})$ (7)

where $T$ is the predefined number of samples. This estimator allows us to
evaluate the uncertainty of our predictions by the definition of variance,
hence called predictive variance and denoted as $\text{Var}_{q}$:

$\displaystyle\text{Var}_{q}\big{(}p(y^{*}|x^{*})\big{)}=\mathbb{E}_{q}[y^{*}y^{*T}]-\mathbb{E}_{q}[y^{*}]\mathbb{E}_{q}[y^{*}]^{T}$
(8)
$\displaystyle=\underbrace{\int\Big{[}\text{diag}\Big{(}\mathbb{E}_{p}[y^{*}]\Big{)}-\mathbb{E}_{p}[y^{*}]\
\mathbb{E}_{p}[y^{*}]^{T}\Big{]}\ q_{\theta}(w|\mathcal{D})\
dw}_{\text{aleatoric}}$
$\displaystyle+\underbrace{\int\Big{(}\mathbb{E}_{p}[y^{*}]-\mathbb{E}_{q}[y^{*}]\Big{)}\
\Big{(}\mathbb{E}_{p}[y^{*}]-\mathbb{E}_{q}[y^{*}]\Big{)}^{T}q_{\theta}(w|\mathcal{D})\
dw}_{\text{epistemic}}$

where $\mathbb{E}_{p}[y^{*}]=\mathbb{E}_{p(y^{*}|x^{*})}[y^{*}]$ and
$\mathbb{E}_{q}[y^{*}]=\mathbb{E}_{q_{\theta}(y^{*}|x^{*})}[y^{*}]$.
It is important to split the uncertainty in form of the predictive variance
into aleatoric and epistemic quantities since it allows the modeler to
evaluate the room for improvements: while aleatoric uncertainty (also known as
statistical uncertainty) is merely a measure for the variation of ("noisy")
data, epistemic uncertainty is caused by the model. Hence, a modeler can see
whether the quality of the data is low (i.e. high aleatoric uncertainty), or
the model itself is the cause for poor performances (i.e. high epistemic
uncertainty). The former can be improved by gathering more data, whereas the
latter requests to refine the model [32].

###  3.1 Related work

As shown in (8), the predictive variance can be decomposed into the aleatoric
and epistemic uncertainty. [33] and [6] have proposed disparate methods to do
so. In this paragraph, we first review these two methods before we explain our
own method and how it overcomes deficiencies of the two aforementioned ones,
especially the usage of Softmax in the output layer by [6].
[33] derived a method how aleatoric and epistemic uncertainties can directly
be estimated by constructing a Bayesian neural network with the last layer
before activation consisting of mean and variance of logits, denoted
$\hat{\mu}_{t}$ and $\hat{\sigma}_{t}^{2}$. In other words, the pre-activated
linear output of the neural network has $2K$ dimensions with $K$ being the
number of output units, i.e. potential classes. They propose an estimator as

$\text{Var}_{q}\big{(}p(y^{*}|x^{*})\big{)}=\frac{1}{T}\sum_{t=1}^{T}\text{diag}(\hat{\sigma}_{t}^{2})+\frac{1}{T}\sum_{t}=1^{T}(\hat{\mu}_{t}-\bar{\mu})\
(\hat{\mu}_{t}-\bar{\mu})^{T}$ (9)

where $\bar{\mu}=\sum_{t=1}^{T}\hat{\mu}_{t}/T$. [6] mention the deficiencies
of this approach: first, it models the variability of the linear predictors
$\hat{\mu}_{t}$ and $\hat{\sigma}_{t}^{2}$ and not the predictive
probabilities; second, it ignores the fact that the covariance matrix of a
multinomial random variable is a function of the mean vector; and third, the
aleatoric uncertainty does not reflect correlations because of the diagonal
matrix.
To overcome these deficiencies, [6] propose

$\text{Var}_{q}\big{(}p(y^{*}|x^{*})\big{)}=\underbrace{\frac{1}{T}\sum_{t=1}^{T}\text{diag}(\hat{p}_{t})-\hat{p}_{t}\
\hat{p}_{t}^{T}}_{\text{aleatoric}}+\underbrace{\frac{1}{T}\sum_{t=1}^{T}(\hat{p}_{t}-\bar{p})(\hat{p}_{t}-\bar{p})^{T}}_{\text{epistemic}}$
(10)

where $\bar{p}=\frac{1}{T}\sum_{t=1}^{T}\hat{p}_{t}$ and
$\hat{p}_{t}=\text{Softmax}\big{(}f_{w_{t}}(x^{*})\big{)}$. By doing so, they
do not need the pre-activated linear outputs $\hat{\mu}_{t}$ and
$\hat{\sigma}_{t}^{2}$. Consequently, they can directly compute the
variability of the predictive probability and do not require additional
sampling steps. With increasing $T$, (10) converges in probability to (8).

###  3.2 Softplus normalization

Deploying the Softmax function assumes we classify according to the logistic
sigmoid function:

$\displaystyle a$ $\displaystyle=\text{log}\
\frac{p(c=1|x)}{p(c=2|x)}=\frac{p(x|c=1)\ p(c=1)}{p(x|c=2)\ p(c=2)}$ (11)
$\displaystyle\sigma(a)$ $\displaystyle=\frac{1}{1+\exp{(-a)}}=p(c=1|x)$ (12)

where the label $c=1$ is assigned to the input $x$ if
$\sigma(a)=p(c=1|x)=p(y=1|x^{*})\geq 0.5$. Assume we have $C$ classes, this
logistic sigmoid function can be generalized to the Softmax function:

$p(c|x)=\frac{p(x|c)\ p(c)}{\sum_{c=1}^{C}p(x|c)\ p(c)}$ (13)

Here, we propose a novel method how aleatoric and epistemic uncertainty
estimations can be computed without having an additional non-linear
transformation by the generalized logistic sigmoid function in the output
layer.
The predictive variance remains in the form as stated in (10), except the
subsequent two amendments compared to [6]: we circumvent the exponential term
in the Softmax function, hence generate consistency of activation functions in
the entire Cnn. Despite not punishing wrong predictions particularly with the
exponential term, we also resolve concerns about robustness by employing not
more than one type of non-linear transformations within the same neural
network while we preserve the desideratum of having probabilities as outputs.
We implement this in two subsequent computations: first, we employ the same
activation function in the output layer as we employ in any other hidden
layer, namely the Softplus function, see (16); second, these outputs are
normalized by dividing the output of each class by the sum of all outputs:

$y^{*}_{norm}=\frac{y^{*}_{c}}{\sum_{C}y^{*}_{c}}$ (14)

This procedure can be summarized as a replacement of
$\hat{p}_{t}=\text{Softmax}\big{(}f_{w_{t}}(x^{*})\big{)}$ with
$\hat{p}_{t}=\text{Softplus}_{\text{n}}\big{(}f_{w_{t}}(x^{*})\big{)}$ where
$\text{Softplus}_{\text{n}}$ is the aforementioned normalization of the
Softplus output. We call this computation Softplus normalization. Our
intuition for this is as follows: we classify according to a categorical
distribution as we outlined previously in section 3. This can be seen as an
approximation of the outputs with one-hot vectors. If those vectors already
contain a lot of zeros, the approximation is going to be more accurate. For a
Softmax, predicting zeros requires a logit of $-\infty$ which is hard to
achieve in practice. For the Softplus normalization, we just need a negative
value roughly smaller than $-4$ to get nearly a zero in the output.
Consequently, Softplus normalization can easily produce vectors that are in
practice zero, whereas Softmax cannot. In sum, our proposed method is as
follows:

$\text{Var}_{q}\big{(}p(y^{*}|x^{*})\big{)}=\underbrace{\frac{1}{T}\sum_{t=1}^{T}\text{diag}(\hat{p}_{t})-\hat{p}_{t}\
\hat{p}_{t}^{T}}_{\text{aleatoric}}+\underbrace{\frac{1}{T}\sum_{t=1}^{T}(\hat{p}_{t}-\bar{p})(\hat{p}_{t}-\bar{p})^{T}}_{\text{epistemic}}$
(15)

where $\bar{p}=\frac{1}{T}\sum_{t=1}^{T}\hat{p}_{t}$ and
$\hat{p}_{t}=\text{Sofplus}_{\text{n}}\big{(}f_{w_{t}}(x^{*})\big{)}$.

##  4 Experiments

For all conducted experiments, we implement the foregoing description of
Bayesian Cnns with variational inference in LeNet-5 [34], AlexNet [3], and VGG
[2]. The exact architecture specifications can be found in the Appendix 6 and
in our GitHub repository111https://github.com/kumar-shridhar/PyTorch-Softplus-
Normalization-Uncertainty-Estimation-Bayesian-CNN. We train the networks with
the MNIST dataset of handwritten digits [34], and with the CIFAR-10 and
CIFAR-100 datasets [35] since these datasets serve widely as benchmarks for
Cnns’ performances. The originally chosen activation functions in all
architectures are ReLU, but we must introduce another, called Softplus, see
(16), because of our method to apply two convolutional or fully-connected
operations. As aforementioned, one of these is determining the mean $\mu$, and
the other the variance $\alpha\mu^{2}$. Specifically, we apply the Softplus
function because we want to ensure that the variance $\alpha\mu^{2}$ never
becomes zero. This would be equivalent to merely calculating the MAP, which
can be interpreted as equivalent to a maximum likelihood estimation (MLE),
which is further equivalent to utilising single point-estimates, hence
frequentist inference. The Softplus activation function is a smooth
approximation of ReLU. Although it is practically not influential, it has the
subtle and analytically important advantage that it never becomes zero for
$x\rightarrow-\infty$, whereas ReLU becomes zero for $x\rightarrow-\infty$.

$\text{Softplus}(x)=\frac{1}{\beta}\cdot\log\big{(}1+\exp(\beta\cdot
x)\big{)}$ (16)

where $\beta$ is by default set to $1$.
All experiments are performed with the same hyper-parameters settings as
stated in the Appendix 6.

###  4.1 Results of Bayesian CNNs with variational inference

| MNIST | CIFAR-10 | CIFAR-100
---|---|---|---
Bayesian VGG (with VI) | 99 | 86 | 45
Frequentist VGG | 99 | 85 | 48
Bayesian AlexNet (with VI) | 99 | 73 | 36
Frequentist AlexNet | 99 | 73 | 38
Bayesian LeNet-5 (with VI) | 98 | 69 | 31
Frequentist LeNet-5 | 98 | 68 | 33
Bayesian LeNet-5 (with Dropout) | 99 | 83 |

Table 1: Comparison of validation accuracies (in percentage) for different
architectures with variational inference (VI), frequentist inference and
Dropout as a Bayesian approximation as proposed by Gal and Ghahramani [21] for
MNIST, CIFAR-10, and CIFAR-100.

We evaluate the performance of our Bayesian Cnns with variational inference.
Table 1 shows a comparison of validation accuracies (in percentage) for
architectures trained by two disparate Bayesian approaches, namely variational
inference, i.e. Bayes by Backprop and Dropout as proposed by Gal and
Ghahramani [21], plus frequentist inference for all three datasets. Bayesian
Cnns trained by variational inference achieve validation accuracies comparable
to their counter-architectures trained by frequentist inference. On MNIST,
validation accuracies of the two disparate Bayesian approaches are comparable,
but a Bayesian LeNet-5 with Dropout achieves a considerable higher validation
accuracy on CIFAR-10, although we were not able to reproduce these reported
results.
In Figure 2, we show how Bayesian networks incorporate naturally effects of
regularization, exemplified on AlexNet. While an AlexNet trained by
frequentist inference without any regularization overfits greatly on
CIFAR-100, an AlexNet trained by Bayesian inference on CIFAR-100 does not. It
performs equivalently to an AlexNet trained by frequentist inference with
three layers of Dropout after the first, fourth, and sixth layers in the
architecture. In initial epochs, Bayesian Cnns trained by variational
inference start with a low validation accuracy compared to architectures
trained by frequentist inference. Initialization for both inference methods is
chosen equivalently: the variational posterior probability distributions
$q_{\theta}(w|\mathcal{D})$ is initially approximated as standard Gaussian
distributions, while initial point-estimates in architectures trained by
frequentist inference are randomly drawn from a standard Gaussian
distribution. These initialization methods ensure weights are neither too
small nor too large in the beginning of training.

Figure 2: AlexNet trained on CIFAR-100 by Bayesian and frequentist inference.
The frequentist AlexNet without Dropout overfits while the Bayesian AlexNet
naturally incorporates an effect of regularization, comparable to a
frequentist AlexNet with three Dropout layers.

###  4.2 Results of uncertainty estimations by Softplus normalization

Our results for the Softplus normalization estimation method of aleatoric and
epistemic uncertainties are combined in Table 2. It compares the means over
epochs of aleatoric and epistemic uncertainties for our Bayesian Cnns LeNet-5,
AlexNet and VGG. We train on MNIST and CIFAR-10, and compute homoscedastic
uncertainties by averaging over all classes, i.e. over (heteroscedastic)
uncertainties of each class. We see a correlating pattern between validation
accuracies and epistemic uncertainty: with increasing validation accuracy,
epistemic uncertainty decreases, observable by the different models. In
contrast, aleatoric uncertainty measures the irreducible variability of the
datasets, hence is only dependent on the datasets and not the models, which
can be seen by the constant aleatoric uncertainties of each dataset across
models.

| Aleatoric uncertainty | Epistemic uncertainty | Validation accuracy
---|---|---|---
Bayesian VGG (MNIST) | 0.00110 | 0.0004 | 99
Bayesian VGG (CIFAR-10) | 0.00099 | 0.0013 | 85
Bayesian AlexNet (MNIST) | 0.00110 | 0.0019 | 99
Bayesian AlexNet (CIFAR-10) | 0.00099 | 0.0002 | 73
Bayesian LeNet-5 (MNIST) | 0.00110 | 0.0026 | 98
Bayesian LeNet-5 (CIFAR-10) | 0.00099 | 0.0404 | 69

Table 2: Aleatoric and epistemic uncertainty for Bayesian VGG, AlexNet and
LeNet-5 calculated for MNIST and CIFAR-10, computed by Softplus normalization
(3.2). Validation accuracy is displayed to demonstrate the negative
correlation between validation accuracy and epistemic uncertainty.

We further investigate influences on aleatoric uncertainty estimations by
additive standard Gaussian noise of different levels. Here, a random sample of
the standard Gaussian distribution is multiplied by a level $\gamma$ and added
to each pixel value. We test for $\gamma=0,0.1,0.2,0.3$ and see that aleatoric
uncertainty is independent of the added noise level (see Figure 3), i.e.
aleatoric uncertainty is constant to six decimal places. Our intuition for
this phenomenon is as follows: by normalizing the output of the Softplus
function among one batch, of which all images have the same added level of
noise, the aleatoric uncertainty captures the variability among the images in
one batch - and not across batches with different levels of noise.

$\gamma=$

$\text{aleatoric}=$

$0$

$0.00099897$

$0.1$

$0.00099885$

$0.2$

$0.00099875$

$0.3$

$0.00099838$

Figure 3: Aleatoric uncertainty is computed with Bayesian VGG on the same
MNIST input image with different noise levels $\gamma$.

##  5 Conclusion

We propose a novel method how aleatoric and epistemic uncertainties can be
estimated in Bayesian Cnns. We call this method Softplus normalization.
Firstly, we discuss briefly the differences of a Bayes by Backprop based
variational inference to the Dropout based approximation of Gal & Ghahramani
[21]. Secondly, we derive the Softplus normalization uncertainty estimation
method with previously published methods for these uncertainties. We evaluate
our approach on three datasets (MNIST, CIFAR-10, CIFAR-100).
As a base, we show that Bayesian Cnns with variational inference achieve
comparable results as those achieved by the same network architectures trained
by frequentist inference, but include naturally a regularization effect and an
uncertainty measure. Building on that, we examine how our proposed Softplus
normalization method to estimate aleatoric and epistemic uncertainties is
derived by previous work in this field, how it differs from those, and explain
why it is more appropriate to use in computer vision than these previously
proposed methods.

## References

  * [1] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [2] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [3] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in neural information processing systems, pages 1097–1105, 2012.
  * [4] Alex Graves.  Practical variational inference for neural networks.  In Advances in Neural Information Processing Systems, pages 2348–2356, 2011.
  * [5] Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra.  Weight uncertainty in neural networks.  arXiv preprint arXiv:1505.05424, 2015.
  * [6] Yongchan Kwon, Joong-Ho Won, Beom Joon Kim, and Myunghee Cho Paik.  Uncertainty quantification using bayesian neural networks in classification: Application to ischemic stroke lesion segmentation.  2018\.
  * [7] Zachary C Lipton, Jianfeng Gao, Lihong Li, Xiujun Li, Faisal Ahmed, and Li Deng.  Efficient exploration for dialogue policy learning with bbq networks & replay buffer spiking.  arXiv preprint arXiv:1608.05081, 2016.
  * [8] Rein Houthooft, Xi Chen, Yan Duan, John Schulman, Filip De Turck, and Pieter Abbeel.  Curiosity-driven exploration in deep reinforcement learning via bayesian neural networks.  arXiv preprint arxiv.1605.09674, 2016.
  * [9] Meire Fortunato, Charles Blundell, and Oriol Vinyals.  Bayesian recurrent neural networks.  arXiv preprint arXiv:1704.02798, 2017.
  * [10] Kumar Shridhar, Felix Laumann, and Marcus Liwicki.  A comprehensive guide to bayesian convolutional neural network with variational inference.  CoRR, abs/1901.02731, 2019.
  * [11] Kirill Neklyudov, Dmitry Molchanov, Arsenii Ashukha, and Dmitry Vetrov.  Variance networks: When expectation does not meet your expectations.  arXiv preprint arXiv:1803.03764, 2018.
  * [12] Wray L Buntine and Andreas S Weigend.  Bayesian back-propagation.  Complex systems, 5(6):603–643, 1991.
  * [13] Geoffrey E Hinton and Drew Van Camp.  Keeping the neural networks simple by minimizing the description length of the weights.  In Proceedings of the sixth annual conference on Computational learning theory, pages 5–13. ACM, 1993.
  * [14] Sepp Hochreiter and Jürgen Schmidhuber.  Simplifying neural nets by discovering flat minima.  In Advances in neural information processing systems, pages 529–536, 1995.
  * [15] John S Denker and Yann LeCu.  Transforming neural-net output levels to probability distributions.  In Advances in neural information processing systems, pages 853–859, 1991.
  * [16] David JC MacKay.  Probable networks and plausible predictions—a review of practical bayesian methods for supervised neural networks.  Network: Computation in Neural Systems, 6(3):469–505, 1995.
  * [17] Radford M Neal.  Bayesian learning for neural networks, volume 118.  Springer Science & Business Media, 2012.
  * [18] Alex Graves.  Stochastic backpropagation through mixture density distributions.  arXiv preprint arXiv:1607.05690, 2016.
  * [19] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  The Journal of Machine Learning Research, 15(1):1929–1958, 2014\.
  * [20] Sida Wang and Christopher Manning.  Fast dropout training.  In international conference on machine learning, pages 118–126, 2013.
  * [21] Yarin Gal and Zoubin Ghahramani.  Bayesian convolutional neural networks with bernoulli approximate variational inference.  arXiv preprint arXiv:1506.02158, 2015.
  * [22] Diederik P Kingma, Tim Salimans, and Max Welling.  Variational dropout and the local reparameterization trick.  In Advances in Neural Information Processing Systems, pages 2575–2583, 2015.
  * [23] Christos Louizos and Max Welling.  Multiplicative normalizing flows for variational bayesian neural networks.  In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 2218–2227. JMLR. org, 2017.
  * [24] Lars Maaløe, Casper Kaae Sønderby, Søren Kaae Sønderby, and Ole Winther.  Auxiliary deep generative models.  arXiv preprint arXiv:1602.05473, 2016.
  * [25] Chunyuan Li, Andrew Stevens, Changyou Chen, Yunchen Pu, Zhe Gan, and Lawrence Carin.  Learning weight uncertainty with stochastic gradient mcmc for shape classification.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5666–5675, 2016.
  * [26] Solomon Kullback and Richard A Leibler.  On information and sufficiency.  The annals of mathematical statistics, 22(1):79–86, 1951.
  * [27] Radford M Neal and Geoffrey E Hinton.  A view of the em algorithm that justifies incremental, sparse, and other variants.  In Learning in graphical models, pages 355–368. Springer, 1998\.
  * [28] Jonathan S Yedidia, William T Freeman, and Yair Weiss.  Constructing free-energy approximations and generalized belief propagation algorithms.  IEEE Transactions on information theory, 51(7):2282–2312, 2005\.
  * [29] Karl Friston, Jérémie Mattout, Nelson Trujillo-Barreto, John Ashburner, and Will Penny.  Variational free energy and the laplace approximation.  Neuroimage, 34(1):220–234, 2007.
  * [30] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [31] Dmitry Molchanov, Arsenii Ashukha, and Dmitry Vetrov.  Variational dropout sparsifies deep neural networks.  arXiv preprint arXiv:1701.05369, 2017.
  * [32] Armen Der Kiureghian and Ove Ditlevsen.  Aleatory or epistemic? does it matter?  Structural Safety, 31(2):105–112, 2009.
  * [33] Alex Kendall and Yarin Gal.  What uncertainties do we need in bayesian deep learning for computer vision?  In Advances in neural information processing systems, pages 5574–5584, 2017.
  * [34] Yann LeCun, Leon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [35] Alex Krizhevsky and Geoffrey Hinton.  Learning multiple layers of features from tiny images.  Technical report, Citeseer, 2009.

##  6 Appendix

###  6.1 Experiment specifications

variable | value
---|---
learning rate | 0.001
epochs | 100
batch size | 128
sample size | 10
$(\alpha\mu^{2})_{init}$ of approximate posterior $q_{\theta}(w|\mathcal{D})$ | -10
optimizer |  Adam [30]
$\lambda$ in $\ell$-2 normalization | 0.0005
$\beta_{i}$ |  $\frac{2^{M-i}}{2^{M}-1}$ [5]

###  6.2 Model architectures

####  6.2.1 LeNet-5

layer type | width | stride | padding | input shape | nonlinearity
---|---|---|---|---|---
convolution ($5\times 5$) | 6 | 1 | 0 | $M\times 1\times 32\times 32$ | Softplus
Mmax-pooling ($2\times 2$) |  | 2 | 0 | $M\times 6\times 28\times 28$ |
convolution ($5\times 5$) | 16 | 1 | 0 | $M\times 1\times 14\times 14$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 16\times 10\times 10$ |
fully-connected | 120 |  |  | $M\times 400$ | Softplus
fully-connected | 84 |  |  | $M\times 120$ | Softplus
fully-connected | 10 |  |  | $M\times 84$ | Softplus normalization

####  6.2.2 AlexNet

layer type | width | stride | padding | input shape | nonlinearity
---|---|---|---|---|---
convolution ($11\times 11$) | 64 | 4 | 5 | $M\times 3\times 32\times 32$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 64\times 32\times 32$ |
convolution ($5\times 5$) | 192 | 1 | 2 | $M\times 64\times 15\times 15$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 192\times 15\times 15$ |
convolution ($3\times 3$) | 384 | 1 | 1 | $M\times 192\times 7\times 7$ | Softplus
convolution ($3\times 3$) | 256 | 1 | 1 | $M\times 384\times 7\times 7$ | Softplus
convolution ($3\times 3$) | 128 | 1 | 1 | $M\times 256\times 7\times 7$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 128\times 7\times 7$ |
fully-connected | 128 |  |  | $M\times 128$ | Softplus normalization

####  6.2.3 VGG

layer type | width | stride | padding | input shape | nonlinearity
---|---|---|---|---|---
convolution ($3\times 3$) | 64 | 1 | 1 | $M\times 3\times 32\times 32$ | Softplus
convolution ($3\times 3$) | 64 | 1 | 1 | $M\times 64\times 32\times 32$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 64\times 32\times 32$ |
convolution ($3\times 3$) | 128 | 1 | 1 | $M\times 64\times 16\times 16$ | Softplus
convolution ($3\times 3$) | 128 | 1 | 1 | $M\times 128\times 16\times 16$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 128\times 16\times 16$ |
convolution ($3\times 3$) | 256 | 1 | 1 | $M\times 128\times 8\times 8$ | Softplus
convolution ($3\times 3$) | 256 | 1 | 1 | $M\times 256\times 8\times 8$ | Softplus
convolution ($3\times 3$) | 256 | 1 | 1 | $M\times 256\times 8\times 8$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 256\times 8\times 8$ |
convolution ($3\times 3$) | 512 | 1 | 1 | $M\times 256\times 4\times 4$ | Softplus
convolution ($3\times 3$) | 512 | 1 | 1 | $M\times 512\times 4\times 4$ | Softplus
convolution ($3\times 3$) | 512 | 1 | 1 | $M\times 512\times 4\times 4$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 512\times 4\times 4$ |
convolution ($3\times 3$) | 512 | 1 | 1 | $M\times 512\times 2\times 2$ | Softplus
convolution ($3\times 3$) | 512 | 1 | 1 | $M\times 512\times 2\times 2$ | Softplus
convolution ($3\times 3$) | 512 | 1 | 1 | $M\times 512\times 2\times 2$ | Softplus
max-pooling ($2\times 2$) |  | 2 | 0 | $M\times 512\times 2\times 2$ |
fully-connected | 512 |  |  | $M\times 512$ | Softplus normalization
  *[Cnns]: Convolutional neural network
  *[Cnn]: Convolutional neural network


masunaga1998MMI lowlyndenbell1976LLB masunaga:opacity.limit.tauMI99

††thanks: $\dagger$NASA Hubble Fellow

# The opacity limit

Michael Y. Grudić1,† Philip F. Hopkins2 1Carnegie Observatories, 813 Santa
Barbara St, Pasadena, CA 91101, USA  2TAPIR, Mailcode 350-17, California
Institute of Technology, Pasadena, CA 91125, USA  Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

The opacity limit is an important concept in star formation: isothermal
collapse cannot proceed without limit, because eventually cooling radiation is
trapped and the temperature rises quasi-adiabatically, setting a minimum Jeans
mass $M_{\rm J}^{\rm min}$. Various works have considered this scenario and
derived expressions for $M_{\rm J}^{\rm min}$, generally $\sim
10^{-3}-10^{-2}M_{\odot}$ in normal star-forming conditions, but with
conflicting results about the scaling with ambient conditions and material
properties. We derive expressions for the thermal evolution of dust-cooled
collapsing gas clumps in various limiting cases, given a general ambient
radiation field ($u_{\rm rad}$, $T_{\rm rad}$) and a general power-law dust
opacity law $\sigma_{\rm d}=A_{\rm d}T^{\beta}$. By accounting for temperature
evolution self-consistently we rule out a previously-proposed regime in which
the adiabatic transition occurs while the core is still optically-thin. If the
radiation field is weak or dust opacity is small, $M_{\rm J}^{\rm min}$ is
insensitive to dust properties/abundance ($\sim A_{\rm
d}^{-\frac{1}{11}}-A_{\rm d}^{-\frac{1}{15}}$), but if the radiation field is
strong and dust is abundant it scales $\propto A_{\rm d}^{1/3}$. This could
make the IMF less bottom-heavy in dust-rich and/or radiation-dense
environments, e.g. galactic centers, starburst galaxies, massive high-$z$
galaxies, and proto-star clusters that are already luminous.

###### keywords:

stars: mass function – galaxies: star clusters: general – stars: formation

##  1 Introduction

We live in a hierarchical universe, consisting of nested structures spanning a
vast range of scales from the cosmic web down to to moons and asteroids. It is
natural for inhabitants of such a universe to ask: where does structure
formation stop? Gravity is scale-free but other processes are not, so surely
at times the laws of physics must step in and say no more, preventing some
structures from spawning substructures (Hoyle, 1953).

This question is particular relevant for understanding star formation and the
initial mass function (IMF). Stars in the present-day Universe form by
fragmentation and accretion within cold, magnetized, supersonically-turbulent
giant molecular clouds (GMCs) that are roughly isothermal over a wide range of
densities. The equations of self-gravitating isothermal magnetohydrodynamics
are scale-free (McKee & Ostriker, 2007), and hence do not unambiguously imply
that GMCs should fragment to objects of stellar mass instead of sub-
fragmenting to smaller objects (Martel et al., 2006; Krumholz, 2014; Guszejnov
et al., 2018). The Jeans mass,

$M_{\rm J}=\frac{\uppi^{5/2}}{6}\left(\frac{k_{\rm B}T}{Gm_{\rm
p}\mu}\right)^{3/2}\rho^{-1/2},$ (1)

can become arbitrarily small if the density $\rho$ continues to increase at
fixed $T$ (assuming $\mu\approx 2.4$ for molecular gas). Why, then, does the
observed IMF almost always have a feature around $\sim 0.1-0.3M_{\odot}$, with
increasing rarity of lower-mass objects?

One possibility is that radiative transfer in the interstellar medium places
restrictions on when collapse and fragmentation can proceed: the opacity
limit. Low & Lynden-Bell (1976) (hereafter ), Rees (1976), and Silk (1977)
solved the balance of heating and cooling in collapsing, opaque, Jeans-mass
clumps, and argued that if optically-thick objects can no longer cool, they
will evolve quasi-adiabatically and the Jeans mass must rise. When essentially
any dust is present, the opacity of dust to its own $\sim 1\,\rm mm$ thermal
emission sets this limit, and the predicted minimum Jeans mass is on the order
of $10^{-3}-10^{-2}\,M_{\odot}$. There is some sensitivity to the precise
hydrodynamics and geometry of the collapse (Boyd & Whitworth, 2005; Whitworth
& Stamatellos, 2006), but weak sensitivity to the opacity: for an opacity law
$\sigma_{\rm d}=A_{\rm d}T^{\beta}$, $M_{\rm J}^{\rm min}\propto A_{\rm
d}^{-\frac{1}{4\beta+7}}$, $\beta\sim 1-2$. This minimum Jeans mass is much
smaller than a star, but the dynamics of accretion after hydrostatic core
formation may have a multiplicative effect on the typical mass that a fragment
eventually accretes (Lee & Hennebelle, 2018a; Hennebelle et al., 2019; Colman
& Teyssier, 2020).

Masunaga & Inutsuka (1999) (hereafter , ) revisited the opacity limit, and
were able to calibrate the uncertain form factors in this theory to a suite of
1D radiation hydrodynamics simulations of protostellar collapse (Masunaga et
al., 1998) (hereafter , ). They corrected the misconception that the condition
of being optically-thick ($\tau\gtrsim 1$) was sufficient to induce adiabatic
evolution, showing that adiabatic evolution begins when the compressional
heating rate exceeds the maximum radiative energy diffusion rate. Their
prediction for $M_{\rm J}^{\rm min}$ is rather more sensitive, with one regime
(their “Case 1”) in which $M_{\rm J}^{\rm min}\propto A_{\rm d}^{-1}$, which
would make it difficult to make low-mass stars in dust-poor environments such
as proto-globular clusters and ultra-faint dwarf galaxies.

We believe the “Case 1” opacity limit of  is unrealistic, and are not aware of
a work that contains the general set of correct expressions for the opacity-
limited density, temperature, and Jeans mass, and related quantities like the
first hydrostatic core (“Larson”) mass. For this reason, we revisit the
subject of dust-cooled protostellar collapse, writing down approximations for
the temperature evolution of collapsing fragments as a function of density in
various asymptotic regimes, and arriving at expressions for the opacity-
limited Jeans mass and the associated transition density $n_{\rm ad}$ that
demarcates the quasi-isothermal and -adiabatic regimes, while correcting some
unrealistic assumptions made in previous works. With these expressions in hand
we then comment briefly on the scaling of the minimum Jeans mass in different
environments, and implications for the IMF.

##  2 Model for dust-cooled collapse

Consider the scenario of a spherical, self-gravitating, thermally-supported
gas clump undergoing approximately self-similar, nearly-isothermal collapse,
as has been studied extensively in previous work (Larson, 1969; Penston, 1969;
Masunaga et al., 1998; Vaytet & Haugbølle, 2017; Bhandare et al., 2018). This
is a highly-idealized scenario, and other works have considered scenarios that
could be more representative of the 3D turbulent collapse process of molecular
clouds (Boyd & Whitworth, 2005; Whitworth & Stamatellos, 2006), but such
details only affect the dimensionless prefactors on the opacity-limited
quantities of interest, which are themselves degenerate with dust properties.
Hence we default to the simple spherical setup, which contains the essential
physics, and which has an exhaustive literature of more-detailed calculations
that can guide our analytic model.

###  2.1 Assumptions

We solve for the thermal evolution of the collapsing core under the following
simplifying assumptions and approximations, all of which follow previous work:

  1. 1.

Dust is the only coolant. This is generally valid at the high ($\gtrsim
10^{5}\,\rm cm^{-3}$) densities characteristic of prestellar collapse, where
dust emission does tend to exceed molecular and fine-structure emission. A
dust-coupled phase should only be absent in extremely metal-poor ($\lesssim
10^{-4}Z_{\odot}$) or metal-pristine star-forming conditions (Sharda &
Krumholz, 2022; Klessen & Glover, 2023).

  2. 2.

The dust temperature $T_{\rm d}$ and gas temperature $T$ are equal. The two
temperatures should asymptotically approach each other because the rate of
dust-gas heat transfer per particle scales $\propto n_{\rm
H}\sqrt{T}\left(T-T_{\rm d}\right)$, and no plausible heating mechanism scales
as steeply with density if $T$ is nearly isothermal, so $T-T_{\rm d}$ will
drop.

  3. 3.

The  assumption that effective shielding length $l=\tau/\sigma n_{\rm H}$ for
determining the optical depth is on the order of the Jeans length
$\lambda_{\rm J}=\sqrt{\uppi c_{\rm s}^{2}/\left(G\rho\right)}$:

$l=f_{\rm J}\lambda_{\rm J},$ (2)

where the numerical prefactor $f_{\rm J}\approx 0.75$ was measured from
simulations by .

  4. 4.

The core has no substructure in velocity, gas density, or dust density. This
is not so clearly justifiable in nature, and is simply our choice of problem
scope. Turbulent/convective motions or clumping of gas and dust could all
conceivably affect the rate of energy transport in general and challenge the
opacity limit scenario, but this is beyond our scope.

  5. 5.

The temperature structure before the adiabatic transition follows from
equating the heating and cooling rate. This is only especially accurate when
$t_{\rm cool}<<t_{\rm dyn}$, but  showed that it is a fair description even
when $t_{\rm cool}\sim t_{\rm dyn}$.

  6. 6.

The core is bathed in an isotropic radiation field specified by an energy
density $u_{\rm rad}$ and an effective absorption-weighted temperature $T_{\rm
rad}$, which need not be in local thermodynamic equilibrium (i.e.
$u_{\mathrm{rad}}\neq aT_{\rm rad}^{4}$) in general. This choice is motivated
by star-forming conditions in the Solar neighborhood, where typically $T_{\rm
rad}\sim 20\rm K$ for the IR component, and $u_{\rm rad}\sim
1\rm\,eV\,cm^{-3}<<aT_{\rm rad}^{4}$. This choice of radiation field and the
dust opacity constitute our parameter space, and the temperature structure
follows from the balance of heating and cooling.

  7. 7.

The dust mass fraction, grain size distribution, and composition remain
constant throughout the collapse. This may be the least-realistic assumption
of current protostellar collapse studies: in dense core conditions many
processes can affect the dust grain properties, with a tendency for grains to
grow and the grain size distribution to become more top-heavy (Weidenschilling
& Ruzmaikina, 1994; Harada et al., 2017; Tu et al., 2022; Lebreuilly et al.,
2023). A transition in the dust properties to a new equilibrium independent of
initial conditions is all but inevitable, because the processes regulating
grain size typically operate on a collisional timescale $\propto
n^{-1}T^{-1/2}$ that will generally overtake the free-fall time $\propto
n^{-1/2}$.

###  2.2 Dust model

Both optically-thin dust emission and absorption depend on the Planck-mean
dust opacity, while optically-thick emission and absorption depend on the
Rosseland mean opacity. Will will ignore the distinction here because they
tend to be close in value, and the key aspect of the problem is how they scale
with temperature, which is similar (Semenov et al., 2003). For dust
temperatures $\lesssim 100\rm K$ the dust cross section per H nucleus may be
approximated as

$\sigma_{\rm d}=A_{\rm d}T^{\beta},$ (3)

where $T$ is either the grain temperature $T_{\rm d}$ that determines the
emission cross-section, or the radiation temperature $T_{\rm rad}$ that
determines the absorption cross-section. The normalization $A_{\rm d}$ rolls
together all possible variation due to the abundance, composition, and grain
size distribution of dust:

$A_{\rm d}=\sigma_{\rm 0}\hat{\sigma}_{\rm d}Z_{\rm
d}\left(10\,\mathrm{K}\right)^{-\beta},$ (4)

where $\sigma_{\rm 0}=4\times 10^{-26}\,\rm cm^{2}\,\rm H^{-1}$ is the dust
cross section per H at $10\,\rm K$ and Solar neighborhood dust abundance,
$Z_{\rm d}$ is the normalized dust mass fraction, and $\hat{\sigma}_{\rm d}$
contains any variations in the intrinsic grain properties relative to the
fiducial model. Dust properties in dense prestellar cores are somewhat
uncertain and may vary considerably, so we will retain arbitrary values in our
expressions. Our fiducial dust model will be $Z_{\rm d}=\hat{\sigma}_{\rm
d}=1$ and $\beta=1.5$, corresponding to model “a” used by , based on a fit to
the results of Pollack et al. (1994). It should be noted that this dust model
lies at the low end of the range of opacity values at $\sim 10\,\rm K$
compared to various other dust models in common use (e.g. Weingartner &
Draine, 2001; Semenov et al., 2003), but on the other hand dust models
calibrated to diffuse ISM diagnostics may not represent the conditions in
dense cores where significant grain growth may have occurred. In general,
$A_{\rm d}$ and $\beta$ must be adjusted here to determine specific prediction
of a given dust model.

###  2.3 Optically-thin cooling regime

Figure 1: Map of the parameter space of density versus temperature (left) or
Jeans mass $M_{\rm J}$ (right) for collapsing cores cooled by dust emission.
Solid curves plot the various simulation results for the central temperature
in the  radiation hydrodynamics simulations, using their fiducial dust opacity
model while varying the initial radiative equilibrium temperature from
$5-30\,\rm K$, and corresponding black dotted lines correspond to the
prediction of the model given in §2. Dashed/dotted lines plot various
thresholds, as labeled, corresponding to Eqs. 10, 8, 14. The shaded regions
are optically-thick ($\tau\gtrsim 1$), assuming our fiducial dust opacity
model (Eq. 3) and a Jeans-scale clump (Eq. 5).

Thermal emission by dust will be optically-thin if the optical depth
$\tau_{\rm out}$ to the emitted radiation is less than unity. For the central
core of the collapsing region, we make the assumption of  that the shielding
length is on the order of the Jeans length , and hence

$\tau_{\rm out}\approx f_{\rm J}\lambda_{\rm J}\sigma_{\rm d}n_{\rm
H}=0.25Z_{\rm d}\hat{\sigma}_{\rm d}\left(\frac{n_{\rm
H}}{10^{10}\,\mathrm{cm}^{-3}}\right)^{1/2}\left(\frac{T_{\rm d}}{10\,\rm
K}\right)^{\frac{2\beta+1}{2}},$ (5)

For a general radiation field where the radiation and dust temperatures are
distinct, it is important to distinguish between $\tau_{\rm out}$ and the
optical depth to ingoing radiation $\tau_{\rm in}$, given by substituting
$T_{\rm rad}$ for $T_{\rm d}$ in the above expression. While $\tau_{\rm
out}<<1$, the dust cooling rate per H nucleus may be written

$\dot{E}_{\mathrm{thin}}=-ac\sigma_{\rm d}T_{\rm d}^{4}$ (6)

where $a$ is the radiation constant and $c$ is the speed of light. While Eq. 6
holds, the gas-dust temperature follows from equating it to the dominant heat
source. Assuming $T\sim T_{\rm d}$, the cooling time in the optically-thin
regime is a function of temperature and dust properties only:

$t_{\rm cool,thin}=\frac{-k_{\rm B}T}{\left(\gamma-1\right)\dot{E}_{\rm
thin}}=700\mathrm{yr}\left(\frac{T}{10\mathrm{K}}\right)^{-\beta-3}\left(Z_{\rm
d}\hat{\sigma}_{\rm d}\right)^{-1}.$ (7)

If at any time $t_{\rm cool}>t_{\rm ff}$, the fragment will virialize and
quasi-isothermal collapse cannot be maintained. The ratio of the two
timescales is

$\frac{t_{\rm cool,thin}}{t_{\rm
ff}}=1.6\left(\frac{T}{10\mathrm{K}}\right)^{-\beta-3}\left(Z_{\rm
d}\hat{\sigma}_{\rm d}\right)^{-1}\left(\frac{n_{\rm H}}{10^{10}\,\rm
cm^{-3}}\right)^{1/2}.$ (8)

Compressional heating (i.e. “$P\mathrm{d}V$ work”) due to gravitational
collapse is the heat source that is always present in this scenario. The
mechanical heating rate per H nucleus of gas collapsing on the free-fall
timescale $\sim\left(G\rho\right)^{-1/2}$ may be written as (; ):

$\dot{E}_{\rm grav}=C_{\rm grav}\sqrt{4\uppi G\rho}\frac{c_{\rm s}^{2}m_{\rm
p}}{X_{\rm H}}\\\ $ (9)

where $C_{\rm grav}$ encodes the specific rate of collapse and  calibrated
$C_{\rm grav}\approx 1$ from their simulations. Note that the condition
$\dot{E}_{\rm grav}=-\dot{E}_{\rm cooling}$ and the condition $t_{\rm
cool}\sim t_{\rm ff}$ are identical, modulo order-unity factors in the
different quantities’ definitions (see Fig. 1). So a collapsing clump heated
only by $P\mathrm{d}V$ work will always trace a track in the $n_{\rm H}-T$
plane corresponding to $t_{\rm cool}/t_{\rm ff}\sim 1$.

Equating Eq. 9 to the optically-thin dust cooling rate gives the dust and gas
temperature where compression heating balances dust cooling:

$\begin{split}T_{\rm grav}&=\left(\frac{2C_{\rm grav}}{acA_{\rm
d}\mu}\sqrt{\frac{\uppi Gm_{\rm p}n_{\rm H}}{X_{\rm
H}^{3}}}\right)^{\frac{1}{3+\beta}}\\\ &\approx 10\mathrm{K}\,\left(Z_{\rm
d}\hat{\sigma}_{\rm d}\right)^{-\frac{1}{3+\beta}}\left(\frac{n_{\rm
H}}{10^{10}\rm cm^{-3}}\right)^{\frac{1}{2\beta+6}}.\end{split}$ (10)

Hence, in the absence of heat sources other than compression, optically-thin
cores do not collapse isothermally: they will rise gradually in temperature as
they collapse, which can be seen in the lowest-temperature calculations in
radiation hydrodynamics simulations (Masunaga et al., 1998; Vaytet &
Haugbølle, 2017).

If the absorption of radiation by dust is the main heat source and the core is
optically-thin to both ingoing and outgoing radiation, then the gas and dust
will approach an equilibrium temperature that is independent of the dust
abundance and the normalization of the opacity, because emission and
absorption are both $\propto A_{\rm d}$. This temperature,

$T_{\mathrm{abs}}=\left(\frac{T_{\rm
rad}^{\beta}u_{\mathrm{rad}}}{a}\right)^{\frac{1}{4+\beta}}\approx
6\mathrm{K}\left(\frac{T_{\rm rad}}{20\,\rm
K}\right)^{\frac{1}{3}}\left(\frac{u_{\rm rad}}{1\,\rm
eV\,\mathrm{cm}^{-3}}\right)^{\frac{1}{6}},$ (11)

depends only on $\beta$, the local density of radiation $u_{\rm rad}$, and
spectral energy distribution represented by the effective radiation
temperature $T_{\rm rad}$. This will typically set the minimum temperature in
well-shielded regions of giant molecular clouds. The radiation may come from
the surrounding galaxy, nearby newly-formed stars, or the cosmic microwave
background. For black-body radiation, $u_{\rm rad}=aT_{\rm rad}^{4}$, so
$T_{\rm abs}=T_{\rm rad}$, and e.g. for the CMB $T_{\rm
CMB}=2.73\mathrm{K}\left(1+z\right)$, which may be important for star
formation at high redshift (see §3.2).

If $T_{\rm rad}>T_{\rm d}$ then the core will become optically-thick to the
background radiation before it becomes optically-thick to its own emission.
This causes a transition to a new regime. Instead of reaching the central
region, incident radiation will be absorbed at larger radii and re-radiated at
longer wavelengths (lower $T_{\rm rad}$), with a limiting optically-thick
value

$T_{\rm bb}=\left(\frac{u_{\rm rad}}{a}\right)^{1/4}\approx
4\,\mathrm{K}\left(\frac{u_{\rm rad}}{1\,\mathrm{eV\,cm^{-3}}}\right)^{1/4}.$
(12)

Strictly speaking this limit corresponds to radiative equilibrium where both
$\tau_{\rm in}>>1$ and $\tau_{\rm out}>>1$, and of course intermediate
temperatures will be reached e.g. if $\tau_{\rm in}>>1$ but $\tau_{\rm
out}<<1$: dust re-radiated from a thin absorption zone at $T_{\rm abs}$ will
reach the center with density equal to $u_{\rm rad}$. But this marginal case
would give a temperature only slightly above Eq. 12 for plausible values of
$T_{\rm rad}$ and $u_{\rm rad}$. The net result is that both the radiation and
gas/dust temperatures at the center will drop from $T_{\rm rad}$ and $T_{\rm
abs}$ respectively, to $\sim T_{\rm bb}=\left(u_{\rm rad}/a\right)^{1/4}$ as
$\tau_{\rm in}$ crosses unity. For the purposes of our present calculation,
this simply means that cases where $T_{\rm rad}>T_{\rm bb}$ initially will map
onto cases where the original incident radiation field was already a black-
body at $T_{\rm bb}$; in effect, the case where $T_{\rm rad}>T_{\rm bb}$ can
effectively be ignored because this condition cannot persist all the way
through the optically-thick and adiabatic transition. An isotropic ambient
radiation field with $T_{\rm rad}<T_{\rm abs}$ (i.e. $u_{\rm rad}>aT_{\rm
rad}^{4}$) is also not thermodynamically possible to create, so although we
have assumed an a priori general radiation field, this does not actually probe
any new parameter space for the opacity limit. The parameter space for the
opacity limit reduces to just $u_{\rm rad}$ and the dust opacity law.

####  2.3.1 General optically-thin behavior

As long as the system is optically-thin to cooling radiation and our other
assumptions hold, $T$ may be approximated by equating compressional plus
radiative heating with radiative cooling, which can only be solved numerically
in general. However we may approximate the overall evolution interpolating
between the temperatures set by the dominant heat source, and accounting for
the case where the clump is optically-thick to the ambient radiation field:

$T_{\rm thin}\approx\begin{cases}\max\left(T_{\rm grav},T_{\rm
abs}\right)&\tau_{\rm in}<1\\\ \max\left(T_{\rm grav},T_{\rm
bb}\right)&\tau_{\rm in}\geq 1\end{cases}$ (13)

This behaviour is seen in the simulation starting at $5\rm K$ in  plotted in
Fig. 1: it starts out in the radiation-domianted regime with $T\sim T_{\rm
bb}\sim const.$, but eventually the temperature starts to rise as we
transition to compression-heated regime $T\sim T_{\rm grav}\propto
n^{\frac{1}{2\beta+6}}$. Note that in any of the above regimes for optically-
thin cooling, the Jeans mass will continue to drop for any plausible value of
$\beta$ – the adiabatic transition cannot occur until $\tau_{\rm out}\gtrsim
1$.

###  2.4 Optically-thick cooling regime

At a given temperature $T$, the clump will be optically-thick to its own
emission if $\tau_{\rm out}>1$, at a density threshold of

$\begin{split}n_{\rm thick}&=\frac{G\mu m_{\rm p}^{2}T^{-2\beta-1}}{\uppi
A_{\rm d}f_{\rm J}^{2}k_{\rm b}X_{\rm H}}\\\ &\approx 1.5\times
10^{11}\,\mathrm{cm}^{-3}\,\left(Z_{\rm d}\hat{\sigma}_{\rm
d}\right)^{-2}\left(\frac{T}{10\,\rm K}\right)^{-2\beta-1},\end{split}$ (14)

which we plot as a dashed line in the $n_{\rm H}-T$ plane in Fig. 1. Above
$n_{\rm thick}$, the maximum total cooling rate per H due to optically-thick
transport of dust-emitted photons is approximated asymptotically by the
radiative diffusion rate (Masunaga & Inutsuka, 1999; Rafikov, 2007):

$\dot{E}_{\rm thick}\approx-\frac{ac\sigma_{\rm
d}T^{4}}{\tau^{2}}=\frac{\dot{E}_{\rm thin}}{\tau^{2}}.$ (15)

When optically-thick, the clump develops a photosphere in radiative
equilibrium with the environment at $T_{\rm bb}$; if radiation were the only
heat source then collapse could continue in radiative equilibrium
indefinitely. This is evident in the simulation with $T_{\rm abs}=30\,\rm K$
shown in Fig. 1: the transition across $\tau=1$ has no significant effect on
the thermal behaviour. Hence the condition $\tau\gtrsim 1$ is not sufficient
to commence adiabatic evolution.

###  2.5 Adiabatic transition and minimum Jeans mass

The adiabatic transition occurs when $\dot{E}_{\rm grav}=-\dot{E}_{\rm
thick}$, or equivalently when the optically-thick cooling time is longer than
the freefall time (Masunaga et al., 1998). Past this point, heat from
gravitational collapse cannot escape as quickly as it is generated. This
criterion is both necessary and sufficient to cause adiabatic evolution. The
demarcation in the $n_{\rm H}-T$ plane is at a density of:

$\begin{split}n_{\rm H}&\left(\dot{E}_{\rm thick}=-\dot{E}_{\rm
grav}\right)\\\ &=5\times 10^{10}\,\mathrm{cm}^{-3}\left(\hat{\sigma}_{\rm
d}Z_{\rm
d}\right)^{-2/3}\left(\frac{T}{10\rm\,K}\right)^{\frac{2\beta-4}{3}},\end{split}$
(16)

which looks like a near-vertical line in Fig. 1 for our fiducial $\beta=1.5$
model, and is exactly vertical for $\beta=2$. Considering heating by
compression and radiation respectively, we consider 2 distinct cases for the
adiabatic transition. The true adiabatic transition density will be roughly
lesser of the two, while the minimum Jeans mass will be the larger.

####  2.5.1 Case 1: Compression heating

If compression is the sole heat source, there is only one characteristic
density $n_{\rm ad}$ at which $\tau\sim 1$, $t_{\rm cool}\sim t_{\rm ff}$, and
$\dot{E}_{\rm grav}\sim-\dot{E}_{\rm thick}\sim-\dot{E}_{\rm thin}$
simultaneously, at the intersection of the lines in Fig. 1. This density is:

$\begin{split}n_{\rm ad,grav}&=\left[\left(\frac{ac}{2C_{\rm
grav}}\right)^{4\beta+2}\left(\frac{\mu}{k_{\rm
B}}\right)^{6\beta+8}\times\right.\\\ &\left.\frac{G^{5}m_{\rm
p}^{2\beta+11}X_{\rm H}^{4\beta-3}}{f_{\rm J}^{4\beta+12}A_{\rm
d}^{10}}\right]^{\frac{1}{4\beta+7}}\\\ &\approx 5.7\times 10^{10}\,\rm
cm^{-3}\left(\hat{\sigma}_{\rm d}Z_{\rm d}\right)^{-\frac{10}{4\beta+7}}.\\\
\end{split}$ (17)

Extrapolating Eq. 10 to $n_{\rm ad}$, the temperature at the adiabatic
transition is

$\begin{split}T_{\rm ad,grav}&=\left(\frac{4C_{\rm grav}^{2}G^{2}k_{\rm
B}m_{\rm p}^{3}}{a^{2}c^{2}A_{\rm d}^{4}f_{\rm J}^{2}X_{\rm
H}^{4}\mu}\right)^{\frac{1}{4\beta+7}}\\\ &\approx
13\,\mathrm{K}\,\left(\hat{\sigma}_{\rm d}Z_{\rm
d}\right)^{-\frac{4}{4\beta+7}},\end{split}$ (18)

and the minimum Jeans mass is Eq. 1 evaluated at $n_{\rm ad,grav}$ and $T_{\rm
ad,grav}$:

$\begin{split}M_{\rm J,grav}^{\rm min}&=\left[\left(\frac{k_{\rm
B}}{\mu}\right)^{9\beta+16}\left(\frac{C_{\rm
grav}}{ac}\right)^{2\beta+4}\times\right.\\\
&\left.\frac{\uppi^{12\beta+21}f_{\rm
J}^{2\beta+3}}{2^{2\beta+3}3^{4\beta+7}G^{6\beta+10}m_{\rm
p}^{9\beta+15}A_{\rm d}X_{\rm H}}\right]^{\frac{1}{4\beta+7}}\\\ &\approx
2.3\times 10^{-3}M_{\odot}\left(\hat{\sigma}_{\rm d}Z_{\rm
d}\right)^{-\frac{1}{4\beta+7}}.\end{split}$ (19)

These expressions correspond to Eqs. 9-11 in Low & Lynden-Bell (1976), and Eq.
19 agrees well with the minimum Jeans mass seen in the lowest-temperature
simulation plotted in Fig. 1. Compression heating is always present, so this
$M_{\rm J}^{\rm min}$ is effectively a lower bound on the dust opacity-limited
minimum Jeans mass, and is indeed remarkably insensitive to dust properties:
$\propto\left(\hat{\sigma}_{\rm d}Z_{\rm
d}\right)^{-\frac{1}{11}}-\left(\hat{\sigma}_{\rm d}Z_{\rm
d}\right)^{-\frac{1}{15}}$ for $\beta\sim 1-2$. However the temperature and
density of the adiabatic transition point may both be significantly higher if
dust opacity is low.

####  2.5.2 Case 2: Radiation heating

If radiative equilibrium is setting the temperature when $\dot{E}_{\rm
grav}=-\dot{E}_{\rm thick}$ is crossed, then the temperature at the adiabatic
transition will be $T_{\rm ad}\sim T_{\rm bb}=\left(u_{\rm
rad}/a\right)^{1/4}$. The density threshold can be obtained by inserting Eq.
11 into Eq. 16:

$\begin{split}n_{\rm ad,rad}&=\left(\frac{a^{2}c^{2}Gm_{\rm p}^{3}T_{\rm
bb}^{4-2\beta}\mu^{4}X_{\rm H}}{4\uppi^{3}A_{\rm d}^{2}C_{\rm grav}^{2}f_{\rm
J}^{4}k_{\rm B}^{4}}\right)^{1/3}\\\ &\approx 5.2\times
10^{10}\,\mathrm{cm}^{-3}\left(\frac{T_{\rm
ad}}{10\,\mathrm{K}}\right)^{\frac{4-2\beta}{3}}\left(\hat{\sigma}_{\rm
d}Z_{\rm d}\right)^{-2/3},\\\ \end{split}$ (20)

and evaluating Eq. 1 at $n_{\rm ad}$ and $T_{\rm ad}$ gives:

$\begin{split}M_{\rm J,rad}^{\rm min}&=\left(\frac{\uppi^{18}A_{\rm
d}^{2}C_{\rm grav}^{2}f_{\rm J}^{4}k_{\rm B}^{13}T_{\rm bb}^{2\beta+5}X_{\rm
H}^{2}}{11664a^{2}c^{2}G^{10}m_{\rm p}^{15}\mu^{13}}\right)^{1/6}\\\ &\approx
1.7\times 10^{-3}M_{\odot}\left(\frac{T_{\rm bb}}{10\rm
K}\right)^{\frac{2\beta+5}{6}}\left(\hat{\sigma}_{\rm d}Z_{\rm
d}\right)^{1/3},\end{split}$ (21)

which corresponds to Eq. 26 in . For the case $T_{\rm bb}=30\,\mathrm{K}$
simulated in , this expression predicts $M_{\rm J}^{\rm min}=7.4\times
10^{-3}M_{\odot}$, again in good agreement with the data in Fig. 1. If the
radiation field is the cosmic microwave background with $T_{\rm
bb}=2.73\left(1+z\right)\,\rm K$, then $M_{\rm J}^{\rm min}=3\times
10^{-4}M_{\odot}\left(1+z\right)^{\frac{2\beta+5}{6}}\left(\hat{\sigma}_{\rm
d}Z_{\rm d}\right)^{1/3}$, which would exceed the compression-dominated value
above $z\sim 3$, assuming $\hat{\sigma}_{\rm d}Z_{\rm d}\sim 1$.

The relevant regime is that which results in the higher temperature, so the
overall minimum Jeans mass will be roughly the larger of Eq. 19 and Eq. 21:

$M_{\rm J}^{\rm min}\approx\max\left(M_{\rm J,rad}^{\rm min},M_{\rm
J,rad}^{\rm min}\right),$ (22)

which we plot in Fig. 2.

###  2.6 Hydrostatic core mass

Another important mass scale that reached the first hydrostatic core before
$\rm H_{2}$ dissociation causes the second collapse at $T_{\rm diss}\sim
1500\rm K$ (Larson, 1969), the “Larson mass” $M_{\rm L}$. Hennebelle et al.
(2019) relate $M_{\rm J}^{\rm min}$ to $M_{\rm L}$ by assuming the entropy of
the core is determined at the opacity limit by $n_{\rm ad}$ and $T_{\rm ad}$,
and the subsequent evolution is isentropic with adiabatic index
$\gamma_{1}=5/3$ below $T_{\rm ex}=150\rm K$, and $\gamma_{2}=7/5$ above it.
The adiabatic scaling of the Jeans mass gives:

$M_{\rm L}=k_{\rm L}M_{\rm J}^{\rm min}\left(\frac{T_{\rm ex}}{T_{\rm
ad}}\right)^{\frac{3\gamma_{1}-4}{2\gamma_{1}-2}}\left(\frac{T_{\rm
diss}}{T_{\rm ex}}\right)^{\frac{3\gamma_{2}-4}{2\gamma_{2}-2}}$ (23)

where $k_{\rm L}$ is a dimensionless factor to account for the new hydrostatic
structure after the adiabatic transition, and $k_{\rm L}\approx 1.3$ recovers
the $\sim 0.03M_{\odot}$ hydrostatic core mass measured in numerical
simulations in the compression-heated limit (Vaytet & Haugbølle, 2017). In the
compression regime we get

$M_{\rm L}=0.03\,M_{\odot}\left(\hat{\sigma}_{\rm d}Z_{\rm
d}\right)^{\frac{2}{4\beta+7}},$ (24)

which is an increasing, but still fairly-shallow function of dust opacity. In
the radiative limit we have

$M_{\rm L}=0.03M_{\odot}\,\left(\frac{T_{\rm
bb}}{10\,\mathrm{K}}\right)^{\frac{4\beta+1}{12}}\left(\hat{\sigma}_{\rm
d}Z_{\rm d}\right)^{\frac{1}{3}},$ (25)

which scales in the same way with dust opacity as the corresponding Jeans
mass, but is somewhat less sensitive to the ambient radiation field.

##  3 Discussion

###  3.1 Comparison with  and

In the previous section we derived approximate expressions for the adiabatic
transition density and minimum Jeans mass under similar assumptions to , but
with one key difference that strongly affects how our predictions generalize
across parameter space. When deriving the minimum Jeans mass,  assumed that
the condition $\dot{E}_{\rm grav}\sim-\dot{E}_{\rm thin}$ was sufficient for
the Jeans mass to reach a minimum (their “Case 1”). This is not a realistic
assumption. From Eqs. 10 and 11 for the optically-thin temperature it is clear
that $T$ cannot scale as steeply as $\propto\rho^{1/3}$ for any plausible dust
model $\beta\sim 1-2$. Hence the Jeans mass $M_{\rm J}\propto
T^{3/2}\rho^{-1/2}$ will always decrease while the clump is optically-thin.
The gas must become optically-thick to undergo an adiabatic transition or
reach a minimum Jeans mass. Although  used the  simulations to guide their
model, crucially there was no example of a simulation that reached its minimum
Jeans mass when $\tau<1$, and we are not aware of such a numerical calculation
in the literature. We do concur that the condition $\tau\gtrsim 1$ is not
sufficient for the adiabatic transition, as was assumed by , but it is a
necessary condition.

The correct choice is to ignore Case 1 and consider only the criterion
$\dot{E}_{\rm grav}\approx-\dot{E}_{\rm thick}$ (their “Case 3”), which is
necessary and sufficient for the adiabatic transition. This eliminates the
predicted regime in which $M_{\rm J}^{\rm min}\propto\left(Z_{\rm
d}\hat{\sigma}_{\rm d}\right)^{-1}$ that would imply strong variations in the
opacity limit with dust abundance or properties. A Jeans-scale core with
$M_{\rm J}^{\rm min}=3.7M_{\odot}$, $Z_{\rm d}\hat{\sigma}_{\rm d}=10^{-2}$,
and $T=10\rm\,K$ would have an optical depth $\tau\sim 10^{-6}-10^{-5}$, and
thus would be very well described by the optically-thin behavior in §2.3, and
would continue to decrease in Jeans mass. Two regimes remain: the compression-
dominated regime (§2.5.1) where $M_{\rm J}$ is nearly independent of dust
properties (Eq. 19), and the radiation-dominated regime where $M_{\rm
J}\propto\left(Z_{\rm d}\hat{\sigma}_{\rm d}\right)^{1/3}$ (Eq. 21), a much
milder sensitivity to dust properties.

###  3.2 Implications for the IMF

Figure 2: Variation of the minimum Jeans mass $M_{\rm J}^{\rm min}$ as a
function of the Solar-normalized dust opacity parameter $Z_{\rm
d}\hat{\sigma}_{\rm d}$ and the radiative equilibrium temperature $T_{\rm bb}$
(Eq. 12). We assume a $\beta=1.5$ dust opacity model and estimate $M_{\rm
J}^{\rm min}$ as the greater of the compression-dominated limit (Eq. 19) and
the radiation-dominated limit (Eq. 21). The dotted line demarcates the
transition between the two regimes. The location of the Solar neighborhood in
this parameter space is indicated, assuming $u_{\rm rad}\sim
1\,\mathrm{eV}\,\mathrm{cm}^{-3}$. Vertical lines plot the lower bound $T_{\rm
bb}\geq 2.73\left(1+z\right)\mathrm{K}$ imposed by the CMB at various
indicated redshifts. Corresponding values of $\log\Sigma_{\rm SFR}$ in
$M_{\odot}\,\mathrm{yr}^{-1}\mathrm{kpc}^{-2}$ are also indicated, assuming
$u_{\rm rad}\propto\Sigma_{\rm SFR}$ and scaling from the Solar neighborhood
value.

The opacity-limited minimum Jeans mass has long been thought to be responsible
for some feature of the IMF. Most naturally, it could shape the low-mass
cutoff and the abundance of brown dwarfs. Bate (2005) performed controlled
numerical experiments in which their varied it explicitly with a barotropic
equation of state, and found that it determined the low-mass IMF cutoff, while
the mean cloud Jeans mass determined the mean stellar mass (e.g. Bate &
Bonnell, 2005). More recently, Lee & Hennebelle (2018a, b, 2019); Hennebelle
et al. (2019); Colman & Teyssier (2020) used simulations and analytic theory
to explore how $M_{\rm J}^{\rm min}$ might set the IMF peak, on the order of
$0.1M_{\odot}$. The jump from Jupiter-like to M dwarf masses has two steps.
First, from the opacity limit to the hydrostatic core $\sim 10$ times more
massive (§2.6). Then, there is a region around this core in which
fragmentation is suppressed due to the core’s tidal field, so the gas
reservoir for a typical star is several times larger than the hydrostatic core
mass.

In ambient radiation fields similar to the Solar neighborhood ($\sim
1\rm\,eV\,cm^{-3}$), $M_{\rm J}^{\rm min}$ is nearly constant, scaling
$\propto A_{\rm d}^{-1/13}-A_{\rm d}^{-1/15}$, leading to factor of $\sim 2$
variations over the entire parameter space where a dust-coupled phase exists.
So any IMF feature set by $M_{\rm J}^{\rm min}$ should not vary strongly in
galactic regions with moderate radiation fields. On the other hand, where
radiation fields are strong, such as in galactic centers, starburst galaxies,
or in the high-redshift Universe where the CMB is $\left(1+z\right)^{4}$ times
more intense, Eq. 21 predicts non-negligible variations in $M_{\rm J}^{\rm
min}$, and by extension the IMF, as a function of the ambient radiation field
and metallicity. In Figure 1 we map out the parameter space for variations in
$M_{\rm J}^{\rm min}$ as a function of $Z_{\rm d}\hat{\sigma}_{\rm d}$ and
$u_{\rm rad}$/$T_{\rm bb}$. The most severely-affected regions would be those
with $Z_{\rm d}\hat{\sigma}_{\rm d}\gtrsim 1$ and $T_{\rm abs}\gtrsim
10\rm\,K$, corresponding to regions of galaxies of Milky Way mass or greater,
with ambient radiation fields $\gtrsim 100\times$ the Solar neighborhood. Such
regions would generally have a significantly less bottom-heavy IMF, or a
greater IMF peak mass.

Bate (2023) presented a suite of star formation simulations with a $z\sim 5$
CMB covering a range of metallicities from $0.01-1$, and found their results
only deviated significantly from their $z\sim 0$ results when $Z_{\rm
d}\hat{\sigma}_{\rm d}\approx 1$, which produced an IMF scaled by a factor of
$\sim 3$ in mass. Their dust model (detailed in Bate 2014) is approximated by
$\beta\sim 2.2$, and for these parameters our model predicts $M_{\rm J}^{\rm
min}=5\times 10^{-3}$. Meanwhile the lower-metallicity $z=5$ models and Bate
(2014) $z=0$ models would all still be in the compression-dominated regime
with $M_{\rm J}^{\rm min}\sim 2\times 10^{-3}M_{\odot}$. Therefore the results
of these simulations are fairly well-fit by model in which the IMF is
determined by the opacity limit.

Once star formation has begun in a cluster, the ambient radiation field may be
of secondary importance to that generated by the star cluster itself,
particularly once the spatial correlation of young stars is taken into account
(Lee & Hopkins, 2020). For example, in the Solar neighborhood GMC simulated in
Grudić et al. (2022) we found that the radiation field seen by a median H atom
in the cloud was elevated by a factor of 100 compared to the ISRF once massive
stars started appearing - and the enhancement in dense regions forming most of
the stars would be significantly greater. In this scenario too the opacity
limit would move to higher masses and the formation of very low-mass stars
would be suppressed, although it would have been only marginally-resolved at
the $10^{-3}M_{\odot}$ mass resolution of our simulation. Note that this
scenario, in which protostellar radiation modifies the opacity limit, is
distinct from the local self-regulation of the Jeans mass around individual
stars due to protostellar radiation, as explored in other works (Bate, 2009;
Krumholz, 2011).

####  3.2.1 Importance of other physics for the IMF

The IMF is probably not determined by the opacity limit alone. This was noted
by Bate & Bonnell (2005), who also found sensitivity to the cloud bulk
properties in hydrodynamic simulations with a barotropic equation of state. It
is also evident from the results of the isothermal MHD star formation
experiments of Haugbølle et al. (2018) and Guszejnov et al. (2020), neither of
which included dust opacity physics in any form. Both works demonstrated a
turnover feature in the mass spectrum from the Salpeter-like
$\mathrm{d}N/\mathrm{d}M\sim M^{-2}$ to the flat $\mathrm{d}N/\mathrm{d}M\sim
M^{-1}$, argued by Haugbølle et al. (2018) to be related to the turbulence-
compressed Jeans mass. The GMC simulations with protostellar jet feedback by
Guszejnov et al. (2021) also included an isothermal run that produced a very
similar IMF to the runs with realistic cooling physics (except at small
masses), indicating that those IMF results were not being determined by the
opacity limit. These results suggest that the thermal structure of GMCs at
densities much lower than the opacity limit can also influence the IMF. This
is qualitatively consistent with the predictions of turbulent fragmentation
theories (Padoan & Nordlund, 2002; Hennebelle & Chabrier, 2008; Hopkins,
2013).

Another important regulator of fragmentation is accretion-powered protostellar
radiation mediated by dust heating, which inevitably dominates the thermal
structure of protostellar disks, suppressing disk fragmentation and brown
dwarf formation (Bate, 2009; Offner et al., 2009). The fundamental physics of
protostellar evolution and radiation have even been proposed to determine the
characteristic mass of the IMF more generally, if the thermal structure of
individual collapsing fragments can be regulated by radiation from the central
star (Krumholz, 2011; Krumholz et al., 2016; Guszejnov et al., 2016; Federrath
et al., 2017). However, in this picture one finds that the self-regulated
fragment mass decreases with increasing ambient ram pressure and/or surface
density (Krumholz, 2011; Tanvir et al., 2022). The detailed numerical
prediction in Fig. 5 of Krumholz (2011) is well-approximated by

$M_{\rm K}\approx 0.03M_{\odot}\left(\frac{P}{10^{10}\rm
cm^{-3}}\right)^{-0.28}\left(\hat{\sigma}_{\rm d}Z_{\rm d}\right)^{-0.22}.$
(26)

This might explain why Hennebelle et al. (2020) did not find that their IMF
peak prediction was strongly affected by the addition of protostellar
radiation: they simulated very dense star-forming clumps with pressures
upwards of $G\Sigma^{2}\sim 10^{8}-10^{10}\,k_{\rm B}\,\rm cm^{-3}\,\rm K$,
where the Krumholz mass is predicted to be at most comparable to $M_{\rm L}$,
making the effects difficult to distinguish.

The natural synthesis of these findings is that initial and boundary
conditions should be taken into account when asking which processes determine
the IMF. There are at least three distinct, physically-plausible fragment mass
scales that could produce corresponding IMF features: one related to the
opacity limit, another related to the physics of turbulent fragmentation, and
another to protostellar heating.

##  4 Conclusion

In this work we have revisited the physics of dust-cooled protostellar
collapse, under a set of analytic asymptotic approximations guided and
validated by numerical simulations. We have adopted a physically-consistent
model in which all relevant quantities are expressed in terms of dust
properties and ambient conditions, and adopted calibrations in quantitative
agreement with the results of radiation hydrodynamics simulations (e.g. Fig.
1). Crucially, we have avoided fallacious assumptions about the conditions for
the adiabatic transition (c.f. §3.1). Lastly, we have explored some
implications of the environmental variation of the minimum Jeans mass in the
context of the results of recent star formation simulations. Nothing we have
done here is fundamentally new, but we have aimed to put a fresh coat of paint
on these important concepts, and to sort out some confusing or conflicting
results.

Our key result is that the opacity-limited minimum Jeans mass is roughly the
larger of a compression-heated limit, Eq. 19 (Eq. 11 in Low & Lynden-Bell
1976), and a radiation-heated limit, Eq. 21 (generalizing Eq. 26 in Masunaga &
Inutsuka 1999). The implied variations with metallicity and the radiative
equilibrium temperature (Fig. 2) may warrant consideration in modern
discussions of the role of the opacity limit or the Larson mass in setting the
IMF’s lower cutoff or peak, particularly in extreme or high-redshift
conditions.

## Acknowledgements

We are grateful for helpful feedback from Stella Offner, Matthew Bate and
Patrick Hennebelle. Support for MYG was provided by NASA through the NASA
Hubble Fellowship grant #HST-HF2-51479 awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555. This work used
computational resources provided by Frontera allocation AST21002. This
research is part of the Frontera computing project at the Texas Advanced
Computing Center. Frontera is made possible by National Science Foundation
award OAC-1818253.

## References

  * Bate (2005) Bate M. R., 2005, MNRAS, 363, 363
  * Bate (2009) Bate M. R., 2009, MNRAS, 392, 1363
  * Bate (2014) Bate M. R., 2014, MNRAS, 442, 285
  * Bate (2023) Bate M. R., 2023, MNRAS, 519, 688
  * Bate & Bonnell (2005) Bate M. R., Bonnell I. A., 2005, MNRAS, 356, 1201
  * Bhandare et al. (2018) Bhandare A., Kuiper R., Henning T., Fendt C., Marleau G.-D., Kölligan A., 2018, A&A, 618, A95
  * Boyd & Whitworth (2005) Boyd D. F. A., Whitworth A. P., 2005, A&A, 430, 1059
  * Colman & Teyssier (2020) Colman T., Teyssier R., 2020, MNRAS, 492, 4727
  * Federrath et al. (2017) Federrath C., Krumholz M., Hopkins P. F., 2017, in Journal of Physics Conference Series. p. 012007, doi:10.1088/1742-6596/837/1/012007
  * Grudić et al. (2022) Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, MNRAS, 512, 216
  * Guszejnov et al. (2016) Guszejnov D., Krumholz M. R., Hopkins P. F., 2016, MNRAS, 458, 673
  * Guszejnov et al. (2018) Guszejnov D., Hopkins P. F., Grudić M. Y., Krumholz M. R., Federrath C., 2018, MNRAS, 480, 182
  * Guszejnov et al. (2020) Guszejnov D., Grudić M. Y., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2020, MNRAS, 496, 5072
  * Guszejnov et al. (2021) Guszejnov D., Grudić M. Y., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 502, 3646
  * Harada et al. (2017) Harada N., Hasegawa Y., Aikawa Y., Hirashita H., Liu H. B., Hirano N., 2017, ApJ, 837, 78
  * Haugbølle et al. (2018) Haugbølle T., Padoan P., Nordlund Å., 2018, ApJ, 854, 35
  * Hennebelle & Chabrier (2008) Hennebelle P., Chabrier G., 2008, ApJ, 684, 395
  * Hennebelle et al. (2019) Hennebelle P., Lee Y.-N., Chabrier G., 2019, ApJ, 883, 140
  * Hennebelle et al. (2020) Hennebelle P., Commerçon B., Lee Y.-N., Chabrier G., 2020, ApJ, 904, 194
  * Hopkins (2013) Hopkins P. F., 2013, MNRAS, 430, 1653
  * Hoyle (1953) Hoyle F., 1953, ApJ, 118, 513
  * Klessen & Glover (2023) Klessen R. S., Glover S. C. O., 2023, arXiv e-prints, p. arXiv:2303.12500
  * Krumholz (2011) Krumholz M. R., 2011, ApJ, 743, 110
  * Krumholz (2014) Krumholz M. R., 2014, Phys. Rep., 539, 49
  * Krumholz et al. (2016) Krumholz M. R., Myers A. T., Klein R. I., McKee C. F., 2016, MNRAS, 460, 3272
  * Larson (1969) Larson R. B., 1969, MNRAS, 145, 271
  * Lebreuilly et al. (2023) Lebreuilly U., Vallucci-Goy V., Guillet V., Lombart M., Marchand P., 2023, MNRAS, 518, 3326
  * Lee & Hennebelle (2018a) Lee Y.-N., Hennebelle P., 2018a, A&A, 611, A88
  * Lee & Hennebelle (2018b) Lee Y.-N., Hennebelle P., 2018b, A&A, 611, A89
  * Lee & Hennebelle (2019) Lee Y.-N., Hennebelle P., 2019, A&A, 622, A125
  * Lee & Hopkins (2020) Lee E. J., Hopkins P. F., 2020, MNRAS, 495, L86
  * Low & Lynden-Bell (1976) Low C., Lynden-Bell D., 1976, MNRAS, 176, 367
  * Martel et al. (2006) Martel H., Evans II N. J., Shapiro P. R., 2006, ApJS, 163, 122
  * Masunaga & Inutsuka (1999) Masunaga H., Inutsuka S.-i., 1999, ApJ, 510, 822
  * Masunaga et al. (1998) Masunaga H., Miyama S. M., Inutsuka S.-i., 1998, ApJ, 495, 346
  * McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
  * Offner et al. (2009) Offner S. S. R., Klein R. I., McKee C. F., Krumholz M. R., 2009, ApJ, 703, 131
  * Padoan & Nordlund (2002) Padoan P., Nordlund Å., 2002, ApJ, 576, 870
  * Penston (1969) Penston M. V., 1969, MNRAS, 144, 425
  * Pollack et al. (1994) Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., Roush T., Fong W., 1994, ApJ, 421, 615
  * Rafikov (2007) Rafikov R. R., 2007, ApJ, 662, 642
  * Rees (1976) Rees M. J., 1976, MNRAS, 176, 483
  * Semenov et al. (2003) Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611
  * Sharda & Krumholz (2022) Sharda P., Krumholz M. R., 2022, MNRAS, 509, 1959
  * Silk (1977) Silk J., 1977, ApJ, 214, 152
  * Tanvir et al. (2022) Tanvir T. S., Krumholz M. R., Federrath C., 2022, MNRAS, 516, 5712
  * Tu et al. (2022) Tu Y., Li Z.-Y., Lam K. H., 2022, MNRAS, 515, 4780
  * Vaytet & Haugbølle (2017) Vaytet N., Haugbølle T., 2017, A&A, 598, A116
  * Weidenschilling & Ruzmaikina (1994) Weidenschilling S. J., Ruzmaikina T. V., 1994, ApJ, 430, 713
  * Weingartner & Draine (2001) Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296
  * Whitworth & Stamatellos (2006) Whitworth A. P., Stamatellos D., 2006, A&A, 458, 817