Physical conditions in spherical prestellar cores with analytical equations

1 Introduction

A major issue for modelling complex processes in star-forming cloud cores is the evolution of the core. General star-formation models can afford focusing on different aspects of macroscopic physics (e.g. Grudić et al. 2021, Lebreuilly et al. 2021, Bate 2022). On the other hand, models addressing microscopic-level processes in great detail, such as chemistry or radiative transfer, often require star-forming cloud core evolution solutions that are simple, general, and computationally cheap, yet sufficiently close to reality. In the age of high-performance computing, such solutions are still useful when new phenomena are investigated, and the results must be understood clearly, without interference of complex background physics.

In the case of astrochemistry, the simple standard is a pseudo-time-dependent model that considers a stable, dark dense cloud core with constant physical conditions. The numerical density of H atoms usually is assumed to be in the interval n H = 1 0 4 … 1 0 5 cm − 3 , temperature T = 10  K, and interstellar extinction A V = 10  mag (see Hasegawa et al. 1992 for a classic example). Such a simple model does not account for cloud macroscopic evolution in any way. However, of great interest are scenarios where the cloud core collapses to form a protostar because radio-astronomical observations of molecules help tracing processes at different stages in star formation (e.g. Cazaux et al. 2003, Sakai et al. 2008, Bacmann et al. 2012, Lefloch et al. 2018, Harju et al. 2020).

Understanding the molecular-level processes is essential for explaining observations of collapsing cores. However, such prestellar cores cannot be represented by a pseudo-time-dependent model, which has been so useful in qualitative investigations of different proposed molecular processes (see Roberts et al. 2007, Acharyya et al. 2011, Reboussin et al. 2014, Hincelin et al. 2015, Rawlings and Williams 2021, Kalvāns and Silsbee 2022, to name a few) because of two reasons. First, the variation of physical conditions during the evolution of the core is too great and T , n H , and A V cannot be approximated to single values during the simulation. Second, the chemical history from previous stages with their specific conditions is partially preserved for further stages in the form of icy mantles accumulated on grain surfaces. The composition of such interstellar ices is less affected by changes in temperature or radiation, and, when the composition does change, it occurs within the ice elemental build-up established at the stage when the ice layer accreted. This is true until the mantle starts to evaporate. Evaporation is a major chemical event (Viti and Williams 1999), during which complex molecules synthesized in ices and their daughter species temporally become observable in the gas (Blake et al. 1987). The type and abundance of evaporated molecules depend directly on the freeze-out of ices during the cold prestellar stage.

