# Unsupervised Domain Adaptation via Calibrating Uncertainties

Ligong Han  Department of Computer Science, Rutgers University  Yang Zou
Department of Electrical and Computer Engineering, Carnegie Mellon University
Ruijiang Gao  McCombs School of Business, The University of Texas at Austin
Lezi Wang  Department of Computer Science, Rutgers University  Dimitris
Metaxas  Department of Computer Science, Rutgers University

###### Abstract

Unsupervised domain adaptation (UDA) aims at inferring class labels for
unlabeled target domain given a related labeled source dataset. Intuitively, a
model trained on source domain normally produces higher uncertainties for
unseen data. In this work, we build on this assumption and propose to adapt
from source to target domain via calibrating their predictive uncertainties.
The uncertainty is quantified as the Rényi entropy, from which we propose a
general Rényi entropy regularization (RER) framework. We further employ
variational Bayes learning for reliable uncertainty estimation. In addition,
calibrating the sample variance of network parameters serves as a plug-in
regularizer for training111This is an extended version of a CVPRW’19 workshop
paper with the same title. In the current version the gradient variance
regularizer is discussed as a new quantification of uncertainty and a way of
regularizing the training dynamics.. We discuss the theoretical properties of
the proposed method and demonstrate its effectiveness on three domain-
adaptation tasks.

##  1 Introduction

The ability to model uncertainty is important in unsupervised domain
adaptation (UDA). For example, self-training [16, 32] often requires the model
to reliably estimate the uncertainty of its prediction on target domain in the
pseudo-label selection phase. However, traditional deep neural networks (DNN)
can easily assign high confidence to a wrong prediction [5, 18], thus are not
able to reliably and quantitatively render the uncertainty given data.

On the one hand, Bayesian neural networks (BNN) [20, 5, 1, 12] tackles this
problem by taking a Bayesian view of the model training. Instead of obtaining
a point estimate of weights, BNN tries to model the distribution over weights.
We leverage BNN as a powerful tool for uncertainty estimation. On the other
hand, one can estimate the empirical uncertainty of the model by the variance
of network parameters, which we call gradient variance regularization (GVR).

Finally, our approach builds on the intuition that a model gives similar
uncertainty estimates on the two domains learns to adapt from source to target
well. Thus, we propose to directly calibrate the estimated uncertainties
between source and target domains during training. This calibration can be
considered in three-folds, from which we listed our contributions as follows:

  * •

We propose a general framework for unsupervised domain adaptation by
calibrating the predictive uncertainty, and discuss its relationship with
entropy regularization [8] and self-training [16].

  * •

We introduce variational Bayes neural networks to provide reliable uncertainty
estimations.

  * •

We propose to calibrate the variance of network parameters as a model-and-
objective-agnostic regularization (GVR) on the optimization dynamics.

##  2 Related Work

Shannon entropy is commonly used to quantify the uncertainty of a given
distribution. Entropy-based UDA has already been proposed in [28]. Unlike
[28], we avoid using adversarial learning which tends to be unstable and hard
to train. Also, entropy regularization is proposed in [8] for semi-supervised
learning and can be directly applied to UDA. However, our framework is more
general since the uncertainty is not necessarily to be the Shannon entropy. In
fact, we formalize the uncertainty as Rényi entropy which is a generalization
of Shannon entropy. Many other methods in UDA can be modeled under this
framework, for example, self-train [16, 32] can be viewed as minimizing the
min-entropy which is a special case of Rényi entropy.

As pointed out by [7], directly optimizing the estimated Shannon entropy given
data requires the classifier to be locally-Lipschitz [19]. Co-DA [14] and
DIRT-T [25] propose to solve this problem by incorporate the locally-Lipschitz
constraint via virtual adversarial training (VAT) [19].

Another complimentary line of research employs self-ensemble and shows
promising results [3]. Indeed, BNN [5] performs Bayesian ensembling by nature.
This is part of the reason why BNN provides a better uncertainty estimation.

In supervised learning, regularization is proposed to avoid overfitting.
Besides weight decay, typical regularization techniques include label
smoothing [6, 26], network output regularization [23], knowledge distillation
[11]. We believe our proposed gradient variance regularizer GVR can also be
used in supervised settings.

##  3 Uncertainty in Deep Neural Networks

Rényi entropy. For a discrete probability distribution
$P=(P_{1},\ldots,P_{K})$, the Rényi entropy [29] of order $\alpha$
($\alpha>0$) is defined as

$\displaystyle\text{H}_{\alpha}(P)=\frac{1}{1-\alpha}\log(\sum_{k}{P_{k}^{\alpha}}).$
(1)

The limiting value of $\text{H}_{\alpha}$ when $\alpha\rightarrow 1$ is the
Shannon entropy, and $\alpha\rightarrow\infty$ corresponds to the min-entropy,
$\text{H}_{\infty}(P)=\min_{k}{-\log(P_{k})}=-\log\max_{k}P_{k}$. A typical
deep neural network for classification usually produces a discrete
distribution over possible classes given the input data. Thus, we quantify the
predictive uncertainty by the Rényi entropy on this probability distribution.

Bayesian neural networks. BNN estimates the posterior over network weights
while optimizing the training objective. Given the dataset
$\mathcal{D}=\\{x^{(i)},y^{(i)}\\}_{i=1}^{N}$, the output of BNN is denoted as
$f(x|w)$ where $x$ is input data and $w$ are the weights. For a classification
task, $f$ is the predicted logits and the resulting probability vector is
given by a softmax function: $P(y|x,w)=\text{softmax}(f(x|w))$. The predictive
distribution over labels given input $x$ is
$P(y|x)=\mathbb{E}_{P(w|\mathcal{D})}{P(y|x,w)}$. We define the uncertainty
evaluated by BNNs as the entropy $\text{H}_{\alpha}(P(y|x))$.

We adopt the method from [12], where aleatoric and epistemic uncertainties are
jointly modeled. In [12], the logits are assumed to be Gaussian and the
reparameterization trick is utilized. The predicted logit is
$\hat{f}(x)=\mu_{\theta}(x)+\sigma_{\theta}(x)\epsilon$ with $\epsilon\sim
N(0,I)$. The final predicted probability vector $P(y|x)$ is approximated by
Monte Carlo sampling (with $M$ samples),

$\displaystyle
P(y|x;\theta)=\frac{1}{M}\sum_{m=1}^{M}{\text{softmax}(\hat{f}^{(m)}(x))}.$
(2)

Variational inference. As estimating the posterior $P(w|\mathcal{D})$ is often
intractable [1, 12], variational inference is commonly adopted, where the
posterior of weights is approximated by $Q_{\theta}(w)\approx
P(w|\mathcal{D})$ with parameter $\theta$. Specifically, in supervised
learning, $Q_{\theta}(w)$ is estimated by maximizing the evidence lower bound
(ELBO) [13, 5]:

$\displaystyle
ELBO=\underbrace{\mathbb{E}_{Q_{\theta}(w)}{\log{P(y|x,w)}}}_{(\text{I})}-\underbrace{\text{D}_{KL}(Q_{\theta}\|P(w))}_{(\text{II})},$
(3)

where $P(w)$ is the prior, and term (I) is the standard cross-entropy loss
evaluated at $x$ with parameter $w$. Gal _et al_. [4, 5] proposes to view
dropout together with weight decay as a Bayesian approximation, where sampling
from $Q_{\theta}$ is equivalent to performing dropout and term (II) becomes
$L_{2}$ regularization (or weight decay) on $\theta$.

Gradient variance. Rather than finding a variational approximation of the
posterior $Q_{\theta}(w)$, one can instead estimate the model-dependent
uncertainty by the sample variance of $\theta$ (or the sample variance of $w$
in the case of non-Bayesian networks). To be precise, sampling mini-batches
$\mathcal{B}_{i}$ from a batch $\mathcal{B}$, one can compute the adapted
parameters by performing one gradient step (at $\theta$):
$\theta^{\prime}_{i}=\theta-\eta\nabla_{\theta}\mathcal{L}_{\mathcal{B}_{i}}(f_{\theta})$,
where $\mathcal{L}$ is the objective and $\eta$ is the inner learning-rate.
Then the variance of $\theta$ can be defined as the trace of the covariance of
vectorized $\theta^{\prime}_{i}$s:

$\displaystyle\text{Var}(\\{\theta^{\prime}_{i}\\})=\text{trace}(\text{Cov}({\\{\text{vec}(\theta^{\prime}_{i})\\}})),$
(4)

where
$\text{Cov}(\\{\text{vec}(\theta^{\prime}_{i})\\})=\eta^{2}\text{Cov}\left(\\{\text{vec}(\nabla_{\theta}\mathcal{L}_{\mathcal{B}_{i}}(f_{\theta}))\\}\right)$
and $\\{\cdot\\}$ denotes a collection or a set. It can be easily seen that
regularizing the variance of parameters is essentially regularizing the
variance of gradients. We will discuss the usage of this gradient variance as
a regularizer as well as its relationship with MAML [2] in the next section.

