A numerical approach to model chemistry of complex organic molecules in a protoplanetary disk

Multiphase astrochemical modeling presents a numerical challenge especially for the simulation of objects with the wide range of physical parameters such as protoplanetary disks. We demonstrate an implementation of the analytical Jacobian for the numerical integration of the system of differential rate equations that govern chemical evolution in star-forming regions. The analytical Jacobian allowed us to greatly improve the stability of the code in protoplanetary disk conditions. We utilize the MONACO code to study the evolution of abundances of chemical species in protoplanetary disks. The chemical model includes 670 species and 6,015 reactions in the gas phase and on interstellar grains. The specific feature of the utilized chemical model is the inclusion of low-temperature chemical processes leading to the formation of complex organic molecules (COMs), included previously in the models of chemistry of COMs in prestellar clouds. To test the impact of analytical Jacobian on the stability of numerical simulations of chemical evolution in protoplanetary disks, we calculated the chemical composition of the disk using a twophase model and four variants of the chemical reaction network, three values of the surface diffusion rates, and two types of the initial chemical composition. We also show a preliminary implementation of the analytical Jacobian to a three-phase model.


Introduction
Studies of the chemical composition of objects in the interstellar medium, especially for the content of complex organic molecules (COM), is an important prerequisite for understanding the origin of life in the Universe. Protoplanetary disks are dust-and gas-rich objects that could possibly form planetary systems. They are ubiquitous around young low-mass stars (e.g., Manara et al. 2016, Kim et al. 2017. A study of the chemical composition in the disks around Sun-like stars will provide an idea of the origin of organic molecules in the early Solar System, which in turn can serve as a key to understanding the early chemical composition of the Earth and other planets. To date, complex organic molecules (which are defined to have six or more atoms, including carbon, Herbst and van Dishoeck 2009) Favre et al. 2018) molecules were also found in the disks. Other molecules representative of the COMs content in the ISM, such as methyl formate HCOOCH 3 and dimethyl ether CH OCH 3 3 , are still not found in the Class II protoplanetary disks, but are widely observed in star-forming regions representing earlier stages of low-mass protostellar evolution, such as prestellar cores (Jiménez-Serra et al. 2016) and hot cores/corinos (Jorgensen andPILS Team 2020, Manigand 2020) and FU Ori type young stars (Lee et al. 2019). Aikawa et al. (1997) were among the first to study the evolution of the molecular composition of protoplanetary disks. They considered stationary minimum mass solar nebula (MMSN) without radial mixing; density and temperature did not change over time. Their chemical model included gas-phase reactions, adsorption onto dust grains, and thermal desorption from dust particles. Ionization and dissociation by interstellar and stellar ultraviolet radiation  were neglected. The chemical network of reactions was based on the UMIST94 database (Millar et al. 1991).
Over the next 20 years, protoplanetary disk models became more sophisticated (e.g., see review by Henning and Semenov 2013). However, the applied chemical models remained mostly two phase, that is, only gas-grain interactions were considered. Ruaud and Gorti (2019) were able to apply the three-phase chemical model to protoplanetary disks for the first time.
In this article, for the first time, we apply a scenario of the formation of complex organic molecules in cold gas of prestellar cores proposed by Vasyunin and Herbst (2013) and further developed by Vasyunin et al. (2017) to a protoplanetary disk around a Sun-like star. The evolution of the chemical composition was calculated for 1 Myr assuming that the disk structure is in a quasi-stationary mode for a given time period (Akimkin et al. 2013).
To numerically solve the system of differential equations that determine chemical evolution, the three-phase MONACO code uses the DVODE integrator (Brown et al. 1989). In the current state, the application of the MONACO code to protoplanetary disks is challenging due to a wide range of physical conditions typical for disks. Numerical integration of a system of ordinary differential equations requires the calculation of the Jacoian matrix of the system. The DVODE can work in two regimes: with internally generated numerical Jacobian and with user-supplied analytical Jacobian. The latter option typically results in much higher numerical stability of integration. On the other hand, it requires additional efforts from researcher aimed at derivation and implementation of the analytical expressions for the Jacobian matrix into the numerical code. To solve this problem, we added to the code the implementation of specifying the analytical Jacobian of the system of differential equations.
In this study, we set the following goals: by supplying the analytical Jacobian, to increase the stability of the MONACO code for efficiently calculating the evolution of the chemical composition of the protoplanetary disk under the wide range of physical parameters and conditions typical of protoplanetary disks. Also, we aim to study the formation of COMs in the disk, especially midplane, using the model suggested by Vasyunin and Herbst (2013) tested on prestellar cores, the conditions that are close to the conditions in midplane.

