The Tri-Agency Cosmological Simulations (TACS) Task Force was formed when Program Managers from the Department of Energy (DOE), the National Aeronautics and Space Administration (NASA), and the National Science Foundation (NSF) expressed an interest in receiving input into the cosmological simulations landscape related to the upcoming DOE/NSF Vera Rubin Observatory (Rubin), NASA/ESA's Euclid, and NASA's Wide Field Infrared Survey Telescope (WFIRST). The Co-Chairs of TACS, Katrin Heitmann and Alina Kiessling, invited community scientists from the USA and Europe who are each subject matter experts and are also members of one or more of the surveys to contribute. The following report represents the input from TACS that was delivered to the Agencies in December 2018.

Cosmological probes pose an inverse problem where the measurement result is obtained through observations, and the objective is to infer values of model parameters which characterize the underlying physical system -- our Universe. Modern cosmological probes increasingly rely on measurements of the small-scale structure, and the only way to accurately model physical behavior on those scales, roughly 65 Mpc/h or smaller, is via expensive numerical simulations. In this paper, we provide a detailed description of a novel statistical framework for obtaining accurate parameter constraints by combining observations with a very limited number of cosmological simulations. The proposed framework utilizes multi-output Gaussian process emulators that are adaptively constructed using Bayesian optimization methods. We compare several approaches for constructing multi-output emulators that enable us to take possible inter-output correlations into account while maintaining the efficiency needed for inference. Using Lyman alpha forest flux power spectrum, we demonstrate that our adaptive approach requires considerably fewer --- by a factor of a few in Lyman alpha P(k) case considered here --- simulations compared to the emulation based on Latin hypercube sampling, and that the method is more robust in reconstructing parameters and their Bayesian credible intervals.