##  4 Domain Adaptation via Calibrating Uncertainties

Rényi entropy regularization. Denote source and target dataset as
$\mathcal{D}_{S}=\\{x^{(s)},y^{(s)}\\}_{s\in\mathcal{S}}$ and
$\mathcal{D}_{T}=\\{x^{(t)}\\}_{t\in\mathcal{T}}$ respectively, where
$x^{s},x^{t}$ indicate the samples and $y^{s}$ is the label in source domain,
and $\mathcal{D}=\mathcal{D}_{S}\cup\mathcal{D}_{T}$. We propose to calibrate
the predictive uncertainty of target dataset with the source domain
uncertainties. Concretely, we minimize the cross-entropy (CE) loss in the
source domain while constraining the predicted entropy in the target domain:

$\displaystyle\min_{\theta}\mathcal{L}_{CE}$
$\displaystyle=\frac{1}{|\mathcal{S}|}\sum_{s\in\mathcal{S}}{\text{H}_{CE}(y^{(s)},P(y|x^{(s)};\theta))}$
s.t.
$\displaystyle\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}{\text{H}_{\alpha}(P(y|x^{(t)};\theta))}\leq
C,$ (5)

where $\text{H}_{CE}(\cdot,\cdot)$ is the cross-entropy and $C$ indicates the
strength of the applied constraint. In practice, the network is first
pretrained on labeled source dataset using ELBO in Equation 3. Then, unlabeled
target data is introduced in the above Equation 5, and $P(y|x;\theta)$ is
computed from Equation 2. Note that the resulting CE loss is no longer the
term (I) in ELBO, since the expectation is inside logarithm. We simply treat
$Q_{\theta}(w)$ “as is” the true posterior and evaluate CE using
$P(y|x;\theta)$. For a non-Bayesian network, $P(y|x,w)$ is used as a
replacement of $P(y|x;\theta)$.

To solve Equation 5, rewrite it as a Lagrangian with a multiplier $\beta$,

$\displaystyle\mathcal{F}=$
$\displaystyle\frac{1}{|\mathcal{S}|}\sum_{s\in\mathcal{S}}{\text{H}_{CE}(y^{(s)},P(y|x^{(s)};\theta))}+$
$\displaystyle\beta\left(\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}{\text{H}_{\alpha}(P(y|x^{(t)};\theta))}-C\right).$
(6)

Since $\beta,C\geq 0$, an upper bound on $\mathcal{F}$ is obtained,

$\displaystyle\mathcal{F}\leq$
$\displaystyle\frac{1}{|\mathcal{S}|}\sum_{s\in\mathcal{S}}{\text{H}_{CE}(y^{(s)},P(y|x^{(s)};\theta))}+$
$\displaystyle\frac{\beta}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}{\text{H}_{\alpha}(P(y|x^{(t)};\theta))}=\mathcal{L}_{\alpha}.$
(7)

Ideally, Equation 6 can be optimized via dual gradient descent and $\beta$ is
jointly updated along with $\theta$. For simplicity, we follow the work of
[10] and fix $\beta$ as a hyper-parameter in the experiment and minimize the
upper bound $\mathcal{L}_{\alpha}$.

Note that letting $\alpha\rightarrow 1$ in Equation 7 is in fact the (Shannon)
entropy regularization as described in [7, 8], except that here we consider a
variational BNN. As pointed out in [8], directly optimizing Equation 7 can be
difficult and expectation maximization (EM) algorithms are often used.
Proposed in [30, 8], deterministic annealing EM anneals the predicted
probabilities as soft-labels and minimizes the resulting cross-entropy. In an
extreme case, soft-labels become one-hot vectors and the algorithm terms out
to be self-training with pseudo-labels [16]. In our Rényi entropy
regularization framework, self-training is essentially optimizing the min-
entropy ($\alpha\rightarrow\infty$). Then the objective reads

$\displaystyle\mathcal{L}_{\infty}=$
$\displaystyle\frac{1}{|\mathcal{S}|}\sum_{s\in\mathcal{S}}{\text{H}_{CE}(y^{(s)},P(y|x^{(s)};\theta))}+$
$\displaystyle\frac{\beta}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}{\text{H}_{CE}(\hat{y}^{(t)},P(y|x^{(t)};\theta))},$
(8)

with
$\hat{y}^{(t)}=\text{onehot}(\operatorname*{arg\,max}_{k\in\\{1,\ldots,K\\}}{P(y_{k}|x^{(t)};\theta)})$
to be pseudo-labels in target domain. Subscript $k$ denotes the $k$-th element
in a given $K$-dim vector. The relationship between $\mathcal{L}_{1}$ and
$\mathcal{L}_{\infty}$ can be immediately realized by noticing that the
Shannon entropy is an upper bound of the min-entropy:

$\displaystyle\text{H}_{1}(P)=$
$\displaystyle-\sum_{k}{P_{k}\log(P_{k})}\geq-\sum_{k}{P_{k}}\log(\max_{k}{P_{k}})=$
$\displaystyle-\log(\max_{k}{P_{k}})=\text{H}_{\infty}(P)=\text{H}_{CE}(\hat{y},P).$
(9)

We build our method on top of class-balanced self-training (CBST) proposed in
[32] and use it as the backbone of RER. CBST seeks to select most confident
predictions pseudo-labels in a self-paced (“easy-to-hard”) scheme, since
jointly learning the model and optimizing pseudo-labels on all unlabeled data
is naturally difficult. The authors also propose to normalize the class-wise
confidence levels in pseudo-label generation to balance the class
distribution. For a detailed formulation, we suggest readers referring Section
4.1 and 4.2 in [32].

Gradient variance regularization. The entropy regularization or self-training
framework as formulated above implicitly encourages cross-domain feature
alignment. However, pseudo-labels can be quite noisy even if BNN is employed
to estimated their reliability. Trusting all selected pseudo-labels as one-
hot-encoded “ground-truth” is overconfident and self-training with noisy
pseudo-labels can lead to incorrect entropy minimization. Indeed, we observe
that the model can quickly converge to its overconfident predictions.
Therefore, the parameter variance evaluated in target domain using pseudo-
labels via Equation 4 can be even smaller than that of the source domain. To
address this problem, we again propose to regularize the self-training by
maximizing the gradient variance. Algorithm 1 illustrates the regularized
self-training procedure on target domain (the training on source and target
domains are preformed alternately, which is omitted in the algorithm box).
$\eta$ and $\eta^{\prime}$ are the inner- and outer-stepsize, and $\lambda$ is
the hyper-parameter weighting the regularization term.

Algorithm 1 Gradient Variance Regularization

1:while not done do

2: Sample mini-batches $\mathcal{B}_{i}\sim\mathcal{D}_{T}$

3: for all $\mathcal{B}_{i}$ do

4: Evaluate $\nabla_{\theta}\mathcal{L}_{\mathcal{B}_{i}}(f_{\theta})$

5: Compute
$\theta^{\prime}_{i}=\theta-\eta\nabla_{\theta}\mathcal{L}_{\mathcal{B}_{i}}(f_{\theta})$

6: end for

7: Collect $\Theta^{\prime}=\\{\theta^{\prime}_{i}\\}$

8: Compute
$\text{Var}(\Theta^{\prime})=\text{trace}(\text{Cov}({\\{\text{vec}(\theta^{\prime}_{i})\\}}))$

9: Update
$\theta\leftarrow\theta-\eta^{\prime}\nabla_{\theta}\left(\sum_{i}{\mathcal{L}_{\mathcal{B}_{i}}(f_{\theta})}-\lambda\text{Var}(\Theta^{\prime})\right)$

10:end while

Notice that the proposed GVR shares similarities with MAML [2], comparing from
a dynamical systems standpoint and despite that MAML samples mini-batches of
different tasks. Taking a first-order Taylor expansion of the loss function
around $\theta$,

$\displaystyle\mathcal{L}_{\mathcal{B}_{i}}(f_{\theta^{\prime}_{i}})\approx\mathcal{L}_{\mathcal{B}_{i}}(f_{\theta})-\eta\|{\nabla_{\theta}\mathcal{L}_{\mathcal{B}_{i}}(f_{\theta})}\|_{2}^{2},$
(10)

we demonstrate that MAML tries to maximize the sensitivity of the loss
functions with respect to the parameters by maximizing the $L_{2}$ norms of
the gradients. On the contrary, GVR maximizes the variance of gradients, which
intuitively encourages the model to escape from bad local minima.

It is worth mentioning that GVR is not only model-agnostic but also objective-
agnostic. This is useful when the regularizer itself is the objective to be
optimized. Moreover, GVR is complementary to VAT [19] since in VAT the
gradient is computed with respect to input data. We conjecture that the data
gradient somewhat captures the aleatoric (data-dependent) uncertainty, which
we leave for future work.