Models
2.1 Physical model of the protoplanetary disk As a physical model of a protoplanetary disk (PPD), we used the model presented by Molyarova et al. (2017). This model is the PPD model around a T Tauri type star with the mass of M 1 ⊙ and utilizes the 1+1D approach to calculate disk density and temperature. The disk is considered quasi-stationary, axisymmetric, and hydrostatic in the vertical direction. The protoplanetary disk model used in this article is a grid of 4,400 points (55 radial and 80 vertical points). The model parameters are presented in Table 1.
Dust temperature in the upper disk is calculated using multifrequency ray tracing (RT) procedure for the stellar and background radiation similar to the study by Molyarova et al. (2018). RT is done in 2D in r φ , ( )-plane and includes four directions (to and from the central star, upward and downward relative to the disk midplane). The corresponding angle-averaged intensity is used to calculate radiation field strength and rates of photoreactions. We assume M 1 ⊙ for the stellar mass, stellar effective temperature (T * ), and radius (R * ) are taken from the evolutionary tracks in the study by Baraffe et al. (2015).
In Figure 1, the distributions of disk grid points (top left), the strength of UV radiation (top right), gas temperature (bottom left), and gas density (bottom right) in the disk (at each grid point) are also presented as a function of radial distance from the central star (r) and vertical distance from disk midplane z r ( ) .

