To improve our understanding of high-z galaxies we study the impact of H$_{2}$ chemistry on their evolution, morphology and observed properties. We compare two zoom-in high-resolution (30 pc) simulations of prototypical $M_{\star}\sim 10^{10} {\rm M}_{\odot}$ galaxies at $z=6$. The first, "Dahlia", adopts an equilibrium model for H$_{2}$ formation, while the second, "Alth{\ae}a", features an improved non-equilibrium chemistry network. The star formation rate (SFR) of the two galaxies is similar (within 50\%), and increases with time reaching values close to 100 ${\rm M}_{\odot}/\rm yr$ at $z=6$. They both have SFR-stellar mass relation consistent with observations, and a specific SFR of $\simeq 5\, {\rm Gyr}^{-1}$. The main differences arise in the gas properties. The non-equilibrium chemistry determines the H$\rightarrow$ H$_{2}$~transition to occur at densities $> 300\,{cm}^{-3}$, i.e. about 10 times larger than predicted by the equilibrium model used for Dahlia. As a result, Alth{\ae}a features a more clumpy and fragmented morphology, in turn making SN feedback more effective. Also, because of the lower density and weaker feedback, Dahlia sits $3\sigma$ away from the Schmidt-Kennicutt relation; Alth{\ae}a, instead nicely agrees with observations. The different gas properties result in widely different observables. Alth{\ae}a outshines Dahlia by a factor of 7 (15) in [CII]~$157.74\,\mu{\rm m}$ (H$_{2}$~$17.03\,\mu{\rm m}$) line emission. Yet, Alth{\ae}a is under-luminous with respect to the locally observed [CII]-SFR relation. Whether this relation does not apply at high-z or the line luminosity is reduced by CMB and metallicity effects remains as an open question.

Dust plays a key role in the evolution of the ISM and its correct modelling in numerical simulations is therefore fundamental. We present a new and self-consistent model that treats grain thermal coupling with the gas, radiation balance, and surface chemistry for molecular hydrogen. This method can be applied to any dust distribution with an arbitrary number of grain types without affecting the overall computational cost. In this paper we describe in detail the physics and the algorithm behind our approach, and in order to test the methodology, we present some examples of astrophysical interest, namely (i) a one-zone collapse with complete gas chemistry and thermochemical processes, (ii) a 3D model of a low-metallicity collapse of a minihalo starting from cosmological initial conditions, and (iii) a turbulent molecular cloud with H-C-O chemistry (277 reactions), together with self-consistent cooling and heating solved on the fly. Although these examples employ the publicly available code KROME, our approach can be easily integrated into any computational framework.

Understanding the formation of the extremely metal poor star SDSS-J102915+172927 is of fundamental importance to improve our knowledge on the transition between the first and second generation of stars in the Universe. In this paper, we perform three-dimensional cosmological hydrodynamical simulations of dust-enriched halos during the early stages of the collapse process including a detailed treatment of the dust physics. We employ the astrochemistry package \krome coupled with the hydrodynamical code \textsc{enzo} assuming grain size distributions produced by the explosion of core-collapse supernovae of 20 and 35 M$_\odot$ primordial stars which are suitable to reproduce the chemical pattern of the SDSS-J102915+172927 star. We find that the dust mass yield produced from Population III supernovae explosions is the most important factor which drives the thermal evolution and the dynamical properties of the halos. Hence, for the specific distributions relevant in this context, the composition, the dust optical properties, and the size-range have only minor effects on the results due to similar cooling functions. We also show that the critical dust mass to enable fragmentation provided by semi-analytical models should be revised, as we obtain values one order of magnitude larger. This determines the transition from disk fragmentation to a more filamentary fragmentation mode, and suggests that likely more than one single supernova event or efficient dust growth should be invoked to get such a high dust content.

Observations of high redshift quasars at $z>6$ indicate that they harbor supermassive black holes (SMBHs) of a billion solar masses. The direct collapse scenario has emerged as the most plausible way to assemble SMBHs. The nurseries for the direct collapse black holes are massive primordial halos illuminated with an intense UV flux emitted by population II (Pop II) stars. In this study, we compute the critical value of such a flux ($J_{21}^{\rm crit}$) for realistic spectra of Pop II stars through three-dimensional cosmological simulations. We derive the dependence of $J_{21}^{\rm crit}$ on the radiation spectra, on variations from halo to halo, and on the impact of X-ray ionization. Our findings show that the value of $J_{21}^{\rm crit}$ is a few times $\rm 10^4$ and only weakly depends on the adopted radiation spectra in the range between $T_{\rm rad}=2 \times 10^4-10^5$ K. For three simulated halos of a few times $\rm 10^{7}$~M$_{\odot}$, $J_{21}^{\rm crit}$ varies from $\rm 2 \times 10^4 - 5 \times 10^4$. The impact of X-ray ionization is almost negligible and within the expected scatter of $J_{21}^{\rm crit}$ for background fluxes of $J_{\rm X,21} \leq 0.1$. The computed estimates of $J_{21}^{\rm crit}$ have profound implications for the quasar abundance at $z=10$ as it lowers the number density of black holes forming through an isothermal direct collapse by a few orders of magnitude below the observed black holes density. However, the sites with moderate amounts of $\rm H_2$ cooling may still form massive objects sufficient to be compatible with observations.