##  5 Experiments

We first show results on three digit datasets MNIST [15], USPS and SVHN [21],
where we consider MNIST$\rightarrow$USPS and SVHN$\rightarrow$MNIST. Then we
present preliminary results on a challenging benchmark: VisDA17
(classification) [22] which contains 12 classes. We follow the standard
protocol in [22, 27, 24]. Classification accuracies on source and target
domains for base models are reported in Table 1. We use DTN [31] as our base
model for MNIST$\rightarrow$USPS and SVHN$\rightarrow$MNIST. To implement its
Bayesian variant (BDTN), we add another classifier to predict the logarithm of
variance.

Model | Source | Target
---|---|---
DTN | 100.00 | 83.94
BDTN-M1 | 100.00 | 83.78
BDTN-M5 | 100.00 | 86.83
BDTN-M10 | 100.00 | 86.28
BDTN-M20 | 100.00 | 86.78
BDTN-M100 | 100.00 | 87.06

(a) MNIST

Model | Source | Target
---|---|---
DTN | 97.42 | 72.91
BDTN-M1 | 95.91 | 65.51
BDTN-M5 | 99.16 | 71.12
BDTN-M10 | 99.42 | 71.38
BDTN-M20 | 99.50 | 73.64
BDTN-M100 | 99.33 | 74.91

(b) SVHN

Table 1: Training base models on MNIST and SVHN. BDTN is a modified Bayesian
DTN [31], with different $M$ values (as defined in Equation 2). Classification
accuracies in source and target domains are reported.

Domain adaptation results on digit datasets are shown in Table 2. Our proposed
Rényi entropy regularization methods with non-Bayesian and Bayesian base
models are listed as RERs and BRERs respectively. We see self-training with
pseudo-labels ((B)RER-$\infty$) are in general more stable than directly
minimizing the Shannon entropy ((B)RER-1). Also, adding GVR in (B)RER-$\infty$
improves the performance. However, we also observe that GVR is not helpful in
(B)RER-1 settings.

Model | Target Acc (%) | Acc Gain (%)
---|---|---
Source-DTN | 83.94 | -
Source-BDTN | 84.89 | -
RER-1 |  91.57$\pm$0.13 | 7.63
RER-1-GVR |  91.97$\pm$0.26 | 8.03
RER-$\infty$ |  93.57$\pm$0.30 | 9.63
RER-$\infty$-GVR |  93.88$\pm$0.14 | 9.94
BRER-1 |  92.78$\pm$0.42 | 7.89
BRER-1-GVR |  93.07$\pm$0.72 | 8.18
BRER-$\infty$ |  94.42$\pm$0.12 | 9.53
BRER-$\infty$-GVR | 94.53$\pm$0.23 | 9.64

(a) MNIST$\rightarrow$USPS

Model | Target Acc (%) | Acc Gain (%)
---|---|---
Source-DTN | 64.48 | -
Source-BDTN | 70.98 | -
RER-1 |  88.46$\pm$0.90 | 23.98
RER-1-GVR |  85.48$\pm$4.71 | 21.00
RER-$\infty$ |  88.16$\pm$1.19 | 23.68
RER-$\infty$-GVR |  90.31$\pm$2.31 | 25.83
BRER-1 |  92.49$\pm$4.73 | 21.51
BRER-1-GVR |  92.37$\pm$4.76 | 21.39
BRER-$\infty$ |  96.06$\pm$0.68 | 25.08
BRER-$\infty$-GVR | 96.38$\pm$0.05 | 25.40

(b) SVHN$\rightarrow$MNIST

Table 2: Results on MNIST$\rightarrow$USPS and SVHN$\rightarrow$MNIST. RER
uses DTN [31] as the base model. BRER-$\infty$ uses BDTN as the base model and
optimizes $\mathcal{L}_{\infty}$, while BRER-1 optimizes $\mathcal{L}_{1}$.
Results are averaged over 4 runs with different random seeds.

Mean accuracies on VisDA17 dataset are reported in Table 3. Following the
protocol in [32], we train a standard ResNet101 [9] as the base model and add
a second classifier (denoted as BRes101) to predict the logarithm of variance
on logits. Results show that BNN improves upon non-Bayesian baselines by a
large margin. GVR has not been tested on VisDA17 with (B)Res101 since the
memory requirement exceeds our GPU capacities.

Model | Target mean-Acc (%) | Acc Gain (%)
---|---|---
Source-Res101 | 48.02 | -
Source-BRes101 | 46.03 | -
MMD [17] | 61.1 | -
GTA-Res152 [24] | 77.1 | -
RER-$\infty$ |  76.81$\pm$2.73 | 28.79
BRER-$\infty$ | 80.59$\pm$1.39 | 34.56

Table 3: Preliminary results on VisDA17 [22] classification benchmark
(validation set). Results are averaged over 5 runs with different random
seeds.

##  6 Conclusion

In this work, we propose to approach unsupervised domain adaptation via
calibrating the predictive uncertainties between source and target domains.
The uncertainty is quantified under a general Rényi entropy regularization
framework, within which we introduce Bayesian neural networks for accurate and
reliable uncertainty estimations. From a frequentist point of view, we in
addition propose to approximate the model uncertainty via the sample variance
of network parameters (or gradients) during training. Results show that the
uncertainty estimation by Bayesian networks and gradient variances is
effective and leads to stable performance in unsupervised domain adaptation.

## References

  * [1] C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra.  Weight uncertainty in neural networks.  arXiv preprint arXiv:1505.05424, 2015.
  * [2] C. Finn, P. Abbeel, and S. Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 1126–1135. JMLR. org, 2017.
  * [3] G. French, M. Mackiewicz, and M. Fisher.  Self-ensembling for visual domain adaptation.  arXiv preprint arXiv:1706.05208, 2017.
  * [4] Y. Gal and Z. Ghahramani.  Bayesian convolutional neural networks with bernoulli approximate variational inference.  arXiv preprint arXiv:1506.02158, 2015.
  * [5] Y. Gal and Z. Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In international conference on machine learning, pages 1050–1059, 2016.
  * [6] I. Goodfellow, Y. Bengio, and A. Courville.  Deep learning.  MIT press, 2016.
  * [7] Y. Grandvalet and Y. Bengio.  Semi-supervised learning by entropy minimization.  In Advances in neural information processing systems, pages 529–536, 2005.
  * [8] Y. Grandvalet and Y. Bengio.  Entropy regularization.  Semi-supervised learning, pages 151–168, 2006.
  * [9] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [10] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner.  beta-vae: Learning basic visual concepts with a constrained variational framework.  In International Conference on Learning Representations, volume 3, 2017.
  * [11] G. Hinton, O. Vinyals, and J. Dean.  Distilling the knowledge in a neural network.  arXiv preprint arXiv:1503.02531, 2015.
  * [12] A. Kendall and Y. Gal.  What uncertainties do we need in bayesian deep learning for computer vision?  In Advances in neural information processing systems, pages 5574–5584, 2017.
  * [13] D. P. Kingma and M. Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [14] A. Kumar, P. Sattigeri, K. Wadhawan, L. Karlinsky, R. Feris, B. Freeman, and G. Wornell.  Co-regularized alignment for unsupervised domain adaptation.  In Advances in Neural Information Processing Systems, pages 9345–9356, 2018.
  * [15] Y. LeCun, L. Bottou, Y. Bengio, P. Haffner, et al.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [16] D.-H. Lee.  Pseudo-label: The simple and efficient semi-supervised learning method for deep neural networks.  In Workshop on Challenges in Representation Learning, ICML, volume 3, page 2, 2013.
  * [17] M. Long, Y. Cao, J. Wang, and M. I. Jordan.  Learning transferable features with deep adaptation networks.  arXiv preprint arXiv:1502.02791, 2015.
  * [18] C. Louizos and M. 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.
  * [19] T. Miyato, S.-i. Maeda, M. Koyama, K. Nakae, and S. Ishii.  Distributional smoothing with virtual adversarial training.  arXiv preprint arXiv:1507.00677, 2015.
  * [20] R. M. Neal.  Bayesian learning for neural networks, volume 118.  Springer Science & Business Media, 2012.
  * [21] Y. Netzer, T. Wang, A. Coates, A. Bissacco, B. Wu, and A. Y. Ng.  Reading digits in natural images with unsupervised feature learning.  2011\.
  * [22] X. Peng, B. Usman, N. Kaushik, D. Wang, J. Hoffman, and K. Saenko.  Visda: A synthetic-to-real benchmark for visual domain adaptation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops, pages 2021–2026, 2018.
  * [23] G. Pereyra, G. Tucker, J. Chorowski, Ł. Kaiser, and G. Hinton.  Regularizing neural networks by penalizing confident output distributions.  arXiv preprint arXiv:1701.06548, 2017.
  * [24] S. Sankaranarayanan, Y. Balaji, C. D. Castillo, and R. Chellappa.  Generate to adapt: Aligning domains using generative adversarial networks.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8503–8512, 2018.
  * [25] R. Shu, H. H. Bui, H. Narui, and S. Ermon.  A dirt-t approach to unsupervised domain adaptation.  arXiv preprint arXiv:1802.08735, 2018.
  * [26] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna.  Rethinking the inception architecture for computer vision.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2818–2826, 2016.
  * [27] E. Tzeng, J. Hoffman, K. Saenko, and T. Darrell.  Adversarial discriminative domain adaptation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 7167–7176, 2017.
  * [28] T.-H. Vu, H. Jain, M. Bucher, M. Cord, and P. Pérez.  Advent: Adversarial entropy minimization for domain adaptation in semantic segmentation.  arXiv preprint arXiv:1811.12833, 2018.
  * [29] Wikipedia contributors.  Rényi entropy — Wikipedia, the free encyclopedia, 2018.  [Online; accessed 13-May-2019].
  * [30] A. L. Yuille, P. Stolorz, and J. Utans.  Statistical physics, mixtures of distributions, and the em algorithm.  Neural Computation, 6(2):334–340, 1994.
  * [31] X. Zhang, F. X. Yu, S.-F. Chang, and S. Wang.  Deep transfer network: Unsupervised domain adaptation.  arXiv preprint arXiv:1503.00591, 2015.
  * [32] Y. Zou, Z. Yu, B. Vijaya Kumar, and J. Wang.  Unsupervised domain adaptation for semantic segmentation via class-balanced self-training.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 289–305, 2018.


