© R. Guadarrama et al. 2022

1 Introduction

Protoplanetary disks form during the gravitational collapse of rotating cloud cores and consist of molecular hydrogen, dust grains of various sizes, and chemical species of varying complexity. While molecular hydrogen constitutes the main mass reservoir of the disk, dust plays an important role in determining the disk’s thermal balance and supplying the planet formation process with building material. Chemical species, in particular, water and carbon-based molecules, are the main prerequisites of life and also play a key role in the process of dust growth (e.g., Okuzumi et al. 2016). According to Qi et al. (2013), the fronts where volatile species condense (snowlines) play an important role in the planetary formation process. They define the location where dust coagulation is favoured because of an increase in the stickiness of icy grains.

The position of the snowlines for the most abundant volatiles may also have an influence on the elemental and molecular composition of planetesimals and planets. The position of the respective snowlines for different chemical species determine the radii at which molecules freeze onto dust grains. This influences the chemical composition of icy mantels of dust grains at specific regions in the disk. As the dust grains are the seeds of planet formation, the chemical composition of their icy mantels also influences the chemical composition of forming planetesimals, and, subsequently, planets. In the case of gas giants, the composition of the atmosphere is also linked to the position of the different snowlines and, thus, to their location on the disk. For instance, the carbon to oxygen ratio, a defining characteristic of the atmosphere, is set by the amount of gas phase oxygen and carbon a gas giant accretes from the disk. The accreted amount of oxygen and carbon-bearing molecules in the gas phase is dependent on the location on the disk where they freeze out (snowline). Therefore, planetary formation efficiencies as well as their core and atmosphere compositions are strongly linked to the snowlines of volatile species.

Until recently, most numerical studies of protoplanetary disks were limited to the solar metallicity case. However, rising numbers of exoplanet detections have revealed that the frequency of the occurrence of giant planets depends on the metallicity of the host star (Santos et al. 2005; Benz et al. 2014). In addition, studies of protostars and pre-main-sequence stars in the Magellanic Clouds indicate that mass accretion rates are inversely proportional to metallicity (De Marchi et al. 2013). This finding seems to be corroborated by Yasui et al. (2010, 2016), who inferred shorter lifetimes to low-metallicity disks in the outer regions of the Galaxy. However, this trend could also be explained by higher rates of disk photoevaporation due to a more diluted medium that allows a stronger radiation field (Nakatani et al. 2018). Another example of these studies is the work performed by Berg et al. (2016). By observing the chemical abundances of HII regions in 12 nearby dwarf galaxies, the study suggests that there are different chemistry cycles in non-solar metallicity environments. This evidence suggests that disk evolution and planet formation in protoplanetary disks may proceed differently in low-metallicity environments.

The chemical evolution and composition of protoplanetary disks have been extensively studied in the context of solar-metallicity disks (see a review by Pontoppidan et al. 2014; Henning & Semenov 2013). Among various compounds found in protoplanetary disks, volatile species with low sublimation temperatures are optimal for studies of this kind because they are central to the process of planet formation. In this work, we use the radiation thermo-chemical code PRODIMO (Woitke et al. 2009) to examine the impact of low metallicity on important molecules (e.g., H₂O, CO) on protoplanetary disks.

When we refer to metallicity, we mean the values of element abundances and the gas-to-dust ratio. To simulate a lower metallicity, we reduce the element abundance of metals together with the amount of dust. Then we analyze the resulting temperature structure, radiation field, and chemical abundances of various molecules in different metallicities. The formation and destruction reactions of chemical species are also analyzed. Except for ionized molecules, we track the position of the respective snowlines. For a deeper analysis of the chemical abundances, we performed a comparison of the vertical column density between the models with lower metallicity and an artificially scaled-down reference model. Although it is expected that the behaviour of the column densities as a function of the metallicity may not be linear, we use the comparison with scaled-down models to be able to define the regions of the disk where the discrepancy is largest and to quantify the effect of non-linearity. We use an evolved, externally irradiated Class II disk for our models.

In Sect. 2, we describe the gas and dust structure of the disk in the models as well as the chemical network and element abundances used for this work. We also explain how we applied the modifications in the disk models to simulate lower metallicities in the disks. Our results are presented in Sect. 3. In Sect. 4, we discuss some aspects of the used model and the results and in Sect. 5, we present our conclusions and describe our next steps.

2 Method

We employed the radiation thermo-chemical code PRODIMO (Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011) to model the chemical composition of young protoplanetary disks with distinct metallicities. PRODIMO includes the radiative transfer of stellar and background radiation, along with the gas and dust thermal balance to provide the temperature structure and the local radiation field of the disk. The chemical abundances are then calculated using steady-state chemistry employing a chemical reaction network used in Kamp et al. (2017). Various heating and cooling processes are also considered to solve consistently for the gas temperature and the chemistry. The implementation of X-rays and X-ray chemistry is taken from Aresu et al. (2011).

2.1 Metallicities and initial abundances

We used different element abundances and dust-to-gas ratios to study the impact of lower metallicities on the chemical composition of the disk, focusing on water and several most abundant carbon-bearing molecules. We considered three model realizations of protoplanetary disks with three distinct metallicities, Z. In all cases, the mass of the central star, stellar effective temperature, stellar luminosity, and X-ray luminosity were fixed and set equal to M_* = 0.7 M_⊙, T_* = 4000 K, L_* = 1.0 L_⊙, and L_X = 2.6 × 10⁻⁴ L_⊙. We assumed a solar metallicity value of Z_⊙ = 0.02, which is in agreement with the value of Z_⊙ = 0.0196 ± 0.0014 from von Steiger & Zurbuchen (2016). The abundances for our reference model are the abundances used for standard TTauri disk in the ProDiMo code (Woitke et al. 2016). As a result, our reference model Z1 has a metallicity value of Z1 = 0.017 and it therefore displays a value of Z1 = 0.85 Z_⊙. The element abundances corresponding to model Z1 are shown in Table 1. The models are named according to their disk metallicity. Model Z01 and Z001 correspond to a metallicity of 0.1 and 0.01 (10% and 1% of the reference metallicity Z1, respectively). Thus, the metallicities of our models lie within the 0.85 × 10⁻²−0.85 Z_⊙ limits. In order to produce models with different metallicities, we performed two modifications to the initial setup of the models: (1) we reduced the initial abundances of 13 of the 15 elements (all except H and He) by a factor of 10 and 100 for the gas and the dust and (2) we reduced the dust-to-gas mass ratio in the input parameters by the same factors (models Z01 and Z001, respectively).

This process resulted in the above-mentioned three models (from herein referred to as Z1, Z01, and Z001) characterized by the reference, 1/10, and 1/100 of reference metallicity, respectively. The disk mass, the metallicity relative to the solar value, and the dust-to-gas ratio are displayed in Table 2.

Further variations of the disk and star properties with metallicity will be considered in follow-up studies. The grid resolution is 150 points in the horizontal direction and 100 in the vertical direction. The details of the initial disk’s physical structure, dust content, and chemical composition are provided below.

2.2 The disk model

The disk model we used for this study assumes certain simplifications in its physical properties such as the star, the disk geometry, dust settling, and so on. The set of assumptions was taken to reproduce the most commonly observed multi-wavelength properties of Class II protoplanetary disks. The basic structural properties of the model are described in detail in Woitke et al. (2016). The chemical network used in this study is introduced and explained in Kamp et al. (2017). In Sects. 2.2.1, 2.2.2, and 2.2.3, we provide a brief overview of this model.

2.2.1 Gas disk structure