Context. The seeds of the first supermassive black holes may have resulted from the direct collapse of hot primordial gas in $\gtrsim 10^4$ K haloes, forming a supermassive or quasistar as an intermediate stage. Aims. We explore the formation of a protostar resulting from the collapse of primordial gas in the presence of a strong Lyman-Werner radiation background. Particularly, we investigate the impact of turbulence and rotation on the fragmentation behaviour of the gas cloud. We accomplish this goal by varying the initial turbulent and rotational velocities. Methods. We performed 3D adaptive mesh refinement simulations with a resolution of 64 cells per Jeans length using the ENZO code, simulating the formation of a protostar up to unprecedentedly high central densities of $10^{21}$ cm$^{-3}$, and spatial scales of a few solar radii. To achieve this goal, we employed the KROME package to improve modelling of the chemical and thermal processes. Results. We find that the physical properties of the simulated gas clouds become similar on small scales, irrespective of the initial amount of turbulence and rotation. After the highest level of refinement was reached, the simulations have been evolved for an additional ~5 freefall times. A single bound clump with a radius of $2 \times 10^{-2}$ AU and a mass of ~$7 \times 10^{-2}$ M$_{\odot}$ is formed at the end of each simulation, marking the onset of protostar formation. No strong fragmentation is observed by the end of the simulations, regardless of the initial amount of turbulence or rotation, and high accretion rates of a few solar masses per year are found. Conclusions. Given such high accretion rates, a quasistar of $10^5$ M$_{\odot}$ is expected to form within $10^5$ years.

The ability of metal free gas to cool by molecular hydrogen in primordial halos is strongly associated with the strength of ultraviolet (UV) flux produced by the stellar populations in the first galaxies. Depending on the stellar spectrum, these UV photons can either dissociate $\rm H_{2}$ molecules directly or indirectly by photo-detachment of $\rm H^{-}$ as the latter provides the main pathway for $\rm H_{2}$ formation in the early universe. In this study, we aim to determine the critical strength of the UV flux above which the formation of molecular hydrogen remains suppressed for a sample of five distinct halos at $z>10$ by employing a higher order chemical solver and a Jeans resolution of 32 cells. We presume that such flux is emitted by PopII stars implying atmospheric temperatures of $\rm 10^{4}$~K. We performed three-dimensional cosmological simulations and varied the strength of the UV flux below the Lyman limit in units of $\rm J_{21}$. Our findings show that the value of $\rm J_{21}^{crit}$ varies from halo to halo and is sensitive to the local thermal conditions of the gas. For the simulated halos it varies from 400-700 with the exception of one halo where $\rm J_{21}^{crit} \geq 1500$. This has important implications for the formation of direct collapse black holes and their estimated population at z > 6. It reduces the number density of direct collapse black holes by almost three orders of magnitude compared to the previous estimates.

Recent discoveries of carbon-enhanced metal-poor stars like SMSS J031300.36-670839.3 provide increasing observational insights into the formation conditions of the first second-generation stars in the Universe, reflecting the chemical conditions after the first supernova explosion. Here, we present the first cosmological simulations with a detailed chemical network including primordial species as well as C, C$^+$, O, O$^+$, Si, Si$^+$, and Si$^{2+}$ following the formation of carbon-enhanced metal poor stars. The presence of background UV flux delays the collapse from $z=21$ to $z=15$ and cool the gas down to the CMB temperature for a metallicity of Z/Z$_\odot$=10$^{-3}$. This can potentially lead to the formation of lower mass stars. Overall, we find that the metals have a stronger effect on the collapse than the radiation, yielding a comparable thermal structure for large variations in the radiative background. We further find that radiative backgrounds are not able to delay the collapse for Z/Z$_\odot$=10$^{-2}$ or a carbon abundance as in SMSS J031300.36-670839.3.