# The New Worlds Simulations: Large-scale Simulations across Three Cosmologies

Katrin Heitmann11affiliationmark:  Thomas Uram22affiliationmark:  Nicholas
Frontiere33affiliationmark:  Salman Habib11affiliationmark:
33affiliationmark:  Adrian Pope33affiliationmark:  Silvio
Rizzi22affiliationmark:  Joe Insley33affiliationmark:  1HEP Division, Argonne
National Laboratory, 9700 S. Cass Ave., Lemont, IL 60439, USA  2Argonne
Leadership Computing Facility, Argonne National Laboratory, Lemont, IL 60439,
USA  3CPS Division, Argonne National Laboratory, 9700 S. Cass Ave., Lemont, IL
60439, USA

###### Abstract

In this paper we describe the set of “New Worlds Simulations”, three very
large cosmology simulations, Qo’noS, Vulcan, and Ferenginar, that were carried
out on the Summit supercomputer with the Hardware/Hybrid Cosmology Code, HACC.
The gravity-only simulations follow the evolution of structure in the Universe
by each employing 12,2883 particles in (3$h^{-1}$Gpc)3 volumes, leading to a
mass resolution of $m_{p}\sim 10^{9}h^{-1}$M⊙. The simulations cover three
different cosmologies, one $\Lambda$CDM model, consistent with measurements
from Planck, one simulation with massive neutrinos, and one simulation with a
varying dark energy equation of state. All simulations have the same phases to
allow a detailed comparison of the results and the investigation of the impact
of different cosmological parameters. We present measurements of some basic
statistics, such as matter power spectra, correlation function, halo mass
function and concentration-mass relation and investigate the differences due
to the varying cosmologies. Given the large volume and high resolution, these
simulations provide excellent bases for creating synthetic skies. A subset of
the data is made publicly available as part of this paper.

###### keywords:

methods: statistical — cosmology: large-scale structure of the universe

Heitmann et al. The New Worlds Simulations

##  1 Introduction

Future cosmological surveys, such as the Legacy Survey of Space and Time
(LSST) to be carried out with the Rubin Observatory (Abell et al., 2009), and
newly developed instruments and missions, such as the Dark Energy
Spectroscopic Instrument (DESI, Aghamousa et al. 2016), Euclid (Laureijs et
al., 2011), and Roman (Doré et al., 2019), aim to shed light on the nature of
the dark Universe and unravel the mysteries of the origin of dark matter and
the accelerated expansion of the Universe. Measurements from these surveys
have reached a level of accuracy that they can also provide constraints on
neutrino masses. To interpret the wealth of data, simulations have become ever
more important. They provide predictions in the nonlinear regime of structure
formation which holds key information about cosmology. They also function as
testbeds for complex analysis and processing tasks as well as for the
exploration of possible new cosmological probes. For many of these tasks,
large simulation volumes ($\sim$ (3$h^{-1}$Gpc)3 or larger) with excellent
mass resolution (particle masses of $m_{p}\sim 10^{9}h^{-1}$M⊙ or better) are
required to capture survey size and depth.

Figure 1: Visualization of the particle density from Rank 0 (out of 24576
ranks) from the Ferenginar simulation at $z=0$.

Large-scale simulations with high mass resolution have become increasingly
important for cosmological surveys to prepare for data arrival. For example,
the LSST Dark Energy Science Collaboration (LSST DESC) (Abate et al., 2012)
has created the second data challenge (DC2) (LSST Dark Energy Science
Collaboration et al., 2021) based on the Outer Rim simulation (Heitmann et
al., 2019a). DC2 covers 400 deg2 of the LSST footprint and five years of
observations. The Outer Rim simulation is based on a best-fit WMAP-7 cosmology
and has been now updated by the Last Journey simulation, featuring a
cosmological model close to the best-fit Planck results (Planck Collaboration,
2020). Another example includes the Euclid Flagship simulation (Potter,
Stadel, & Teyssier, 2017) to model the ongoing Euclid survey. In addition,
major simulation suites have been created to support survey science, including
the Buzzard and MICE simulations for the Dark Energy Survey (DES) (DeRose et
al., 2019; Crocce et al., 2015), and the AbacusSummit suite (Maksimova et al.,
2021) for DESI science.

In this paper, we introduce three new simulations, similar in size to the
Outer Rim, Last Journey and Euclid Flagship simulations. The Qo’noS, Vulcan,
and Ferenginar simulations were carried out with the Hardware/Hybrid
Accelerated Cosmology Code (HACC), described in Habib et al. (2016). Each
simulation evolves 12,2883 particles in a (3$h^{-1}$Gpc)3 volume. The models
feature three different cosmologies, a $\Lambda$CDM model with the current
best-fit Planck values (Qo’noS), the same cosmology with a massive neutrino
component added (Vulcan), and a dynamical dark energy model, with a dark
energy equation of state parameterized by $(w_{0},w_{a})$ (Ferenginar). Figure
1 shows the visualization of the data from one rank from the Ferenginar
simulation at $z=0$. All simulations have the same initial phases, enabling
detailed comparisons of subtle differences that are not washed out due to
uncertainties from cosmic variance. The set of simulations also lends itself
well to a blind analysis challenge to investigate how sensitive cosmological
probes are to small changes in cosmology.

The creation of such large simulations poses a challenge with regard to code
scalability, as both the computational and memory requirements can only be
satisfied by scaling to large supercomputers. We have developed HACC for more
than a decade and have optimized the code to run on the fastest supercomputers
in the world available at any given time. The architectures have changed
considerably during that time span, including hardware-accelerated machines
like Roadrunner111https://en.wikipedia.org/wiki/Roadrunner_(supercomputer),
Titan222https://en.wikipedia.org/wiki/Titan_(supercomputer) and
Summit333https://en.wikipedia.org/wiki/Summit_(supercomputer) as well as many-
core architectures like Mira444https://www.alcf.anl.gov/alcf-resources/mira,
Sequoia555https://en.wikipedia.org/wiki/Sequoia_(supercomputer) and
Theta666https://www.alcf.anl.gov/alcf-resources/theta. HACC has run
successfully on all of these machines and has been recently prepared to run on
the exascale platforms Frontier and Aurora at the Leadership Computing
Facilities at Oak Ridge and Argonne, respectively. For each simulation
described in this paper, we utilized 4096 nodes on Summit, close to the full
machine. A set of first level analysis tasks was carried out on the fly. These
tasks include halo finding with different halo definitions, power spectrum
estimates, and halo core finding to identify substructures. Results from
analyzing these first level data products are presented in this paper and the
differences due to the cosmological parameter choices are highlighted.

The paper is organized as follows. In Section 2 we provide an overview of the
three simulations and the outputs that were stored. Next, in Section 3 we
present measurements for a range of statistics and show the deviations due to
the differences in the cosmological models simulated. We briefly discuss
publicly available data products in Section 4, and conclude in Section 5.

##  2 Simulation Set-up and Data Products

In this section we provide details about the set-up of the simulations and
data products generated from the simulation runs.

###  2.1 Cosmological Models

The three simulations introduced in this paper, Qo’noS, Vulcan, and Ferenginar
each evolve 12,2883 particles in a (3$h^{-1}$Gpc)3 volume. We choose the same
total matter density value, baryon density value, Hubble parameter,
$\sigma_{8}$ and spectral index for all three cosmologies as follows:

$\displaystyle\Omega_{m}$ $\displaystyle=$ $\displaystyle 0.3096,$ (1)
$\displaystyle\omega_{b}$ $\displaystyle=$ $\displaystyle 0.02242,$ (2)
$\displaystyle h$ $\displaystyle=$ $\displaystyle 0.6766,$ (3)
$\displaystyle\sigma_{8}$ $\displaystyle=$ $\displaystyle 0.8102,$ (4)
$\displaystyle n_{s}$ $\displaystyle=$ $\displaystyle 0.9665.$ (5)

In addition, we restrict the simulations to flat universes with $\Omega_{k}=0$
and hence $\Omega_{\Lambda}=1-\Omega_{m}$. These parameters lead to a particle
mass of $m_{p}=1.26\cdot 10^{9}h^{-1}$M⊙, the same for all simulations. The
force resolution employed is 2.4$h^{-1}$kpc. The simulations were started at
$z=200$ using the Zel’dovich approximation (Zel’dovich, 1970) to assign the
initial particle positions and velocities. The input transfer functions were
generated using CAMB (Lewis, Challinor & Lasenby, 2000). Each simulation has
the same initial phases, and were individually run on 4096 nodes on Summit
using forty two cores and all six GPUs per node. This enables detailed
comparisons of the results and investigations of the differences due to
changes of cosmological parameters.

in addition to the five parameters given in Eqs. (1) – (5), for the Qo’noS
simulation we have chosen a constant dark energy of state and massless
neutrinos:

$\displaystyle w_{0}$ $\displaystyle=$ $\displaystyle-1,$ (6) $\displaystyle
w_{a}$ $\displaystyle=$ $\displaystyle 0,$ (7) $\displaystyle\Omega_{\nu}$
$\displaystyle=$ $\displaystyle 0.$ (8)

The Vulcan simulation is similar to Qo’noS with the addition of a small
massive neutrino component, but with $\Omega_{m}$ unchanged, so that
$\Omega_{cdm}=0.25833$ (a slightly lower value than for Qo’noS). We treat
massive neutrinos at the background level instead of simulating them as a
separate particle species. This approach avoids numerical inaccuracy induced
by simulating particle populations with very different masses. A detailed
description of our approach is given in Upadhye et al. Upadhye et al. (2014)
and Heitmann et al. (2016). Its validity with regard to power spectrum
measurements on large to quasi-linear scales is discussed extensively in
Upadhye et al. (2014) and the effect of neutrinos on the mass function is
investigated for two models in Biswas et al. (2019), including one for high-
mass neutrinos. We chose a model model the neutrino mass of three species via
their sum to be $\sum m_{\nu}=0.1$eV. This leads to the following parameter
choices (where we have used the relationship $\Omega_{\nu}=\sum
m_{\nu}/93.1{\rm eV}h^{2}$):

$\displaystyle w_{0}$ $\displaystyle=$ $\displaystyle-1,$ (9) $\displaystyle
w_{a}$ $\displaystyle=$ $\displaystyle 0,$ (10) $\displaystyle\Omega_{\nu}$
$\displaystyle=$ $\displaystyle 0.00234.$ (11)

Finally, the Ferenginar simulation employs the same cosmology as the Qo’noS
simulation but where we replace the cosmological constant with a dynamical
dark energy model, parameterized via

$w(a)=w_{0}+w_{a}(1-a).$ (12)

The three additional parameters to Eqs. (1) – (5) are given by:

$\displaystyle w_{0}$ $\displaystyle=$ $\displaystyle-0.9,$ (13)
$\displaystyle w_{a}$ $\displaystyle=$ $\displaystyle 0.7,$ (14)
$\displaystyle\Omega_{\nu}$ $\displaystyle=$ $\displaystyle 0.$ (15)

To gain an intuitive understanding about the three different models and how we
expect structure-formation results to differ, we briefly discuss the resulting
Hubble expansion rates $H(z)$ for the Qo’noS and Vulcan simulations compared
to the Ferenginar simulation. The Hubble expansion rate for a dark energy
model with an equation of state parameterized via Eq. (12) is given by

$\begin{split}H(z)&=H_{0}\biggl{\\{}\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}(1+z)^{3(1+w_{0}+w_{a})}\\\
&\qquad\quad\times\exp\left\\{-\frac{3w_{a}z}{1+z}\right\\}\biggr{\\}}^{\\!1/2},\end{split}$
(16)

which simplifies for $\Lambda$CDM cosmologies (Qo’noS and Vulcan) to

$H(z)=H_{0}\biggl{\\{}\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}\biggr{\\}}^{\\!1/2}.$
(17)

Figure 2 shows the evolution of the Hubble parameter $H(z)$ as a function of
redshift for the three models considered in this paper. For comparison, we
also show a model with no dark energy, i.e. $\Omega_{m}=1$ and 31 measurement
points based on cosmic chronometers (CCH), provided in Table 1 in Gómez-Valent
& Amendola (2018). We chose this data set due to its large redshift coverage,
out to $z\sim 2$. The measurements illustrate the uncertainties at higher
redshifts that still hold clues about the exact nature of the cause for the
accelerated expansion of the Universe. The evolution of the Hubble parameter
in our dynamical dark energy model is steeper at early times than for a
$\Lambda$CDM model but more comparable at lower redshifts due to choice of the
same value for $H(0)$ for the models considered.

Figure 2: Redshift evolution of the Hubble parameter for the different
cosmologies considered in the paper and a model without dark energy for
comparison (purple: $\Lambda$CDM, green: dynamical dark energy, blue: no dark
energy). Additionally, we show data from cosmic chronometers (CCH), provided
in Gómez-Valent & Amendola (2018), Table 1.

Another measurement that provides insights into the differences of the models
considered is the linear power spectrum. Figure 3 shows the ratio of the
initial power spectra from the three simulations, scaled to $z=0$ using the
growth functions corresponding to each model. Since for all three models we
chose the same value for $\sigma_{8}$, the power spectra cross at the scale
corresponding to the normalization. Both the dynamical dark energy model and
the neutrino model have more power on large scales compared to the
$\Lambda$CDM cosmology. The differences on small scales are minimal. The
figure demonstrates the major challenge that we will face when trying to
disentangle subtle physics effects such as small neutrino masses or deviations
from $w=-1$ from next generation data. Simulations as those presented here
will be important to demonstrate that our cosmological probes are indeed
sensitive enough for this endeavour.

Finally, we provide some visual insights into the different cosmologies in
Fig. 4. We show a zoom-in for the outputs of one of the ranks from the
simulations. The camera position in these ParaView777https://www.paraview.org/
images is the same for all three models and the different redshifts shown. The
center column shows the result for the $\Lambda$CDM cosmology (Qo’noS), on the
left, we show the varying dark energy equation of state simulation
(Ferenginar) and on the right the cosmology that includes neutrinos (Vulcan).
From top to bottom we show three redshifts, $z=2,1$ and $z=0$. At earlier
times (first row), differences between Ferenginar and Qo’noS/Vulcan
(center/right) are more apparent. For example, the large nodal structure in
the upper right corner is tighter in the Ferenginar simulation as is the
structure in the lower half of the figure close to the right edge of the
figure. Moving forward in time (second row at $z=1$ and third row at $z=0$)
the differences are less pronounced. The visualizations also show how the void
structures grow over time and the filaments become more defined as expected.
Throughout Section 3 we will quantitatively show how the three simulations
display different clustering properties early on and how the differences
reduce over time.

Figure 3: Ratios of the power spectra from the initial conditions of the three
simulations, all scaled to $z=0$ via their linear growth functions. We use the
$\Lambda$CDM model as base. For Qo’noS and Ferenginar, the exact, scale-
independent growth function was used and the full matter power spectrum is
shown. For Vulcan, only the $bc$ component of the power spectrum is shown and
the scale dependence induced by the neutrinos was not taken into account when
scaling to $z=0$. Figure 4: Zoom-in visualization at three redshifts $z=2,1,0$
(top to bottom) for the three different simulations. Particles are overlayed
on the density representation. Differences are clearly visible at early times,
the Ferenginar simulation shows tighter and larger structures in the high-
density regions. Over time, the structures appear more and more similar.

###  2.2 Outputs

The data sizes for the simulations are very large. For example, a single time
snapshot, containing positions and velocities, a particle ID and a potential
value for each particle, has a size of close to 70TB (after lossless
compression). This makes it prohibitive to save the full raw particle data
output across all the time steps. To reduce the data size to a manageable
amount, we ran all first level analysis on the fly (halo finding, power
spectrum measurements, downsampling of particle snapshots) and did not save
full raw data snapshots. We have analyzed 101 time snapshots between $z=10$
and $z=0$, evenly spaced in $\log_{10}(a)$ and saved a diverse set of analysis
products at the following output values in $z$:

$\displaystyle z$ $\displaystyle=$
$\displaystyle\left\\{10.04,9.81,9.56,9.36,9.15,8.76,8.57,8.39,8.05,\right.$
(18) $\displaystyle 7.89,7.74,7.45,7.31,7.04,6.91,6.67,6.56,6.34,$
$\displaystyle 6.13,6.03,5.84,5.66,5.48,5.32,5.24,5.09,4.95,$ $\displaystyle
4.74,4.61,4.49,4.37,4.26,4.10,4.00,3.86,3.76,$ $\displaystyle
3.63.3.55,3.43,3.31,3.21,3.10,3.04,2.94,2.85,$ $\displaystyle
2.74,2.65,2.58,2.48,2.41,2.32,2.25,2.17,2.09,$ $\displaystyle
2.02,1.95,1.88,1.80,1.74,1.68,1.61,1.54,1.49,$ $\displaystyle
1.43,1.38,1.32,1.26,1.21,1.15,1.11,1.06,1.01,$ $\displaystyle
0.96,0.91,0.86,0.82,0.78,0.74,0.69,0.66,0.62,$ $\displaystyle
0.58,0.54,0.50,0.47,0.43,0.40,0.36,0.33,0.30,$ $\displaystyle
0.27,0.24,0.21,0.18,0.15,0.13,0.10,0.07,$
$\displaystyle\left.0.05,0.02,0.00\right\\}.$

For these redshifts we save halo catalogs based on the Friends of Friends
(FOF) definition (Davis et al., 1985) using a linking length $b=0.168$. For
FOF halos with more than 500 particles, we measure the mass overdensity in
spheres around the potential minimum center until it falls below 200 times
$\rho_{c}=3H^{2}/8\pi G$, defining the radius of the SO halos. Then we
determine the enclosed overdensity mass M200c.

We also saved all FOF halo particles. We store FOF halos down to 20 particle
halo size, while SO halos have at least 500 particles per halo. In addition to
a range of halo properties, we also saved SO profiles. To enable substructure
tracking, we stored halo core particle data (for more details about the core
concept, see Rangel et al. 2020 and Korytov et al. 2023 for an implementation
example concerning galaxies in clusters). Finally, we saved down-sampled
particle at 1% snapshots for 59 steps, starting at $z=3.04$, following the
list above.

##  3 Results in the Nonlinear Regime

Next, we display measurements comparing the simulations in the nonlinear
regime, showing results for both particle- and halo-based statistics. We
compare some of the predictions with recently released emulators to verify the
accuracy of the simulations.

###  3.1 Power Spectra

We begin with measurements of the matter power spectra shown in the upper
panel of Fig. 5 at three redshifts, $z=0,1,2$. Similar to the linear power
spectra shown in Fig. 3, we plot the $bc$ component only for the Vulcan run
and the full matter power spectra for the Qo’noS and Ferenginar simulations.
The nonlinear evolution for the $\Lambda$CDM (Qo’noS) and the neutrino
simulation (Vulcan) are very similar, while the dynamical dark energy model in
the Ferenginar run exhibits larger nonlinear power at small scales compared to
the other two simulations. All three cosmological models are normalized with
respect to the same value for $\sigma_{8}$ which leads to closer agreement at
$z=0$ compared to the higher redshifts. The middle panel of Fig. 5 shows the
ratio of the Vulcan and Qo’noS simulation for the three redshifts and the
lower panel shows the ratio of the Ferenginar and Qo’noS simulations. The
difference between the $bc$ component of the neutrino simulation and the
$\Lambda$CDM simulation is larger on large scales (small $k$ values) while for
the dynamical dark energy, there is also a difference at small scales. The
dynamical dark energy model exhibits more power than the other two models. The
excess power for the Ferenginar run at early redshifts manifests itself in
other structure formation measurements as we will show later.

Figure 5: Upper panel: Power spectra at three redshifts, $z=0,1$ and $z=2$.
For Qo’noS and Ferenginar, the full matter power spectrum is shown while for
Vulcan, we show only the $bc$ component of the power spectrum. The difference
between Qo’noS and Vulcan is very small. The difference with respect to
Ferenginar decreases at lower redshift. Middle and lower panel: Ratio of power
spectra for the same three redshifts. The $\Lambda$CDM simulation (Qo’noS) is
used as base in all cases.

In Fig. 6 we compare the simulation results at three redshifts with the
CosmicEmu, an emulator released in Moran et al.
(2023)888https://github.com/lanl/CosmicEmu and based on the Mira-Titan
Universe simulations (Heitmann et al., 2016). Using emulators to predict
cosmological measurements from simulations over a wide range of cosmologies
has become a commonly used approach after it was introduced in the
cosmological context in Heitmann et al. (2006); Habib et al. (2007). The
emulator schemes include a wide range of predictions, for e.g. matter power
spectra (Lawrence et al., 2010; Heitmann et al., 2014; Euclid Collaboration et
al., 2019; Angulo et al., 2021; Moran et al., 2023), mass functions
(McClintock et al., 2019; Bocquet et al., 2020), and galaxy correlation
functions (Zhai et al., 2019; Kwan et al., 2023). The major advantage of
emulators over fitting functions such as Halofit (Smith et al., 2003) for the
power spectrum or the Tinker mass function prediction (Tinker et al., 2008) is
their high accuracy, even away from $\Lambda$CDM. In our comparison here, the
accuracy of the power spectrum prediction is at the $1\%$ level at $z=0$ for
all three simulations, covering $\Lambda$CDM, neutrinos and a varying dark
energy equation of state. Overall, the simulations and emulator predictions
across models and redshifts are in agreement at better than 5% out to $k\sim
5h$Mpc-1, consistent with the results reported in Moran et al. (2023). The
error bars shown in Fig. 6 are derived by assuming a Gaussian Random Field
with the measured power spectrum where the standard deviation of power in the
$i$th wavenumber bin is given by $\sigma_{P(k_{i})}={P(k_{i})}/{\sqrt{N_{\rm
modes}(k_{i})}}$.

Figure 6: Comparison of the simulation results to the CosmicEmu introduced in
Moran et al. (2023). The upper panel shows the comparison of the emulator to
the Qo’nos simulation, the middle panel the Vulcan results and the last panel
shows the Ferenginar results. In all three cases, three redshifts are shown
(purple: $z=2$, light blue: $z=1$ and light green: $z=0$). The best agreement
is achieved for $z=0$ at better than 1% throughout the $k$-range considered.
The agreement degrades with higher redshifts (4% at $z=2$) which is consistent
with the findings in Moran et al. (2023). The error bars represent the
uncertainties due to cosmic variance.

###  3.2 Matter Correlation Function

Following on, we show results for the two-point matter correlation function,
$\xi_{mm}(r)$, at $z=0$ in Fig. 7, the equivalent of the power spectrum
measurement in real space. We follow the same strategy to measure
$\xi_{mm}(r)$ as described in Heitmann et al. (2021) and refer the reader to
that publication for details. To make the evaluation of the correlation
function from the particles computationally feasible, we only use 0.01% of the
full particle data set. We focus on the region of the baryon acoustic
oscillation peak as this is the major target of ongoing spectroscopic surveys
such as DESI (Fig. 7 shows in the inset the BAO peak region). As for the power
spectra, at $z=0$ the results for the different models are very similar and
the measurement for the Vulcan simulation only includes the CDM and baryon
components of the simulation. On all scales, the Qo’noS simulation has the
lowest amplitude, followed by the Ferenginar simulation. The Vulcan results
are the highest. On smaller scales before the BAO peak ($20h^{-1}$Mpc$<r<$
80$h^{-1}$Mpc) the Qon’oS and Ferenginar results are very similar, while the
amplitude of the Vulcan result is higher at the 1-2% level.

Figure 7: Matter correlation functions for all three models at $z=0$. The
small inset image shows the difference for the three models at the BAO peak
location. Overall, at this redshift, the differences between the three models
are minimal.

###  3.3 Halo-Matter Correlation Function

The halo-matter correlation functions for all three models are shown in Fig. 8
for five different halo mass bins at redshift $z=0$. We use the FOF halo
sample for this measurement. The ratios for Vulcan and Ferenginar with respect
to the Qo’noS result are shown in Fig. 9. The difference between the two
$\Lambda$CDM simulations with (Vulcan) and without neutrinos (Qo’noS) shown in
the upper panel, is very small, at the percent level. In the largest halo mass
bin there is some noise on small scales due to the limited number of halos.
The differences between the $\Lambda$CDM and the dynamical dark energy model
(Ferenginar) are small on large scales (down to $\sim 1h^{-1}$Mpc) and
increase on small scales up to 30%. Again, the mass bin with the cluster sized
objects is noisy due to the small number of halos. On larger scales, beyond
$r=2h^{-1}$Mpc the differences between the halo-matter correlation functions
are minimal for the different cosmologies.