We used a fixed parameterized density structure for the disk based on viscous evolution models. We refer to Woitke et al. (2016) for details of this step. The gas density structure is an axisymmetric flared (2D) function of the radius and the vertical height of the disk (r and h, respectively) and it is given by: $$\rho \left(r,z\right)=\frac{\text{Σ}\left(r\right)}{\sqrt{2\pi }·h\left(r\right)}\exp \left(-\frac{z^2}{2h{\left(r\right)}^2}\right)\left[\text{g\hspace{0.17em}}{\text{cm}}^{-3}\right].$$(1)

Here, ∑(r) is the radial surface density profile of the disk. We assume a simple power-law distribution with a tapered outer edge $$\text{Σ}\left(r\right)={\text{Σ}}_0{\left(\frac{r}{R_{\text{in}}}\right)}^{-\lambda }\exp \left(-{\left(\frac{r}{R_{\text{tap}}}\right)}^{2-\lambda }\right)\left[\text{g\hspace{0.17em}\hspace{0.17em}}{\text{cm}}^{-2}\right].$$(2)

Here, R_in is the inner disk radius and R_tap is the characteristic radius. The constant ∑₀ is determined from the specified disk mass, M_disk = 2π ∫ ∑(r)r dr, and Eq. (2). The vertical scale height h(r) is given by a radial power law: $$h\left(r\right)=H_{100}{\left(\frac{r}{100\text{\hspace{0.17em}\hspace{0.17em}au}}\right)}^{\beta },$$(3)

where H₁₀₀ is the disk scale height at r = 100 au and has a value of 10 au and β is the flaring power index. Figure 1 shows the density structure and the radial surface density. The parameters for the disk density structure are listed in Table 3.

Fig. 1 Gas density structure (top) and radial surface density profile (bottom) of the reference model Z1.

2.2.2 Dust disk structure

For the fiducial model Z1 with the reference metallicity of Z1 = 0.85 Z_⊙, we assumed a dust-to-gas mass ratio of 0.01 in the disk. The respective gaseous element abundances are listed in Table 1 and are taken from Kamp et al. (2017). We take dust growth into account by using a dust size distribution with a minimum and maximum dust grain sizes of a_min = 0.05 µm and a_max = 3000 µm, respectively. We used a simple power-law for the dust size distribution f(a) ∝ a^−p with the canonical value for interstellar grains of p = 3.5 (Mathis et al. 1977). The dust composition consists of a mixture of 60% amorphous laboratory silicate, 15% amorphous carbon, and 25% vacuum. The vacuum represents the porosity of the dust grains. We used a distribution of hollow spheres with a maximum hollow volume ratio, $V_{\text{hollow}}^{\max }=0.8$. All the relevant dust properties are listed in Table 3 and are taken from Woitke et al. (2016).

In this work, we applied dust settling based on the method of Dubrulle et al. (1995). The dust composition and dust size distribution are constant throughout the entire disk. Furthermore, the dust composition, its maximum and minimum sizes, and the slope of the dust size distribution are assumed to be identical for all considered metallicities. This may be an oversimplification, but the process of dust growth in low-metallicity disks is still poorly known.

2.2.3 Chemical model

After modeling the radial and vertical physical structure of the disk, PRODIMO calculates a detailed continuum and line radiative transfer as well as the heating and cooling balance and the gas and surface chemistry of the disk. The net formation rate of a chemical species, i, is calculated on the basis of the following equation: $$\begin{array}{lll}\frac{\text{d}{n}_i}{\text{d}t}\hfill & =\hfill & {\displaystyle \sum _{jkl}{R}_{jk\to il}\left(T_g\right)n_j{n}_k+{\displaystyle \sum _{jl}\left(R_{j\to il}^{\text{ph}}+R_{j\to il}^{\text{cr}}\right)n_j+\cdots }}\hfill \\ \hfill & \hfill & -n_i\left({\displaystyle \sum _{jkl}{R}_{il\to jk}}{n}_l+{\displaystyle \sum _{jk}\left(R_{i\to jk}^{\text{ph}}+R_{i\to jk}^{\text{cr}}\right)+\cdots }\right),\hfill \end{array}$$(4)

Here, the first row represents all the formation reactions that produce species i and the second row all the destruction reactions that destroy species, i. Here, R_jk→il indicates the two-body gas-phase reaction rate between two reactants (j and k), forming two products (i and l), while $R_{j\to il}^{\text{ph}}$ stands for the photo-reaction rate. This rate depends on the local strength of the UV radiation field. $R_{j\to il}^{\text{cr}}$ indicates a cosmic ray-induced reaction rate. We only show three kinds of reactions in Eq. (4). The dots indicate that there are other kinds of reactions involved (freeze-out, X-ray chemistry, etc.). The technical details of the chemistry are fully described in Woitke et al. (2009). We assume kinetic chemical equilibrium (steady-state chemistry) for this work. This means we have $\frac{\text{d}{n}_i}{\text{d}t}=0$ for the left side of Eq. (4) for the steady-state model. We obtained i = 1… N_sp non-linear equations with j = 1… N_sp unknown particle densities n_j, where N_sp is the number of considered species. The chemistry of the disk is computed by solving Eq. (4) using a globally convergent Newton–Raphson method as a steady-state solver Woitke et al. (2009) The steady-state chemistry approach has also been used and compared with time-dependent chemistry in Woitke et al. (2016).

The chemical network used in this work is based on the one used in Kamp et al. (2017). We use additionally several surface chemistry reactions described in Thi et al. (2020, and in prep.). The surface chemistry network (Thi et al. 2020) is based on Hasegawa & Herbst (1993). The total amount of species is 250 and 5790 chemical reactions, respectively. The modifications to the element abundances done for this work will be explained in detail in the next section.

Fig. 2 Surface density in the midplane for dust (black) and gas (blue) as a function of radius. The models Z1 (solid line), Z01 (dashed line), and Z001 (dotted line) are shown.

Fig. 3 Spectral energy distribution of the different disk models and stellar spectrum. The models Z1 (solid line), Z01 (dashed line), and Z001 (dotted line) are shown. A lower metallicity leads to a decrease in the flux density for the SED.

2.3 Surface density

Figure 2 shows the resulting gas and dust surface density after the mentioned changes in the models have been applied. It shows a decrease in the dust surface density for lower metallicities. We note that we keep the gas surface density and disk mass constant in all our models.

3 Results

3.1 Spectral energy distribution (SED)

Figure 3 shows the impact of the metallicity in the models on the spectral energy distribution (SED). For the calculation of the SED, a distance to the source of 140 pc is assumed. It is clear that for reduced dust-to-gas ratio, the flux density of the disk also decreases. This is a consequence of the lesser amount of dust suggested by a model with reduced metallicity, which is shown clearly in Fig. 2 in terms of the dust surface density. At shorter wavelengths, the emission of the disk comes from the part where it is optically thick. Therefore, it does not vary significantly if a fraction of dust is removed. At λ < 1 µm, the SED flux density is dominated by the stellar emission, which remains unchanged for the three models as the star parameters do not change. This suggests that low-metallicity disks are more difficult to observe at longer wavelengths. At shorter wavelengths (e.g., near-infrared), they are as easy to observe as disks with higher metallicity. The impact that a lower metallicity has on the stellar parameters is a subject that will be covered in future studies.

Fig. 4 Contour plots of the models with different metallicities. The left, middle and right columns of the plot represent the models Z1, Z01, and Z001; the latter two exhibit 1/10 and 1/100 of the reference metallicity, respectively. The top, middle, and bottom rows correspond to the dust temperature, gas temperature, and the difference between the gas and dust temperature relative to the gas temperature. The decreasing metallicity leads to an increase in both dust and gas temperature but with the gas temperature increasing stronger (bottom row). The black dashed line represents the surface corresponding to a visual extinction of Av = 1.