The aforementioned means that chemistry in star-forming clouds has to be simulated with time-dependent physical models. Comprehensive 1D and 2D models are able to trace chemistry in different places and times of a star-forming core (Wakelam et al. 2011, Tassis et al. 2012, Furuya et al. 2017, Kalvāns 2018, Aikawa et al. 2020, Flores-Rivera et al. 2021, Jin et al. 2022, Murillo et al. 2022). Some 3D models have also been employed in astrochemistry (Hocuk and Cazaux 2015, Hincelin et al. 2016, Quénard et al. 2018). However, new chemical processes and molecules often have to be first investigated with the help of simpler 0D models that trace the evolution of a single gas parcel in a spherically symmetric star-forming core undergoing free-fall collapse followed by heat-up of the remaining envelope by the protostar (Garrod et al. 2022, Taniguchi et al. 2019, Zhao et al. 2021, Entekhabi et al. 2022). The simplified approach is rather straightforward to implement (e.g. Nejad et al. 1990, Taquet et al. 2014). The chemistry is typically calculated only for a single gas parcel at the centre of the star-forming core, where most of its mass resides. This approach offers clearly interpretable results but presents a problem because the centre of the core is where the star forms. Moreover, astrochemical observations of prestellar cores, protostars, and very low luminosity objects (VeLLOs, early star-forming cores, some possibly containing the first hydrostatic cores (Young et al. 2004, Belloche et al. 2006) trace sight-lines along the central objects as well as off-centre regions (e.g. Tobin and Megeath 2019, Di Francesco et al. 2020, Imai et al. 2022, Redaelli et al. 2022, van Gelder et al. 2022). Therefore, it is important to have physical conditions for 0D models at locations near and far of the protostar, and even for prestellar gas parcels that fall into the protostar and the initial and eventual final locations of the parcel have to be carefully chosen.

For addressing and easing the above problems with 0D models, and expanding their scope of applicability, we aim to create a series of precalculated datasets tracing the conditions in gas parcels at different initial radii within prestellar cores with a range of masses. The goal is to provide astrochemists with simple-to-use physical data for prestellar cores. These data would be more reliable and of higher quality than the relatively simple free-fall collapse formulas used in the references above and would be applicable for specific cases, such as calculating chemistry in gas parcels that end up in the protostar or parcels that remain at the margin of the envelope. Therefore, the novelty of this study will be equations that provide an accessible and relatively accurate framework for studies of microscopic processes in prestellar cores.

The aforementioned aim – providing a physics package of prestellar cores for 0D astrochemical calculations – presents several tasks. The data have to represent a range of core masses and a range of histories for infalling gas parcels, including parcels that end up in the protostar and in the remaining circumstellar envelope. Nevertheless, the data must be sufficiently compact, with limited choices to be attractive and not to complicate the concise nature of investigative 0D models. Finally, the data must be presented in a format that can be easily implemented in astrochemical codes.

2 Methods

The diffuse interstellar medium, molecular clouds, and dense star forming cloud cores differ by their density, temperature, motions, and evolutionary time scales. Denser than average regions of the ISM are considered clouds or nebulae. Molecular clouds are dense nebulas, shielded from the interstellar radiation, permitting the formation of molecules. Only the densest molecular clouds foster star formation. They have typical sizes between 1 and 200 pc and temperatures between 10 and 20 K. Smaller clumps are nested within molecular clouds, with size < 10  pc and mass of 10...100   M ⊙ . Within these clumps, low-mass cores can be found with masses below 10   M ⊙ and sizes below 0.1 pc (Enoch et al. 2008, Ballesteros-Paredes et al. 2020). This data-generating study is concerned with dense cores starting runaway gravitational collapse, up to the point of initial formation of the protostar, that is, the prestellar core stage.

For the creation of dataset for prestellar cores, we employ a 1D code for the dynamical evolution of the thermal structure of a spherically symmetric star-forming core (Pavlyuchenkov et al. 2015). It is a two-component model for computing the thermal structure of protostellar clouds. A key feature of this model is the separate descriptions of dust and gas temperatures, which can differ during the earliest stages of clouds evolution. For radiative-transfer calculations, the overall frequency range has been split into low (IR) and high (UV) intervals, representing the dust emission and interstellar radiation. Mean UV intensity has been approximated via direct integration of the radiative-transfer equation along specified directions. The adopted temperature of background UV radiation is 1 0 4 K with a dilution of background UV radiation of 1 0 − 14 . IR intensity is based on solving a system of momentum equations realizing diffusion approximation. The model includes four interacting components - gas, dust, IR, and UV radiation. While star-forming cores come in a variety of spheroid shapes (Myers et al. 1991, Jones et al. 2001, Anathpindika and Di Francesco 2022), a sphere allows the generality necessary for 0D models aiming to reproduce observed chemical abundances in a multitude of sources. Significant deviation from spheroid geometry occurs only in the protostellar stage.

The Pavlyuchenkov et al. (2015) model has been derived from the 2D code of Zhilkin et al. (2009), which stems from the works of Dudorov et al. (2003), Pavlyuchenkov and Shustov (2004), Semenov et al. (2005), Pavlyuchenkov et al. (2006), and Dudorov and Zhilkin (2008). Pavlyuchenkov and Zhilkin (2013) converted the Zhilkin et al. (2009) code to 1D, while adding a more complex and precise description of radiation transfer in a prestellar cloud core. In the transition from 2D to 1D, cloud core rotation and magnetic fields were excluded from the model. This bodes well with our aim to provide standardised data for simple 0D chemical simulations, where magnetic fields and rotation are ignored to reduce the number of variables and allow easy comparison between different models. Pavlyuchenkov et al. (2015) improved the model by adding separate gas and dust temperatures, among other changes. Model results have a sufficiently high accuracy for comparison with observational data (Pavlyuchenkov and Zhilkin 2013).

Initially, the model calculates the structure of a core in thermal and hydrostatic equilibrium. After this point, we considered the collapse and evolution of the dense core, up to the formation of the first hydrostatic core.

Technically, the model is using a one-dimensional, constant sink-cell approach. Each spherical shell corresponds to 1/40,000 of total prestellar core mass. The sphere shell count choice (in this case, 40,000) depends on the required data accuracy, preferred simulation length and total data size. The core initial conditions are described by radius values, which correspond to each shell, and a constant mass value that is the same for all shells. In essence these are mass coordinates specifically formatted for the model. First, this structure is used to model a cloud at a hydrostatic equilibrium, based on joint solutions of the equations for the thermal structure and hydrostatic equilibrium. The hydrostatic equilibrium structure of the cloud is computed by iteratively calculating radiation transport (UV, IR) and dust and gas temperatures for each shell from heating and cooling balance, until there are no significant temperature, pressure or density changes from one iteration to the next. Afterward, the collapse is initiated starting at the centre. The compression of the cloud is initiated by destabilizing the core in equilibrium by increasing the density throughout the cloud by a factor of two. When running the simulation, the cloud is evolved a few thousand years, until all cells have a positive in-fall velocity and the core can be considered to undergo a full gravitational collapse. The time for full core collapse to be initiated depends on the core parameters, mainly on the density and total core mass. This is also the time in cloud evolution when we start to use the data for the polynomial equation modelling, with end being the time of a formed protostar. In Figure 1, shell radius, density, dust temperature and gas temperature are displayed for a modelled 5.0   M ⊙ cloud at the start of a full gravitational collapse.

Figure 1

State of a 5.0  M ⊙ prestellar cloud at a full gravitational collapse. The parameters displayed are radius, density, dust, and gas temperatures for all shells.

In the original model, the core parameters were hardcoded in the source code. A method was developed to generate specifically formatted input data and read parameters from a file. The model updating was achieved by interpolating a cubic spline curve created from mass coordinates, which correspond to shell distance from centre (in cm) and density (g  cm − 3 ). The curve is integrated to find the column density σ , g  cm − 2 , from the centre to the edge of the cloud core. The total mass of the core is determined by a simple equation for a sphere, r 3 × π × σ r . From the total mass, single shell mass can be obtained, after which the integration points, or radius values were determined by using the Monte Carlo method with a low convergence tolerance. These data – shell mass and radius for each shell, are then used as initial conditions for the model. The cubic spline curve fitting and integration was done by using the ALGLIB C++ library.^[1]

For generating the new data, we picked seven cores with masses of 0.5, 1, 3, 5, 10, 20, 30  M ⊙ , and seven shells with numbers 500, 1,000, 3,000, 5,000, 10,000, 20,000, and 30,000, in which 4 physical values were parametrized: the shell’s density, column density, as well as gas and dust temperatures. The shells have different distances from the centre of the core, and this distance is also affected by the total core mass. These specific shells were picked to sufficiently probe the entire prestellar cloud. The shell number is a good indication of how deep each shell peneterates into the core, with smaller numbers being closer to the centre. All of the modelled cores have the same number of total shells (40,000), and each shell resides near the same relative positions in all cores at the beginning of simulation.

For creating the fits, the curve fitting was done by using a widely used python library NUMPY and determining the necessary amount of coefficients for a low error tolerance and convenient data usage. For applications in astrochemistry, a smooth evolution without artificial curves is essential. Thus, a target maximum error of 10% was deemed more than acceptable to ensure that the data quality is sufficient for use in studies of prestellar cores. In Figure 2, the fits are displayed for a 5.0  M ⊙ prestellar cloud core. In these plots, if the original data were plotted, they would be indistinguishable from the fitted data because of low errors.

Figure 2

Graphical representation of polynomial fit curves for a 5.0  M ⊙ prestellar cloud core (data from Tables 5 and 9). Darker and thinner curves represent conditions in infalling gas parcels closer to the centre, starting with shell No. 500.

3 Results

Four physical parameters were chosen to be published – gas density ρ , g  cm − 3 , column density σ , g  cm − 2 , gas temperature T g  K, and dust grain temperature T d , K. These macroscopic variables provide the essential input for 0D studies of microscopic processes in prestellar cores. While other parameters were also available, such as infall velocity, they were not deemed crucial. The values of ρ and σ can be easily converted to hydrogen atom numerical densities by multiplying with a factor 1 ⁄ ( m p M ¯ η ) , where m p is the mass of H atom, M ¯ ≈ 1.43 is average atomic mass (amu) of interstellar gas, and η ≈ 1.1 is the total number of atoms (of all elements) per one H atom. Column density σ considers the column of matter between the gas parcel in consideration toward the nearest outer margin of the modelled prestellar core. The shortest column is what primarily determines the amount of external irradiation by the interstellar UV radiation field. While interstellar grains have a certain size distribution (Mathis et al. 1977), parameter T d gives temperature for the standard grain size of 0.1  μ m, often considered in astrochemical studies.

The main results are the coefficients obtained from curve fitting the four parameters for seven shells for each of the seven cloud core masses, which are summarized in Tables 2–9, each for a different stellar mass. Curves were fitted from the parameter and the elapsed time in years. Modeled prestellar cores have different evolution times, from 122290 to 179359 years, since we are only modelling the collapsing core evolution stage with prestellar core lifetimes t end summarized in Table 1. These evolutionary times start at infall and end when a central protostar is formed, which is also our stopping criteria all modelled clouds. The polynomial parameter fits (Tables 2–9) are valid for prestellar core evolution times t up to t end . Extrapolation for up to a few 1 0 4  yr after t end can be useful in some applications. It can be a reasonable approach for the low-mass cores, where the appearance of the sink particle may not mean formation of a luminous protostar with immediate effects on its environment. Extrapolation is not verified by our error assessment and the necessity and outcome for each such data use must be cautiously reviewed as it could become less physically accurate if a hot core phase ensues.

To use the data and get the value of a parameter at a certain evolution time, replace the variable X in the polynomial with the time of evolution. For example, for a 1  M ⊙ cloud core shell No. 5,000, the gas temperature will be given by values in Table 3 and give the equation:

This equation at a time X = 81,527 years gives 16.3388 K, which compared to modelled value of 16.3653K is with an error of 0.0265 K or 0.16%. In Tables A1–A7, the maximum and mean error is given for each curve compared to modelled values. Tables A1–A7 also contain the starting and ending distances for each shell from the centre of core. The core’s evolution time in years and the total radius of the core (i.e., the outermost shell No. 40,000) are given in Table 1. When fitting the data, we adopted the following criteria: the number of coefficients must be between 5 and 15, and the maximum error should not exceed 10%. These criteria were arrived at during the modelling, as most of the parameters provided stable fits at five coefficients. Some of the fits with more than five coefficients were too unstable to be left at only five coefficients, having the average error of up to 1,000%.

4 Summary

In total, there are 196 fits, of which 22 that have more than 5 coefficients are given in a separate longer Table 9. Unfortunately, the low mass core shells that start the closest to the centre of the core are the most unstable for polynomial fitting, which means a maximum error of up to 31% even with 15 coefficients. There are only five fits that have a maximum error of above 10% with 15 coefficients, but they are still included in the tables. Given our spherical-core assumption and the uncertainties in prestellar core observations, even the erroneous fits can be considered sufficiently accurate for application in astrochemical or other studies.

The model used in our work is significantly more advanced and detailed than in some previous works, for example (Rawlings et al. 1992), which opens an opportunity to update and extend their results using the equations presented in this article. There are similarities in used models, when collapse is initiated both models exhibit a characteristic density profile of r − 3 ⁄ 2 and similar cloud sizes. While Rawlings et al. (1992) are also discussing the phase after a protostar has been formed, which is out of scope of this article, our employed model can be used to evolve clouds for a time after formation of the central protostar.

The polynomial fits allow a simple calculation of physical conditions in a variety of evolutionary scenarios relevant to prestellar cores. While central parcels close to the would-be protostar are most interesting, many observations, such as those of background stars, sample gas in the outer envelope, which is also well represented in our datasets. The relatively narrow time-span of a prestellar core may require additional input of conditions relevant to dense clouds and starless cores before the prestellar stage, as well as conditions for the protostellar envelope and, eventually, protoplanetary disk after the prestellar stage.

5 Data tables

Tables 2, 3, 4, 5, 6, 7, 8, 9,