Figure 8: Halo-matter correlation function for five mass bins, separately
shown for each model at $z=0$.

Figure 9: Ratios of the halo-matter correlation function for five mass bins at
$z=0$. Upper plot: Ratio between Vulcan and Qo’noS, lower plot: ratio between
Ferenginar and Qo’nos. The effect of the neutrinos in the Vulcan simulation is
minimal while the effect of the dynamical dark energy in Ferenginar changes
the halo-matter correlation by up to 30% in the largest halo mass bin at small
distances.

###  3.4 Halo Mass Function

Before we show results for the halo mass function, we compare the number of
particles in halos as a function of redshift $z$ in Fig. 10. The halo model
assumes that all particles reside in halos. In a simulation, due to the finite
mass resolution, this is not the case. It is interesting, however, to
investigate what fraction of particles do reside in halos and how this
fraction depends on cosmology. This measurement is very similar to the halo
mass function in spirit but has in addition information about the time
evolution. However, this measurement does not contain information about the
number of halos in different mass bins. The upper panel in Fig. 10 shows the
direct comparison of the three simulations. The difference between Qo’noS
($\Lambda$CDM) and Vulcan ($\Lambda$CDM + massive neutrinos) is very small,
the points are basically on top of each other. For the Vulcan simulation, we
do not include neutrino particles since our simulation approach only evolves
the $bc$ component of the matter in the simulation. The dynamical dark energy
model (Ferenginar, light blue) exhibits more particles in halos at earlier
redshift – this is consistent with the excess in the nonlinear power spectrum
shown in Fig. 5. At later stages the difference decreases. The lower panel
measures the ratio of the simulations, showing that indeed the Ferenginar
simulation has early on (at $z=4$) 40% more particles in halos than the
$\Lambda$CDM models, with the ratio declining steadily until reaching $z=0$.

Figure 10: Upper panel: Evolution of the number of particles in halos. Shown
is the fraction of particles in halos as a function of redshift out to $z=4$.
The Qo’nos and Vulcan simulation have very similar fractions, while the
Ferenginar simulation has a larger fraction of particles reside in halos, in
particular at earlier times. Lower panel: Ratio of the number of particles in
halos with respect to the $\Lambda$CDM model, Qo’noS.

Figure 11: Halo mass function for $M_{200c}$ at three different redshifts for
the three cosmologies. The Qo’noS and Vulcan results are basically
indistinguishable at all redshifts. At $z=0$ the difference between the
cosmologies is very small, at $z=2$ the dynamical dark energy model shows an
excess of halos at all masses. This is consistent with the results shown in
Fig. 10.

Next, we show the results for the mass function itself. We focus on the
measurements of SO masses, M200c. Our halo finder first identifies FOF halos
using a linking length of $b=0.168$ and determines the potential minimum for
each halo as its center. For FOF halos with more than 500 particles, we
measure M200c. In Fig. 11 we show the results for the SO halo mass function at
three redshifts, $z=0,1,2$. The error bars provided are simply Poisson error
bars. The two $\Lambda$CDM cosmologies, with and without neutrinos, are
basically identical at all three redshifts and over the full halo mass range.
The Ferenginar simulation has more halos in particular at early redshifts. At
$z=0$, the halo mass functions of all three simulations are very similar. The
finding for the Fereginar simulation, a higher halo mass function at early
times and less of a difference at $z=0$ compared to the $\Lambda$CDM
simulations, is consistent with the large-scale structure measurements we have
shown so far. The slowing-down of structure formation over time is due to a
dark energy component that sets in later for the Ferenginar simulation.

Figure 12: Comparison of the simulation results to the CosmicEmu introduced in
Bocquet et al. (2020). The solid lines indicate 5% uncertainty. We only show
mass bins with at least 100 halos. The upper panel shows the comparison of the
emulator to the Qo’nos simulation, the middle panel the Vulcan results and the
last panel shows the Ferenginar results. In all three cases, three redshifts
are shown (purple: $z=2$, light blue: $z=1$ and light green: $z=0$). The best
agreement is achieved for $z=0$ at better than 2-3% for most of the mass range
considered. The agreement degrades with higher redshifts ($\sim$10% at $z=2$)
which is in agreement with the findings in Bocquet et al. (2020). We use
Poisson error bars.

Similar to the power spectrum analysis, we show a comparison of our results
with a mass function emulator from Bocquet et al. (2020) in Fig. 12. The
emulator provides predictions for SO masses, $M_{200c}$, for group and
cluster-sized objects, starting at a mass of $M_{200c}=10^{13}h^{-1}$M⊙.
Accordingly, the $x$-axis in Fig. 12 is reduced from the full mass function
plot shown in Fig. 11. Results for three redshifts are presented for the
comparison and for all three models. We restrict the upper mass end to bins
that have at least 100 halos. For $z=0$ we find very good agreement between
the simulations and the emulator, at the 2-3% level. At redshift $z=1$, the
agreement is still over most of the mass range better than 5% (indicated by
the solid lines in the figure). The accuracy degrades somewhat for $z=2$.
These findings are fully consisent with the accuracy reported in Bocquet et
al. (2020).

###  3.5 Concentration-Mass Relation

The final example measurement we provide is for the concentration-mass (c-M)
relation at three redshifts, $z=0,1,2$. The concentration is obtained by
fitting a Navarro-Frenk-White (NFW) profile (Navarro et al., 1996, 1997) to
each halo via

$\rho(r)=\frac{\delta\rho_{\rm{crit}}}{(r/r_{s})(1+r/r_{s})^{2}}.$ (19)

Here $\delta$ is a characteristic dimensionless density, and $r_{s}$ is the
scale radius of the NFW profile. The concentration of a halo is then simply
given by $c_{\Delta}=r_{\Delta}/r_{s}$, where $\Delta$ is the overdensity
measured with respect to the critical density $\rho_{c}$, and $r_{\Delta}$ is
the radius at which the enclosed mass is equal to the volume of the sphere
times the density. As for the mass function, we choose $\Delta=200$, leading
to $c_{200c}=R_{200c}/r_{s}$. Details about our fitting approach can be found
in Bhattacharya et al. (2013) and Child et al. (2018). The error bars shown
combine the error contributions from the individual concentration measurements
and the Poisson error due to the finite number of halo in each bin. For a
detailed discussion, the reader is referred to Bhattacharya et al. (2013).

The results are shown in Fig. 13 for the three simulations and three
redshifts. The concentration-mass relation for the Qo’noS and Vulcan
simulations are basically indistinguishable at all three redshifts. This is to
be expected given our assumption that the neutrino mass considered here is
small enough to allow us to neglect their nonlinear impact on structure
formation (for a more extensive discussion about the effect of neutrinos on
group and cluster sized halos see Biswas et al. (2019)).

The concentration-mass relation for the Ferenginar simulation shows a
different behavior than the $\Lambda$CDM simulations. While at high redshift,
$z=2$, all three simulations have a very similar c-M relation, at late times,
the dynamical dark energy model exhibits larger concentrations, a $\sim 20\%$
increase, than the other two simulations. The mass function measurements (Fig.
11) show that the halos in the Ferenginar simulation formed earlier. This
early formation then translates into tighter, more concentrated, halos at
later times. This result is also consistent with the matter power spectrum
measurements (Fig. 5) which has more power on small scales in the Ferenginar
simulation compared to the other two simulations.

Figure 13: Concentration-mass relation at three different redshifts, $z=0,1,2$
for the three different cosmologies. The difference between the Qo’noS and
Vulcan simulations is minimal. At high redshift ($z=2$, left panel), the
differences between the three models are small. The Ferenginar simulation
shows a more interesting behavior: While at early redshifts, the
concentration-mass relation is almost the same as for the $\Lambda$CDM cases,
at lower redshifts the concentrations are considerably higher.

##  4 Data Release

Several HACC simulation outputs have been released via the HACC Simulation
Data Portal999https://cosmology.alcf.anl.gov/, described in Heitmann et al.
(2019b). The HACC Simulation Data Portal is hosted at the Argonne Leadership
Computing Facility (ALCF) and data transfer from the portal is incorporated
via Globus101010https://www.globus.org/. As part of this paper, we release a
selected set of outputs from the three
simulations111111https://cosmology.alcf.anl.gov/transfer/newworlds. The data
selection is the similar as for previous large simulation releases. We focus
on FOF halo properties and downsampled particle files:

$z=\\{0.0,0.05,0.21,0.50,0.54,0.78,0.86,1.43,1.49\\}.$ (20)