Fig. 5 Contour plots that show the difference between the temperatures of the model Z1 and the Z01 model (top panel) and the Z001 model (bottom panel). The right-hand column shows the difference in dust temperature and the left-hand column the difference in the gas temperature. The red color represents the regions where the model with lower metallicity has a higher temperature. The white and blue colors represent the regions where both temperatures are similar and where the Z1 model has a higher temperature, respectively.

3.2 Dust and gas temperature

The first result we analyzed in the produced models was the influence of a lower metallicity on the gas and dust temperatures. For a description of the properties of model Z1, we separated the structure of the disk into three radial zones and three vertical layers following Henning & Semenov (2013). The radial structure consists of three zones: (1) inner zone: it is located closest to the star and extends up to ≈ 1 au and has a temperature above 100 K. (2) The middle zone: it is located between 1 and 30 au. The temperatures in this zone are in the range of 100 (middle and upper parts of the disk) ≥ T ≥ 20 K (midplane). 3) The outer zone: it starts at 30 au and extends to the edge of the disk. The temperature in this zone is below 20 K.

The vertical structure from bottom to top is as follows: (1) The midplane: the coldest vertical layer of the disk. Here the chemical reaction is predominantly gas grain interaction and dust surface chemistry. (2) The rich molecular layer: the temperature in this layer lies between 20 and 100 K. It is warm enough to keep some species in the gas phase (e.g., CO). It is also deep enough that the molecules are mostly shielded from photodissociation. The conditions in this layer allow rich molecular chemistry to take place. (3) The photon-dominated layer: the hottest layer is characterized by temperatures between 100 and 5000 K. X-rays and UV radiation are the main drivers of the chemistry. Most of the species here are in atomic form and very often ionized.

Figure 4 displays the dust and gas temperatures and their differences in all the models. The figure shows that the decrease in metallicity produces higher gas and dust temperatures. As already mentioned, the top layer of the disk (photon-dominated layer) is characterized by gas temperatures of T_gas > 1000 K. Both gas and dust increase in temperature as the metallicity decreases, mostly because the radiation is able to penetrate the disk more deeply with less of a shielding effect from the dust mass in the environment. This aspect also affects the thermal accommodation process between dust and gas and hinders the dust from reaching the same temperatures as the gas.

The middle layers (rich molecular layers) exhibit a greater temperature for the gas for a small planetary radius but a greater dust temperature for radii greater than 100 au. In this part of the disk, the collisional exchange of energy is not enough to make both dust and gas temperatures equal. Gas and dust are thermally decoupled and depending on the position in this layer, the dust or gas will be warmer. For a radius smaller than 100 au, this region is being directly hit by the stellar high energy radiation and the dominant heating processes are chemical heating followed by heating by the formation of H₂ on dust and heating by collisional de-excitation of H₂. These UV-driven heating processes elevate the gas temperature. For radii beyond 100 au, we have regions where the dust is warmer than the gas. These regions of the disk are shielded more effectively from the stellar radiation so that the heating of the gas by direct radiation (mainly X-rays) is not efficient. Additionally, OI and other line cooling processes of the gas dominate this region, leading to the dust becoming warmer than gas. As the metallicity decreases, the region where dust is warmer than the gas shrinks because the radiation is not being shielded as strongly and the gas heats up again.

Finally, for all three models with different metallicity, we observe that the gas and dust temperatures are in equilibrium in the midplane of all models. This section of the disk has a relatively low temperature and high density. This means that the energy exchange rate by collisions between dust and gas particles is frequent enough to reach thermalization. Thus, the gas and dust temperature become the same through a process called thermal accommodation (Woitke et al. 2009). As the metallicity decreases, the dust density of these deep layers also decreases and the thermal accommodation becomes less efficient. This leads to a shrinking of the region (white area) where the dust and gas have the same temperature.

We note that a general trend of increasing temperature with decreasing metallicity is true for the type of models considered in this work, which are often referred to as passive disks. For passive disks, the only heating source is stellar radiation and cosmic rays. In the case of active disks, the hydrodynamical processes provide additional heating via viscous and compressional heating, which operate predominantly in the disk’s midplane. The active disk tends to decrease its temperature with a lower metal-licity as the heat in a more diluted disk (characterized by lower opacity) can escape more easily and thus cool the disk (Vorobyov et al. 2020).

In order to take a closer look at the impact of the different metallicity on the gas and dust temperature of the disk, we compared the temperature of the lower metallicity models with model Z1. Figure 5 shows the influence of a lower metallicity over the gas and dust temperature by displaying 2D plots of the relative temperature difference between the model Z1 and those with the reduced metallicity. A decreased metallicity leads to an increase in the gas and dust temperature in the disk. This effect is stronger in the middle layers of the disk but it also takes place in the midplane beyond a radius of 20 au for the Z01 model and a radius of 4 au for the Z001 model.

Fig. 6 UV radiation field for the Z1 model (left panel), model Z01 (middle panel), and Z001 (right panel), in units of the Draine field. The dashed contour lines correspond to the tick values in the color bar. The decreasing metallicity enables the UV radiation to get to the midplane at a smaller radius.

3.3 Radiation field

In Fig. 6, the UV radiation (91.2 nm < λ < 250 nm) field in Draine units is shown. The contour lines in the plots show that for the reduced metallicity models the radiation has higher values along the midplane of the disk. For the Z1 model, the radiation value of 10⁻¹⁵ reaches the midplane at a radius greater than 100 au. For the Z01 model, the same radiation value is reached at a radius around 10 and 20 au and for the Z001 model, this value is reached at a radius around 1 au. This confirms that models with lower metallicity allow the radiation of the star to penetrate deeper regions of the disk. In this case, this is simply the effect of reducing the dust mass. The same effect would be achieved by using unchanged elemental abundances. The effect of a stronger radiation field in the disk on the respective disk chemistry is presented in Sects. 3.4, 3.6, and 3.5.

Fig. 7 Contour plots of the abundances. From left to right and top to bottom: H2O, CH4, CO, CO2, HCN, and CN are represented relative to the total hydrogen nuclei number density. The dashed contour lines describe the tick values on the colorbar. The left column corresponds to model Z1. The middle column corresponds to the model Z01 and the left column to the model Z001. The upper row represents the gas species, the lower row the ice species, and the bottom row the total abundances. With decreasing metallicity the cavity in the midplane for all gas-phase species shown here shrinks (Z01) and vanishes (Z001). In the case of CH4, a second cavity in the upper layers of the disk grows as metallicity decreases.

3.4 Spatial distribution of abundances

In this section, we present the impact of metallicity on the abundances of various molecules relative to hydrogen. The ice and gas abundances of H₂O, CH₄, CO, CO₂, HCN, CN, HCO⁺, and N₂H⁺ are displayed as a function of the radius and the height of the disk in Figs. 7 and 8. These species were selected because they are a product of reactions involving the two most abundant molecules in the gas phase (H₂ and CO) and because their molecular lines are known tracers of physical and chemical disk properties. For instance, lines of H₂O, CO in the IR regime can give us information about the temperature in the disk (Henning & Semenov 2013). Then, the CN, HCN lines in the sub-mm and mm regime serve as photochemistry tracers (Henning & Semenov 2013) and HCO⁺ and N₂H⁺ lines in the same regime are used to trace the ionization (Henning & Semenov 2013). In the study performed by Hobbs et al. (2022), one of the main findings was that the abundances of CH₄ on cooler hot Jupiters can be used to link the composition of a planet’s atmosphere to its formation location. In addition, even though CO is still considered to be the best gas mass tracer, CO₂ and H₂O can potentially be used as alternative mass tracers as well (Molyarova et al. 2017). We chose to take a deep look into these specific molecules also because they play an important role as the building blocks of complex carbon-based molecules – that are also the building blocks for organic life.

