Chemical evolution during the formation of a protoplanetary disk

(Abridged) The aim of this study is to investigate the chemical evolution from the prestellar phase to the formation of the disk, and to determine the impact that the chemical composition of the cold and dense core has on the final composition of the disk. We performed 3D nonideal magneto-hydrodynamic (MHD) simulations of a dense core collapse using the adaptive-mesh-refinement RAMSES code. For each particle ending in the young rotationally supported disk, we ran chemical simulations with the three-phase gas-grain chemistry code Nautilus. Two different sets of initial abundances, which are characteristic of cold cores, were considered. The final distributions of the abundances of common species were compared to each other, as well as with the initial abundances of the cold core. We find that the spatial distributions of molecules reflect their sensitivity to the temperature distribution. The main carriers of the chemical elements in the disk are usually the same as the ones in the cold core, except for the S-bearing species, where HS is replaced by H$_2$S$_3$, and the P-bearing species, where atomic P leads to the formation of PO, PN, HCP, and CP. However, the abundances of less abundant species change over time. This is especially the case for"large"complex organic molecules (COMs) such as CH$_3$CHO, CH$_3$NH$_2$, CH$_3$OCH$_3$, and HCOOCH$_3$ which see their abundances significantly increase during the collapse. These COMs often present similar abundances in the disk despite significantly different abundances in the cold core. In contrast, the abundances of many radicals decrease with time. A significant number of species still show the same abundances in the cold core and the disk, which indicates efficient formation of these molecules in the cold core. This includes H$_2$O, H$_2$CO, HNCO, and"small"COMs such as CH$_3$OH, CH$_3$CN, and NH$_2$CHO.