For each of the redshifts listed, we provide access to 1% of the particle
data. For each particle, a particle ID, comoving positions $(x,y,z)_{\rm COM}$
in $h^{-1}$Mpc, and comoving peculiar velocities $(vx,vy,vz)_{\rm COM}$ in
km/s are given. The particles were randomly selected for the first HACC
output, but for all following outputs, the same particles are chosen. In
addition, halo information is made available. The halos are identified with an
FOF halo finder using a linking length of $b=0.168$. For each halo, its center
position and velocity are provided. For the center definition, both the
gravitational potential minimum and the center of mass measurement are
available. For the FOF halo definition, these two quantities can be
significantly different, in particular for large, unrelaxed halos. The halo
mass is given in $h^{-1}$M⊙.

##  5 Summary and Outlook

In this paper we introduced the “New Worlds Simulations”, three new gravity-
only, extreme-scale simulations. Each simulation covers a volume of
(3$h^{-1}$Gpc)3 and evolves 12,2883 particles, making them some of the largest
simulations of structure formation currently available. A good value for the
mass resolution, $m_{p}\sim 10^{9}h^{-1}$M⊙, in combination with the large
volume make these simulations particularly useful for investigations related
to ongoing and upcoming cosmological surveys. The three simulations have the
same initial phases and span three cosmological models, a $\Lambda$CDM
cosmology informed by Planck constraints, a $\Lambda$CDM cosmology including
massive neutrinos, and a dynamical dark energy model. The dark energy equation
of state for the dynamical dark energy model is expressed via the simple
$(w_{0},w_{a})$ parameterization, commonly used by ongoing surveys.

The simulations were carried out with HACC on the Summit supercomputer at
OLCF. HACC provides a set of analysis capabilities that are carried out during
simulation runtime. In this paper, we showed selected results from this in
situ analysis for both particle and halo statistics. The results are
consistent with previous findings as demonstrated via comparisons to publicly
available emulators. Several interesting differences between the results from
the three simulations were discussed in the light of the different physics
included in the simulations. Most notably, the varying dark energy equation of
state led to stronger clustering at early times, translating for example into
higher concentration values and halo counts at higher redshifts. All
simulations were normalized to the same value of $\sigma_{8}$ at $z=0$ and
results at the final time are very similar. The addition of a massive neutrino
component did not have a large effect on the different statistics we measured
due to the small neutrino mass chosen. Overall, the small differences between
the measurements we presented in this paper hint at the difficult tasks laying
ahead of us to find new fundamental physics, clearly disentangled from
systematic effects. The availability of detailed data at different epochs in
the evolution of the Universe will be crucial to solidly confirm or reject
physics beyond the Standard Model of Cosmology. We publicly release some of
our simulation results to enable the creation of synthetic sky catalogs to
help the community to further explore the subtle differences between the
cosmologies considered in this paper.

###### Acknowledgements.

We thank Kelly Moran and Sebastian Bocquet for discussions about the emulator
comparisons. Argonne National Laboratory’s work was supported under the U.S.
Department of Energy contract DE-AC02-06CH11357. This research used resources
of the Oak Ridge Leadership Computing Facility at the Oak Ridge National
Laboratory, which is supported by the Office of Science of the U.S. Department
of Energy under Contract No. DE-AC05-00OR22725. This research used resources
of the Argonne Leadership Computing Facility, a U.S. Department of Energy
(DOE) Office of Science user facility at Argonne National Laboratory and is
based on research supported by the U.S. DOE Office of Science-Advanced
Scientific Computing Research Program, under Contract No. DE-AC02-06CH11357.

## References

  * Abate et al. (2012) Abate A. Aldering, G., Allan, S.W. et al. [LSST DESC], arXiv:12110310 [astro-ph.CO]
  * Abell et al. (2009) Abell, P.A., Allison, J., Anderson, S.F. et al. [LSST Science Collaborations and LSST Project Collaboration], arXiv:0912.0201 [astro-ph.IM].
  * Planck Collaboration (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Aghamousa et al. (2016) Aghamousa, A., Aguilar, J., Ahlen, S. et al. [DESI Collaboration], arXiv:1611.00036 [astro-ph.IM]
  * Angulo et al. (2021) Angulo, R. E., Zennaro, M., Contreras, S., et al. 2021, Mon. Not. Roy. Astron. Soc. 507, 5869
  * Bhattacharya et al. (2013) Bhattacharya, S., Habib, S., Heitmann, K., & Vikhlinin, A. 2013, Astrophys. J. 766, 32
  * Biswas et al. (2019) Biswas, R., Heitmann, K., Habib, S., et al. 2019, arXiv:1901.10690
  * Bocquet et al. (2020) Bocquet, S., Heitmann, K., Habib, S., et al. 2020, ApJ 901, 5
  * Child et al. (2018) Child, H.L., Habib, S., Heitmann, K., et al. 2018, Astrophys. J. 859, 55
  * Crocce et al. (2015) Crocce, M., Castander, F. J., Gaztañaga, E., et al. 2015, MNRAS, 453, 1513
  * Davis et al. (1985) Davis, M., Efstathiou, G., Frenk, C.S., White, S.D.M. 1985, ApJ, 292, 371
  * DeRose et al. (2019) DeRose, J., Wechsler, R. H., Becker, M. R., et al. 2019, arXiv:1901.02401
  * Doré et al. (2019) Doré, O., Hirata C, Wang Y. et al, 2019, arXiv:1904.01174 [astro-ph.CO]
  * Euclid Collaboration et al. (2019) Euclid Collaboration, Knabenhans, M., Stadel, J., et al. 2019, MNRAS, 484, 5509
  * Gómez-Valent & Amendola (2018) Gómez-Valent, A. & Amendola, L. 2018, J. Cosmology Astropart. Phys, 2018, 051
  * Habib et al. (2007) Habib, S., Heitmann, K., Higdon, D., et al. 2007, Phys. Rev. D, 76, 083503
  * Habib et al. (2016) Habib, S., Pope, A., Finkel, H., Frontiere, N., et al. 2016, New Astronomy, 42, 49
  * Heitmann et al. (2006) Heitmann, K., Higdon, D., Nakhleh, C., et al. 2006, ApJ, 646, L1
  * Heitmann et al. (2014) Heitmann, K., Lawrence, E., Kwan, J., et al. 2014, ApJ, 780, 111
  * Heitmann et al. (2016) Heitmann, K., Bingham, D., Lawrence, E., et al. 2016, ApJ, 820, 108
  * Heitmann et al. (2019a) Heitmann K., Finkel, H., Pope, A. et al, 2019, ApJS, 245, 16
  * Heitmann et al. (2019b) Heitmann, K., Uram, T., Finkel, H. et al., 2019, ApJS, 244, 17
  * Heitmann et al. (2021) Heitmann, K., Frontiere, N., Rangel, E., et al. 2021, ApJS, 252, 19
  * Korytov et al. (2023) Korytov, D., Rangel, E., Bleem, L., et al. 2023, The Open Journal of Astrophysics, 6, 24. doi:10.21105/astro.2302.04194
  * Kwan et al. (2023) Kwan, J., Saito, S., Leauthaud, A., et al. 2023, ApJ, 952, 80
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193 [astro-ph.CO]
  * Lawrence et al. (2010) Lawrence, E., Heitmann, K., White, M., et al. 2010, ApJ, 713, 1322
  * Lewis, Challinor & Lasenby (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ538, 473
  * LSST Dark Energy Science Collaboration et al. (2021) LSST Dark Energy Science Collaboration (LSST DESC), Abolfathi, B., Alonso, D., et al. 2021, ApJS, 253, 31
  * Maksimova et al. (2021) Maksimova, N. A., Garrison, L. H., Eisenstein, D. J., et al. 2021, MNRAS, 508, 4017
  * McClintock et al. (2019) McClintock, T., Rozo, E., Becker, M. R., et al. 2019, ApJ, 872, 53
  * Moran et al. (2023) K. R. Moran, K. Heitmann, E. Lawrence, S. Habib, D. Bingham, A. Upadhye, J. Kwan, D. Higdon and R. Payne 2023 MNRAS, 520, 3443
  * Navarro et al. (1996) Navarro, J.F., Frenk, C.S., & White, S.D.M. 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro, J.F., Frenk, C.S., & White, S.D.M. 1997, ApJ, 490, 493
  * Potter, Stadel, & Teyssier (2017) Potter, D., Stadel, J., & Teyssier, R. 2017, Computational Astrophysics and Cosmology 4, 2
  * Rangel et al. (2020) Rangel, E., Frontiere, N., Habib, S., et al. 2020, arXiv:2008.08519. doi:10.48550/arXiv.2008.08519
  * Smith et al. (2003) Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709
  * Upadhye et al. (2014) Upadhye, A., Biswas, R., Pope, A., et al. 2014, Phys. Rev. D, 89, 103515
  * Zel’dovich (1970) Zel’dovich, Y.B. 1970, A&A, 5, 84
  * Zhai et al. (2019) Zhai, Z., Tinker, J. L., Becker, M. R., et al. 2019, ApJ, 874, 95