The impact of metallicity on the abundances can be described as follows. On the one hand, as the metallicity decreases the relative abundances of ice and gas species decrease the entire disk. This general decrease is the result of the reduction of the elemental abundances mentioned in Sect. 2.1. On the other hand, the dust reduction applied to obtain lower metallicities affects specific regions of the disk. We focus on the impact on specific regions of the disk:

First, the midplane shows a cavity for the gaseous species that shrinks with decreasing metallicity. The ice species are also affected by this and exhibit a shrinking of their spatial distribution. This is a consequence of the increase of the input energy by a stronger radiation field in the disk as the metallicity decreases. The stronger radiation field desorbs or dissociates ice-phase species and enhances the formation reaction of some gas-phase species. This leads to a depletion of icy species in the location of the midplane cavity and an increase in the abundance of the respective gaseous species.

Second, the total abundances displayed in the bottom row in the panes of Fig. 7 also exhibit a cavity in the midplane that shrinks with decreasing metallicity. In the case of H₂O, the total abundance does not show a midplane cavity for any metallicity, but a general decrease in all the disks is noticeable. Apart from H₂O, the behavior of the midplane cavity for the total abundances mirrors the replenishment of gaseous species as the metallicity decreases.

We decided to carry out a simple comparison of the Z01 and Z001 models with models where we artificially reduced the chemical abundances of the Z1 model by a factor of 10 and 100, respectively. We did this comparison to show the non-linear behavior of the chemistry when reducing the metallicity. We refer to the models with an artificially reduced column density as “scaled-down” models.

3.5 Vertical column densities

In this section, we offer a detailed analysis of the vertical column density of the above-mentioned species, along with the chemical reactions responsible for the formation and destruction of the molecules studied in this work. We analyze the behavior of the total vertical column density of the aforementioned species. The ice species are followed by a # symbol.

Figures 9–12 show the total vertical column density of the mentioned chemical species. We again made a simple comparison of the Z01 and Z001 models with models, where we artificially reduced the vertical column density of the Z1 model by a factor of 10 and 100, respectively. The red shaded area in the figures highlights the cases where the column density is higher than the simple scaled-down of the Z1 model abundances and the green shaded area where it is lower. For each molecule, we describe the general trends of the column density. We also give in Table 4 the respective changes in the formation and destruction reactions for the molecules for different metallicities. A more detailed display of the behavior of the molecules in a frozen and gaseous state is given in Appendix A. The relevant reactions are given in a detailed form in Table A.1 as well. For completeness, we also show the impact of the dust reduction and the elemental reduction of the metals separately in Appendix B.

Fig. 8 Contour plots of the abundances for HCO+ and N2H+. We only show the gaseous abundances because HCO+ and N2H+ do not freeze-out in our models.

3.5.1 H₂O

The top panel of Fig. 9 shows a total column density decrease with radius for all three models. However, there are small local drops that are located in different radii for each model. The total column densities (gas + ice) shown in the middle and bottom panels have equal or slightly lower values than the scaled-down values. The reason for the lower than the scaled-down values lies in the stronger radiation field and efficient photodissociation of H₂O# (Eq. (A.1)) in the deeper parts of the disk. It is the consequence of the radiation of the star being able to enter deeper regions of the disk as the metallicity of the model decreases. For Z1, the main destruction reaction for H₂O# in the disk region of underabundance is thermal desorption, which does not affect the total column density. Table 4 shows that for lower metallicities (Z01 and Z001), photodissociation takes over as the main destruction reaction because the radiation field is stronger in that part of the disk. Thus, H₂O# is not only being sublimated by the rising temperatures, but is also dissociated by photons that are able to reach the disk midplane (Andersson et al. 2006). In particular for Z001, the small drop beyond 2 au starts at the same radius where a great part of the UV radiation is not being strongly shielded and it reaches the midplane for that model (see Fig. 6). The resulting effect is a more efficient depletion of all H₂O as the metallicity decreases. This is shown clearly in Fig. 13, where the total ratios $R_{{\text{H}}_2\text{O}}<1$.

Fig. 9 Total vertical column density of H2O for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red and green areas show the cases where the vertical column density is higher or lower than the simple scaled-down values. The bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

3.5.2 CH₄

All three models show a decrease that starts at 1 au and continues until a radius of 20 au, where a strong increase follows (top row of Fig. 10). The decrease from 1 au to 20au is due to the cavity for the total CH₄ that is visible in Fig. 7. The middle and bottom rows show an underabundance inside 0.4 au for the lower metallicity models. This is mainly due to the enhanced depletion by photodissociation of gaseous CH₄ caused by a stronger radiation field. It is also shown that roughly between 0.8 and 8 au, the lower metallicity models have higher values than scaled-down values. The overabundance of Z01 and Z001 between 0.8 and 8 au is the effect of the larger formation rates that the lower metallicity models have near the midplane and around a radius of 1.5 au. For instance, the CH₄ formation rates in that region for the Z1, Z01, and Z001 model are 3.44 × 10⁻¹³, 3.84 × 10⁻¹¹, and 1.21 × 10⁻⁷ cm⁻³ s⁻¹, respectively. The main destruction rates of CH₄ also increase with decreasing metallicity but not enough to counteract the enhanced formation of CH₄ (see Table 4). Beyond 20 au, the three models follow the scaling-down by one and two orders of magnitude very closely and have either equal or slightly lower values than the scaled-down values. The total molecular amount (Fig. 13) shows that the Z01 and Z001 models have a lower than 1 ratio in comparison to the scaled-down values. Thus, the overabundance is present in regions of the disk where the total abundance of CH₄ is very low and it does not compensate for the underabundance in the parts of the disk with a radius greater than 25 au.

3.5.3 CO

We note that CO shows a very similar trend to H₂O, with the very important exception: the magnitudes of the differences to the scaled-down values are greater than for H₂O. For all three models, the total column density decreases with increasing radius. The Z1 model has a strong bump for a radius between 0.1 and 1 au. This feature is absent or much less pronounced in the case of lower metallicity models. The reason for the existence of the mentioned bump is the different formation reactions of CO for the three different models. In that part of the disk, the reactant for the most important formation reaction of CO is gaseous H₂O. This reactant has an abundance that is orders of magnitude lower for models Z01 and Z001, so it is to be expected that CO will be much less abundant in that region for models Z01 and Z001. This is explained in greater detail below (Sect. 3.5.7). As a result, CO is less abundant in the lower metallicity models than the scaled-down values in the vicinity of the bump. Apart from that bump, models Z01 and Z001 show a CO overabundance along the complete disk. This overabundance increases as the metallicity decreases. Table 4 shows that the formation reaction for CO in that region changes with decreasing metallicity and also that the respective reaction rate increases. For example, the most important formation reaction of CO in the Z1 model is dissociative recombination of H₂CO⁺, while for the Z01 and Z001 models the formation reaction is the neutral-neutral reaction of H and HCO and photodissociation of CO₂, respectively. We can see in Fig. 13 that CO is overabundant relative to the scaled-down values as a result of the enhanced formation reactions. Table 5 confirms this overabundance with ratios for the gaseous state of 4.78 and 16.92 for models Z01 and Z001, respectively. This is an important finding as CO is often used for observations. The resulting column density seen in Fig. 11 for models Z01 and Z001 does not show differences of order of magnitude compared to model Z1. Thus, if observations were performed on a sub-solar metallicity region, it would be incorrect to expect lower CO abundances than in the solar metallicity case. This assumption would then result in an underestimation of CO column density by an order of magnitude.