Introduction
Planets are known to form in the disk of gas and dust that arises during the star formation process following the accretion of material onto the protostar. The chemical composition of the disk should consequently determine the composition of the planetesimals that will lead to the future planets (e.g., Thiabaud et al. 2015;Kamp 2020). Characterizing the disk composition in young stellar objects is therefore of significant interest for our understanding of planet formation.
To make an inventory of the molecules present in disks, observational studies have been carried out towards several protoplanetary disks (e.g., Dutrey et al. 2014;Loomis et al. 2020). Although some complex organic molecules (COMs) such as methanol (CH 3 OH) and methyl cyanide (CH 3 CN) have been detected in the gas phase with the Atacama Large Millimeter/submillimeter Array (ALMA), their detected emission does not stem from the midplane, but from higher regions of the disk (Öberg et al. 2015;Walsh et al. 2016). Many species remain undetected compared to what is found in much younger protostars, as these species are expected to be frozen out on grains due to the low temperatures in the midplane of disks. The composition of ices was investigated with the Spitzer Space Telescope and the Infrared Space Observatory, but the detections were limited to the most abundant species (H 2 O, CO, CO 2 , CH 4 , CH 3 OH, NH 3 , e.g., Boogert et al. 2008;Öberg et al. 2008) and associated to protostellar envelopes and cold cores instead of disks. It is consequently necessary to wait for the James Webb Space Telescope to hopefully detect a large variety of species in the ices of disks.
Another efficient way to investigate the chemical composition of disks consists in running chemical simulations for realistic evolutions of the physical conditions from the molecular cloud phase up to formation of the disk. Two approaches are usually considered for the physical evolution: two dimensional (2D) semi-analytic models (e.g., Visser et al. 2009Visser et al. , 2011Drozdovskaya et al. 2014Drozdovskaya et al. , 2016 and multi-dimensional (2D or 3D) hydrodynamic or magneto-hydrodynamic (MHD) simulations (e.g., van Weeren et al. 2009;Furuya et al. 2012;Hincelin et al. 2013Hincelin et al. , 2016. The advantage of chemical models is that they also allow us to explore the chemical diversity in more detail and in particular to characterize the factors that impact chemical composition. Physical changes are known to affect the chem-Article number, page 1 of 44 arXiv:2010.05108v2 [astro-ph.SR] 14 Oct 2020 istry. It remains to be known whether the chemical composition of disks is mainly determined by their physical evolution or if the composition of cold cores can play a major role. Cold cores have been shown to present a diversity of chemical compositions (Lefloch et al. 2018;Ruaud et al. 2018), which begs the question of whether or not this early diverse composition can be inherited by the disk.
In this study, we used 3D MHD simulations of a dense core collapse to investigate the chemical evolution that occurs from the cold core phase to the formation of a young disk. We tested two different sets of initial abundances and compared the abundances of the particles ending in the disk to characterize their impact on the composition of the disk. The paper is organized as follows. The physical and chemical simulations are described in Sections 2 and 3, respectively. The results are presented in Section 4, and then discussed in Section 5. Finally, we conclude in Section 6.

Physical model
We performed 3D dense-core-collapse calculations using the adaptive-mesh-refinement (AMR) RAMSES code (Teyssier 2002), which integrates the equation of resistive MHD Fromang et al. 2006;Masson et al. 2012). RAMSES uses a finite-volume second-order Godunov scheme to integrate the conservative MHD equations and a constrained transport algorithm for the induction equation (Evans & Hawley 1988). In this study, our models only account for ambipolar diffusion, which has been found to be the most efficient resistive effect regulating the disk formation Masson et al. 2016).
We use Boss & Bodenheimer (1979) initial conditions, which consist of a uniform temperature and density sphere of mass M 0 = 1 M . The initial temperature, density, and radius are T 0 = 10 K, ρ 0 = 2.3 × 10 −18 g cm −3 (∼ 6 × 10 5 cm −3 ), and R 0 = 3 940 au. The initial ratio of thermal energy to gravitational energy is α = 0.4. We do not apply density perturbation or solid body rotation, but we introduce initial turbulent velocity perturbations following a Kolmogorov power spectrum. More details on the initial setup of the turbulent velocity field can be found in Hennebelle et al. (2011), Commerçon et al. (2011b), and Joos et al. (2013). The amplitude of the velocity perturbations are scaled such that the initial turbulent Mach number is trans-sonic: M turb,0 = v rms,0 /c s,0 = 0.9, with v rms,0 being the initial rms turbulent velocity and c s,0 the initial isothermal sound speed at 10 K. To mimic the thermal behavior of the collapsing dense core during the initial collapse and first Larson core formation (Larson 1969;Commerçon et al. 2011a;Vaytet et al. 2013), we use a barotropic equation of state (EOS) with P the thermal pressure, ρ the gas density, and γ = 5/3 the ratio of specific heats. ρ c = 10 −13 g cm −3 represents the critical density at which the gas turns adiabatic. The gas temperature is given by the barotropic EOS (eq. 1), which relates the density to the temperature. Last, the initial magnetic field (B 0 ∼ 75 µG) is set to be uniform in the z-direction, with a mass-to-flux over critical mass-to-flux ratio of µ = 5, where µ = (M 0 /φ)/(M/φ) crit , φ = πB 0 R 2 0 is the magnetic flux, and (M/φ) crit = 0.53/3π × (5/G) 1/2 (Mouschovias & Spitzer 1976).
To estimate the ambipolar diffusion resistivity, we use the equilibrium abundances tables from Marchand et al. (2016).
The initial grid has a resolution of 2 6 cells per direction, and the simulations box length is 8R 0 . Outside the initial dense core, the gas is in thermal equilibrium (i.e., T = 10 K) but has a density of ρ 0 /100. The grid is refined according to the Jeans length criterion, with at least 12 cells per Jeans length. The maximum resolution is 2 16 , which is 0.5 au. This type of initial conditions and numerical setup has been extensively used in previous studies (e.g., Commerçon et al. 2008Commerçon et al. , 2012Hincelin et al. 2016;Masson et al. 2016;Vaytet et al. 2018).
We introduced 10 6 tracer particles with an initial uniform distribution within the initial sphere, which allowed us to reconstruct their trajectories and extract the evolution of their density, temperature, and visual extinction with time. This information was subsequently used as input parameters for the chemical modeling (see Section 3). The dust temperature is assumed to be equal to the gas temperature. The time steps are determined by the stability conditions to integrate the Euler equations for the gas dynamics. In order to compute the visual extinction, we used the column density estimator presented in Valdivia & Hennebelle (2014). A similar strategy for the tracer particles was used in Hincelin et al. (2013Hincelin et al. ( , 2016. Our physical model is well adapted to the early phases of star formation, but has a few limitations, which are detailed below. First, the barotropic EOS is a crude approximation of the radiative transfer modeling of the star-disk system, but remains efficient as long as no central object with an internal luminosity (the protostar) has formed. Second, we do not employ sink particles either to describe the central collapsing part, as is done in Hennebelle et al. (2020) for instance, because we restrict our study to the first Larson core stage. Sink particles and a more accurate radiation transport scheme (e.g., Mignon-Risse et al. 2020) will be used in future works investigating a longer time evolution.

Characteristics of the disk at the final time of the simulation
We define t = 0 as the starting time of the MHD simulations. The collapsing dense core reaches the first hydrostatic core (FHSC) stage at t FHSC = 5.01 × 10 4 yr. The calculations are continued for an additional 8.2 × 10 3 yr. At the final time of the simulation (i.e., 5.83 × 10 4 yr), the particles reach different components: the collapsing envelope, the outflow, the magnetic pseudo-disk (Galli & Shu 1993), the rotationally supported disk, and the central core (i.e., the FHSC). In this paper, we focus on the rotationally supported disk only. The rotationally supported disk is the precursor of the protoplanetary disk. Hincelin et al. (2016) explained in their Section 3.1.1 the different criteria to determine the component the particles belong to. In this work, we use the coordinate system shifted to the FHSC center of mass with a vertical axis aligned with the angular momentum of the inner 100 au. For the rotationally supported disk, four criteria have to be fulfilled. The azimuthal velocity v φ has to be at least twice larger than the radial velocity v r , so that we do not consider particles that collapse too quickly. The azimuthal velocity should also be twice larger than the vertical velocity v z . Particles with a smaller azimuthal velocity would rather come from the outflow cavity wall. The rotational support is at least twice the thermal support and the density is above 10 9 cm −3 . In the end, 14 948 out of the 10 6 tracer particles introduced in the simulation end in the disk. Figure 1 shows the spatial distribution of the temperature and the density of the particles in the disk at the final time of the simulation. The disk presents spiral arms and has a diameter of about 100 au. From the outer region to the inner region of the disk, the temperature goes from 11 K to 802 K, while the density ranges between 1 × 10 9 cm −3 and 2 × 10 13 cm −3 . The upper panels of Figure B.1 show the distribution of final density and temperature of the disk particles with histograms. The particles are relatively well spread radially, while they are vertically concentrated on the midplane (see Figure B.2). We count 81% of particles with |z| < 5 au, 15% with 5 ≤ |z| < 10 au, and 4% with |z| ≥ 10 au. For the same radius, no strong variation in density or temperature is seen vertically. The particles with the highest |z| present low densities and temperatures, but mainly because their radii are also high. The vertical structure of the disk is affected by the low resolution of the AMR grid within the disk, which may smooth the density and temperature distribution to the grid scale. There are typically 5-20 cells in the vertical extent of the disk depending on the radius. Describing the vertical structure of the disk would require a much higher resolution (and physical modeling) and is beyond the scope of the present study.
2.3. Evolution of the physical conditions for the particles ending in the disk Figure 2 shows the evolution of density and temperature as a function of time of all the particles ending in the disk. The temperature remains low (10 K) for all particles until ∼5.4 × 10 4 yr. The particles follow different trajectories and encounter different physical conditions. For example, among the 1 632 particles with a final temperature < 20 K, 329 (20%) remain at low temperature along the entire trajectory (see blue line in Figure 2), while 145 (9%) reach a temperature above 100 K at some point. The remaining 1 158 particles (71%) experience a maximum temperature of between 20 and 100 K (see Figure 3). The distribution of the maximum density and temperature reached by the disk particles during their trajectories is presented with histograms in Figure B.1 (lower panels).

Chemical model
The evolution of the physical conditions of the 14 948 particles was used as input parameters in the Nautilus three-phase gas-grain chemical model to calculate the chemical evolution of the particles ending in the disk. This model is fully described in Ruaud et al. (2016). It includes gas-phase, grain-surface, and grain-mantle chemistry, along with adsorption of gas-phase species onto grain surfaces, thermal and nonthermal desorption of species from the grain surface into the gas phase, and surfacemantle and mantle-surface exchange of species. The diffusion energies of the species are equal to 0.4 times the binding energies for the surface and 0.8 times the binding energies for the mantle. The cosmic-ray ionization rate is fixed at 1.3 × 10 −17 s −1 . Direct photo-processes are negligible because the visual extinction is always higher than 30 mag. Chemical desorption (see Garrod et al. 2007) is included with an efficiency of 1% for all species. The gas-phase network is derived from kida.uva.2014 and the grain network is described in Wakelam et al. (2016). The gas and grain networks have been updated according to the following works: Wakelam et al. (2015), Loison et al. (2016), Hickson  Vidal et al. (2017), andLoison et al. (2017). In total, the network consists of 589 gas-phase species and 540 solid-phase species (either on the surface or in the ice mantle) with 13 384 reactions. It should be noted that this version of the chemical network is not complete regarding complex organic chemistry. Some COMs detected in several young solar-type protostars, such as for example glycolaldehyde and ethylene glycol (e.g., Coutens et al. 2015), are not included in this version. Molecules recently detected for the first time in such objects (e.g., CH 3 Cl, HONO, Fayolle et al. 2017;Coutens et al. 2019) are not included either. A chemical network is currently being developed to include more complex species (up to threecarbon species) and grain-surface reactions (e.g., Manigand et al. 2020a). The results presented here for the most complex species should consequently be seen more qualitatively than quantitatively.
The abundances are expressed with respect to the total density of hydrogen nuclei (n H = n(H) + 2 n(H 2 )). The abundances of molecules in ices presented in this paper always correspond to the sum of abundances obtained for the grain mantle and the grain surface. When only referring to the solid-phase species, "s-" precedes the species. Fig. 3: Temperature evolution of six particles ending in the outer region of the disk (beyond a radius of 30 au) with a temperature lower than 20 K. Two particles reach a maximum temperature above 100 K (red), two below 20 K (blue), and two between 20 and 100 K (green).

Initial abundances
To determine the impact of the initial conditions on the final chemical composition of the protoplanetary disk, we consider two different sets of initial abundances that we call A and B. These sets were obtained by Ruaud et al. (2018) with the chemical model described above for smoothed-particle hydrodynamics (SPH) simulations of dense core formation (Bonnell et al. 2013). The abundances were extracted at the time where the density has sufficiently increased to form a dense core. The two clouds that produced the two cores considered in this study underwent different physical histories. The physical conditions at the time we extracted the initial abundances for cores (maximum of density) are similar however. For A (B), the density is 4.0 × 10 5 (2.1 × 10 5 ) cm −3 , the temperature is 10 (11) K, and the visual extinction A V is 22 (18) mag. In case A, the cloud has experienced a relatively quiescent phase, while the history of B is more complex with significantly warmer and colder phases (see Figure B.3). Cloud A has also experienced a moderately dense and cold phase (at 4.3 ×10 7 yr) before reaching its maximum "cold core" density that significantly changed its composition. During this phase, a large fraction of the electron donors and the atomic oxygen have depleted onto the cold dust grains, which led to the formation of solid water (see Figure 4). Even when the density decreases (between the two density peaks), water remains on grains. Less atomic oxygen is available in the gas phase. The main consequence is that water and gas-phase carbon chains are more abundant in cloud A, while oxygen-rich organics and O 2 are more abundant in cloud B (see Figures B.4-B.7). As explained in Ruaud et al. (2018), the SPH simulations follow the gas dynamics of the interstellar medium at galactic scales. Consequently, the timescales shown in Figures 4 and B.3 do not represent the age of dense and cold cores. If we estimate the ages of the cold cores by calculating the time spent between the moment where the density reached 10 3 cm −3 and the final density, the core formed in Cloud A is older than the one in Cloud B with respective ages of 2.1 × 10 6 yr and 4 × 10 5 yr. The main carriers of the different elements are found on grains, and are summarized in Table 1.. Indeed, at the initial time, the temperature is low and the density is relatively high, which leads to the depletion of many species. For oxygen, water is dominant (91%) in case A, while in case B water only represents 48% of the O elemental abundance and the rest is in the form of organic molecules (H 2 CO, CH 3 OH, HCOOH, and HCO). The main carbon carriers are O-bearing organic molecules in case B, while HCN accounts for 27% in case A. The contributions of CH 4 are relatively similar in both cases, and CH 3 OH is more abundant in case B. In Case A, a large range of species with very small contributions (less than 5% of the C atomic abundance) explains the total amount of carbon. Regarding nitrogen, the main carrier is HCN in case A and NH 3 in case B. The sulfur reservoirs are dominated by HS and H 2 S in both cases, but differ by their slightly less abundant contributors, CH 3 SH and CH 3 S in case A, and NS and OCS in case B. SiH 4 is the most abundant Si carrier followed by SiC 8 H in case A, while SiC 4 H and SiO dominate in case B. Atomic P is the major carrier of phosphorus. CP and HCP also contribute in case A, while PN is the second-most abundant carrier in case B. The carriers of Cl and F are the same in both cases: HCl and HF.
Regarding the ionization, the most abundant ions are H + 3 , H + , and S + in case A, while they are Fe + , H + 3 , H + , S + , and Mg + in case B. In both cases, the abundances are about a few 10 −9 .

Results
The gas-grain chemical model Nautilus was run for the 14 948 disk particles and the two sets of initial abundances A and B, which are the chemical compositions at the maximum density of two clouds (Ruaud et al. 2018). The analysis of the chemical results focuses on three different aspects. First, we analyze the spatial distribution of molecules in the disk. We then study the chemical evolution during the formation of the disk. In particular, we investigate the impact of the initial abundances on the disk abundances. Finally, we explore the ionization and its role in the disk resistivity. B.8-B.20 show the spatial distribution of various molecules with non-negligible abundances in the x-y plane of the disk (stacked z axis) at the end of the simulations. Although the abundances of the same molecules can vary between cases A and B, their spatial distribution is usually very similar. Exceptions are found for HCO, HNO, OH, SO, SO 2 , and CH 3 . The spatial distribution does vary strongly from molecule to molecule. While the abundances of some species peak towards the center, others are more abundant in the colder regions that correspond to the spiral arms. Some molecules are also found to peak in the gas phase in an intermediary region. We classified in Table 2 the molecules depending on the region(s) where their gas-phase and total (gas + solid) abundances peak (within about one order of magnitude). We considered the three regions to be: the inner region of the disk (n 10 12 cm −3 , T 100 K), the spiral arms (n 10 11 cm −3 , T 40 K), and the transition region that takes the shape of a ring surrounding the inner region and includes the beginning of the spiral arms (10 11 cm −3 n 10 12 cm −3 , 40 K T 100 K, see Figure 5). It should be noted that molecules can present variations of their abundances inside these regions and consequently molecules in the same category do not necessarily present the same spatial distribution. For example, the spatial extent of the molecule can vary in the inner regions or the molecule peaks in different regions of the spiral arms. In regard to this classification, we observe different trends:

Figures
• Except for a few ions, most of the molecules that peak in the spiral arms are also abundant in the transition region. • The molecules with their highest gas-phase abundance in the spiral arms + transition region or the transition region only have their total abundance peak in the spiral arms + transition region. These molecules are only present in the cold regions either in the gas phase or in the ices. Their absence from the inner regions can be explained by efficient destruction reactions at relatively high temperatures ( 100 K). • No molecules appear in the transition region only when we consider total abundance peaks. The transition region observed for some molecules in the gas phase can be explained by the release of these species trapped on grains below a certain temperature, but they are easily destroyed in the gas phase at warmer temperatures. Indeed, this category of molecules includes a lot of radicals. • Most of the molecules with their highest gas-phase abundance in the inner regions present their total abundance peak in all the regions. This means that these species are also abundant on grains in the region where the temperature is low. Thermal desorption explains their high gas-phase abundances in the inner regions. Many COMs and water belong to this category.
Vertically, we do not see strong variations in abundances. Figure 6 shows the distribution as a function of the radius and height (z) of the gas-phase abundances of six species that belong to the different components discussed above (H + 3 and CH + 5 for s-H 2 CO (14%), s-CH 3 OH (13%), s-HCOOH (11%), s-CH 4 (10%), s-HCO (10%), s-CO (8%), s-CO 2 (5%), s-CH 2 OH (5%) Note: s-means that the species are on grains. The gas-phase abundances are negligible here.  Note: We consider that the variation of abundance between two components is significant when the abundances change by at least an order of magnitude. ‡ These species appear in two categories, as their spatial distribution is sensitive to the chosen set of initial abundances.
the spiral arms only, CCH for the spiral arms + transition region, CN for the transition region only, H 2 O for the inner disk, and CH 4 for all regions). As explained in Section 2.2, the vertical structure of the disk is affected by the low resolution of the AMR grid within the disk. The disk is also embedded, so no chemical stratification due to UV effects is expected contrary to evolved disks. Figure 6 clearly shows that the spatial distributions of the molecules are correlated with the temperature and density.

Final reservoirs
The final reservoirs of the different elements are summarized in Table 3. To do so, we averaged the abundances found for all the disk particles. The reservoirs of oxygen are very similar to the initial ones. Only in case B, does HCO decrease from 5% to 1%. The same is seen for the carbon. In case B, only HCO and CH 2 OH become less abundant. As we see later, these two rad-  : the inner regions (n 10 12 cm −3 , T 100 K, red), the spiral arms (n 10 11 cm −3 , T 40 K, blue), and the transition region (10 11 cm −3 n 10 12 cm −3 , 40 K T 100 K, green).
icals are consumed to form more COMs. The main carriers of nitrogen, silicon, chlorine, and fluorine do not change either. A change is seen for phosphorus, but only in case B. Atomic P is still the main carrier, but the contributions of PO, HCP, and CP significantly increase at the final time. Sulfur is the only element for which major changes are observed in both cases. At the end of the simulations, HS is no longer the main carrier, while the contribution of H 2 S 3 has strongly increased. However, the other carriers are the same as in the initial conditions. H 2 S, CH 3 SH, and CH 3 S still contribute to more than 5% in case A, while H 2 S, NS, and OCS remain similarly abundant in case B. Based on these results, we can conclude that the main carriers of the different elements in the disk are generally inherited from the molecular cloud, but this does not necessarily mean that the chemical content does not evolve during the formation of the disk, as we see in the following section.

Evolution of the chemical composition
Here we investigate the evolution of the chemical content during the formation of the disk more thoroughly. The histograms presented in Figures B.21-B.27 show the abundance distribution for all the disk particles at the end of the simulations, while the initial abundances and the mean, median, minimum, and maximum abundances at the final time are summarized in Table A.1 and shown in Figures B.4-B.7. The distributions obtained for the set of initial abundances A and B are compared with each other as well as with their respective initial abundances to determine the impact of the initial conditions on the chemical composition of the disk. The distributions of the final abundances are considered similar if they overlap by more than 75% over bins of one order of magnitude. The initial abundances are considered similar if they do not differ by more than one order of magnitude. The final abundances are considered similar to the initial abundance if the median of the final abundances and the initial abundance differ by less than one order of magnitude.
Based on these criteria, we classified the molecules into different categories (see Table 4). To highlight the change in chemistry during the formation of the disk, we first separated the molecules that have similar initial and final abundances (third column of Table 4) and the ones for which the abundance significantly changes between the initial and final times (fourth column of Table 4). In addition, we separated the molecules with similar initial abundances (rows 1 and 2) and the ones with different initial abundances (rows 3 and 4). Finally, we also separated the species with similar final distributions (rows 1 and 3) and the ones that show significantly different distributions (rows 2 and 4). A visual inspection was made to check the classification. NH 3 and HNC (indicated in Table 4 with brackets) have been moved to a different category, as their initial abundances differ by almost an order of magnitude.
From this classification, we can infer: ) as a function of radius ( x 2 + y 2 ) and the height z at the final time of the simulations (5.83 × 10 4 yr) for the set of initial abundances A. For the molecules, the color scales represent log 10 (n(X)/n H ).
• All the molecules in the third column of Table 4 present similar initial and final abundances. Although we cannot exclude the destruction and reformation of these species during the collapse, it is more likely that these species are essentially formed in the molecular cloud or prestellar core.
This list includes various types of molecules. Water is one of them, which is not so surprising since both observations and models are in favor of formation at a relatively early stage of the star formation process (e.g., Jones & Williams 1984;Coutens et al. 2012;Furuya et al. 2016Furuya et al. , 2017. Some "small"   Note: [X] i,x and [X] f,x correspond, respectively, to the initial and final abundances for the set of initial abundances x (A or B). We consider the total gas and ice abundances. On the last row, some molecules are found with both [X] i = [X] f and [X] i [X] f . The category depends on the set of initial abundances. To differentiate them, we used a for set A and b for set B. NH 3 and HNC (within brackets) have been moved to a different category, as their initial abundances differ by almost an order of magnitude.
COMs such as methanol, methyl cyanide, formamide, and their potential precursors (HNCO, H 2 CO) also belong to this category. Although deuteration is not included in our network, it is interesting to note that these molecules are the organic species with the lowest deuterium fractionation towards the Class 0 protostar IRAS 16293-2422 B, while the most complex species (absent from this category) present a deuteration level at least a factor of two higher Jørgensen et al. 2016Jørgensen et al. , 2018Persson et al. 2018). The lower deuteration of the "small" COMs would consequently be explained by their formation at an earlier stage of the star formation process in the molecular cloud.
• All the molecules in the fourth column of Table 4 show different initial and final abundances for at least one of the chosen sets of initial abundances. The physical changes produced during the dense core collapse consequently play an important role in the chemical evolution of these species.
In Table 5, we separate the molecules that show a clear increase in abundance from the ones showing a clear decrease.
The molecules that present an increase in abundances in both cases are "large" COMs (CH 3 CHO, CH 3 NH 2 , CH 3 OCH 3 , HCOOCH 3 ). This emphasizes the importance of the collapse in the formation of these molecules typically found in "hot cores" and "hot corinos". Although these molecules are detected towards prestellar cores (Bacmann et al. 2012;Vastel et al. 2014), our model shows that they only acquire their highest (gas + solid) abundances during the collapse. In contrast, the molecules that show decreasing abundances in both cases are relatively small molecules and radicals. In particular, they include radicals such as HCO and CH 3 O, which are consumed to form more COMs. The radicals CCH, C 4 H, CH 3 , CN, NH, NH 2 , NO, and OH also show this decrease in abundances. The abundances of several ions also decrease significantly. O 2 is also destroyed during the collapse. • For the molecules in row 1, the impact of the initial conditions on the final abundances cannot be properly tested because the initial abundances of these species are similar. These species seem to have little sensitivity to the history of the physical conditions from the diffuse medium to the disk formation.
• The molecules in row 2 present different final abundances despite having similar initial abundances. These species are sensitive to the initial abundances of other species that differ in cases A and B; they include three radicals (CH 3 O, HCO, and OH) and two ions (HCO + and PCH 3 + ). • The molecules in row 3 present similar final abundances despite having different initial abundances. These molecules are particularly interesting, as they are not sensitive to the initial conditions. Their abundances are determined during the collapse. Here, we again find some "large" COMs: CH 3 CCH, CH 3 CHO, CH 3 OCH 3 , H 2 CCO, and HCOOCH 3 . • Finally, a large number of species present different initial abundances, but also different final abundances (row 4). The species with similar initial and final abundances ([X] i = [X] f ) are likely only sensitive to the evolution of the physical conditions leading to the cold and dense core. Once the collapse starts, they do not evolve further. In contrast, for the few species with different initial and final abundances ([X] i [X] f ), both the initial prestellar phase and the collapse play a role in their different abundances. This is the case of some large COMs such as CH 3 COCH 3 and CH 3 NH 2 , as well as some small molecules such as CO, NO, O 2 , NH, NH 2 , and PO.

Ionization and resistivity
The chemistry module we used to estimate the ambipolar diffusion resistivity in the models is not consistent with the full gas-grain chemistry we follow on the tracer particles. Currently, performing full MHD calculations with full gas-grain chemistry is not achievable and models can follow the chemical evolution of networks with about 50 species (e.g., Dzyurkevich et al. 2017). Nevertheless, we can compute the resistivities from the full gas-grain chemistry and compare it with the ones we used from Marchand et al. (2016).
At the end of the simulations, the ionization level has significantly decreased compared to its value at the initial time (initial abundances of a few 10 −9 for the most abundant ions; see Sect. 3.2). The most abundant ions are CH + 5 and H + 3 in both Notes: As the abundances vary from particle to particle, we consider that the abundance increases or decreases if it concerns more than two thirds of the particles (10 000). The molecules in the column "others" are the ones that cannot be classified. For example, they can present a spread distribution around the value of the initial abundance or they have similar initial and final abundances in the given case. cases with similar average abundances of a few 10 −12 . In case A, C 2 H + 7 also presents an average abundance of 1 × 10 −12 . In case B, H 5 C 2 O + 2 , PCH + 3 , HNS + , and C 3 H + 9 have abundances of (4-7) × 10 −13 . However, these different ions peak in different regions of the disk (see Section 4.1) and their abundances vary significantly among the different particles. For example, in Figure 7 we see that H + 3 peaks in the regions of lower densities, while CH + 5 is particularly abundant in the regions with a density between 10 10 and 10 12 cm −3 . Several molecular ions including H 5 C 2 O + 2 contribute to the positive charge in the densest regions. Electrons and grains mainly contribute to the negative charge (see Figure 7).
Ionization has a strong effect on the resistivity of the disk. The nonideal magneto-hydrodynamic resistivities are of three types: ambipolar diffusion, Ohmic diffusion, and Hall effect. To estimate them, we used the equations 2-8 from Marchand et al. (2016) and calculated them separately for each particle. Figure 8 shows the spatial distribution of the three resistivities as well as the ratios of the resistivities. The initial abundances do not appear to have any significant effect on the resistivities of the disk. The Hall effect is the dominant mechanism, except for densities of about 10 12 cm −3 , where the ambipolar resistivity is sometimes higher than the Hall resistivity. The Ohmic resistivity is always lower. Our results show that electrons, negatively charged grains, and a large variety of molecular ions (H + 3 , CH + 5 , H 5 C 2 O + 2 , C 2 H + 7 , NH + 4 , HCNH + , C 3 H + 9 , H 2 COH + , etc.) need to be included in the chemical network to properly calculate the resistivity of the disk. While at low densities ( 10 10 cm −3 ) H + 3 is the most important molecular ion to take into account, more complex molecular ions contribute to the resistivity at the highest densities (> 10 12 cm −3 , see Figure 7).
Although we consider a unique size of grains (0.1 µm) in the chemical modeling, we examined whether or not grain growth, which is known to occur in disks, has an impact on the resistiv-y (au) y (au) y (au) x (au) x (au) x (au) ities. Tests were carried out for two particles. The first particle ends in a dense and warm region (n = 1.7 × 10 12 cm −3 , T = 170 K), the other one in a colder and less dense region (n = 7.0 × 10 9 cm −3 , T = 14 K). We find that assuming a grain size ten times larger does not change the previous results. The values of the resistivities are different by less than a factor of 2, except for the ambipolar resistivity of the warmest particle, which differs by a factor of 2.5. In both cases, the Hall resistivity dominates and the species that contribute to the resistivities are also the same ones as for a grain size of 0.1 µm.
The resistivity we obtain from the full gas-grain chemistry is in qualitative agreement with that obtained in the literature, where MRN grain size distributions are used Dzyurkevich et al. 2017;Zhao et al. 2018). The ambipolar and Hall resistivities dominate the Ohmic diffusion coefficient by orders of magnitude, and have similar amplitude, within a factor of ten in the disk. We can conclude that, although full gasgrain chemistry is required to follow the chemical history and budget of the forming protoplanetary disk, it is not necessarily required for the MHD resistivity. Even though we do not see any significant effect on the resistivities for grain sizes of 0.1 and 1 µm, other studies show that dust size could be a key physical quantity (Zhao et al. 2016;Dzyurkevich et al. 2017). More work needs to be done in this field, in particular to quantify the effect of dust growth and size distribution on ionization.