Radiative feedback from populations II stars played a vital role in early structure formation. Particularly, photons below the Lyman limit can escape the star forming regions and produce a background ultraviolet (UV) flux which consequently may influence the pristine halos far away from the radiation sources. These photons can quench the formation of molecular hydrogen by photo-detachment of $\rm H^{-}$. In this study, we explore the impact of such UV radiation on fragmentation in massive primordial halos of a few times $\rm 10^{7}$~M${_\odot}$. To accomplish this goal, we perform high resolution cosmological simulations for two distinct halos and vary the strength of the impinging background UV field in units of $\rm J_{21}$. We further make use of sink particles to follow the evolution for 10,000 years after reaching the maximum refinement level. No vigorous fragmentation is observed in UV illuminated halos while the accretion rate changes according to the thermal properties. Our findings show that a few 100-10, 000 solar mass protostars are formed when halos are irradiated by $\rm J_{21}=10-500$ at $\rm z>10$ and suggest a strong relation between the strength of UV flux and mass of a protostar. This mode of star formation is quite different from minihalos, as higher accretion rates of about $\rm 0.01-0.1$ M$_{\odot}$/yr are observed by the end of our simulations. The resulting massive stars are the potential cradles for the formation of intermediate mass black holes at earlier cosmic times and contribute to the formation of a global X-ray background.

While Population III stars are typically thought to be massive, pathways towards lower-mass Pop III stars may exist when the cooling of the gas is particularly enhanced. A possible route is enhanced HD cooling during the merging of dark-matter halos. The mergers can lead to a high ionization degree catalysing the formation of HD molecules and may cool the gas down to the cosmic microwave background (CMB) temperature. In this paper, we investigate the merging of mini-halos with masses of a few 10$^5$ M$_\odot$ and explore the feasibility of this scenario. We have performed three-dimensional cosmological hydrodynamics calculations with the ENZO code, solving the thermal and chemical evolution of the gas by employing the astrochemistry package KROME. Our results show that the HD abundance is increased by two orders of magnitude compared to the no-merging case and the halo cools down to $\sim$60 K triggering fragmentation. Based on Jeans estimates the expected stellar masses are about 10 M$_\odot$. Our findings show that the merging scenario is a potential pathway for the formation of low-mass stars.

Chemistry plays a key role in many astrophysical situations regulating the cooling and the thermal properties of the gas, which are relevant during gravitational collapse, the evolution of disks and the fragmentation process. In order to simplify the usage of chemical networks in large numerical simulations, we present the chemistry package KROME, consisting of a Python pre-processor which generates a subroutine for the solution of chemical networks which can be embedded in any numerical code. For the solution of the rate equations, we make use of the high-order solver DLSODES, which was shown to be both accurate and efficient for sparse networks, which are typical in astrophysical applications. KROME also provides a large set of physical processes connected to chemistry, including photochemistry, cooling, heating, dust treatment, and reverse kinetics. The package presented here already contains a network for primordial chemistry, a small metal network appropriate for the modelling of low metallicities environments, a detailed network for the modelling of molecular clouds, a network for planetary atmospheres, as well as a framework for the modelling of the dust grain population. In this paper, we present an extended test suite ranging from one-zone and 1D-models to first applications including cosmological simulations with ENZO and RAMSES and 3D collapse simulations with the FLASH code. The package presented here is publicly available at http://kromepackage.org/ and https://bitbucket.org/krome/krome_stable

A detailed quantum analysis of a ionic reaction with a crucial role in the ISM is carried out to generate ab initio reactive cross sections with a quantum method. From them we obtain the corresponding CH+ depletion rates over a broad range of temperatures. The new rates are further linked to a complex chemical network that shows the evolution in time of the CH+ abundance in photodissociation region (PDR) and molecular cloud (MC) environments. The evolutionary abundances of CH+ are given by numerical solutions of a large set of coupled, first-order kinetics equations by employing the new chemical package KROME. The differences found between all existing calculations from low-T experiments are explained via a simple numerical model that links the low-T cross section reductions to collinear approaches where nonadiabatic crossings dominate. The analysis of evolutionary abundance of CH+ reveals that the important region for the depletion reaction of this study is that above 100 K, hence showing that, at least for this reaction, the differences with the existing low-temperature experiments are of essentially no importance within the astrochemical environments. A detailed analysis of the chemical network involving CH+ also shows that a slight decrease in the initial oxygen abundance might lead to higher CH+ abundance since the main chemical carbon ion depletion channel is reduced in efficiency. This simplified observation might provide an alternative starting point to understand the problem of astrochemical models in matching the observed CH+ abundances.

Population III stars are the first stars in the Universe to form at z=20-30 out of a pure hydrogen and helium gas in minihalos of 10^5-10^6 M$_\odot$ . Cooling and fragmentation is thus regulated via molecular hydrogen. At densities above 10^8 cm$^{-3}$, the three-body H2 formation rates are particularly important for making the gas fully molecular. These rates were considered to be uncertain by at least a few orders of magnitude. We explore the impact of new accurate three-body H2 formation rates derived by Forrey (2013) for three different minihalos, and compare to the results obtained with three-body rates employed in previous studies. The calculations are performed with the cosmological hydrodynamics code ENZO (release 2.2) coupled with the chemistry package KROME (including a network for primordial chemistry), which was previously shown to be accurate in high resolution simulations. While the new rates can shift the point where the gas becomes fully molecular, leading to a different thermal evolution, there is no trivial trend in how this occurs. While one might naively expect the results to be inbetween the calculations based on Palla et al. (1983) and Abel et al. (2002), the behavior can be close to the former or the latter depending on the dark matter halo that is explored. We conclude that employing the correct three-body rates is about as equally important as the use of appropriate initial conditions, and that the resulting thermal evolution needs to be calculated for every halo individually.