Fig. 10 Total vertical column density of CH4 for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red and green areas show the cases where the vertical column density is higher or lower than the simple scaled-down values. The bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

Fig. 11 Total vertical column density of CO, CO2, HCN, and CN for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red and green areas show the cases where the vertical column density is higher or lower than the simple scaled-down values. Bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

3.5.4 CO₂

A similar bump to the CO case for small radii is also present here for model Z1. This bump is a consequence of the bump in the CO case. The origin of this bump will also be explained in more detail in Sect. 3.5.7. Apart from that, all three models show a drop that reaches its minimum value at around 2 au. This drop is consistent with the location of the midplane cavity for this molecule. From 2 au on all models show an increase of the column density with radius until it starts to drop again towards the outer disk radius. All the models show a similar increase that drops again at the end of the disk’s radius. This decrease and rise of the total column density around 2 au follow the depletion of gas-phase CO₂ due to ion-neutral reactions (model Z1 and Z01) and photodissociation (model Z001). On the other hand, CO₂# is also being depleted in that part of the disk due to photodissociation (model Z1) and photodesorption (models Z01 and Z001). It is only beyond that radius where CO₂# builds up and is responsible for the increase of the total column density. The comparison between the Z01 and Z001 case, where middle and bottom panel for CO₂ in Fig. 11 shows an underabundance in the region where the bump is placed; whereas between the aforementioned bump and the outer parts of the disk, CO₂ shows an overabundance. This overabundance present in the Z01 model becomes smaller in model Z001. In the part of the disk with overabundance, we again see a change in the formation reaction type and in the reaction rate as metallicity decreases. For instance, the formation reaction changes from thermal desorption to photon desorption to the neutral-neutral reaction of OH and CO for models Z1, Z01, and Z001, respectively. There is an underabundance in the outer parts of the disk that becomes more noticeable and extends into the disk with decreasing metallicity. This increase in size and magnitude in the underabundance as the metallicity decreases leads to the resulting ratios in Fig. 13 for CO₂ in the sense that while the Z01 model shows a net overabundance, the Z001 model has a net underabundance.

Fig. 12 Total vertical column density of HCO+ and N2H+ for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red show the cases where the vertical column density is higher than the simple scaled-down values. The bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

3.5.5 HCN

The trend shown in the top panel for HCN in Fig. 11 exhibits a bump in the Z1 model for small radii. Aside from this bump, all models have a depression that reaches a minimum value at around 6 au and then show an increase beyond that radius. The values for a radius smaller than 2 au are caused by the amount of gaseous HCN, while the values beyond 9 au by the amount of HCN#. In the middle and bottom panel, models Z01 and Z001 show an overabundance for different parts of the disk (for very small radii, for a radius between 0.9 and 9, and for a region around 35 au). The region with the most overabundance is located around 1.8 au and Table 4 confirms the rising of the respective reaction rates as the metallicity is reduced. The formation reaction also changes. For models Z1 and Z01, it is the dissociative recombination of HCNH⁺ and for model Z001 it is the neutral-neutral reaction of CN and C₂H₂. The outer regions of the disk show an underabundance that increases for the Z001 model. Then Z001 model also shows an underabundance between 9 and 20 au. This underabundance is the result of the depletion of HCN# due to the greater values for the respective destruction reaction as metallicity decreases (green shaded rows for HCN# in Table 4. Although the underabundance seems to increase as the metallicity is reduced, as we see in Fig. 13, the ratios to the scaled-down values are still greater than 1 for the complete disk.

3.5.6 CN

The top panel displays the total column density for CN in Fig. 11, showing that none of the models exhibit a bump for a small radius. However, an oscillating behavior among the three models (Z1, Z01, and Z001) is noticeable. This is because most CN at a small radius is located in a very thin layer of the disk (see Fig. 7). The disk model is not able to resolve this thin layer properly and the result is the oscillating behavior of the total column density. We note two important aspects regarding this oscillatory behavior. First, the values for the CN column density are orders of magnitude less than the other molecules presented in this study. Second, we produced models with higher resolution to check how this would influence the oscillations. We found out that the produced column densities of CN are almost the same as in Fig. 11, but with weaker oscillations. Beyond 1 au, CN behaves in a unique way compared to the other molecules. For a radius between 1 and 5 au, the total column density increases by a couple of orders of magnitude as the metallicity decreases. This is the opposite impact that metallicity has on the total column density of the other molecules. It mirrors the clear increase of the CN abundance in the midplane (between a radius of 1 and 5 au) as metallicity decreases (see bottom row of the CN panels in Fig. 7). For the rest of the disk models, Z1 and Z01 have similar values, while model Z001 has lower values with a clear drop at 90 au. Apart from the drop at 90 au for the Z001 model, both models (Z01 and Z001) exhibit an overabundance throughout the complete disk. The region around 2.5 au shows the greatest overabundance and the reaction rates in Table 4 show a clear increase as the metallicity decreases. Models Z1 and Z01 have dissociative recombination of HCNH⁺ as the main formation reaction and for model Z001, the respective reaction is photodissociation of the HCN molecule (Eq. (A.12)). The clear dominant overabundance along the disk confirms the ratios for CN, as shown in Fig. 13.

3.5.7 Bump at small radii for CO, CO₂, and HCN

The aforementioned bump at small radii for CO is of the following origin: in the location of the bump, the rate of the most important formation reaction of CO for model Z1 is greater than the main formation reaction for the Z01 and Z001 models by more than three orders of magnitude. This is mainly the result of the differences in the abundances of gaseous H₂O between the different models. The main formation reaction of CO in that disk region is the ion-neutral reaction of H₂O and HCO⁺ (Eq. (A.4)). Gas-phase H₂O is the reactant of the formation reaction for CO and is much more abundant in the Z1 model in that particular region. Therefore, the formation of CO with the mentioned ion-neutral reaction is much more efficient for the Z1 model relative to models Z01 and Z001. Similarly, in the case of CO₂, the bump is the result of the greater CO abundance of the Z1 model; CO is the reactant for the main formation reaction of CO₂, therefore, a greater CO abundance leads to a greater CO₂ abundance. For HCN, all three models have the same main formation reaction. It is the ion-neutral reaction of NH₃ with HCNH⁺ producing HCN and ${\text{NH}}_4^{+}$. Although all the models have the same main formation reaction, the bump is the result of the different reaction rates characterizing each model.

3.5.8 HCO⁺ and N₂H⁺