Comparison to 2D semi-analytical models
Two-dimensional semi-analytical models were previously used to investigate the chemical history of the material present in the disk at the end of the collapse. In particular, Visser et al. (2009) investigated H 2 O and CO ices. No chemical reaction was taken into account in their model apart from desorption and adsorption of H 2 O and CO. However, their study predicted that the infalling material spent enough time at warm temperatures during the collapse to abundantly form COMs on grains, which agrees with our findings.
Later, Visser et al. (2011) extended the chemical analysis with a gas-phase network, including adsorption onto and desorption from dust grains, and a basic grain-surface network with only hydrogenation reactions. The different physical conditions encountered along the computed trajectories merited division of the disk into several zones with different chemical histories. In our simulations, the situation is somewhat more complex. We see spatial variations in abundances correlated with the temperature, but the trajectories can also vary significantly for particles ending in the same area of the disk (cold or warm regions; see Sect. 4.2.2), which can lead to a large range of abundances for some species in the same region.
Following the same approach, Drozdovskaya et al. (2014) investigated the evolution of methanol with a larger and more comprehensive gas-grain chemical network. These latter authors obtained an abundance of methanol ice of 8 × 10 −7 at the end of the prestellar phase and found that the ice abundance in the disk either increases until a value slightly above 10 −6 or significantly decreases in some cases depending on the region of the disk and the considered dominant disk-growth mechanism (viscous spreading or continuous infall of matter). In our model, the initial abundance is higher, (1-3) × 10 −5 . The reaction-diffusion competition can explain part of the difference. The formation of methanol relies on two activation barriers, which are lowered due to the reaction-diffusion competition. Drozdovskaya et al. (2016) found a higher value of 3.7 × 10 −6 at the end of the prestellar phase, when taking into account reaction-diffusion competition. Different physical conditions and timescales during the precollapse phase could also explain the different initial abundances obtained for methanol (Kulterer et al. 2020). The total methanol abundance in our model does not evolve significantly during the disk formation and cannot be enhanced as methanol is already formed abundantly in the cold and dense core. The disk obtained with the 3D MHD simulations is also still young (∼ 5.8 × 10 4 yr), while the one in Drozdovskaya et al. (2014) reaches 2.5 × 10 5 yr. More than 6 × 10 5 tracer particles are still in the envelope at the end of our simulations. No powerful outflow is seen. Consequently, we do not consider UV irradiation from the star contrary to Drozdovskaya et al. (2014). Once the protostellar object has sufficiently evolved, high stellar UV fluxes should photodissociate both solid and gaseous methanol in the innermost region of the disk as well as on the disk surface, which explains the decreasing abundance of methanol seen in Drozdovskaya et al. (2014).
With the same model, Drozdovskaya et al. (2016) studied the evolution of various COMs. Their abundances are found to be high in the outer part of the disk and decrease towards the center because of the high UV stellar irradiation. To compare the results of these latter authors with those from the present study, it is consequently better to focus on the abundances in the outer midplane of the disk. The abundance of HCOOCH 3 found by these latter authors increases significantly during the collapse similar to our results. They find that the average abundance of methanol is relatively similar at the end of the initial prestellar phase and at the final time in the disk, in agreement with our findings. For the other complex molecules common to the two studies (H 2 CCO, HCOOH, CH 3 CHO, CH 3 OCH 3 ), the situation is different. Their abundances of H 2 CCO at the end of the initial prestellar phase and in the disk are relatively similar, while we see an increase in H 2 CCO. HCOOH is more abundant at the end in their model, while no significant change is observed in our study. CH 3 CHO and CH 3 OCH 3 become more abundant during the formation of the disk in our study, while the average abundance of these molecules decreases more or less depending on the chosen dominant disk growth mechanism in Drozdovskaya et al. (2016). In spite of these different behaviors, the final abundances of the COMs we obtained are in good agreement with those of these latter authors for an infall dominated disk (except for H 2 CCO which is significantly less abundant in our study). In their spread-dominated disk, the abundances are generally much lower. Regarding the smaller molecules, the results are also in agreement, as Drozdovskaya et al. (2016) found that most of them (H 2 O, NH 3 , CH 4 , H 2 S, N 2 ) show relatively similar initial and final average abundances. Differences are seen for CO 2 . In our case, CO 2 does not evolve significantly, while the abundance of CO 2 is found to increase by these latter authors and to reach higher values than in our study. The stellar FUV irradiation would explain this difference as well. The abundance of H 2 S is also two orders of magnitude higher in our study than in Drozdovskaya et al. (2016). Indeed, the initial elemental sulfur abundance is lower in their study. A depleted sulfur abundance is imposed, because their chemical network is not designed to reproduce the observed depletion of sulfur in dark clouds (Vidal et al. 2017). Aikawa et al. (2008) initiated studies of chemical evolution of collapsing cores with radiation hydrodynamic simulations. However, the 1D symmetry adopted in that study prevented any investigation of the molecular evolution during the formation of the disk. van Weeren et al. (2009) first studied molecular evolution in a forming disk with a 2D hydrodynamical model. This demonstrated that the employed method could be extended to 3D modeling and higher resolution as long as the chemical modeling of the particles is parallelized. Furuya et al. (2012) were the first ones to investigate the chemical evolution from a cold core to a first hydrostatic core using a 3D radiation hydrodynamic simulation. Contrary to our study, they did not take into account any magnetic field. The presence of a magnetic field leads to the formation of more realistic disks, which are less inclined to early fragmentation (Commerçon et al. 2011b;Hincelin et al. 2013). As Furuya et al. (2012) mainly focused on the central core (< 10 au), no direct chemical comparison can be made with this work, but their main conclusion was that the total ice and gas abundances of most species were not altered until a temperature of 500 K, as the warm-up phase was too short to form large organic species. In our model, COMs such as HCOOCH 3 , CH 3 CHO, and CH 3 OCH 3 have time to form efficiently.

