Following a cosmological gaseous collapse over many orders of magnitude in length makes such a simulation technically difficult, but with improvements in algorithms and physical models, several groups have been able to make substantial progress. In the late 1990's and early 2000's, two independent groups used three-dimensional simulations to study Population III star formation. Using smoothed particle hydrodynamics (SPH) simulations of isolated and virialized dark matter halos with masses 2 × 106 M at z = 30, one group found that the object cooled to 300 K through H2 formation, using the chemical network of , and fragmented into a filamentary structure with a Jeans mass of 103 M [13, 14]. In the second paper, they followed the collapse to higher gas densities of 108 cm-3 and studied the continued formation of dense clumps with the same 103 M characteristic mass. The other group used cosmological adaptive mesh refinement (AMR) simulations to focus on the formation of a molecular cloud hosted by a 7 × 105 M DM halo [1, 3, 4]. In each successive paper, the H2 chemistry model [2, 8] was improved to include more processes, such as the three-body H2 formation process, to follow the central collapse to gas densities up to 3 × 1013 cm-3. Both groups came to the conclusion that further fragmentation was suppressed because of a lack of cooling below 300 K and that Population III stars were very massive in the range 30-300 M.
These simulations only represented a limited sample of collapses and could not provide much insight to the Population III initial mass function (IMF). Twelve additional AMR simulations were conducted to look at any variations in primordial gas collapses . They found that the collapses occurred in DM halos with masses in the range 1.5-7 × 105 M with the scatter caused by differing halo formation histories. The mass accretion rate onto the central molecular cloud was higher at 10-4 M yr-1 at z ~ 30, and it can increase by two orders of magnitude at z ~ 20 in some halos, agreeing with .
Following the collapse to densities higher than 1013 cm-3 required the inclusion of collisionally induced emission, chemical heating from H2 formation, and gas opacity above 1018 cm-3.  found with SPH simulations that the initial collapsing region did not fragment as it condensed to protostellar densities n ~ 1021 cm-3, forming a protostellar shock in the process. The inner 10 M had an accretion rate varying between 0.01 and 0.1 M yr-1, possibly growing to 10 M within 1000 yr.
In the past few years, multiple groups have been focusing on the subsequent growth of these protostars over several dynamical times, improving upon the earlier works that stopped at the first collapse. This has proved to be challenging because of the ever-decreasing Courant factors at higher densities. One workaround is the creation of "sink particles" that accrete nearby gravitationally-bound gas, allowing the simulation to progress past the first collapse; however, one loses all hydrodynamical information above some density threshold. In one out of five realizations in AMR calculations without sink particles,  found that the collapsing core fragmented at a density of 1011 cm-3 into two clumps that are separated by 800 AU with 100 M of gas within a sphere with radius twice their separation. At the same time,  also found that disk instabilities causes fragmentation into a binary system with a 40 M and 10 M. This was later confirmed by simulations of a collapse of an isolated Jeans-unstable primordial gas cloud that fragmented into many multiple systems with some very tightly bound to separations less than an AU . Utilizing a new moving mesh code,  studied the collapse in five different primordial DM halos, and they evolved them for 1000 yr after the first protostar forms. By evolving these protostars further, they included the effects of protostellar radiative feedback in the infrared in the optically-thin limit. In all cases, the molecular cloud fragments into ~ 10 sink particles, some of which later merge to form more massive protostars. The mass function from these simulations is relatively flat, i.e. a top-heavy IMF.
After the protostar has reached ~ 10 M, radiative feedback from ionizing radiation will begin to suppress further accretion. Only recently has this been incorporated into numerical simulations of Population III star formation. Starting from initial conditions extracted from a cosmological simulation , two-dimensional axi-symmetric radiation hydrodynamical simulations showed that an accretion disk forms around a new protostar with the ionizing radiation preferentially escaping through the polar regions . The disk itself is slowly photo-evaporated, halting accretion after 70,000 yr. At this point, the final mass of the Population III star is 43 M. Without any radiative feedback, the protostar would have continued to grow to ~ 100 M. In a cosmological setting,  found a binary system still forms in the presence of radiative feedback. Without feedback, the primary star grows to 28 M over 5,000 yr. With feedback, the primary and secondary stars only grow to 19 and 10 M, respectively. An extrapolation of the mass accretion history shows that both stellar masses will asymptote to 30 M, creating an equal-mass binary. Once the stars have entered the main sequence, they will start to ionize and heat their cosmic neighborhood, which I will review next.