The formation of the first stars in the Universe is regulated by a sensitive interplay of chemistry and cooling with the dynamics of a self-gravitating system. As the outcome of the collapse and the final stellar masses depend sensitively on the thermal evolution, it is necessary to accurately model the thermal evolution in high resolution simulations. As previous investigations raised doubts regarding the convergence of the temperature at high resolution, we investigate the role of the numerical method employed to model the chemistry and the thermodynamics. Here we compare the standard implementation in the adaptive-mesh refinement code \verb|ENZO|, employing a first order backward differentiation formula (BDF), with the 5th order accurate BDF solver \verb|DLSODES|. While the standard implementation in \verb|ENZO| shows a strong dependence on the employed resolution, the results obtained with \verb|DLSODES| are considerably more robust, both with respect to the chemistry and thermodynamics, but also for dynamical quantities such as density, total energy or the accretion rate. We conclude that an accurate modeling of the chemistry and thermodynamics is central for primordial star formation.

Chemistry has a key role in the evolution of the interstellar medium (ISM), so it is highly desirable to follow its evolution in numerical simulations. However, it may easily dominate the computational cost when applied to large systems. In this paper we discuss two approaches to reduce these costs: (i) based on computational strategies, and (ii) based on the properties and on the topology of the chemical network. The first methods are more robust, while the second are meant to be giving important information on the structure of large, complex networks. To this aim we first discuss the numerical solvers for integrating the system of ordinary differential equations (ODE) associated with the chemical network. We then propose a buffer method that decreases the computational time spent in solving the ODE system. We further discuss a flux-based method that allows one to determine and then cut on the fly the less active reactions. In addition we also present a topological approach for selecting the most probable species that will be active during the chemical evolution, thus gaining information on the chemical network that otherwise would be difficult to retrieve. This topological technique can also be used as an a priori reduction method for any size network. We implemented these methods into a 1D Lagrangian hydrodynamical code to test their effects: both classes lead to large computational speed-ups, ranging from x2 to x5. We have also tested some hybrid approaches finding that coupling the flux method with a buffer strategy gives the best trade-off between robustness and speed-up of calculations.

We present a new computational scheme aimed at reducing the complexity of the chemical networks in astrophysical models, one which is shown to markedly improve their computational efficiency. It contains a flux-reduction scheme that permits to deal with both large and small systems. This procedure is shown to yield a large speed-up of the corresponding numerical codes and provides good accord with the full network results. We analyse and discuss two examples involving chemistry networks of the interstellar medium and show that the results from the present reduction technique reproduce very well the results from fuller calculations.

We present a quantum-mechanical study of the exothermic 7LiH reaction with H. Accurate reactive probabilities and rate coefficients are obtained by solving the Schrodinger equation for the motion of the three nuclei on a single Born-Oppenheimer potential energy surface (PES) and using a coupled-channel hyperspherical coordinate method. Our new rates indeed confirm earlier, qualitative predictions and some previous theoretical calculations, as discussed in the main text. In the astrophysical domain we find that the depletion process largely dominates for redshift (z) between 400 and 100, a range significant for early Universe models. This new result from first-principle calculations leads us to definitively surmise that LiH should be already destroyed when the survival processes become important. Because of this very rapid depletion reaction, the fractional abundance of LiH is found to be drastically reduced, so that it should be very difficult to manage to observe it as an imprinted species in the cosmic background radiation (CBR). The present findings appear to settle the question of LiH observability in the early Universe. We further report several state-to-state computed reaction rates in the same range of temperatures of interest for the present problem.

Calculations have been carried out for the vibrational quenching of excited H$_2$ molecules which collide with Li$^+$ ions at ultralow energies. The dynamics has been treated exactly using the well known quantum coupled-channel expansions over different initial vibrational levels. The overall interaction potential has been obtained from the calculations carried out earlier in our group using highly correlated ab initio methods. The results indicate that specific features of the scattering observables, e.g. the appearance of Ramsauer-Townsend minima in elastic channel cross sections and the marked increase of the cooling rates from specific initial states, can be linked to potential properties at vanishing energies (sign and size of scattering lengths) and to the presence of either virtual states or bound states. The suggestion is made that by selecting the initial state preparation of the molecular partners, the ionic interactions would be amenable to controlling quenching efficiency at ultralow energies.