Comparison to (magneto-)hydrodynamic models
Later, Yoneda et al. (2016) used a radiation-hydrodynamics model to investigate the evolution of the chemical composition of 1.5 × 10 3 SPH particles from the cold prestellar core to the rotationally supported gravitationally unstable disk. In agreement with our results, they found that molecules such as H 2 O, CH 4 , NH 3 , CO 2 , H 2 CO, and CH 3 OH are already abundant at the onset of gravitational collapse and are later sublimated as the SPH particles migrate inside the water snow line, while COMs such as HCOOCH 3 and CH 3 OCH 3 mainly form during the collapse. However, contrary to Yoneda et al. (2016), we found that HCOOH, H 2 CS, and SO 2 mainly form before the collapse. Hincelin et al. (2013) explored the survival of molecules during the collapse of a prestellar core up to the first Larson core using 3D MHD models. Compared to our study, their simulations assumed ideal magnetohydrodynamics and stopped at an earlier time (3.8 × 10 4 yr). The magnetic fields are also lower (µ=10 and µ=200). The physical structure of the obtained disks differ depending on the intensity and inclination Θ of the magnetic field. Most of the particles in their simulated disks have low temperatures (T < 100 K), whereas in our longer simulation the temperature has time to significantly increase in the inner re-gions (because of the use of the barotropic equation of state, this is due to the density increase in the center). These latter authors concluded that the abundances of the most abundant species do not significantly evolve during the collapse. This is in agreement with our findings relative to the main carriers of the C, N, and O elements. Furthermore, Hincelin et al. (2013) only see strong variations for HNC, CO 2 , and HNO. In our study, we do not see such changes for these species, but these variations are not particularly important for their model with µ=10, Θ=0 which is the closest to our model (µ=5, Θ=0). It should be noted that, in our study, a significant decrease is observed for CO when the set of initial abundances A is chosen. The abundance of CO can consequently vary when its initial abundance is not particularly high. These latter authors also noted a similar trend for this molecule when testing the impact that an older molecular cloud would have on the chemistry of a particle arriving in the warm region (∼ 30 K). Hincelin et al. (2016) presented an additional simulation with µ=2, Θ=0, but the strong magnetic field in this case prevents the formation of a rotationally supported disk.