We decided to include HCO⁺ and N₂H⁺ in this section together because they are electrically charged and are only present in the gas phase. The column density of the two ions N₂H⁺ and HCO⁺ is displayed in Fig. 12. The column density of N₂H⁺ present an oscillating behavior similar to CN for a small radius. As with CN, higher resolution models show that the oscillations become weaker as the resolution of the model increases. This effect is also the result of a very thin layer in the inner parts of the disk which cannot be properly resolved in the model. In the top panel, HCO⁺ shows similar values for the three models until a radius of 2 au is reached. Beyond that radius, models Z01 and Z001 have clearly higher total column densities than model Z1. This is expected because of the stronger radiation field that results in a more efficient ionization in the lower metallicity models; N₂H⁺ shows a strong bump between 1 and 5 au and a following increase for models Z1 and Z01. Model Z001 does not have these strong features but it increases slightly as the radius increases. In addition, HCO⁺ and N₂H⁺ are more abundant than the scaled-down values for all radii. In fact, the overabundance is almost everywhere a couple of orders of magnitudes in both lower metallicity models (Z01 and Z001). Table 4 shows that for HCO⁺ the respective ion-neutral formation reaction changes from H⁺ + H₂CO for model Z1 to ${\text{H}}_3^{+}+\text{CO}$ for models Z01 and Z001. For N₂H+, the ion-neutral formation reaction ${\text{H}}_3^{+}+{\text{N}}_2$ does not change with metallicity, but the rate of the reaction has a dramatic rise. The destruction reaction for the two ions is dissociative recombination (Table 4). The electrons needed for this reaction usually come from the ionization of metals. However, as metallicity decreases, the number of available metals is also reduced. The reduction of metals has a stronger impact on the destruction rate of HCO⁺ and N₂H⁺ than a stronger radiation field. This is shown in Appendix B, where the reduction of the metals and the dust reduction are handled separately. Eventually, this leads to an overabundance of the two ions. This general overabundance of both ions is displayed in a very clear way in the respective ratios shown in Fig. 13. All the chemical reactions mentioned in this section are shown in detail in Table A.1.

3.5.9 General trend for column densities and quantitative analysis

In most cases, the decreased metallicity leads to a decrease in the column density in comparison to the Z1 model. However, the actual reduction does not follow the simple reduction of the abundances by factors 10 and 100. It is important to reiterate that the over- or underabundances of the species are not only the result of different rates of adsorption and desorption of each species. The different trends of the column density for the lower metallicity models are also the result of surface reactions on the dust grains, photodissociation in the gaseous and ice phase, and recombination processes. The different values for the formation and destruction rates of the species are triggered by the higher temperatures and the stronger radiation field as the metallicity decreases.

In Table 5, we summarize the quantitative analysis of this comparison for each of the considered species. For each species, the vertical column density, N_sp, is defined as: $$N_{\text{sp}}\left(r\right)={\displaystyle {\int }_0^{z_{\max }\left(r\right)}{n}_{\text{sp}}\left(r,z\right)\text{d}z,}$$(5)

with z as the disk height, z_max as the maximal disk height, r as the radius, and n_sp as the number density of a chemical species at each point of the disk. For each species, we integrated the vertical column density N_sp over the radius r of models Z01 and Z001 and their scaled-down counterparts. We used a simple integration from the inner radius R_in to the outer radius R_out of the disk for this step (Eq. (6)) as follows: $$n_{\text{tot}}=2\pi {\displaystyle {\int }_{R_{\text{in}}}^{R_{\text{out}}}{N}_{\text{sp}}r\text{\hspace{0.17em}d}r},$$(6) $$R=n_{\text{tot}}/n_{\text{tot}\left(\text{sd}\right)}.$$(7)

Next, we calculated the ratio of the total molecular amount of models Z01 and Z001 (n_tot) relative to the scaled-down results (n_tot(sd)). The ratio is calculated with Eq. (7). These results are also presented in Fig. 13. The ratio is higher than 1 in almost all gaseous species. It is particularly high for the ion-species HCO⁺ and N₂H⁺. The exception of gaseous CO₂ is due to the fact that the underabundance in the disk zone closest to the star is also the zone where most part of CO₂ is stored. However, all gaseous species show that the value for the ratio has a positive trend as the metallicity decreases.

The ice species show a smaller than 1 ratio for H₂O# and CH₄# which decreases as the metallicity decreases. The rest of the ice species show a higher than 1 ratio for the Z01 model which decreases for the lowest metallicity model Z001. Here, CO₂# and CN# even exhibit a ratio that is lower than 1 for the Z001 model. All ice species show a negative trend as the metallicity decreases. This indicates that although there is still chemical production of certain ice species underway, the destruction due to photodesorption or thermal desorption of these species becomes more dominant with a decreasing metallicity. The total (gas + ice) values show a smaller than 1 ratio for H₂O and CH₄ with a negative trend with decreasing metallicity. Also, CO and CN show a ratio that is higher than 1 with about three times higher ratios for the lowest metallicity Z001. CO₂ and HCN have greater than 1 ratios for the Z01 model but exhibit a negative trend as the metallicity decreases. CO₂ even shows a smaller than 1 ratio for the Z001 model. In summary, the strongest deviation from the scaled-down values for gaseous species is found in the ions HCO⁺ and N₂H⁺, followed by CN, HCN, CO, and CH₄. Also, H₂O shows a slight deviation and CO₂ even shows smaller particle numbers than the scaled-down values. For the frozen species CN, CO, and CO₂ exhibit the strongest deviation from scaled-down values. Regarding the total (gas + ice) numbers, CO and CN are the molecules that have the strongest deviation.

Fig. 13 Ratio, R, (see Eq. (7)) of the total molecular amount of the scaled-down Z1 and the low metallicity models Z01 and Z001 for each species. The red, blue, and black circles correspond to the ratios of the gaseous, ice, and total (gas + ice) species, respectively. The dashed horizontal line shows the value for the ratio equal to 1.

3.6 Snowline for each molecule

Snowlines are regions in protoplanetary disks where the temperature reaches the sublimation point for a particular volatile species. The position of the snowline is of great importance to the formation and basic architecture of a planetary system. It plays an important role in the formation and chemical composition of planetesimals and, therefore, in the formation and composition of planets (e.g., Öberg et al. 2011; Okuzumi et al. 2016; Banzatti et al. 2015). The temperature in a protoplanetary disk has a radial and a vertical gradient. Therefore the snowline is not a single radius but, rather, a two-dimensional (2D) surface. We show the location of this surface by displaying the regions in the disk where the chemical species are more abundant in the gas phase (red area) or the ice phase (blue area), as shown in Fig. 14. The white regions represent a similar abundance of ice and gas and the black dashed line represents the location in the disk where gas and ice abundances are equal. Each row in the figure represents a different molecule and each column stands for a different metallicity model. As a result of the radial and vertical temperature gradient in the disk, a snow region is enclosed spatially by the 2D surface where the ice and gas abundances are equal. We further define the intersection of this 2D surface with the midplane as the snowline. The innermost intersection is defined as the inner snowline and the next intersection as the outer snowline. For most of the molecules studied here, the outer snowline is not present. However, there are some cases where the outer snowline lies inside the disk and we refer to an ice ring inside the disk.

Based on Fig. 14, it is clear that as the metallicity decreases, the snow region shrinks in its vertical and radial extension for all molecules, and in some cases, the snow region almost completely vanishes (e.g., CO₂ and CN). This is the result of the general increase in the disk’s temperature caused by the stronger radiation field in the disk. Figure 15 summarizes the extension of the snow region and the location of the respective snowline at the midplane for each molecule for all models. The comparison between the models shows an influence of the metallicity on the location of the snowlines for all molecules. All the inner snowlines of each species are pushed outwards with decreasing metallicity. In some cases, the decrease in metallicity moves the outer snowline inside the disk so that an ice ring inside the disk is present. Two of the species (CO₂ and CN) used for this section exhibit the presence of ice rings inside the disk. Figures 14 and 15 show that these ice rings shrink with decreasing metallicity (CO₂ and CN in particular).

The effect of metallicity on the inner (outwards shift) and outer (inwards shift) snowlines is an indication that reducing the metallicity leads to a higher dust temperature (which is the temperature that determines if thermal desorption of frozen species from dust grains takes place) and an enhanced photo-desorption of ice species due to a stronger radiation field in the disk. Additionally, the fact that lower metallicity models have less amount of dust means that the gas phase molecules have less surface available to freeze out. All this favors the depletion of ice species and, consequently, causes the ice rings to shrink as the metallicity decreases.

