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

2.2.1 Supplying the analytical Jacobian

When using a two-phase model (gas–grain) and a chemical network containing ∼ 660 reactants, including ∼ 200 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 contains ∼ 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.

2.2.2 Simulation setups

So far we have applied the following simulation setups:

1. Two phase (gas–grain) with reactive desorption taken into account, without specifying the Jacobian;

2. Two phase with reactive desorption taken into account, with specifying the Jacobian;

3. Three phase (gas–surface–mantle) with reactive desorption taken into account, without specifying the Jacobian;

4. “Incomplete” three phase with reactive desorption taken into account, with specifying the Jacobian.

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 3 OCH 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: E B E D = 0.3 (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.

3.1.1 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 t final = 1 0 6 years with using the network from Vasyunin et al. (2017), E B E D = 0.3 and atomic initial chemical composition. 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 1 0 − 15 were taken into account.

As shown in Figure 2(a) and (b), the relative hydrogen abundance of the hydroxyl group OH reaches ∼ 1 0 − 7 in the gas phase in the outer region of the disk and ∼ 1 0 − 9 on the grain surfaces in the outer disk. Figures 2(c)–3(f) shows that organic molecules such as formaldehyde H 2 CO , methanol CH 3 OH , methyl formate HCOOCH 3 , dimethyl ether CH 3 OCH 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 ∼ 1 0 − 8 to ∼ 1 0 − 5 , whereas the maximum fractional abundances of these molecules in the gas phase are several orders of magnitude lower (from ∼ 1 0 − 13 to ∼ 1 0 − 7 ).

Figure 2

The two-dimensional distribution of the fractional abundance of OH (top row), H 2 CO (middle row), and CH 3 OH (bottom row) with respect to hydrogen in the disk in the gas phase (left column) and on the grain surfaces (right column) at the time t final = 1 0 6 years color scale is lg n mol n H + 2 n H 2 .

The maximum fractional abundance on the grain surfaces of H 2 CO is reached within the radial distances from ∼ 60 to ∼ 200 au (see Figure 2(d)), CH 3 OH – from ∼ 500 to ∼ 600 au (Figure 2(f)), HCOOCH 3 – from ∼ 40 to ∼ 50 au (Figure 3(b)), CH 3 OCH 3 – from ∼ 20 to ∼ 200 au (Figure 3(d)), and HCOOH – from ∼ 80 to ∼ 100 au (Figure 3(f)).

Figure 3

The two-dimensional distribution of the fractional abundance of HCOOCH 3 (top row), CH 3 OCH 3 (middle row) and HCOOH (bottom row) with respect to hydrogen in the disk in the gas phase (left column), and on the grain surfaces (right column) at the time t final = 1 0 6 years color scale is lg n mol n H + 2 n H 2 .

As follows from Figures 2(c)–3(f), the values of the maximum fractional abundances on the grain surfaces of H 2 CO and HCOOH are ∼ 1 0 − 6 , CH 3 OH – ∼ 1 0 − 5 , HCOOCH 3 – ∼ 1 0 − 7 , and CH 3 OCH 3 – ∼ 1 0 − 8 . In the gas phase, the maximum values of abundances H 2 CO and CH 3 OH are ∼ 1 0 − 7 , CH 3 OCH 3 and HCOOCH 3 – ∼ 1 0 − 13 , and HCOOH – ∼ 1 0 − 9 . These results are consistent with the results obtained in the study by Walsh et al. (2014), where complex organic molecules formed more efficiently on the surface of dust particles in the middle plane of the disk and also reached maximum values from ∼ 1 0 − 6 to ∼ 1 0 − 5 (see Figures 6 and 7 in the study by Walsh et al. 2014).

The main channel for the formation of COMs, in particular methanol, in this model consists of the sequential hydrogenation of CO on grains:

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 → photon C + O ) 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 2 CO and gCH 3 OH 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)).

3.1.2 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 using the network from Vasyunin et al. (2017), E B E D = 0.3 and atomic (blue line) and molecular (orange line) initial chemical composition are presented.

Figure 4

Total column densities in ( cm − 2 ) of (a) OH, (b) H 2 CO , (c) CH 3 OH , (d) HCOOCH 3 , (e) CH 3 OCH 3 , and (f) HCOOH in the gas phase as a function of disk radius at the time t final = 1 0 6 years for atomic and molecular initial chemical composition.

In the work by Podio et al. (2019), estimates of the column densities for formaldehyde H 2 CO molecule in the protoplanetary disk of DG Tauri were obtained. Their estimates are at the level of ∼ 1 0 14 cm − 2 . As shown in Figure 4(b), in our model, the total column densities for H 2 CO in the gas reach values ∼ 6 × 1 0 13 cm − 2 for atomic initial composition and ∼ 7 × 1 0 14 cm − 2 for molecular initial composition in the inner disk at ∼ 3 au. Walsh et al. (2016) estimate that the peak of the column density of methanol CH 3 OH molecule in the protoplanetary disk of TW Hya is ∼ 1 0 12 cm − 2 at 30 au. 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 2 CO and CH 3 OH were estimated as ∼ 1 0 12 − ∼ 1 0 13 cm − 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 ⊙ (Kama et al. 2016).

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.