Comparison to observations
The disk obtained at the end of the simulations shows high densities (10 9 -10 13 cm −3 ) and covers a large range of temperatures between 11 and 802 K. The region where the temperature is higher than 100 K has a diameter of about 50 au. In terms of physical conditions, this region is reminiscent of hot corinos, the warm inner regions of young low-mass protostars, which present a rich complex organic chemistry.
IRAS 16293-2422 (hereafter IRAS16293) is certainly the most studied of these objects. It has been the target of many observational projects with ALMA. The kinematics were studied at small scales towards the two components A and B of the system. Observations of rotation motions around source A suggested the presence of an almost edge-on disk (Pineda et al. 2012;Oya et al. 2016). More recently, Maureira et al. (2020) showed with very high-spatial-resolution observations that component A consists of two compact sources, A1 and A2. A dust disk with a FWHM size of ∼12 au is seen towards A2, while an upper limit of ∼4 au is derived for the disk around A1. Infall motions are seen towards source B (Pineda et al. 2012;Zapata et al. 2013). A face-on disk has been suggested for Source B (Zapata et al. 2013;Oya et al. 2018). The two components A and B are characterized by a hot corino chemistry (e.g., Bottinelli et al. 2004;Bisschop et al. 2008;Caux et al. 2011;Jørgensen et al. 2016;Manigand et al. 2020b). Many COMs were detected for the first time in low-mass protostars towards this source, especially in the framework of the ALMA Protostellar Interferometric Line Survey (PILS, Jørgensen et al. 2016). The emission of COMs is quite compact around the sources, similarly to what we see in our simulations. Indeed, these species are expected to be released in the gas phase through thermal desorption once the temperature is sufficiently high (T 100 K). The relative abundances of the COMs were measured towards both components of the binary (e.g., Manigand et al. 2020b). Given the incompleteness of our network regarding complex organic chemistry, a comparison is only possible for a few COMs. Figure 9 shows a comparison of the relative abundances with respect to CH 3 OH. While some species show good agreement with our model (e.g., CH 3 OCH 3 , HCOOH, HNCO), it is clear that other species are underproduced (CH 3 OCH 3 , H 2 CCO) or overproduced (H 2 CO), which suggests that some key reactions are missing in this network. A chemical network is currently being developed to include all COMs detected in low-mass star-forming regions. Sim-ulations using this future network would be necessary to model the chemical richness of this source in detail.
As stated in Section 4.2.2, even if our network does not include isotopic fractionation, an interesting correlation is seen between the level of deuteration found for several species in IRAS16293 Persson et al. 2018;Jørgensen et al. 2018;Calcutt et al. 2018;Manigand et al. 2020b) and their period of formation in our simulations. On one hand, the species showing a low deuteration of 1-3% in IRAS16293 B (HNCO, NH 2 CHO, HCOOH, CH 3 OH, H 2 CO, and CH 3 CN) are the ones presenting similar initial and final abundances in our model. On the other hand, the species with the highest deuteration, 4-8% (HCOOCH 3 , CH 3 OCH 3 , CH 3 CHO) are the ones that become significantly more abundant during the collapse. It should be noted that the abundance of H 2 CCO, which shows a low deuteration in IRAS16293, increases for the set of initial abundances A, but remains similar for the set B. Jørgensen et al. (2018) proposed that this variation in deuteration would reflect the formation time of each species in the ices before or during warm-up and/or infall of material through the envelope. The results obtained here confirm this interpretation. The species formed during the cold core phase in our simulations also show low deuteration towards IRAS16293 A (Manigand et al. 2020b). However, some of the species efficiently formed during the collapse (e.g., HCOOCH 3 ) do not show D/H ratios as high as in source B (Manigand et al. 2020b), which would reveal differences stemming from the collapse phases of the components A and B.
Chemical differentiation was spatially observed on small scales towards IRAS16293 as well as other low-mass protostars (e.g., Sakai et al. 2014b;Oya et al. 2016;Okoda et al. 2018). These variations were interpreted as a change of chemistry at the centrifugal barrier separating the infalling and rotating envelope from the disk. In our simulations, molecules also show different spatial distributions (see Section 4.1). These spatial changes in gas-phase abundances are explained by the increase of temperature and density towards the center. Similarities are seen with these observational studies. CH 3 OH, HCOOCH 3 , and SO are concentrated near the protostar in both cases (Sakai et al. 2014a(Sakai et al. ,b, 2016Oya et al. 2016Oya et al. , 2018, while CCH and c-C 3 H 2 come from a more distant and colder region: the spiral arms in our simulations and the infalling and rotating envelope in the observational studies by Sakai et al. (2014a,b) and Okoda et al. (2018). Only sulfur-bearing species, CS, H 2 CS, and OCS, show a different behavior. CS and OCS are observed in the infalling and rotating envelope (Sakai et al. 2014a;Oya et al. 2016Oya et al. , 2018, while H 2 CS is seen both in the envelope and inside the centrifugal barrier . In our simulations, these species are only present in the inner regions. The absence of OCS and CS in the inner regions compared to our simulations suggests missing destruction reactions at high temperature in our network. Their presence in more extended and colder regions in the observations, opposite to our predictions, is certainly due to a difference in the physical conditions. Indeed, in our simulations, we only modeled the disk, which shows very high densities (> 10 9 cm −3 ) and low temperatures. In such conditions, these three species are consequently abundantly depleted on the grains. In the envelope, which we did not model in this study, the densities of the particles are lower and these molecules are consequently less depleted onto the grains. This was shown by Vidal & Wakelam (2018), for example, who modeled the chemistry in the envelope of hot cores and hot corinos and found that the abundances of these three sulfur-bearing species are relatively high in the envelope.   Manigand et al. (2020b). The error bars correspond to an uncertainty of 30% and 10% for the components A and B, respectively. For the simulations, the bars cover the range between the minimum and maximum abundance ratios.

Conclusions
In this study, we performed 3D MHD simulations of a dense core collapse to investigate the chemical evolution that occurs from the cold core phase to the formation of a young rotationally supported disk. Among the 10 6 tracer particles introduced in the simulations, ∼1.5 × 10 4 end in the disk. We ran the three-phase gas-grain chemistry code Nautilus to follow the evolution of all the disk particles. Two different sets of initial abundances were considered, taken from SPH simulations of formation of dense cores and carefully chosen to be representative of two prestellar cores with different histories. Comparisons were made between the two final distributions of abundances in the disk as well as with their respective initial abundances in the cold and dense core. Our main findings can be summarized as follows: • The molecules were classified into different categories according to their spatial distribution in the disk. The initial abundances can affect the level of abundances of molecules, but have no significant impact on the distribution of most molecules, apart for HCO, HNO, OH, SO, SO 2 , and CH 3 . The spatial distribution of the molecules reflects their sensitivity to the temperature distribution. Radicals are usually destroyed in the warm inner region, while molecules such as water and complex organics are present in all regions of the disk, in ices in the cold outer regions, and in the gas phase in the warm inner regions.
• The main carriers of the different chemical elements at the final time of the simulations are usually the same as the initial ones. Significant changes are only seen for phosphorus and sulfur. For one set of initial abundances, atomic P decreases in abundance while the abundances of PO, PN, HCP, and CP increase. For both cases of initial abundances, the main contribution of HS is replaced by that of H 2 S 3 . The contribution of radicals such as HCO and CH 2 OH to the oxygen and carbon budgets is also smaller with time. • Despite similar carriers at the initial and final times, tens of molecules show a significant change in abundances during the formation of the disk. In particular, the abundances of "large" COMs such as CH 3 CHO, CH 3 NH 2 , CH 3 OCH 3 , and HCOOCH 3 increase significantly. Despite being present in cold and dense cores, their highest abundances are only reached during the collapse. In contrast, small molecules, especially radicals, see their abundances decrease (CCH, C 4 H, CH 3 , CN, NO, NH, NH 2 , OH, HCO, CH 3 O, O 2 ). This confirms that the collapse with its warm-up phase promotes molecular complexity. The abundances of some complex species (CH 3 CHO, CH 3 OCH 3 , H 2 CCO, HCOOCH 3 ) seem to mainly depend on the physical evolution during the collapse, as their final abundances are similar while the initial abundances differ significantly. • A high number of species show a similar abundance in the dense core and in the disk. These species are expected to mainly form in the molecular cloud or prestellar core, and include molecules such as H 2 O, H 2 CO, HNCO, and "small" COMs (CH 3 OH, CH 3 CN, NH 2 CHO) and their possible precursors. It is interesting to note that these COMs are also the ones presenting the lowest deuteration towards the Class 0 protostar IRAS 16293-2422 Jørgensen et al. 2018;Persson et al. 2018;Calcutt et al. 2018). • In spite of similar initial abundances, some species show different final abundances (CH 3 O, HCO, OH, HCO + , PCH + 3 ). These species are consequently sensitive to the initial abundances of other species.
• The MHD resistivities we obtained from the full gas-grain chemistry are in qualitative agreement with the one we used to estimate the ambipolar diffusion resistivity in the 3D simulations.
In conclusion, the chemical content of prestellar cores has an impact on the final composition of disks, even if it does not affect all species. For some species, the physical evolution during the collapse plays a major role, and for other species both the collapse and the initial abundances of the cold core contribute to the final abundances. We classified the species according to their behavior. This classification should help observers to interpret the origin of molecules observed in young disks. The results of this study should also be of interest for the chemical modelers of more evolved disks. Indeed, the initial abundances of standalone chemical models have to be chosen carefully. Using initial abundances within the range of abundances we obtain here for a young disk should make such models more realistic. Table A.1: Initial and final mean, median, minimum, and maximum abundances (a(b) = a × 10 b ) of the species presented in this study for the two sets of initial abundances of all the computed disk particles. The abundances correspond to the total ice + gas abundances.