We simulate the growth of a Population III stellar system, starting from cosmological initial conditions at z = 100. We follow the formation of a minihalo and the subsequent collapse of its central gas to high densities, resolving scales as small as ∼ 1 AU. Using sink particles to represent the growing protostars, we model the growth of the photodissociating and ionizing region around the first sink, continuing the simulation for ∼ 5000 yr after initial protostar formation. Along with the first-forming sink, several tens of secondary sinks form before an ionization front develops around the most massive star. The resulting cluster has high rates of sink formation, ejections from the stellar disc, and sink mergers during the first ∼ 2000 yr, before the onset of radiative feedback. By this time a warm ∼ 5000 K phase of neutral gas has expanded to roughly the disc radius of 2000 AU, slowing mass flow onto the disc and sinks. By 5000 yr the most massive star grows to 20 M ⊙ , while the total stellar mass approaches 75 M ⊙ . Out of the ∼ 40 sinks, approximately 30 are low-mass (M * < 1 M ⊙ ), and if the simulation had resolved smaller scales an even greater number of sinks might have formed. Thus, protostellar radiative feedback is insufficient to prevent rapid disc fragmentation and the formation of a high-member Pop III cluster before an ionization front emerges. Throughout the simulation, the majority of stellar mass is contained within the most massive stars, further implying that the Pop III initial mass function is top-heavy.

Stars rarely form in isolation. Nearly half of the stars in the Milky Way have a companion, and this fraction increases in starforming regions. However, why some dense cores and filaments form bound pairs while others form single stars remains unclear. We present a set of three-dimensional, gravo-magnetohydrodynamic simulations of turbulent star-forming clouds, aimed at understanding the formation and evolution of multiple-star systems formed through large scale ( 10 3 AU) turbulent fragmentation. We investigate three global magnetic field strengths, with global mass-to-flux ratios of µ φ = 2, 8, and 32. The initial separations of protostars in multiples depends on the global magnetic field strength, with stronger magnetic fields (e.g., µ φ = 2) suppressing fragmentation on smaller scales. The overall multiplicity fraction (MF) is between 0.4 − 0.6 for our strong and intermediate magnetic field strengths, which is in agreement with observations. The weak field case has a lower fraction. The MF is relatively constant throughout the simulations, even though stellar densities increase as collapse continues. While the MF rarely exceeds 60% in all three simulations, over 80% of all protostars are part of a binary system at some point. We additionally find that the distribution of binary spin mis-alignment angles is consistent with a randomized distribution. In all three simulations, several binaries originate with wide separations and dynamically evolve to 10 2 AU separations. We show that a simple model of mass accretion and dynamical friction with the gas can explain this orbital evolution.

We present a highly-parallel multi-frequency hybrid radiation hydrodynamics algorithm that combines a spatially-adaptive long characteristics method for the radiation field from point sources with a moment method that handles the diffuse radiation field produced by a volumefilling fluid. Our Hybrid Adaptive Ray-Moment Method (HARM 2 ) operates on patch-based adaptive grids, is compatible with asynchronous time stepping, and works with any moment method. In comparison to previous long characteristics methods, we have greatly improved the parallel performance of the adaptive long-characteristics method by developing a new completely asynchronous and non-blocking communication algorithm. As a result of this improvement, our implementation achieves near-perfect scaling up to O(10 3 ) processors on distributed memory machines. We present a series of tests to demonstrate the accuracy and performance of the method.

Gravitational instability (GI) of a dust-rich layer at the midplane of a gaseous circumstellar disk is one proposed mechanism to form planetesimals, the building blocks of rocky planets and gas giant cores. Self-gravity competes against the Kelvin-Helmholtz instability (KHI): gradients in dust content drive a vertical shear which risks overturning the dusty subdisk and forestalling GI. To understand the conditions under which the disk can resist the KHI, we perform three-dimensional simulations of stratified subdisks in the limit that dust particles are small and aerodynamically well coupled to gas, thereby screening out the streaming instability and isolating the KHI. Each subdisk is assumed to have a vertical density profile given by a spatially constant Richardson number Ri. We vary Ri and the midplane dust-to-gas ratio µ 0 and find that the critical Richardson number dividing KHunstable from KH-stable flows is not unique; rather Ri crit grows nearly linearly with µ 0 for µ 0 = 0.3-10. Plausibly a linear dependence arises for µ 0 ≪ 1 because in this regime the radial Kepler shear replaces vertical buoyancy as the dominant stabilizing influence. Why this dependence should persist at µ 0 > 1 is a new puzzle. The bulk (height-integrated) metallicity is uniquely determined by Ri and µ 0 . Only for disks of bulk solar metallicity is Ri crit ≈ 0.2, close to the classical value. Our empirical stability boundary is such that a dusty sublayer can gravitationally fragment and presumably spawn planetesimals if embedded within a solar metallicity gas disk ∼4× more massive than the minimum-mass solar nebula; or a minimum-mass disk having ∼3× solar metallicity; or some intermediate combination of these two possibilities. Gravitational instability seems possible without resorting to the streaming instability or to turbulent concentration of particles.

In regions of star formation, protostars and newborn stars will accrete mass from their natal clouds. These clouds are threaded by magnetic fields with a strength characterized by the plasma β-the ratio of thermal and magnetic pressures. Observations show that molecular clouds have β < ∼ 1, so magnetic fields have the potential to play a significant role in the accretion process. We have carried out a numerical study of the effect of large-scale magnetic fields on the rate of accretion onto a uniformly moving point particle from a uniform, non-self-gravitating, isothermal gas. We consider gas moving with sonic Mach numbers of up M ≈ 45, magnetic fields that are either parallel, perpendicular, or oriented 45 • to the flow, and β as low as 0.01. Our simulations utilize adaptive mesh refinement in order to obtain high spatial resolution where it is needed; this also allows the boundaries to be far from the accreting object to avoid unphysical effects arising from boundary conditions. Additionally, we show our results are independent of our exact prescription for accreting mass in the sink particle. We give simple expressions for the steady-state accretion rate as a function of β and M for the parallel and perpendicular orientations. Using typical molecular cloud values of M ∼ 5 and β ∼ 0.04 from the literature, our fits suggest a 0.4 M ⊙ star accretes ∼ 4 × 10 −9 M ⊙ /year, almost a factor of two less than accretion rates predicted by hydrodynamic models. This disparity can grow to orders of magnitude for stronger fields and lower Mach numbers. We also discuss the applicability of these accretion rates versus accretion rates expected from gravitational collapse, and under what conditions a steady state is possible. The reduction in the accretion rate in a magnetized medium leads to an increase in the time required to form stars in competitive accretion models, making such models less efficient than predicted by Bondi-Hoyle rates. Our results should find application in numerical codes, enabling accurate subgrid models of sink particles accreting from magnetized media.