We present a novel implementation of supermassive black hole (SMBH) formation, dynamics, and accretion in the massively parallel tree+SPH code, ChaNGa. This approach improves the modeling of SMBHs in fully cosmological simulations, allowing for a more detailed analysis of SMBH-galaxy co-evolution throughout cosmic time. Our scheme includes novel, physically motivated models for SMBH formation, dynamics and sinking timescales within galaxies, and SMBH accretion of rotationally supported gas. The sub-grid parameters that regulate star formation (SF) and feedback from SMBHs and SNe are optimized against a comprehensive set of z = 0 galaxy scaling relations using a novel, multi-dimensional parameter search. We have incorporated our new SMBH implementation and parameter optimization into a new set of high resolution, large-scale cosmological simulations called Romulus. We present initial results from our flagship simulation, Romulus25, showing that our SMBH model results in SF efficiency, SMBH masses, and global SF and SMBH accretion histories at high redshift that are consistent with observations. We discuss the importance of SMBH physics in shaping the evolution of massive galaxies and show how SMBH feedback is much more effective at regulating star formation compared to SNe feedback in this regime. Further, we show how each aspect of our SMBH model impacts this evolution compared to more common approaches. Finally, we present a science application of this scheme studying the properties and time evolution of an example dual AGN system, highlighting how our approach allows simulations to better study galaxy interactions and SMBH mergers in the context of galaxy-BH co-evolution.

Phylodynamics seeks to estimate effective population size fluctuations from molecular sequences of individuals sampled from a population of interest. One way to accomplish this task formulates an observed sequence data likelihood exploiting a coalescent model for the sampled individuals’ genealogy and then integrating over all possible genealogies via Monte Carlo or, less efficiently, by conditioning on one genealogy estimated from the sequence data. However, when analyzing sequences sampled serially through time, current methods implicitly assume either that sampling times are fixed deterministically by the data collection protocol or that their distribution does not depend on the size of the population. Through simulation, we first show that, when sampling times do probabilistically depend on effective population size, estimation methods may be systematically biased. To correct for this deficiency, we propose a new model that explicitly accounts for preferential sampling by modeling the sampling times as an inhomogeneous Poisson process dependent on effective population size. We demonstrate that in the presence of preferential sampling our new model not only reduces bias, but also improves estimation precision. Finally, we compare the performance of the currently used phylodynamic methods with our proposed model through clinically-relevant, seasonal human influenza examples.

We introduce , an package for phylodynamic analysis based on gene genealogies. The package main functionality is Bayesian nonparametric estimation of effective population size fluctuations over time. Our implementation includes several Markov chain Monte Carlo-based methods and an integrated nested Laplace approximation-based approach for phylodynamic inference that have been developed in recent years. Genealogical data describe the timed ancestral relationships of individuals sampled from a population of interest. Here, individuals are assumed to be sampled at the same point in time (isochronous sampling) or at different points in time (heterochronous sampling); in addition, sampling events can be modeled with preferential sampling, which means that the intensity of sampling events is allowed to depend on the effective population size trajectory. We assume the coalescent and the sequentially Markov coalescent processes as generative models of genealogies. We include several coalescent simulation functions that are useful for testing our phylodynamics methods via simulation studies. We compare the performance and outputs of various methods implemented in and outline their strengths and weaknesses. package is available at https://github.com/mdkarcher/phylodyn.

Supplementary data are available at Bioinformatics online.

The sources that reionized the universe are still unknown, but likely candidates are faint but numerous galaxies. In this paper we present results from running a high resolution, uniform volume simulation, the Vulcan , to predict the number densities of undetectable, faint galaxies and their escape fractions of ionizing radiation, f esc , during reionization. Our approach combines a high spatial resolution, a realistic treatment of feedback and hydro processes, a strict threshold for minimum number of resolution elements per galaxy, and a converged measurement of f esc . We calibrate our physical model using a novel approach to create realistic galaxies at z = 0, so the simulation is predictive at high redshifts. With this approach we can (1) robustly predict the evolution of the galaxy UV luminosity function at faint magnitudes down to M UV ∼ −15, two magnitudes fainter than observations, and (2) estimate f esc over a large range of galaxy masses based on the detailed stellar and gas distributions in resolved galaxies. We find steep faint end slopes, implying high number densities of faint galaxies, and the dependence of f esc on the UV magnitude of a galaxy, given by the powerlaw: log f esc = (0.51 ± 0.04)M UV + 7.3 ± 0.8, with the faint population having f esc ∼ 35%. Convolving the UV luminosity function with f esc (M UV ), we find an ionizing emissivity that is (1) dominated by the faintest galaxies and (2) reionizes the universe at the appropriate rate, consistent with observational constraints of the ionizing emissivity and the optical depth to the decoupling surface τ es , without the need for additional sources of ionizing radiation.