The Last Journey is a large-volume, gravity-only, cosmological N-body simulation evolving more than 1.24 trillion particles in a periodic box with a side-length of 5.025Gpc. It was implemented using the HACC simulation and analysis framework on the BG/Q system, Mira. The cosmological parameters are chosen to be consistent with the results from the Planck satellite. A range of analysis tools have been run in situ to enable a diverse set of science projects, and at the same time, to keep the resulting data amount manageable. Analysis outputs have been generated starting at redshift z~10 to allow for construction of synthetic galaxy catalogs using a semi-analytic modeling approach in post-processing. As part of our in situ analysis pipeline we employ a new method for tracking halo sub-structures, introducing the concept of subhalo cores. The production of multi-wavelength synthetic sky maps is facilitated by generating particle lightcones in situ, also beginning at z~10. We provide an overview of the simulation set-up and the generated data products; a first set of analysis results is presented. A subset of the data is publicly available.

We construct an emulator for the halo mass function over group and cluster mass scales for a range of cosmologies, including the effects of dynamical dark energy and massive neutrinos. The emulator is based on the recently completed Mira-Titan Universe suite of cosmological $N$-body simulations. The main set of simulations spans 111 cosmological models with 2.1 Gpc boxes. We extract halo catalogs in the redshift range $z=[0.0, 2.0]$ and for masses $M_{200\mathrm{c}}\geq 10^{13}M_\odot/h$. The emulator covers an 8-dimensional hypercube spanned by {$\Omega_\mathrm{m}h^2$, $\Omega_\mathrm{b}h^2$, $\Omega_\nu h^2$, $\sigma_8$, $h$, $n_s$, $w_0$, $w_a$}; spatial flatness is assumed. We obtain smooth halo mass functions by fitting piecewise second-order polynomials to the halo catalogs and employ Gaussian process regression to construct the emulator while keeping track of the statistical noise in the input halo catalogs and uncertainties in the regression process. For redshifts $z\lesssim1$, the typical emulator precision is better than $2\%$ for $10^{13}-10^{14} M_\odot/h$ and $<10\%$ for $M\simeq 10^{15}M_\odot/h$. For comparison, fitting functions using the traditional universal form for the halo mass function can be biased at up to 30\% at $M\simeq 10^{14}M_\odot/h$ for $z=0$. Our emulator is publicly available at \url{https://github.com/SebastianBocquet/MiraTitanHMFemulator}.

We describe the first major public data release from cosmological simulations carried out with Argonne's HACC code. This initial release covers a range of datasets from large gravity-only simulations. The data products include halo information for multiple redshifts, down-sampled particles, and lightcone outputs. We provide data from two very large LCDM simulations as well as beyond-LCDM simulations spanning eleven w0-wa cosmologies. Our release platform uses Petrel, a research data service, located at the Argonne Leadership Computing Facility. Petrel offers fast data transfer mechanisms and authentication via Globus, enabling simple and efficient access to stored datasets. Easy browsing of the available data products is provided via a web portal that allows the user to navigate simulation products efficiently. The data hub will be extended by adding more types of data products and by enabling computational capabilities to allow direct interactions with simulation results.

This paper introduces cosmoDC2, a large synthetic galaxy catalog designed to support precision dark energy science with the Large Synoptic Survey Telescope (LSST). CosmoDC2 is the starting point for the second data challenge (DC2) carried out by the LSST Dark Energy Science Collaboration (LSST DESC). The catalog is based on a trillion-particle, 4.225 Gpc^3 box cosmological N-body simulation, the `Outer Rim' run. It covers 440 deg^2 of sky area to a redshift of z=3 and is complete to a magnitude depth of 28 in the r-band. Each galaxy is characterized by a multitude of properties including stellar mass, morphology, spectral energy distributions, broadband filter magnitudes, host halo information and weak lensing shear. The size and complexity of cosmoDC2 requires an efficient catalog generation methodology; our approach is based on a new hybrid technique that combines data-driven empirical approaches with semi-analytic galaxy modeling. A wide range of observation-based validation tests has been implemented to ensure that cosmoDC2 enables the science goals of the planned LSST DESC DC2 analyses. This paper also represents the official release of the cosmoDC2 data set, including an efficient reader that facilitates interaction with the data.

A challenging requirement posed by next-generation observations is a firm theoretical grasp of the impact of baryons on structure formation. Cosmological hydrodynamic simulations modeling gas physics are vital in this regard. A high degree of modeling flexibility exists in this space making it important to explore a range of methods in order to gauge the accuracy of simulation predictions. We present results from the first cosmological simulation using Conservative Reproducing Kernel Smoothed Particle Hydrodynamics (CRK-SPH). We employ two simulations: one evolved purely under gravity and the other with non-radiative hydrodynamics. Each contains 2x2304^3 cold dark matter plus baryon particles in an 800 Mpc/h box. We compare statistics to previous non-radiative simulations including power spectra, mass functions, baryon fractions, and concentration. We find self-similar radial profiles of gas temperature, entropy, and pressure and show that a simple analytic model recovers these results to better than 40% over two orders of magnitude in mass. We quantify the level of non-thermal pressure support in halos and demonstrate that hydrostatic mass estimates are biased low by 24% (10%) for halos of mass 10^15 (10^13) Msun/h. We compute angular power spectra for the thermal and kinematic Sunyaev-Zel'dovich effects and find good agreement with the low-l Planck measurements. Finally, artificial scattering between particles of unequal mass is shown to have a large impact on the gravity-only run and we highlight the importance of better understanding this issue in hydrodynamic applications. This is the first in a simulation campaign using CRK-SPH with future work including subresolution gas treatments.

Cosmological cluster-scale strong gravitational lensing probes the mass distribution of the dense cores of massive dark matter halos and the structures along the line of sight from background sources to the observer. It is frequently assumed that the primary lens mass dominates the lensing, with the contribution of secondary masses along the line of sight being neglected. Secondary mass structures may, however, affect both the detectability of strong lensing in a given survey and modify the properties of the lensing that is detected. In this paper, we utilize a large cosmological N-body simulation and a multiple lens plane (and many source planes) ray-tracing technique to quantify the influence of line of sight halos on the detectability of cluster-scale strong lensing in a cluster sample with a mass limit that encompasses current cluster catalogs from the South Pole Telescope. We extract both primary and secondary halos from the "Outer Rim" simulation and consider two strong lensing realizations: one with only the primary halos included, and the other contains all secondary halos down to a mass limit. In both cases, we use the same source information extracted from the Hubble Ultra Deep Field, and create realistic lensed images consistent with moderately deep ground-based imaging. The results demonstrate that down to the mass limit considered the total number of lenses is boosted by about 13-21% when considering the complete multi-halo lightcone. The increment in strong lens counts peaks at lens redshifts of 0.6 approximately with no significant effect at z<0.3. The strongest trends are observed relative to the primary halo mass, with no significant impact in the most massive quintile of the halo sample, but increasingly boosting the observed lens counts toward small primary halo masses, with an enhancement greater than 50% in the least massive quintile of the halo masses considered.

We describe the Outer Rim cosmological simulation, one of the largest high-resolution N-body simulations performed to date, aimed at promoting science to be carried out with large-scale structure surveys. The simulation covers a volume of (4.225Gpc)^3 and evolves more than one trillion particles. It was executed on Mira, a BlueGene/Q system at the Argonne Leadership Computing Facility. We discuss some of the computational challenges posed by a system like Mira, a many-core supercomputer, and how the simulation code, HACC, has been designed to overcome these challenges. We have carried out a large range of analyses on the simulation data and we report on the results as well as the data products that have been generated. The full data set generated by the simulation totals more than 5PB of data, making data curation and data handling a large challenge in of itself. The simulation results have been used to generate synthetic catalogs for large-scale structure surveys, including DESI and eBOSS, as well as CMB experiments. A detailed catalog for the LSST DESC data challenges has been created as well. We publicly release some of the Outer Rim halo catalogs, downsampled particle information, and lightcone data.

Profiles of dark matter-dominated halos at the group and cluster scales play an important role in modern cosmology. Using results from two very large cosmological $N$-body simulations, which increase the available volume at their mass resolution by roughly two orders of magnitude, we robustly determine the halo concentration-mass $(c-M)$ relation over a wide range of masses, employing multiple methods of concentration measurement. We characterize individual halo profiles, as well as stacked profiles, relevant for galaxy-galaxy lensing and next-generation cluster surveys; the redshift range covered is $0\leq z \leq 4$, with a minimum halo mass of $M_{200c}\sim2\times10^{11} M_\odot$. Despite the complexity of a proper description of a halo (environmental effects, merger history, nonsphericity, relaxation state), when the mass is scaled by the nonlinear mass scale $M_\star(z)$, we find that a simple non-power-law form for the $c-M/M_\star$ relation provides an excellent description of our simulation results across eight decades in $M/M_{\star}$ and for $0\leq z \leq 4$. Over the mass range covered, the $c-M$ relation has two asymptotic forms: an approximate power law below a mass threshold $M/M_\star\sim 500-1000$, transitioning to a constant value, $c_0\sim 3$ at higher masses. The relaxed halo fraction decreases with mass, transitioning to a constant value of $\sim 0.5$ above the same mass threshold. We compare Navarro-Frenk-White (NFW) and Einasto fits to stacked profiles in narrow mass bins at different redshifts; as expected, the Einasto profile provides a better description of the simulation results. At cluster scales at low redshift, however, both NFW and Einasto profiles are in very good agreement with the simulation results, consistent with recent weak lensing observations.

We introduce a new cosmic emulator for the matter power spectrum covering eight cosmological parameters. Targeted at optical surveys, the emulator provides accurate predictions out to a wavenumber k~5/Mpc and redshift z<=2. Besides covering the standard set of LCDM parameters, massive neutrinos and a dynamical dark energy of state are included. The emulator is built on a sample set of 36 cosmological models, carefully chosen to provide accurate predictions over the wide and large parameter space. For each model, we have performed a high-resolution simulation, augmented with sixteen medium-resolution simulations and TimeRG perturbation theory results to provide accurate coverage of a wide k-range; the dataset generated as part of this project is more than 1.2Pbyte. With the current set of simulated models, we achieve an accuracy of approximately 4%. Because the sampling approach used here has established convergence and error-control properties, follow-on results with more than a hundred cosmological models will soon achieve ~1% accuracy. We compare our approach with other prediction schemes that are based on halo model ideas and remapping approaches. The new emulator code is publicly available.

The pairwise kinematic Sunyaev-Zel'dovich (kSZ) signal from galaxy clusters is a probe of their line-of-sight momenta, and thus a potentially valuable source of cosmological information. In addition to the momenta, the amplitude of the measured signal depends on the properties of the intra-cluster gas and observational limitations such as errors in determining cluster centers and redshifts. In this work we simulate the pairwise kSZ signal of clusters at z<1, using the output from a cosmological N-body simulation and including the properties of the intra-cluster gas via a model that can be varied in post-processing. We find that modifications to the gas profile due to star formation and feedback reduce the pairwise kSZ amplitude of clusters by ~50%, relative to the naive 'gas traces mass' assumption. We demonstrate that mis-centering can reduce the overall amplitude of the pairwise kSZ signal by up to 10%, while redshift errors can lead to an almost complete suppression of the signal at small separations. We confirm that a high-significance detection is expected from the combination of data from current-generation, high-resolution CMB experiments, such as the South Pole Telescope, and cluster samples from optical photometric surveys, such as the Dark Energy Survey. Furthermore, we forecast that future experiments such as Advanced ACTPol in conjunction with data from the Dark Energy Spectroscopic Instrument will yield detection significances of at least 20{\sigma}, and up to 57{\sigma} in an optimistic scenario. Our simulated maps are publicly available at: http://www.hep.anl.gov/cosmology/ksz.html

This draft report summarizes and details the findings, results, and recommendations derived from the ASCR/HEP Exascale Requirements Review meeting held in June, 2015. The main conclusions are as follows. 1) Larger, more capable computing and data facilities are needed to support HEP science goals in all three frontiers: Energy, Intensity, and Cosmic. The expected scale of the demand at the 2025 timescale is at least two orders of magnitude -- and in some cases greater -- than that available currently. 2) The growth rate of data produced by simulations is overwhelming the current ability, of both facilities and researchers, to store and analyze it. Additional resources and new techniques for data analysis are urgently needed. 3) Data rates and volumes from HEP experimental facilities are also straining the ability to store and analyze large and complex data volumes. Appropriately configured leadership-class facilities can play a transformational role in enabling scientific discovery from these datasets. 4) A close integration of HPC simulation and data analysis will aid greatly in interpreting results from HEP experiments. Such an integration will minimize data movement and facilitate interdependent workflows. 5) Long-range planning between HEP and ASCR will be required to meet HEP's research needs. To best use ASCR HPC resources the experimental HEP program needs a) an established long-term plan for access to ASCR computational and data resources, b) an ability to map workflows onto HPC resources, c) the ability for ASCR facilities to accommodate workflows run by collaborations that can have thousands of individual members, d) to transition codes to the next-generation HPC platforms that will be available at ASCR facilities, e) to build up and train a workforce capable of developing and using simulations and analysis to support HEP scientific research on next-generation systems.

Large-scale structure surveys in the coming years will measure the redshift-space power spectrum to unprecedented accuracy, allowing for powerful new tests of the LambdaCDM picture as well as measurements of particle physics parameters such as the neutrino masses. We extend the Time-RG perturbative framework to redshift space, computing the power spectrum P_s(k,mu) in massive neutrino cosmologies with time-dependent dark energy equations of state w(z). Time-RG is uniquely capable of incorporating scale-dependent growth into the P_s(k,mu) computation, which is important for massive neutrinos as well as modified gravity models. Although changes to w(z) and the neutrino mass fraction both affect the late-time scale-dependence of the non-linear power spectrum, we find that the two effects depend differently on the line-of-sight angle mu. Finally, we use the HACC N-body code to quantify errors in the perturbative calculations. For a LambdaCDM model at redshift z=1, our procedure predicts the monopole~(quadrupole) to 1% accuracy up to a wave number 0.19h/Mpc (0.28h/Mpc), compared to 0.08h/Mpc (0.07h/Mpc) for the Kaiser approximation and 0.19h/Mpc (0.16h/Mpc) for the current state-of-the-art perturbation scheme. Our calculation agrees with the simulated redshift-space power spectrum even for neutrino masses above the current bound, and for rapidly-evolving dark energy equations of state, |dw/dz| ~ 1. Along with this article, we make our redshift-space Time-RG implementation publicly available as the code redTime.

The halo occupation distribution (HOD) approach has proven to be an effective method for modeling galaxy clustering and bias. In this approach, galaxies of a given type are probabilistically assigned to individual halos in N-body simulations. In this paper, we present a fast emulator for predicting the fully nonlinear galaxy power spectrum over a range of freely specifiable HOD modeling parameters. The emulator is constructed using results from 100 HOD models run on a large LCDM N-body simulation, with Gaussian Process interpolation applied to a PCA-based representation of the galaxy power spectrum. The total error is currently ~3% (~2% in the simulation and ~1% in the emulation process) from z=1 to z=0, over the considered parameter range. We use the emulator to investigate parametric dependencies in the HOD model, as well as the behavior of galaxy bias as a function of HOD parameters. The emulator is publicly available at http://www.hep.anl.gov/cosmology/CosmicEmu/emu.html.

Ground and space-based sky surveys enable powerful cosmological probes based on measurements of galaxy properties and the distribution of galaxies in the Universe. These probes include weak lensing, baryon acoustic oscillations, abundance of galaxy clusters, and redshift space distortions; they are essential to improving our knowledge of the nature of dark energy. On the theory and modeling front, large-scale simulations of cosmic structure formation play an important role in interpreting the observations and in the challenging task of extracting cosmological physics at the needed precision. These simulations must cover a parameter range beyond the standard six cosmological parameters and need to be run at high mass and force resolution. One key simulation-based task is the generation of accurate theoretical predictions for observables, via the method of emulation. Using a new sampling technique, we explore an 8-dimensional parameter space including massive neutrinos and a variable dark energy equation of state. We construct trial emulators using two surrogate models (the linear power spectrum and an approximate halo mass function). The new sampling method allows us to build precision emulators from just 26 cosmological models and to increase the emulator accuracy by adding new sets of simulations in a prescribed way. This allows emulator fidelity to be systematically improved as new observational data becomes available and higher accuracy is required. Finally, using one LCDM cosmology as an example, we study the demands imposed on a simulation campaign to achieve the required statistics and accuracy when building emulators for dark energy investigations.

Modeling large-scale sky survey observations is a key driver for the continuing development of high resolution, large-volume, cosmological simulations. We report the first results from the 'Q Continuum' cosmological N-body simulation run carried out on the GPU-accelerated supercomputer Titan. The simulation encompasses a volume of (1300 Mpc)^3 and evolves more than half a trillion particles, leading to a particle mass resolution of ~1.5 X 10^8 M_sun. At this mass resolution, the Q Continuum run is currently the largest cosmology simulation available. It enables the construction of detailed synthetic sky catalogs, encompassing different modeling methodologies, including semi-analytic modeling and sub-halo abundance matching in a large, cosmological volume. Here we describe the simulation and outputs in detail and present first results for a range of cosmological statistics, such as mass power spectra, halo mass functions, and halo mass-concentration relations for different epochs. We also provide details on challenges connected to running a simulation on almost 90% of Titan, one of the fastest supercomputers in the world, including our usage of Titan's GPU accelerators.

Current and future surveys of large-scale cosmic structure are associated with a massive and complex datastream to study, characterize, and ultimately understand the physics behind the two major components of the 'Dark Universe', dark energy and dark matter. In addition, the surveys also probe primordial perturbations and carry out fundamental measurements, such as determining the sum of neutrino masses. Large-scale simulations of structure formation in the Universe play a critical role in the interpretation of the data and extraction of the physics of interest. Just as survey instruments continue to grow in size and complexity, so do the supercomputers that enable these simulations. Here we report on HACC (Hardware/Hybrid Accelerated Cosmology Code), a recently developed and evolving cosmology N-body code framework, designed to run efficiently on diverse computing architectures and to scale to millions of cores and beyond. HACC can run on all current supercomputer architectures and supports a variety of programming models and algorithms. It has been demonstrated at scale on Cell- and GPU-accelerated systems, standard multi-core node clusters, and Blue Gene systems. HACC's design allows for ease of portability, and at the same time, high levels of sustained performance on the fastest supercomputers available. We present a description of the design philosophy of HACC, the underlying algorithms and code structure, and outline implementation details for several specific architectures. We show selected accuracy and performance results from some of the largest high resolution cosmological simulations so far performed, including benchmarks evolving more than 3.6 trillion particles.

Over the next decade, cosmological measurements of the large-scale structure of the Universe will be sensitive to the combined effects of dynamical dark energy and massive neutrinos. The matter power spectrum is a key repository of this information. We extend higher-order perturbative methods for computing the power spectrum to investigate these effects over quasi-linear scales. Through comparison with N-body simulations we establish the regime of validity of a Time-Renormalization Group (Time-RG) perturbative treatment that includes dynamical dark energy and massive neutrinos. We also quantify the accuracy of Standard (SPT), Renormalized (RPT) and Lagrangian Resummation (LPT) perturbation theories without massive neutrinos. We find that an approximation that neglects neutrino clustering as a source for nonlinear matter clustering predicts the Baryon Acoustic Oscillation (BAO) peak position to 0.25% accuracy for redshifts 1 < z < 3, justifying the use of LPT for BAO reconstruction in upcoming surveys. We release a modified version of the public Copter code which includes the additional physics discussed in the paper.

Oscillons are long-lived, localized excitations of nonlinear scalar fields which may be copiously produced during preheating after inflation, leading to a possible oscillon-dominated phase in the early Universe. For example, this can happen after axion monodromy inflation, on which we run our simulations. We investigate the stochastic gravitational wave background associated with an oscillon-dominated phase. An isolated oscillon is spherically symmetric and does not radiate gravitational waves, and we show that the flux of gravitational radiation generated between oscillons is also small. However, a significant stochastic gravitational wave background may be generated during preheating itself (i.e, when oscillons are forming), and in this case the characteristic size of the oscillons is imprinted on the gravitational wave power spectrum, which has multiple, distinct peaks.

Remarkable observational advances have established a compelling cross-validated model of the Universe. Yet, two key pillars of this model -- dark matter and dark energy -- remain mysterious. Sky surveys that map billions of galaxies to explore the `Dark Universe', demand a corresponding extreme-scale simulation capability; the HACC (Hybrid/Hardware Accelerated Cosmology Code) framework has been designed to deliver this level of performance now, and into the future. With its novel algorithmic structure, HACC allows flexible tuning across diverse architectures, including accelerated and multi-core systems. On the IBM BG/Q, HACC attains unprecedented scalable performance -- currently 13.94 PFlops at 69.2% of peak and 90% parallel efficiency on 1,572,864 cores with an equal number of MPI ranks, and a concurrency of 6.3 million. This level of performance was achieved at extreme problem sizes, including a benchmark run with more than 3.6 trillion particles, significantly larger than any cosmological simulation yet performed.

Oscillons are massive, long-lived, localized excitations of a scalar field. We show that in a large class of well-motivated single-field models, inflation is followed by self-resonance, leading to copious oscillon generation and a lengthy period of oscillon domination. These models are characterized by an inflaton potential which has a quadratic minimum and is shallower than quadratic away from the minimum. This set includes both string monodromy models and a class of supergravity inspired scenarios, and is in good agreement with the current central values of the concordance cosmology parameters. We assume that the inflaton is weakly coupled to other fields, so as not to quickly drain energy from the oscillons or prevent them from forming. An oscillon-dominated universe has a greatly enhanced primordial power spectrum on very small scales relative to that seen with a quadratic potential, possibly leading to novel gravitational effects in the early universe.

Analytical arguments suggest that a large class of scalar field potentials permit the existence of oscillons -- pseudo-stable, non-topological solitons -- in three spatial dimensions. In this paper we numerically explore oscillon solutions in three dimensions. We confirm the existence of these field configurations as solutions to the Klein-Gorden equation in an expanding background, and verify the predictions of Amin and Shirokoff for the characteristics of individual oscillons for their model. Further, we demonstrate that significant numbers of oscillons can be generated via fragmentation of the inflaton condensate, consistent with the analysis of Amin. These emergent oscillons can easily dominate the post-inflationary universe. Finally, both analytic and numerical results suggest that oscillons are stable on timescales longer than the post-inflationary Hubble time. Consequently, the post-inflationary universe can contain an effective matter-dominated phase, during which it is dominated by localized concentrations of scalar field matter.

PSpectRe is a C++ program that uses Fourier-space pseudo-spectral methods to evolve interacting scalar fields in an expanding universe. PSpectRe is optimized for the analysis of parametric resonance in the post-inflationary universe, and provides an alternative to finite differencing codes, such as Defrost and LatticeEasy. PSpectRe has both second- (Velocity-Verlet) and fourth-order (Runge-Kutta) time integrators. Given the same number of spatial points and/or momentum modes, PSpectRe is not significantly slower than finite differencing codes, despite the need for multiple Fourier transforms at each timestep, and exhibits excellent energy conservation. Further, by computing the post-resonance equation of state, we show that in some circumstances PSpectRe obtains reliable results while using substantially fewer points than a finite differencing code. PSpectRe is designed to be easily extended to other problems in early-universe cosmology, including the generation of gravitational waves during phase transitions and pre-inflationary bubble collisions. Specific applications of this code will be pursued in future work.