The formation of the firsts stars some 100-300 Myr after the Big Bang marked the end of the cosmic darks ages and created the elemental building blocks of not only rocky planets but eventually us. Understanding their formation, lifetimes, and contributions to the evolution of our universe is one of the current frontiers in astronomy and astrophysics.

In this work I present an improved model for following the formation of Pop III stars, their effects on early galaxy evolution, and how we might search for them. I make use of a new subgrid model of turbulent mixing to accurately follow the time scales required to mix supernova (SN) ejecta -- enriched with heavy elements -- into the pristine gas. I implement this model within a large-scale cosmological simulation and follow the fraction of gas with metallicity below a critical value marking the boundary between Pop III and metal enriched Population II (Pop II) star formation. I demonstrate that accounting for subgrid mixing results in a Pop III stars formation rate that is 2-3 times higher than standard models with the same physical resolution.

I also implement and track a new "Primordial metals" (PM) scalar that tracks the metals generated by Pop III SNe. These metals are taken up by second generation stars and likely result in a subclass of carbon-enhanced, metal-poor (CEMP) stars. By tracking both regular metals and PM, I can model, in post-processing, the elemental abundances of simulation stars. I find good agreement between observations of CEMP-no Milky Way halo stars and second generation stars within the simulation when assuming the first stars had a typical mass of 60 M☉, providing clues as to the Pop III initial mass function.