The intense radiation emitted by luminous quasars dramatically alters the ionization state of their surrounding IGM. This so-called proximity effect extends out to tens of Mpc, and manifests as large coherent regions of enhanced Lyman-$\alpha$ (Ly$\alpha$) forest transmission in absorption spectra of background sightlines. Here we present a novel method based on Ly$\alpha$ forest tomography, which is capable of mapping these quasar `light echoes' in three dimensions. Using a dense grid (10-100) of faint ($m_r\approx24.7\,\mathrm{mag}$) background galaxies as absorption probes, one can measure the ionization state of the IGM in the vicinity of a foreground quasar, yielding detailed information about the quasar's radiative history and emission geometry. An end-to-end analysis - combining cosmological hydrodynamical simulations post-processed with a quasar emission model, realistic estimates of galaxy number densities, and instrument + telescope throughput - is conducted to explore the feasibility of detecting quasar light echoes. We present a new fully Bayesian statistical method that allows one to reconstruct quasar light echoes from thousands of individual low S/N transmission measurements. Armed with this machinery, we undertake an exhaustive parameter study and show that light echoes can be convincingly detected for luminous ($M_{1450} < -27.5\,\mathrm{mag}$ corresponding to $m_{1450} < 18.4\,\mathrm{mag}$ at $z\simeq 3.6$) quasars at redshifts $3<z_\mathrm{QSO}<5$, and that a relative precision better than $20\,\%$ on the quasar age can be achieved for individual objects, for the expected range of ages between 1 Myr and 100 Myr. The observational requirements are relatively modest - moderate resolution ($R\gtrsim750$) multi object spectroscopy at low $\rm{}S/N > 5$ is sufficient, requiring three hour integrations using existing instruments on 8m class telescopes.

Many galaxy formation models predict alignments between galaxy spin and the cosmic web (i.e. the directions of filaments and sheets), leading to intrinsic alignment between galaxies that creates a systematic error in weak lensing measurements. These effects are often predicted to be stronger at high-redshifts ($z\gtrsim1$) that are inaccessible to massive galaxy surveys on foreseeable instrumentation, but IGM tomography of the Ly$\alpha$ forest from closely-spaced quasars and galaxies is starting to measure the $z\sim2-3$ cosmic web with the requisite fidelity. Using mock surveys from hydrodynamical simulations, we examine the utility of this technique, in conjunction with coeval galaxy samples, to measure alignment between galaxies and the cosmic web at $z\sim2.5$. We show that IGM tomography surveys with $\lesssim5$ $h^{-1}$ Mpc sightline spacing can accurately recover the eigenvectors of the tidal tensor, which we use to define the directions of the cosmic web. For galaxy spins and shapes, we use a model parametrized by the alignment strength, $\Delta\langle\cos\theta\rangle$, with respect to the tidal tensor eigenvectors from the underlying density field, and also consider observational effects such as errors in the galaxy position angle, inclination, and redshift. Measurements using the upcoming $\sim1\,\mathrm{deg}^2$ CLAMATO tomographic survey and 600 coeval zCOSMOS-Deep galaxies should place $3\sigma$ limits on extreme alignment models with $\Delta\langle\cos\theta\rangle\sim0.1$, but much larger surveys encompassing $>10,000$ galaxies, such as Subaru PFS, will be required to constrain models with $\Delta\langle\cos\theta\rangle\sim0.03$. These measurements will constrain models of galaxy-cosmic web alignment and test tidal torque theory at $z\sim2$, improving our understanding of the redshift dependence of galaxy-cosmic web alignment and the physics of intrinsic alignments.

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.

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.

The ever increasing size and complexity of data coming from simulations of cosmic structure formation demands equally sophisticated tools for their analysis. During the past decade, the art of object finding in these simulations has hence developed into an important discipline itself. A multitude of codes based upon a huge variety of methods and techniques have been spawned yet the question remained as to whether or not they will provide the same (physical) information about the structures of interest. Here we summarize and extent previous work of the "halo finder comparison project": we investigate in detail the (possible) origin of any deviations across finders. To this extent we decipher and discuss differences in halo finding methods, clearly separating them from the disparity in definitions of halo properties. We observe that different codes not only find different numbers of objects leading to a scatter of up to 20 per cent in the halo mass and Vmax function, but also that the particulars of those objects that are identified by all finders differ. The strength of the variation, however, depends on the property studied, e.g. the scatter in position, bulk velocity, mass, and the peak value of the rotation curve is practically below a few per cent, whereas derived quantities such as spin and shape show larger deviations. Our study indicates that the prime contribution to differences in halo properties across codes stems from the distinct particle collection methods and -- to a minor extent -- the particular aspects of how the procedure for removing unbound particles is implemented. We close with a discussion of the relevance and implications of the scatter across different codes for other fields such as semi-analytical galaxy formation models, gravitational lensing, and observables in general.

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.

[abridged] We present a detailed comparison of fundamental dark matter halo properties retrieved by a substantial number of different halo finders. These codes span a wide range of techniques including friends-of-friends (FOF), spherical-overdensity (SO) and phase-space based algorithms. We further introduce a robust (and publicly available) suite of test scenarios that allows halo finder developers to compare the performance of their codes against those presented here. This set includes mock haloes containing various levels and distributions of substructure at a range of resolutions as well as a cosmological simulation of the large-scale structure of the universe. All the halo finding codes tested could successfully recover the spatial location of our mock haloes. They further returned lists of particles (potentially) belonging to the object that led to coinciding values for the maximum of the circular velocity profile and the radius where it is reached. All the finders based in configuration space struggled to recover substructure that was located close to the centre of the host halo and the radial dependence of the mass recovered varies from finder to finder. Those finders based in phase space could resolve central substructure although they found difficulties in accurately recovering its properties. Via a resolution study we found that most of the finders could not reliably recover substructure containing fewer than 30-40 particles. However, also here the phase space finders excelled by resolving substructure down to 10-20 particles. By comparing the halo finders using a high resolution cosmological volume we found that they agree remarkably well on fundamental properties of astrophysical significance (e.g. mass, position, velocity, and peak of the rotation curve).

We present a set of ultra-large particle-mesh simulations of the LyA forest targeted at understanding the imprint of baryon acoustic oscillations (BAO) in the inter-galactic medium. We use 9 dark matter only simulations which can, for the first time, simultaneously resolve the Jeans scale of the intergalactic gas while covering the large volumes required to adequately sample the acoustic feature. Mock absorption spectra are generated using the fluctuating Gunn-Peterson approximation which have approximately correct flux probability density functions (PDFs) and small-scale power spectra. On larger scales there is clear evidence in the redshift space correlation function for an acoustic feature, which matches a linear theory template with constant bias. These spectra, which we make publicly available, can be used to test pipelines, plan future experiments and model various physical effects. As an illustration we discuss the basic properties of the acoustic signal in the forest, the scaling of errors with noise and source number density, modified statistics to treat mean flux evolution and misestimation, and non-gravitational sources such as fluctuations in the photo-ionizing background and temperature fluctuations due to HeII reionization.

We study the formation of dark matter halos in the concordance LCDM model over a wide range of redshifts, from z=20 to the present. Our primary focus is the halo mass function, a key probe of cosmology. By performing a large suite of nested-box N-body simulations with careful convergence and error controls (60 simulations with box sizes from 4 to 256 Mpc/h, we determine the mass function and its evolution with excellent statistical and systematic errors, reaching a few percent over most of the considered redshift and mass range. Across the studied redshifts, the halo mass is probed over 6 orders of magnitude (10^7 - 10^13.5 M_sun/h). Historically, there has been considerable variation in the high redshift mass function as obtained by different groups. We have made a concerted effort to identify and correct possible systematic errors in computing the mass function at high redshift and to explain the discrepancies between some of the previous results. We discuss convergence criteria for the required force resolution, simulation box size, halo mass range, initial and final redshift, and time stepping. Because of conservative cuts on the mass range probed by individual boxes, our results are relatively insensitive to simulation volume, the remaining sensitivity being consistent with extended Press-Schechter theory. Previously obtained mass function fits near z=0, when scaled by linear theory, are in good agreement with our results at all redshifts, although a mild redshift dependence consistent with that found by Reed and collaborators exists at low redshifts.

Current and upcoming cosmological observations allow us to probe structures on smaller and smaller scales, entering highly nonlinear regimes. In order to obtain theoretical predictions in these regimes, large cosmological simulations have to be carried out. The promised high accuracy from observations make the simulation task very demanding: the simulations have to be at least as accurate as the observations. This requirement can only be fulfilled by carrying out an extensive code validation program. The first step of such a program is the comparison of different cosmology codes including gravitation interactions only. In this paper we extend a recently carried out code comparison project to include five more simulation codes. We restrict our analysis to a small cosmological volume which allows us to investigate properties of halos. For the matter power spectrum and the mass function, the previous results hold, with the codes agreeing at the 10% level over wide dynamic ranges. We extend our analysis to the comparison of halo profiles and investigate the halo count as a function of local density. We introduce and discuss ParaView as a flexible analysis tool for cosmological simulations, the use of which immensely simplifies the code comparison task.

We study the evolution of the mass function of dark matter halos in the concordance LCDM model at high redshift. We employ overlapping (multiple-realization) numerical simulations to cover a wide range of halo masses, 10^7-10^{15}h^{-1}M_\odot, with redshift coverage beginning at z=20. The Press-Schechter mass function is significantly discrepant from the simulation results at high redshifts. Of the more recently proposed mass functions, our results are in best agreement with Warren et al. (2005). The statistics of the simulations -- along with good control over systematics -- allow for fits accurate to the level of 20% at all redshifts. We provide a concise discussion of various issues in defining and computing the halo mass function, and how these are addressed in our simulations.