As a relevant example for planet formation, the snowline for H₂O (0.39, 0.45, and 0.59 au for models Zl, Z01, and Z001, respectively) does not change dramatically with metallicity. Thus, the position of the H₂O snowline on planet formation would not be affected by metallicity changes. Another snowline that does not exhibit a strong change with metallicity is the CO₂ snowline (1.22, 1.49, and 1.70 au). In the case of the CO snowline, the impact of metallicity is more noticeable (11.64, 14.21, and 85.45 au). In particular, for the Z001 model, the CO snowline is pushed outwards to 85 au. The location of the snowlines for these three molecules has special significance in the formation and composition of planets, asteroids, and comets in a protoplanetary disk as these molecules are among the main carriers of C and O. Whether the atmosphere of a planet or a comet is carbon-rich or oxygen-rich, and whether it will be rich in more complex molecules, will depend on radial variations in disk midplane composition, and specifically the locations of the snowlines (Walsh et al. 2016). As CO is considered to be the starting point for the formation of many complex organic molecules (COMs), the metallicity will certainly impact the location in the disk where COMs are more likely to be found and therefore affect the chemical composition of the planets, asteroids, or comets that are formed at those radii.

Fig. 14 Gas-to-ice ratio of the molecules. The left, middle, and right columns represent the Z1, Z01, and Z001 models, respectively. Red and blue regions represent the regions where gas and ice are dominant respectively. The black dashed line shows where the ice and gas abundances are equal (snowline). The locations of the snowlines at the midplane are summarized in Fig. 15.

Fig. 15 Snow regions at the midplane for the Zl (red), Z01 (green), and Z001 (blue) models for all neutral species. The horizontal bars show the region in the midplane where the abundances of the ice species are larger than the abundances of the gaseous species. These regions are enclosed by the inner and outer snowline.

4 Limitations of the model and outlook for future studies

In this section, we briefly discuss certain limitations of our model and their possible consequences. In Sect. 2.2.1, we introduce the physical structure of our model. We note that the column density power index λ of our models shown in Eq. (2) depends on the mass transport mechanisms operating in the disk, which in turn may depend on metallicity (De Marchi et al. 2013). This could lead to different timescales for accretion processes and the dissipation of the disk. In this study, we neglect these possible complications and assume that the disk structure is identical in the considered metallicity range.

We used a fixed stellar spectrum for all the models. In order to provide a more realistic model, we should take into account that the stellar spectrum changes with a different metallicity. This limitation is due to the fact that the model uses stellar spectra with solar metallicity by default. A consequence of this simplification is that the effective temperature of the star is underestimated, as stars with lower metallicity tend to be hotter.

For this work, we focused on the behavior of certain molecules. We chose this set of molecules because of their relevance in dust growth, planetary atmospheric composition, and importance for observations. The inclusion of more complex organic molecules in the analysis will be addressed in future papers.

For this study, we assumed a steady-state chemistry. Although this does not affect the consistency of our results, time-dependent chemistry becomes relevant when studying the impact of processes (e.g., accretion bursts) that have a strong influence on the disk. In a future study, we intend to include the impact of metallicity and accretion bursts in the chemistry of the disk.

We also note the particular case of CO. This molecule is very relevant because of its importance as a gas mass tracer (Molyarova et al. 2017). Two very important conditions that allow the usage of CO as a gas mass tracer are that its emission has to be optically thin and it cannot be too strongly depleted in the gaseous state. Our results suggest that changing the metallicity does not strongly affect the total abundance (gas + ice) of CO (relative to the response of other molecules and excluding the bump inside 1 au). This would mean that CO-based observations would not be able to differentiate disks with different metallicities. This could be linked to the lower water abundance for lower metallicity disks and would suggest that low metallicity disks have a higher CO/H₂O ratio compared to “normal” disks (see Fig. 16). To investigate the impact on observables requires proper modeling of line emission of the species in question, which is out of the scope of this paper but is planned in a future study.

In future studies, we will be looking more deeply at different aspects of the impact of metallicity on protoplanetary disks. In particular, the effects of time-dependent chemistry in the context of an accretion burst, the implementation of a more realistic stellar spectrum, and dust size distribution in the disk. We will also study behavior of more complex organic species and the modeling of their observation.

Fig. 16 CO/H2O ratio for models with different metallicity.

5 Summary and conclusions

We started this study to look into the impact of lower metallicities on certain chemical species in the protoplanetary disk. We used the radiation thermo-chemical code for protoplanetary disks ProDiMo to create a model with a reference metallicity (Z1) and two models with lower metallicities by a factor of 10 (Z01) and 100 (Z001). We simulated the lower metallicities by decreasing the dust-to-gas ratio in the disk by one and two orders of magnitude and by decreasing the initial element abundances by the same factor in the initial setup of the models respectively. We analyzed the impact of metallicity on the chemical species H₂O, CO, CO₂, CH₄, CN, HCN, HCO⁺, and N₂H⁺ because they are often used in observations as useful tracers of different disk properties. The findings in this work could therefore provide useful information for future observations in low-metallicity environments. One of the most important aspects of this study is the comparison of the lower metallicity model’s vertical column density of chemical species with the reference metallicity vertical column density scaled-down by factors of 10 and 100. We defined the notion of over-and underabundance of the models whenever the lower metallicity models show greater or smaller values for the vertical column density than the scaled-down vertical column density. We analyzed the chemical formation and destruction reactions in the midplane wherever the over-or underabundance is strongest. Our main findings are as follows:

• The effects shown in the lower metallicity models cannot be explained by simple scaled-down values of the abundances of the chemical species. The obtained values for the abundances of the Z01 and Z001 models differ from the scaled-down values almost everywhere in the disk. This is because the lower metallicity models allow the radiation field to penetrate deeper regions of the disk. The temperature also increases in the regions where the radiation is stronger than in the Z1 model. This affects the formation and destruction reactions of the species;

• The main effect, namely, a stronger radiation field appears to be due to the enhancement of the different processes responsible for the destruction of ice species and the formation of gas species and the shifts in the type of destruction and formation reactions. Thus, the stronger radiation field causes a substantial increase in the formation efficiencies of gas-phase species. This is clearly shown in Fig. 13 where the ratios between the scaled-down values and the lower metallicity values of all gaseous species (red dots) and the ice species (blue dots) are displayed. Furthermore, the impact of metallicity on the chemical reactions is different for each individual species. In some cases, the reaction type changes with metallicity, while in other cases, the type of reaction does not change with metallicity, but the rate of the reaction is clearly affected. This is shown in detail in Table 4, where the different main reactions and the respective rates of the three models are compared;

• Our results show that CO shows a greater abundance than the scaled-down values for lower metallicities. This molecule shows a clear non-linear response to the change in metallicity. Beyond the bump at small radii, metallicity does not seem to have a strong influence on the column density of CO. This finding has very relevant implications for observations as CO would not be a useful indicator for the metallicity of the disk;

• The impact of metallicity on the position of the snowlines is also present. With decreasing metallicity the snowline of the molecules analyzed in this study is “pushed” outwards. Additionally, decreasing the metallicity leads to the creation of ice rings for some species (e.g., CO₂) and to the shrinking of already existing ice rings for other species (e.g., CN).

The inclusion of other metallicity-dependent features will be performed in future studies (star spectrum and dust size distribution among others). The analysis of other scenarios (e.g., accretion bursts) will also help us gain more insight into the impact of metallicity in protoplanetary disks.

Acknowledgements