Chemical model
In this study, we utilized a chemical model with the network of gas phase and surface chemical reactions, which was used in the study by Vasyunin and Herbst (2013) with  The chemical reaction network used in this model contains 670 species and 6,015 gas-phase and surface reactions, as well as 198 species and 880 reactions in the ice mantle of dust particles depending on the simulation mode (see details in Section 2.2.2). Following Vasyunin et al. (2017), we utilize five types of desorptions in the model: thermal evaporation, photodesorption, desorption by cosmic ray particles (CRP), CRP-driven photodesorption, and reactive desorption. We do not consider CRP attenuation inside the disk. Cosmic-ray ionization rate in our model is ζ 1.3 10 s CR 17 1 = × − − (Glassgold and Langer 1974). Ionization by short-living radionuclides is ζ 6.5 10 s RN 19 1 = × − − (Umebayashi and Nakano 2009). We also consider thermal hopping across the grain surface for species and quantum tunneling through potential barriers for light species (atomic and molecular hydrogen), depending on the simulation mode (see details in Section 2.2.2).
In Table 2, the atomic initial fractional abundances of elements with respect to the total number of hydrogen nuclei used in the model are presented (according to Wakelam and Herbst 2008). The molecular initial composition is calculated as a result of the chemical evolution of a cold dark cloud at 10 6 years by MONACO with next parameters: density of gas n 10 cm We modified the chemical model to take into account the radiation fields from the central star and interstellar radiation according to the protoplanetary disk model. The rate constants of photoionization reactions K i are calculated as follows: where α is the photoreaction rate in the field of unshielded ultraviolet radiation; G UV and A V are factors characterizing the radiation field and attenuation of the radiation field, respectively (they are the parameters of the protoplanetary disk model, G UV is utilized in units of the Draine's field Draine 1978), and γ is a parameter that takes into account the increased field attenuation in the ultraviolet range compared to the visible one (McElroy et al. 2013).
The grid points of the protoplanetary disk model are not independent in terms of calculating the evolution of the chemical composition because the self-shielding of the molecules H 2 and CO has to be taken into account when calculating the photodissociation of those molecules. In our case of 1+1D physical disk model, chemistry in each radial column of the model grid must be calculated in a specific order, namely, starting from the outermost point in disk atmosphere and toward the disk midplane. This is needed because one has to calculate self-shielding factors for CO and H 2 molecules along the radiation propagation path, that is, 1+1D model is in the vertical direction at each radial column, and for that, it is necessary to know the corresponding vertical column densities of the molecules. H 2 and CO self-shielding are calculated based on the study by Visser et al. (2009).
For the top points in each radial column (points with a maximum height above the disk midplane), the column densities H 2 and CO (N 80 col ) are assumed to be zero. The column densities at the other points N i col are calculated as follows: In the case of a three-phase model, the evolution of the molecular composition is determined by the following system of differential equations change in the abundances in the gas phase, on the surface of dust grains and in the mantle of dust particles, respectively. In the case of a two-phase model, this system of equations is noticeably simplified: the terms of the equations describing abundances in the mantle and the transition of molecules between the mantle and the surface are not considered. For the numerical integration of the system of differential rate Eq. (3), the Adams method is used, implemented in the DVODE integrator. The integration time is 1 Myr.

Supplying the analytical Jacobian
When using a two-phase model (gas-grain) and a chemical network containing~660 reactants, including~200 The fractional abundances are given with respect to the total number of hydrogen nuclei. species on the surface and~6,000 chemical reactions, the Jacobian contains~435,000 elements, among which 6,000 are nonzero, which is about 1.4%. In the case of a three-phase model (gas-surface-mantle), the Jacobian con-tains~739,600 elements, and the number of nonzero terms is rather difficult to determine in advance, but it is very large, much more than in the case of a two-phase model.
The two-phase Jacobian is well described analytically with the exception of reactive desorption (R i des term in Eq. (3)) that violates this harmony of the analytical Jacobian. In different models, we have considered three types of reactive desorption. The first type is based on the probability that a molecule formed as a result of a reaction on the grain surface has energy exceeding its binding energy and on the assumption of a rapid loss of this energy, for example, the transfer of energy to a dust particle more massive than the molecule itself (Garrod et al. 2007). This type of desorption is used for reactions with one product, which adds~700 more nonzero elements of the Jacobian ( 0.16% + ). In the second type of reactive desorption, a fixed value of the desorption efficiency is considered, but for all products of surface reactions (Vasyunin and Herbst 2013), which adds~1,000 nonzero Jacobian elements ( 0.23% + ). The third type of reactive desorption takes into account both the features of the reaction and the properties of the dust grain surface, namely, the effective mass of the surface element from which the molecule is desorbed into the gas phase (Minissale et al. 2016). This most complex type of reactive desorption considered in our model adds~25,000 nonzero elements to the Jacobian ( 5.7% + ). The two-phase Jacobian, modified to take into account reactive desorption, is still sparse, but not as simple.
To obtain the symbolic Jacobian of such a system of ordinary differential equations, the SymPy symbolic computation package for the Python language was chosen. The process of numerically solving the system of equations is as follows. On the basis of the chemical reaction network, a system of differential equations in Fortran is formed, since the numerical solution is performed by means of this language. A system of differential equations in Python is also formed, but in a symbolic form. Then symbolic partial derivatives are calculated using SymPy. Then symbolic expressions are simplified. After that, the Python code forms the Fortran source code, and the system of differential equations is solved numerically.

Simulation setups
So far we have applied the following simulation setups: -Two phase (gas-grain) with reactive desorption taken into account, without specifying the Jacobian; -Two phase with reactive desorption taken into account, with specifying the Jacobian; -Three phase (gas-surface-mantle) with reactive desorption taken into account, without specifying the Jacobian; -"Incomplete" three phase with reactive desorption taken into account, with specifying the Jacobian.
In two-phase approaches, no distinction is made between the surface and the ice mantle of dust particles, which is a noticeable simplification in comparison with the structure of real multilayer ice. In threephase approaches, we allow surface molecules to bury deep into the ice mantle. Thus, such molecules become inaccessible for direct desorption into the gas phase. In the three-phase approach, the change of the mantle abundances occurs due to adsorption to the surface ("tran" terms in Eq. (3)), chemical reactions ("chem" terms), and physical diffusion between the inert bulk and the surface ("diff" terms). In the "incomplete" three-phase approach, physical diffusion is not implemented yet. In all setups, we used the type of reactive desorption according to Minissale et al. (2016). Each simulation approach we tested using four different chemical network options: network used in the study by Vasyunin and Herbst (2013); network with refined binding energies for some species; network that additionally includes reactant CH OCH 3 2 and chemical reactions with its participation; updated network used in the study by Vasyunin et al. (2017). Also, for each simulation approach, we used three values of the diffusion/desorption surface ratio: (with tunneling for light species) and 0.5 and 0.8 (with no tunneling), as well as two variants of the initial chemical composition: atomic and molecular.

Two-phase setups
In total, in each two-phase approach (with and without Jacobian), we calculated 24 models of protoplanetary disks (4 variants of chemical networks, 3 different values of the surface chemistry parameter, and 2 variants of the initial chemical composition). Thus, 105,600 runs of the numerical integration of the differential equations system were performed in each of the two-phase modes under different physical conditions and different parameters of chemistry.
When using the Jacobian, the calculation speed increased by four to five times on average. The number of unsuccessfully completed calculations has decreased significantly, namely, by 150 times. Now that number is eight unsuccessful calculations out of 105,600 runs. Thus, we consider it reasonable to use the symbolic Jacobian in two-phase models with reactive desorption.

Distribution of organic molecules in the disk
In this section, we present two-dimensional distributions of the fractional abundances of selected organic molecules and some chemically related species with respect to the total number of hydrogen nuclei in the disk in the gas phase and on the surface of dust particles at time sition. In all figures in this section, the abscissa shows the radial distance in logarithmic scale, and the ordinate shows the ratio of vertical and radial distances from midplane and central star, respectively. In those coordinates, the grid of model points is uniform spatially (see Figure 1(a)). The color scale displays the decimal logarithm of the fractional abundances of the molecule. The g prefix in front of the name of the molecule denotes molecule on the grain surface. Species without "g", in contrast, reside in the gas phase. To plot the two-dimensional distribution, only the fractional abundances of the molecules greater than 10 15 − were taken into account. As shown in Figure 2(a) and (b), the relative hydrogen abundance of the hydroxyl group OH reaches~10 7 − in the gas phase in the outer region of the disk and~10 9 − on the grain surfaces in the outer disk. Figures 2(c)-3(f) shows that organic molecules such as formaldehyde H CO 2 , methanol CH OH 3 , methyl formate HCOOCH 3 , dimethyl ether CH OCH 3 3 , and formic acid HCOOH, are more abundant on grain surfaces than in the gas phase in the midplane. The maximum fractional abundances of these molecules on the grain surfaces are in the midplane in intermediate and outer regions of the disk and reach values from~10 8 − to~10 5 − , whereas the maximum fractional abundances of these molecules in the gas phase are several orders of magnitude lower (from~10 13 − to~10 7 − ). The maximum fractional abundance on the grain surfaces of H CO 2 is reached within the radial distances from 60 to~200 au (see Figure 2( In our model, COMs are more efficiently formed in the midplane on grains. At 5-30 K and weak UV field, due to attenuation in the depth of the disk, CO in sufficiently large quantities can freeze out to grains. At high altitudes, CO is destroyed under the stronger UV field (CO C O photon → + ) and COMs are no longer formed in such quantities as in the midplane. However, as shown in Figure  2(d) and (f), in the midplane at 300-400 au, there is a region with a reduced abundances of gH CO 2 and gCH OH 3 on grain surfaces. We attribute this to a slight increase in the UV field in this region in the physical model of the protoplanetary disk (see Figure 1(b)).

Total column densities of organic molecules
We also present the column densities profiles for selected molecules as a function of the disk radius (total number of molecules in the column) in cm 2 ( ) − . In Figure 4, total column densities for molecules in the gas phase with  Figure 4(c) shows that in our model gas-phase methanol reaches the same value between~0.6 and~1.5 au for atomic initials.
In the work by Walsh et al. (2014), the model column densities of H CO 2 and CH OH 3 were estimated as 10~10 cm 12 13 2 − − (see Figure 8 in the study by Walsh et al. 2014). We understand that our profiles of column densities do not reproduce the observational data and previous models completely; however, we believe that there are sections of the profiles that satisfy the existing data. We do not pay detailed attention to the reasons for the peaks or troughs of our profiles since the priority for us is the modification of the three-phase model for stable calculations of the protoplanetary disk. We started with our two-phase model to work out the technicalities of the Jacobian assignment and are aiming to apply the Jacobian to our full three-phase model. Note that among the possible reasons for such results may indeed be a low value of the disk mass. In our model, we use the value of 0.01 M ⊙ . While estimates of disk mass, for example, TW Hydra, referred to by Walsh et al. (2016), have an upper limit of 0.1 M ⊙ and a best-fit of 0.023 M ⊙ .
As shown in Figure 4, the use of different variants of the initial chemical compositions has a significant impact on the simulation results, especially with regard to complex organic molecules. The approach in which the protoplanetary disk inherits its initial composition from the previous evolutionary stages of the protostar is more fair.

Three-phase setups
Work on adding the Jacobian to the three-phase model is ongoing. Currently, we have implemented the supplying of analytical Jacobian for the processes of instantaneous surface redefinition due to the adsorption/desorption of molecules from/to the gas phase from the surface. This added another 230,104 (31%) non-zero Jacobian elements. Such Jacobian obviously ceases to be sparse.
The use of the Jacobian of the differential equations system significantly affected the operation of the threephase mode. In the three-phase mode without the Jacobian, the calculations for most of the grid points of the protoplanetary disk physical model do not complete successfully in a reasonable CPU time. The addition of the Jacobian to the "incomplete" three-phase regime significantly improved the stability and speed of calculations, but still does not allow us to calculate the chemical evolution in the inner disk. At the current moment in this mode, we were able to calculate the chemical composition for the outer regions of the disk with the radial distance from the central star of more than 25 au. In the three-phase model, the molecules from the deep of ice mantle cannot be directly desorbed into the gas phase. "Incomplete" three-phase mode is an intermediate step toward "full" three-phase mode. In the "incomplete" mode, there is no physical transfer of molecules between the mantle and the surface. Mantle molecules can become surface molecules only by emptying the upper layers of the surface. Therefore, the results obtained in this mode must be treated with caution. Nevertheless, using the example of methanol ( Figure 5), it can be seen that the abundances in the gas phase and on the grain surfaces in the area accessible to us (especially in midplane) have become several orders of magnitude lower compared to the two-phase regime (Figure 2(e) and (f)). Obviously, this is achieved due to the settling of species deep into the bulk. Hence, the presence of the third solid phase in the chemical model has a significant effect on the abundance of molecules in the gas phase.

Summary
We have implemented an option to supply the analytical Jacobian in the MONACO code for the chemical evolution of interstellar objects, which was previously applied to cold dark clouds and applied it for the first time to the physical model of the protoplanetary disk around a Sunlike star. Specifying the Jacobian for the two-phase (gasgrain) model made it possible to reduce the time of numerical integration of the differential equations system, as well as to increase the stability and accuracy of calculations as applied to the protoplanetary disk.
We also present preliminary results on adding the Jacobian to the three-phase (gas-surface-bulk) model, the intermediate results of which also demonstrated the justification and the necessity of using the Jacobian in modeling the formation of complex organic molecules in protoplanetary disks using the MONACO code. This update is crucial to allow numerically effective modeling of threephase chemistry in protoplanetary disk conditions.
Funding information: MYK and AIV acknowledge the support Ministry of Science and Higher Education of the Russian Federation via the State Assignment contract FEUZ-2020-0038. VVA is grateful to the Foundation for the Advancement of Theoretical Physics and Mathematics "BASIS" for financial support (20-1-2-20-1).
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.

Conflict of interest:
The authors state no conflict of interest.
Data availability statement: The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.