We would like to express our gratitude to the anonymous referee for the valuable comments that helped us improve the presentation of this work. The computational results presented have been achieved using the Vienna Scientific Cluster (VSC). This work was supported by the Austrian Science Fund (FWF) under research grant P31635-N27. Ch. Rab is grateful for support from the Max Planck Society and acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 325594231.

Appendix A Detailed chemical reactions and display of the ice and gaseous column densities

In this appendix, we offer a more complete and detailed analysis of the behavior of the species in the sense that one can see if the ice or gas phase is more affected by the metallicity change. Figure A.1 adds the behavior of the vertical column density for the gas and ice species as presented for the total column density in Figs. 9 to 12. Table A.1 shows the specific chemical reactions referred to in Table 4 that take place at the positions where over-or underabundance is largest.

Fig. A.1 Gas (left column), ice (middle column), and total vertical column density (right column) of all examined molecules (except for HCO+ and N2H+) for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red and green areas show the cases where the vertical column density is higher or lower than the simple scaled-down values. The bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

Appendix B Reducing the element abundance and the dust-to-gas ratio separately

To study the effect of the respective modifications we performed to simulate a lower metallicity we produced models where only the dust-to-gas ratio or the element abundances are lowered by a factor of 10. The resulting vertical column densities are shown in Fig. B.1. In most of the cases, the models where only the elemental abundance was reduced show a profile that resembles a simple scaled-down value from the Z1 model. The models where only the dust-to-gas ratio was reduced exhibit different features for certain radii that match the profile of the lower metallicity model. A good example of this is water-ice shown in the top right panel of Fig. B.1. The model where only the abundances were reduced shows a shift downwards that mostly matches the Z001 model. However, the drop in column density between 10 and 20 au in the Z001 model is absent for the abundance reduced model. The model with only dust reduction is the one that exhibits the drop at the same radius. This suggests that on the one hand, the reduction in the abundances is responsible for a general scaling down of the column density in the entire disk. On the other hand, dust reduction is responsible for the specific features (drop or rise) in certain parts of the disk.

Fig. B.1 Vertical column density of chemical species. With the exception of the bottom row, the left column shows the gaseous species and the left column the frozen species. the last row shows the ions HCO+ (left) and N2H+ (right). The reference model Z1 (solid black line), the model with only an elemental abundance reduction (dashed red line), the model with only dust-to-gas reduction (dotted blue line), and the Z01 (solid green line) show the effect of the two modifications done in this study separately.

References

All Tables

All Figures

	Fig. 1 Gas density structure (top) and radial surface density profile (bottom) of the reference model Z1.
In the text

	Fig. 2 Surface density in the midplane for dust (black) and gas (blue) as a function of radius. The models Z1 (solid line), Z01 (dashed line), and Z001 (dotted line) are shown.
In the text

	Fig. 3 Spectral energy distribution of the different disk models and stellar spectrum. The models Z1 (solid line), Z01 (dashed line), and Z001 (dotted line) are shown. A lower metallicity leads to a decrease in the flux density for the SED.
In the text

Fig. 4 Contour plots of the models with different metallicities. The left, middle and right columns of the plot represent the models Z1, Z01, and Z001; the latter two exhibit 1/10 and 1/100 of the reference metallicity, respectively. The top, middle, and bottom rows correspond to the dust temperature, gas temperature, and the difference between the gas and dust temperature relative to the gas temperature. The decreasing metallicity leads to an increase in both dust and gas temperature but with the gas temperature increasing stronger (bottom row). The black dashed line represents the surface corresponding to a visual extinction of Av = 1.

In the text

Fig. 5 Contour plots that show the difference between the temperatures of the model Z1 and the Z01 model (top panel) and the Z001 model (bottom panel). The right-hand column shows the difference in dust temperature and the left-hand column the difference in the gas temperature. The red color represents the regions where the model with lower metallicity has a higher temperature. The white and blue colors represent the regions where both temperatures are similar and where the Z1 model has a higher temperature, respectively.

In the text

	Fig. 6 UV radiation field for the Z1 model (left panel), model Z01 (middle panel), and Z001 (right panel), in units of the Draine field. The dashed contour lines correspond to the tick values in the color bar. The decreasing metallicity enables the UV radiation to get to the midplane at a smaller radius.
In the text

Fig. 7 Contour plots of the abundances. From left to right and top to bottom: H2O, CH4, CO, CO2, HCN, and CN are represented relative to the total hydrogen nuclei number density. The dashed contour lines describe the tick values on the colorbar. The left column corresponds to model Z1. The middle column corresponds to the model Z01 and the left column to the model Z001. The upper row represents the gas species, the lower row the ice species, and the bottom row the total abundances. With decreasing metallicity the cavity in the midplane for all gas-phase species shown here shrinks (Z01) and vanishes (Z001). In the case of CH4, a second cavity in the upper layers of the disk grows as metallicity decreases.

In the text

	Fig. 8 Contour plots of the abundances for HCO+ and N2H+. We only show the gaseous abundances because HCO+ and N2H+ do not freeze-out in our models.
In the text

Fig. 9 Total vertical column density of H2O for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red and green areas show the cases where the vertical column density is higher or lower than the simple scaled-down values. The bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

In the text

Fig. 10 Total vertical column density of CH4 for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red and green areas show the cases where the vertical column density is higher or lower than the simple scaled-down values. The bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

In the text

Fig. 11 Total vertical column density of CO, CO2, HCN, and CN for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red and green areas show the cases where the vertical column density is higher or lower than the simple scaled-down values. Bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

In the text

Fig. 12 Total vertical column density of HCO+ and N2H+ for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red show the cases where the vertical column density is higher than the simple scaled-down values. The bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

In the text

	Fig. 13 Ratio, R, (see Eq. (7)) of the total molecular amount of the scaled-down Z1 and the low metallicity models Z01 and Z001 for each species. The red, blue, and black circles correspond to the ratios of the gaseous, ice, and total (gas + ice) species, respectively. The dashed horizontal line shows the value for the ratio equal to 1.
In the text

	Fig. 14 Gas-to-ice ratio of the molecules. The left, middle, and right columns represent the Z1, Z01, and Z001 models, respectively. Red and blue regions represent the regions where gas and ice are dominant respectively. The black dashed line shows where the ice and gas abundances are equal (snowline). The locations of the snowlines at the midplane are summarized in Fig. 15.
In the text

	Fig. 15 Snow regions at the midplane for the Zl (red), Z01 (green), and Z001 (blue) models for all neutral species. The horizontal bars show the region in the midplane where the abundances of the ice species are larger than the abundances of the gaseous species. These regions are enclosed by the inner and outer snowline.
In the text

	Fig. 16 CO/H2O ratio for models with different metallicity.
In the text

Fig. A.1 Gas (left column), ice (middle column), and total vertical column density (right column) of all examined molecules (except for HCO+ and N2H+) for the models. Top row shows models Z1, Z01, and Z001 (solid, dashed, and dotted lines, respectively). The middle row shows the difference between the Z01 model (black dashed line) and the values by reducing the vertical column density of Z1 manually by a factor of 10 (red dashed line). The red and green areas show the cases where the vertical column density is higher or lower than the simple scaled-down values. The bottom row shows the same plots as the middle row but for the Z001 model and the Z1 model reduced by a factor of 100.

In the text

Fig. B.1 Vertical column density of chemical species. With the exception of the bottom row, the left column shows the gaseous species and the left column the frozen species. the last row shows the ions HCO+ (left) and N2H+ (right). The reference model Z1 (solid black line), the model with only an elemental abundance reduction (dashed red line), the model with only dust-to-gas reduction (dotted blue line), and the Z01 (solid green line) show the effect of the two modifications done in this study separately.

In the text