© ESO 2020

1 Introduction

Protoplanetary disks provide their host accreting star with material and regulate the formation, growth, and migration of planets. The global evolution and dispersal of disks around nascent stars is regulated by the transport of angular momentum and mass-loss processes (Armitage 2011). Fast jets¹ are ubiquitously observed in accreting young stars of all ages, with universal collimation and connection between accretion and ejection, which is probably of magnetic origin (Cabrit 2002, 2007b). However, the region of the disk that is actually involved in mass ejection and the associated angular momentum extraction remain a topic of hot debate.

The main observational method that has been proposed thus far to locate the launching region of jets relies on the joint measurement of rotation and axial velocities (Anderson et al. 2003; Ferreira et al. 2006). While studies of atomic jets are currently limited by spectral resolution in the optical range (De Colle et al. 2016), the unique combination of high spectral (< 0.5 km s⁻¹) and spatial (50 mas) resolution ofALMA now allows one to conduct similar tests on jets from the youngest so-called Class 0 protostars, which are much brighter in molecules (e.g., Tafalla et al. 2010) than jets from more evolved protostars and pre-main sequence stars (Class I and II), which are mainly atomic. In this context, high angular resolution observations (≃ 8 au) of the fastest part of the HH 212 jet have unveiled rotation signatures in SiO emission, which is suggestive of a disk wind that is launched within 0.3 au (Lee et al. 2017; Tabone et al. 2017). Owing to the high bolometric luminosity of the central protostar (L_bol ≃ 9 L_⊙, Zinnecker et al. 1992), dust is expected to be sublimated within 0.3 au, and the HH212 SiO-rich jet would thus trace a dust-free magnetohydrodynamic (MHD) disk wind.

To test the likelihood of this interpretation, it is now of paramount importance to check if the presence of SiO molecules in the jet of HH 212 is indeed compatible with a dust-free disk wind origin, as suggested by the rotation kinematics. Detailed astrochemical modeling of dusty magnetized disk winds show that molecules can survive the acceleration and the far-ultraviolet (FUV) field emitted by the accreting protostars (Panoglou et al. 2012). However, molecule formation in dust-free disk winds remains an open question. In the absence of dust, the FUV field can more easily penetrate the unscreened flow and photodissociate molecules, whereas H₂ formation on grains, which is the starting point of molecule synthesis, is severely reduced. So far, astrochemical models have only investigated dust-free winds launched from the stellar surface. Pioneering studies have shown that they are hostile to the formation of H₂ due to photodetachment of the key intermediate H⁻ by visible photons from the star (Rawlings et al. 1988; Glassgold et al. 1989; Ruden et al. 1990), whereas other molecules, such as CO, SiO, or H₂O, are destroyed when a strong UV excess is included (Glassgold et al. 1991). However, thus far, there have not been any similar investigations in a disk wind geometry and with a fully self-consistent FUV field.

The origin of the observed molecules in Class 0 jets is actually still debated, as high-velocity molecular emission is not necessarily tracing a pristine wind. Instead of assuming that molecules are material that is ejected from the vicinity of the star (“wind” scenario), another class of scenarios proposes that molecules are “entrained” by a fast and unseen atomic jet through bow-shocks or turbulence (Raga & Cabrit 1993; Canto & Raga 1991). Even if observations of small-scale molecular micro-jets (≃ 10−400 au, Cabrit et al. 2007; Lee et al. 2017) may contradict the entrainment of envelope material, a slow dusty disk wind surrounding the fast jet could still bring fresh molecular material close to the jet axis and explain the observed collimation properties of molecular jets (White et al. 2016; Tabone et al. 2018). Observational diagnostics of the wind launch radius based on rotation signatures would not be reliable anymore in such time-variable jets (e.g., Fendt 2011) or in jets prone to turbulent mixing.

In contrast with an inner disk wind scenario, the entrainment scenario implies that molecular material is rich in dust. Chemistry can then be used as a powerful diagnostic. Early astrochemical models of stellar winds have already pointed-out unique features of dust-free chemistry, such as a low CO/H₂ ratio (Glassgold et al. 1989). It suggests that chemistry is a promising diagnostic of the dust content and, as such, could distinguish between dust-free disk wind and an entrainment scenario. However, because sublimation temperature depends on the composition and sizes of grains, jets launched from the dust sublimation region of silicate and carbonaceous grains may contain a small fraction of surviving dust, such as aluminum oxide grains (Al₂O₃), for which the sublimation temperature is higher, ≃1700 K (Lenzuni et al. 1995). Despite representing a small mass-fraction of the total interstellar dust (≃ 2%, assuming solar elemental abundances, Asplund et al. 2009), aluminum oxide grains could have a strong impact on the chemistry and thus bias the proposed test. Understanding the precise impact of a nonvanishing fraction of dust on the chemistry is thus an important step in distinguishing between the “wind” and “entrainment” scenario.

Determining if molecular jets are indeed launched within the dust sublimation radius is a crucial question, as they would then bring new clues as to planet formation theories. Recent studies suggest that the first steps of planet formation may occur in the protostellar phase (e.g., Greaves & Rice 2010; Manara et al. 2018; Harsono et al. 2018). Probing the bulk elemental compositionand inner depletion pattern of the inner gaseous regions of protoplanetary disks within the dust sublimation radius is a powerful tool to uncover key disk processes related to planet formation, such as dust trapping in the outer parts of the disks (McClure & Dominik 2019; McClure 2019). However, deeply embedded Class 0 disks are too extincted to be probed by optical and near-infrared atomic lines. Hence, if protostellar molecular jets are tracing a pristine dust-free disk wind, they would offer a unique opportunity to have access to the elemental composition of the inner region of Class 0 disks, and thus reveal elusive disk processes. In this perspective, observations of high velocity molecular bullets toward active protostars show abundant oxygen-bearing species (SiO, SO) but a drop in carbon-bearing species, such as HCN or CS, which was interpreted as a low C/O ratio (Tafalla et al. 2010; Tychoniec et al. 2019). However, it remains unclear if abundance ratios of these molecular tracers are indeed directly related to a change in elemental abundances or if they are unique features of dust-poor chemistry or shocked gas. A fine understanding of the chemistry operating in dust-free or dust-poor jets is required when using molecular jets as a probe of elemental abundances of the inner disks.

In the present paper, we revisit pioneering astrochemical models of stellar winds by investigating if molecules can be formed in a disk wind launched within the dust sublimation radius. We focus our analysis on H₂ and on the mostabundant oxygen-bearing molecules observed toward protostellar jets, namely CO, SiO, OH, and H₂O (e.g., Tafalla et al. 2010; Kristensen et al. 2011). Sulfur and nitrogen chemistry, as well as the dependency of molecular abundances on elemental ratios, is beyond the scope of the present paper. Throughout this work, special attention is paid to the impact of a small fraction of dust on the chemistry to determine if molecular abundances can be used as a discriminant diagnostic even when the wind contains surviving refractory dust. The thermal balance, together with shock models, will be presented in the next paper of this series.

In Sect. 2, we present the basic physical and chemical ingredients of the astrochemical model. For the sake of generality, we then study molecule formation by using single point models (Sect. 3). This allowed us to derive simple criteria for molecules to be abundant as well as to specify their formation routes. More realistic models assuming a specific flow geometry are explored in Sect. 4. In contrast with single-point models, they include the effect of the dilution of the radiation field and density, time-dependent chemistry, and shielding of the radiation field. Limitations of the model as well as observational diagnostics of dust-free and dust-poor winds are discussed in Sect. 5. Our findings are summarized in Sect. 6.

2 Numerical method

Models presented throughout this work are based on the publicly available Paris-Durham shock code, which was initially designed to compute the dynamical, thermal, and chemical evolution of a plane-parallel steady-state shock wave. The code includes detailed microphysical processes and a comprehensive time-dependent chemistry (Flower & Pineau des Forêts 2003; Lesaffre et al. 2013). The versatility of this code also allows for computing the thermal and chemical evolution of any 1D stationary flow in a slab approach. In this work, the recent version developed by Godard et al. (2019), which includes key processes of photon-dominated region (PDR) physics, has been further upgraded to ensure the proper treatment of dust-free chemistry². Two types of models have been computed: single spatial point models (Sect. 3) and specific wind models (Sect. 4). Here, we present the basic numerical method used in this work. Details on the specific wind model, including the prescribed density structure, are given in Sect. 4.1.

Fig. 1 Schematic view of the 1D geometry used in this work. The code solves the chemical evolution of a slab of gas flowing at a constantvelocity Vgas and irradiated from the left (upstream) along a direction parallel to the flow. The proton density profile nH (z) is prescribedand the temperature TK is constant. The radiation field is the sum of a FUV part modeled by an ISRF scaled by G0 and a visible part modeled by a black-body radiation field at Tvis = 4000 K, which is diluted by a factor W (see Appendix A). The attenuation of the radiation field by gas-phase photoprocesses and by dust, if any, along z was consistently computed.

2.1 Geometry and parameters

Figure 1 gives a schematic view of the adopted geometry from which specific models can be built. The model is a slab of gas flowing at a velocity V_gas along the direction z and irradiated upstream (z = 0 in Fig. 1). The proton density n_H(z) = n(H) + 2n(H₂) + n(H⁺) along the slab is also prescribed. Throughout this work, the gas is assumed to be isothermal with the kinetic temperature T_K. This allows us to study the influence of the temperature on the chemistry in a parametric way and independent of uncertainties in the thermal balance.

The impinging radiation field is parametrized as a sum of a FUV component modeled as a standard interstellar radiation field (ISRF, Mathis et al. 1983) scaled by a parameter G₀ plus a visible component modeled as a black-body radiation field at T_vis = 4000 K diluted by a factor³ W. The adopted functional form of the FUV field as well as the shape of the resulting unshielded radiation field are given in Appendix A.

2.2 Radiative transfer and photodestruction processes

The attenuation of the radiation field between λ = 91 nm and 1.5 μm through the slab is computed by considering absorption by continuous photoprocesses along rays perpendicular to the slab and parallel to the flow. The attenuation coefficient due to gas-phase processes is then $${\kappa }_G(z,\lambda )={\sum }_i{\sigma }_i(\lambda )n_i(z),$$(1)

where σ_i(λ) is the cross-section of the photoprocess involving the species X_i and n_i is its particle number density. When dust is included, we assume that the grain size follows a power-law distribution of index − 3.5 with a minimum and maximum size of 0.01 and 0.3 μm, respectively (Mathis et al. 1977). Following Godard et al. (2019, see their Appendix B) and for simplicity purposes, absorption of UV photons by grains is calculated assuming the absorption coefficient of single size spherical graphite grains of radius a_g = 0.02 μm derived by Draine & Lee (1984) and Laor & Draine (1993), where $a_{\text{g}}\equiv \sqrt{\langle r_c^2\rangle }$ is calculated from the mean square radius of the grains.

Photodissociation and photoionization rates are computed following two different approaches. If the photoprocess leads to a significant attenuation of the FUV field or if it is a key destruction or formation pathway for a major species, then its rate is consistently computed by integrating the cross-section over the local radiation field including its UV and visible components (see Appendix A). Accordingly, the model includes absorption of photons by photoionization of C, S, Si, Mg, Fe, H⁻, H₂O, SO, O₂, and CH, and photodissociation of OH, H₂O, SiO, CN, HCN, H₂O, H${}_2^{+}$, SO, SO₂, O₂, CH⁺, and CH. Cross-sections are taken from Heays et al. (2017) and then they were subsequently resampled on a coarser grid of irregular sampling (about 100 points) to reduce computing time.

Alternatively, rates associated with other continuous photoprocesses are assumed to depend linearly on the integrated FUV radiation field as $$k_{\text{photo}}=\alpha \frac{F(z)}{F_{\text{ISRF}}},$$(2)

where F_ISRF is the FUV photon flux from 911 to 2400 Å, which is associated with the isotropic standard interstellar radiation field (ISRF), F(z) is the local FUV photon flux computed over the same range of wavelength, and α is the photodissociation rate for an unshielded ISRF. We note that even if this method seems to be crude, the nontrivial FUV attenuation by the dust-free gas prevents us from relying on more sophisticated fits as a function of N_H as used in typical dusty PDR models.

In order to avoid a prohibitive fine sampling of the radiation field around each UV line, CO and H₂ photodissociation are not treated according to the latter method. For CO photodissociation, we include self-shielding as well as shielding by H₂ and by a continuous process by expressing the photodissociation rate k_CO as $$k_{\text{CO}}=2.06\times {10}^{-10}{\text{s}}^{-1}{\theta }_1(N(\text{CO})){\theta }_2(N({\text{H}}_2))\frac{\chi }{1.23},$$(3)

where θ₁(N(CO)) and θ₂ (N(H₂)) are self-shielding and cross-shielding factors tabulated by Lee et al. (1996) and χ is the ratio at 100 nm of the local FUV flux to the mean interstellar radiation field of Draine (1978) (2  × 10⁻⁶ erg s⁻¹ cm⁻² Å⁻¹). The adopted shielding function gives a good approximation of the photodissociation rate of CO; however, it neglects the effect of the excitation of CO (Visser et al. 2009). It is important to note that the factor 1.23 stands for the ratio between the Mathis et al. (1983) and Draine (1978) UV flux at 100 nm. The photodissociation rate of H₂ was computed under the FGK approximation (Federman et al. 1979).

2.3 The excitation of H₂

Time-dependent populations of the first 50 ro-vibrational levels of H₂ in the electronic ground state are computed, including de-excitation by collision with H, He, H₂, and H⁺ (Flower et al. 2003); by the UV radiative pumping of electronic lines followed by fluorescence (Godard et al. 2019); and by the excitation of H₂ due to formation on grain surfaces. For H₂ formation on grain surface, we assume that the levels are populated following a Boltzmann distribution at a temperature of 1/3 of the dissociation energy of H₂ (Black & van Dishoeck 1987). For H₂ gas-phase formation routes, levels are populated in proportion to their local population densities.

2.4 Elemental abundances

Total fractional elemental abundances are constructed from Table 1 of Flower & Pineau des Forêts (2003). For dust-free models, all elements areplaced in the gas phase. For models with nonvanishing dust fraction, the dust content is quantified by the dust-to-gas mass ratio Q. The relative abundances between elements that are locked in the grains are assumed to be constant for any dust-to-gas ratio and equal to those of Flower & Pineau des Forêts (2003). PAHs are expected to be photodissociated by multiphoton events in the inner disk atmospheres (< 0.5 au, Visser et al. 2007). PAHs are consequently not included in the models and all carbon that is locked in PAHs is assumed to be released in the gas phase. The resulting fractional elemental abundances are given in Table 1 for dust-free models and for dusty models with Q = 6 × 10⁻³. This value, which corresponds to a fractional abundance of grain of 6.9 × 10⁻¹¹ (no sublimation of the grains), is taken as the reference for dusty models and we define Q_ref ≡ 6 × 10⁻³.

2.5 Chemical network

The chemical network is constructed from Flower & Pineau des Forêts (2015) and Godard et al. (2019). It includes 140 species and about a thousand gas-phase reactions. Details on the chemical reactions that were added to the former network, together with adopted rate coefficients, are given in Appendix B.

Regarding the formation of H₂, three gas-phase formation routes are included (see Fig. 2). The first route is the electron catalysis through the intermediate anion H⁻, via radiative attachment followed by fast associative detachment: $$\begin{array}{l}\text{H}+{\text{e}}^{-}\to {\text{H}}^{-}+h\nu ,\\ {\text{H}}^{-}+\text{H}\to {\text{H}}_2+{\text{e}}^{-}.\end{array}$$

This route has been shown to be quenched due to the photodetachment of the fragile H⁻ by visible and NIR fields in T Tauri stellar winds (Rawlings et al. 1988; Glassgold et al. 1989) but efficient in the early Universe at z < 100 (Galli & Palla 2013).We reconsider the role of this route in Sect. 3. The second route is the ionic catalysis by any ion noted here as X⁺ (with X = H, C, S, and Si), through the intermediate ion XH⁺ via radiative association followed by an ion-neutral reaction, and can be expressed as follows: $$\begin{array}{l}{\text{X}}^{+}+\text{H}\to {\text{XH}}^{+}+h\nu ,\\ {\text{XH}}^{+}+\text{H}\to {\text{H}}_2+{\text{X}}^{+}.\end{array}$$

This route is very similar to the former though it is built from X⁺ instead of e⁻. In contrast with H⁻, XH⁺ ions are photodestroyed by UV photons and can thus survive the strong visible radiation fields. Surprisingly, previous models of dust-free stellar winds have never discussed formation routes via CH⁺, SiH⁺, or SH⁺, focusing only on the formation route via H${}_2^{+}$. We show below and in Appendix C that the latter is inefficient compared to the formation by CH⁺ or SiH⁺. Lastly, H₂ is also formed through the three-body reaction: $$\text{H}+\text{H}+\text{H}\to {\text{H}}_2+\text{H},$$(8)

which is relevant at a high density. Complementary reactions involved in the chemistry of the intermediates H${}_2^{+}$ and H⁻ have also been included (see Appendix B).

Regarding H₂ formation on dust, we adopt the formation rate of Hollenbach & McKee (1979) assuming a single grain size distribution of radius a_g. An important caveat is the gas temperature dependence of the probability S(T_K) that a H atom sticks when it collides with a grain. In the high temperature regime that is relevant for jets, large discrepancies exist in the literature regarding sticking probabilities, especially when including chemisorption or specific substrates (see Flower & Pineau des Forêts 2013, Appendix A). Our adopted expression for S(T_K) (see Eq. (C.24)) gives a lower limit on the formation rate of H₂ at a high temperature.

In dust-poor gas, atomic species can dominate the opacity of the gas in the UV. When ionized, they can also contribute to the synthesis of the key anion H⁻ by increasing the electron fraction. Thus, a reduced chemical network for Mg/Mg⁺, incorporating photoionization and charge exchange reactions, has been added. Other rarer elements, such as Li, Al, or Na, are not expected to contribute significantly in either the attenuation of the radiation field or the ionization balance of the gas and are consequently not included in the model.

Fig. 2 Dust-free formation routes for H2 where X stands for H, C, S, and Si. Blue, green, and red, arrows represent two-body reactions, the three-body reaction, and the photodestruction of key intermediates, respectively.

3 Chemistry

To examine the formation and destruction routes of the main molecules observed in protostellar jets, single spatial point models have been computed over a range of physical conditions that is representative of protostellar jets. The evolution of the gas is assumed to be isothermal and isochoric. Parameters of the models presented in this section are proton density n_H, kinetic temperature T_K, unshielded FUV radiation field G₀, and dilution factor of the (visible) black-body radiation field at 4000 K noted by W. The explored parameter space is summarized in Table 2. In this section, we discuss steady-state chemical abundances. This approach, though simple, allows us to identify relevant chemical reactions and capture the essential features of the chemistry operating in jets. The results of this section are summarized in Sect. 3.3.

3.1 Formation of H₂

Formation of H₂ constitutes the first step of molecular synthesis. The dominant gas-phase formation route of H₂ and its efficiency depends mainly on the ability of H⁻ to survive to photodetachment (see Fig. 2). The influence of dust on the chemistry then depends critically on the efficiency of gas-phase routes. Our results are summarized in Fig. 3 for specific values of T_K = 1000 K, x_e = 4.8 × 10⁻⁴, and x(C⁺) = 3.6 × 10⁻⁴. We first study H₂ formation inthe absence of dust (bottom part of the Fig. 3 with boundaries ①, ②, and ③) and then the influence of a nonvanishing dust fraction (bulk of the diagram with boundaries ④, ⑤, ⑥, and ⑦). To generalize our results obtained here for a single set of values of T_K, G₀, and W, we also propose in Appendix C a numerically validated analytical approach that provides expressions for each boundary and the associated formation rates of H₂ for any physical condition.

3.1.1 Dust-free

Figure 4 shows that in the absence of dust and for a radiation field characterized by G₀ = 10⁴ and W = 5 × 10⁻⁷, the gas remains atomic up to n_H = 3 × 10¹² cm⁻³. The abundance of H₂ increases from ≃10⁻⁹ up to ≃ 0.2. Over the explored range of density, H₂ is preferentially photodissociated by the unshielded radiation field. Since the photodissociation rate does not depend on the density, the global trend seen in H₂ is mostly due to the formation routes.

We also plotted an analytical model of the abundance of H₂ assuming formation by H⁻ only and destruction by photodissociation (see Appendix C, Eq. (C.11)). The analytic expression reproduces the global increase of H₂ from n_H = 3 × 10⁶ to 3 × 10¹¹ cm⁻³ very well. In this regime, formation by H⁻ is the dominant formation route of H₂ (boundaries ③ to ①, Fig. 3). The route via H⁻ being a catalytic process by electrons, its efficiency depends linearly on the electron fraction. The recombination of ions at high density reduces the electron fraction and thus, the efficiency of this route. Below n_H ≃ 2 × 10⁸ cm⁻³ (boundary ②), photodetachment of H⁻ takes over from the reaction H⁻ + H → H₂ + e⁻, leading to a decrease of H⁻ and limiting the formation rate of H₂. Our analytical approach generalizes this result to any density and radiation field and shows that for a diluted black-body at 4000 K, this transition appears for (boundary ② and Appendix C, Eq. (C.18)) $$\frac{n_{\text{H}}}{W}=4.6\times {10}^{14} {\text{cm}}^{-3}.$$(9)

Despite the photodestruction of H⁻, formation by H⁻ remains the dominant formation route of H₂ down to n_H = 3 × 10⁶ cm⁻³.

Below this value (boundary ③), the analytical model underestimates the H₂ abundance. In this regime, formation via CH⁺ takes over from formation by H⁻. This is due to a quenching of the H⁻ route caused by a rapid photodetachment of H⁻. For example, at the boundary ③, only ~1% of H⁻ that formed by radiative attachment is actually converted in H₂. More generally, this transition appears for (boundary ③, Appendix C, Eq. (C.20)) $$\begin{array}{ll}{n}_{\text{H}}/W=\hfill & 6.7\times {10}^{12}{\text{cm}}^{-3}\left(\frac{x({\text{C}}^{+})}{3.6\times {10}^{-4}}\right){\left(\frac{x_e}{4.8\times {10}^{-4}}\right)}^{-1}\hfill \\ \hfill & \times {\left(\frac{T_{\text{K}}}{1000 \text{K}}\right)}^{-1.32},\hfill \end{array}$$(10)

where x(C⁺) is the abundance of C⁺. Interestingly, the very low radiative association rate of S⁺ with H prevents this already rare element in participating significantly in the formation of H₂ via SH⁺. Formation via H${}_2^{+}$ is found to be negligible over the explored parameter range. As seen in Fig. 4, H${}_2^{+}$ is always at least two orders of magnitude less abundant than H⁻, CH⁺, or SiH⁺ and it does not form H₂ at similar levels.

Above n_H = 3 × 10¹¹ cm⁻³, the analytic model also underestimates the abundance of H₂. In this regime, the three-body reaction route takes over from the formation by H⁻. When H⁻ is not photodetached, this transition appears for (boundary ① and Appendix C, Eq. (C.22)) $$n_{\text{H}}=1.9\times {10}^{13} {\text{cm}}^{-3} \left(\frac{x_e}{4.8\times {10}^{-4}}\right){\left(\frac{T_{\text{K}}}{1000 \text{K}}\right)}^{1.24}.$$(11)

Fig. 3 Schematic view of the dominant H2 formation routes depending on the density, the dust fraction, and visible radiation field W summarizing our results presented in Sect. 3.1 and in Appendix C. The location of the boundaries are given for a temperature of TK = 1000 K, xe = 4.8 × 10−4, and x(C+) = 3.6 × 10−4. The schematic view remains valid from ≃100 K up to ≃ 5000 K. Dependencies on TK, xe, and x(C+) are given inAppendix C. We note that some limits depend on the visible flux W and others do not. Depending on the visible flux, boundaries ①, ②, and ③ can merge.

Fig. 4 Steady-state abundances relative to total H nuclei for H2 and chemical species involved in its formation for dust-free single point models with G0 = 104, W = 5 × 10−7, TK = 1000 K, and nH ranging from105 to 2 × 1012 cm−3. An analytical expression of the steady state abundance of H2, assuming destruction by photodissociation and formation by H− only, is alsoplotted in black dotted line (see Appendix C). Boundaries defined in Fig. 3 are also indicated on the upper axis. It is important to note that because of the decrease of the electron fraction following recombination at high density, a three-body reaction takes over from the formation via H− (boundary ①) at a lower density than indicated in Fig. 3.

3.1.2 Dust-poor

Since these local models assume no attenuating material to the source, the inclusion of dust does not significantly affect the efficiency of gas-phase formation routes for H₂. Consequently, our previous results on dust-free chemistry still hold. The effect of the increasing dust fraction Q∕Q_ref is to add a new formation route that can compete with gas-phase formation routes. The critical amount of dust above which formation on a dust grain takes over depends on the efficiency of the dust-free formation route. Figure 5 shows the variation of the abundance of H₂ as a function of the dust fraction Q∕Q_ref for the fiducial radiation field as well as for n_H = 10¹⁰ cm⁻³ and n_H = 10⁵ cm⁻³. As shown above, gas-phase formation routes are dominated by H⁻ in the first case and by CH⁺ in the latter.

Below Q∕Q_ref ≃ 10⁻³ and for n_H = 10¹⁰ cm⁻³, the H₂ abundance is independent of the gas-to-dust ratio Q and is equal to its dust-free value. In this regime, H₂ is formed through H⁻ with a maximal efficiency (H⁻ is not photodetached) and the formation on dust is negligible. At about Q∕Q_ref ≃ 2 × 10⁻³, formation on dust takes over from gas-phase formation, driving up the H₂ abundance. A more detailed analysis (see Appendix C) shows that when H⁻ is the main formation route and not photodetached, the critical dust fraction below which the H₂ gas-phase formation takes over from formation on grains is (Fig. 3, boundary ⑥) $$Q/Q_{\text{ref}}=9.3\times {10}^{-3}{\left(\frac{T_{\text{K}}}{1000 \text{K}}\right)}^{0.4}{\left(\frac{S(T_{\text{K}})}{S(1000 \text{K})}\right)}^{-1}\left(\frac{x_e}{4.8 {10}^{-4}}\right),$$(12)

where S(T_K) is the sticking coefficient of H on grains adopted from Hollenbach & McKee (1979). When the formation through H⁻ is reduced by photodetachment, this critical dust fraction is proportional to n_H ∕W (Fig. 3, boundary ⑤ and Eq. (C.25)).

Figure 5 shows that for a lower density-to-visible field ratio (n_H ∕W = 2 × 10¹¹ cm⁻³), dust has a significant impact at a much smaller dust fraction. As seen in the previous section, the gas-phase formation of H₂ is then dominated by CH⁺. Formation by CH⁺ being about two orders of magnitude less efficient than electron catalysis, the critical dust fraction above which formation on grains takes over from gas-phase formation is lowered by a similar factor accordingly (Fig. 3, boundary ④ and Eq. (C.26)).

Fig. 5 Steady-state abundances of H2 normalized to its value at Q∕Qref = 1 as a functionof Q∕Qref for two densities: nH =1010 cm−3 (solid red line) and nH = 105 cm−3 (solid blue line). Other parameters are constant and equal to their fiducial values (see Table 2). Gas-phase formation of H2 is dominatedby H− for nH = 1010  cm−3 (nH∕W = 2 × 1016 cm−3) and by CH+ for nH = 105  cm−3 (nH∕W = 2 × 1011 cm−3). Dashed linesindicate H2 abundance in the absence of dust. Abundances for each set of the model are normalized to their value at Q∕Qref = 1. We note that because of the low electron fraction at nH = 1010 cm−3 (xe ≃ 10−4), the formation on grains takes over from the formation via H− (boundary ⑥) at a lower dust fraction than indicated in Fig. 3.

3.2 Other molecules

H₂, even when scarce, constitutes the precursor of other molecules, such as CO, SiO, OH, and H₂O, that they are observed in protostellar jets. The molecular richness then depends on the abundance of H₂. As such, the inclusion of dust increases molecular abundances only by increasing the fraction of H₂. However, other physical variables can regulate molecular abundances, such as the temperature, the FUV radiation field, and the density. Since molecules are essentially formed by two-body reactions and destroyed by photodissociation in the UV domain, molecular abundances mostly depend on the ratio n_H ∕G₀. In this section, results on molecular abundances obtained by varying n_H for fixed G₀ = 10⁴ can thus be generalized for any G₀ by translating them to the ratio n_H∕G₀.

Fig. 6 Steady-state abundances relative to total H nuclei for relevant molecular species from single point models in the absence of dust. Panel a: abundances as a function of temperature for G0 = 104, W = 5 × 10−7, and nH = 109 cm−3, panel b: abundances as a function of nH (lower axis) and nH∕G0 (upper axis) for TK = 1000 K, G0 = 104, and W = 5 × 10−7.

3.2.1 Dust-free

Figure 6a shows that steady-state abundances of OH, CO, H₂O, and SiO increase by several orders of magnitude with temperature and reach maximum abundances above ≃ 1000 K. This trend is due to the activation of endothermic gas-phase formation routes at a high temperature. The formation of CO, SiO, and H₂O is indeed initiated by the formation of OH through the neutral-neutral reaction $$\text{O}+{\text{H}}_2\to \text{OH}+\text{H} \Delta E=+2980 \text{K}.$$(13)

This warm route involving H₂ is fundamental for the formation of all of the considered molecules, even if the H₂ fraction is low.

Figure 6b shows steady-state abundances as a function of density for a temperature of T_K = 1000 K, which is sufficient for efficient molecule formation. As H₂ and OH increase with n_H, the abundances of the other species increase as well. For n_H ≳ 3 × 10¹⁰ cm⁻³, the gas is essentially atomic but rich in CO, SiO, and H₂O. This is one of the most fundamentaland unique characteristics of dust-free chemistry. Still, due to complex formation and destruction pathways, each molecular species behaves differently as a function of n_H, revealing their formation and destruction routes in H₂ poor gas. Figure 7 summarizes the dominant reactions contributing to the formation and destruction of CO, H₂O, and SiO, which depend on the ionization state of C and Si.

The CO formation pathway is essentially regulated by the ionization state of the carbon. Below n_H = 2 × 10⁹ cm⁻³, carbon is ionized and CO is produced via the ion-neutral reaction C⁺ + OH, producing either CO⁺ or CO directly. In the former, CO⁺ is neutralized by a fast charge exchange with H. Above n_H ≃ 2 × 10⁹ cm⁻³, CO is produced directly by the neutral-neutral reaction C+OH → CO. Destruction occurs mostly via photodissociation.

The H₂O abundance exhibits a stiff increase with n_H. Over the explored parameter range, H₂O is produced through the neutral-neutral reaction OH + H₂. Below n_H ≃ 10¹⁰ cm⁻³, destruction occurs via photodissociation and via a reaction with C⁺. At a higherdensity, the main destruction route is via the reverse reaction H₂O + H → OH + H₂.

The SiO abundance exhibits a stiff increase around n_H ≃ 3 × 10¹⁰ cm⁻³, rising by more than four orders of magnitudes in one decade of n_H. This feature is due to a change in both destruction and formation routes. Below n_H ≃ 3 × 10¹⁰ cm⁻³, Si⁺ is the main silicon carrier and SiO synthesis is initiated by the ion-neutral reaction Si⁺ + OH → SiO⁺. However, in contrast to the analogous reaction with CO⁺, SiO⁺ cannot be neutralized through a charge exchange with the main collider, namely H. Alternatively, SiO⁺ decays toward SiO through SiO⁺ + H₂ → SiOH⁺, eventually leading to SiO by dissociative recombination. H₂ is rare in the absence of dust and thus this formation route is much less efficient than the analogous reaction that forms CO. In addition, at low n_H ∕G₀, the abundant C⁺ efficiently destroys SiO to form CO. Consequently when carbon and silicon are ionized, the medium is hostile to the formation and the survival of SiO. On the contrary, above n_H ≃ 3 × 10¹⁰ cm⁻³, SiO is directly formed through the neutral-neutral reaction Si + OH, and destruction by C⁺ is quenched by the recombination of carbon. In contrast to the analogous reaction with H₂O, the reverse reaction SiO + H has a very high endothermicity (3.84 eV) that prevents any destruction by H. Given these favorable factors, SiO becomes the main silicon reservoir above n_H ≃ 3 × 10¹⁰ cm⁻³.

Fig. 7 Dominant reactions controlling the abundance of CO, H2O, and SiO under warm (TK ≥ 800 K) and irradiated conditions. OH appears to be a key intermediate for the three species. The ionization state of carbon controls the destruction of SiO and H2O as well as the formation of CO. The ionization state of silicon controls the formation of SiO. For an unshielded ISRF FUV radiation field, carbon and silicon are ionized for nH ∕G0 ≲ 105 cm−3 and nH ∕G0 ≲ 3 × 106 cm−3, respectively.

Fig. 8 Wind model adopted in this work. Panel a: schematic view of the geometry of the model. Streamlines are assumed to be straight lines launched from the disk (see Fig. 1a). The wind velocity Vj is constant and equal to 100 km s−1. The wind launching region extends from Rin = 0.05 au out to Rout = 0.3 au. We focus on the chemical evolution of a representative streamline launched from 0.15 au in the disk and reduce the problem to 1D (see Sect. 4.1). Panel b: prescribed density nH and unshielded FUV flux G0 profile along the representative streamline launched from 0.15 au for wind solution ⓐ (see Table 3). For other models, nH and G0 are simply rescaled according to the Eqs. (19) and (17), respectively. Panel c: prescribed nH ∕G0 and nH ∕W ratio along the same streamline. We note that due to the collimation of the flow (z0 ≫ R0), these ratios increase with distance by a factor of ${(z_0/R_0)}^2=625$. In the absence of dust, the opacity of the gas in the visible is negligible so that the nH ∕W ratio is expected to be the true local ratio between the density and visible field.

3.2.2 Dust-poor

The inclusion of dust increases the molecular abundances by increasing the H₂ abundance. Thus, the minimal amount of dust required to affect the chemistry of CO, OH, SiO, and H₂O is similar to that determined in the Sect. 3.1.2. Molecular abundances are then increased accordingly, but specific formation and destruction routes remain the same.

3.3 Summary

In this section, the chemistry of dust-free and dust-poor gas is studied. Regarding H₂, we show that the dominant gas-phase formation route and its efficiency critically depends on the ratiobetween the density and the visible radiation field, which is quantified here by n_H ∕W. Above n_H ∕W ≃ 5 × 10¹⁴ cm⁻³, H₂ is efficiently formed via H⁻; whereas below this value, photodetachment reduces its efficiency. When optimally formed via H⁻, a minimum fraction of dust of about Q∕Q_ref ≃ 2 × 10⁻³ is required to have a significant impact on the chemistry.

Regarding CO, OH, SiO, and H₂O, we find that high abundances are reached for n_H ∕G₀ ≥ 10⁶ cm⁻³, despite low abundances of H₂. Efficient formation routes are initiated by OH and require a warm environment (T_K ≥ 800 K). The inclusion of dust increases molecular abundances by increasing the H₂ abundance accordingly. We also find that the abundance of SiO is very sensitive to the ionization state of carbon and silicon. When both species are singly ionized, the SiO abundance is very low due to both destruction by C⁺ and a very low efficiency of the formation by Si⁺ in an H₂-poor environment.

4 1D wind models

In this section, the chemistry studied in the previous section is incorporated in a more realistic model of 1D wind streamline. In addition to the quantities found to control the chemistry in unattenuated static environments (namely n_H ∕G₀, n_H ∕W, T_K, Q∕Q_ref, and T_vis), we include three ingredients: the attenuation of the radiation field with the distance from the source, the time-dependent chemistry, and the differential geometrical dilutions of the density and the radiation field. Below, we present our simple 1D model where these effects are implemented with a minimal number of free parameters, allowing one to investigate a wide range of source evolutionary phases.

4.1 The streamline model

The disk wind model, which is illustrated in Fig. 8a, is built from a simple flow geometry that captures the essentialproperties of MHD disk wind models in a parametric approach without relying on a peculiar wind solution. Following Kurosawa et al. (2006), we assume that the wind is launched from a region of the disk between R_in and R_out, and it propagates along straight streamlines diverging from a point located at a distance − z₀ below the central object (Fig. 8a).

We follow the evolution of only one representative streamline launched from R₀ by using the astrochemical model presented in Sect. 2. To reduce the number of free parameters, the wind velocity, noted V_j, is assumed to be constant with distance. The conservation of mass then yields a geometrical dilution of density along the streamline anchored at R₀ of $$n_{\text{H}}(z)=n_H^0\frac{1}{{\left(1+\frac{z}{z_0}\right)}^2},$$(14)

we note as $n_H^0$ the density at the base (z = 0).

The radiation field is assumed to be emitted isotropically from the star position. Along a given streamline, it is reduced by a geometrical dilution factor and attenuated by gas-phase species and dust, if any. To keep the problem tractable in 1D, the attenuation of the radiation field is assumed to proceed along each streamline. The geometrically diluted, unattenuated FUV field at position z is given by $$G_0(z)=G_0^0\frac{1}{{\left(1+\frac{z}{z_0}\right)}^2+{\left(\frac{z}{R_0}\right)}^2},$$(15)

we note as $G_0^0$ the unattenuated FUV field at the base (z = 0). The same geometrical dilution applies to the unattenuated visible radiation field.

Figures 8b,c show that under this simple wind geometry, the radiation field is diluted on a spatial scale ≃ R₀, while the density field is diluted on a scale z₀. Due to the collimation of the flow (z₀ ≫ R₀), the ratios n_H∕G₀ and n_H∕W, which would determine the flow chemistry in the absence of attenuation and time-dependent effects, increase with distance up to a factor of $1+{(z_0/R_0)}^2$ from their initial values at z = 0. Hence, the differential dilutions of the radiation and density fields in our model is regulated by the wind collimation angle (z₀ ∕R₀).

Our 1D chemical wind-model is thus controlled, in principle, by nine free parameters; six of which are the same as in Sect. 3, namely $n_H^0$, $G_0^0$, W⁰, T_K, Q∕Q_ref, and T_vis. The remaining three parameters that control the wind attenuation, dilution, and nonequilibrium effects are as follows: the opening angle of the streamline z₀ ∕R₀, the density dilution scale z₀ (or alternatively the anchor radius R₀), and the velocity of the wind V_j. Since we cannot explore the full parameter space in the present study, here, we chose to fix most of them to representative values observed in molecularjets and young protostars, and we focus on varying the source evolutionary status. Namely, we take a fixed launch radius R₀ = 0.15 au (as estimated for the SiO jet in HH 212 by Tabone et al. 2017). The ratio z₀ ∕R₀ is fixed to 25, leading to an opening angle of the computed streamline of 5°, which is in line with the observed universal collimation properties of jets across ages (Cabrit 2007a; Cabrit et al. 2007). The wind velocity V_j is taken equal to 100 km s⁻¹, which is in line with typical velocities measured from proper motions toward Class 0 molecular jets (e.g., Lee et al. 2015). The flow temperature is taken as T_K = 1000 K (see discussion in Sect. 5.2).

The radiation field emitted by the accreting protostar is modeled as a black body of photospheric origin with T_vis = 4000 K, plus a FUV component coming from the accretion shock onto the stellar surface. At the base of the streamline, the dilution factor of the stellar black body, defined as $W\equiv \frac{J_{\nu }}{B_{\nu }}$, is (assuming R₀ ≫ R_*) $$W^0=\frac{1}{4}{\left(\frac{R_{*}}{R_0}\right)}^2=2.2\times {10}^{-3}{\left(\frac{R_0}{0.15 \text{au}}\right)}^{-2}{\left(\frac{R_{\star }}{3R_{\odot }}\right)}^2.$$(16)

The radius of the protostar is fixed to R_⋆ = 3 R_⊙, where R_⊙ is the solar radius. Regarding UV excess, FUV observations of BP Tau and TW Hya show that an ISRF provides a good proxy for the shape of the radiation field, though neglecting the line contribution to the FUV flux (> 35%, Bergin et al. 2003). Here, we assume that the FUV excess follows a Mathis radiation field (Mathis et al. 1983, see Appendix A) and we neglect the UV line emission. For BP Tau, Bergin et al. (2003) found that a FUV flux of G₀ = 560 at 100 au is required to match the FUV continuum level. Assuming that the FUV flux scales with the accretion luminosity $L_{\text{acc}}=\frac{G{M}_{*}{\dot{M}}_{\text{acc}}}{R_{*}}$, which is 0.24 L_⊙ for BP Tau, the scaling factor at the base of the streamline anchored at R₀ is $$G_0^0=1.1\times {10}^{10}\left(\frac{{\dot{M}}_{\text{acc}}}{{10}^{-5}{M}_{\odot } {\text{yr}}^{-1}}\right)\left(\frac{M_{\star }}{0.1M_{\odot }}\right){\left(\frac{R_{\star }}{3R_{\odot }}\right)}^{-1}{\left(\frac{R_0}{0.15 \text{au}}\right)}^{-2},$$(17)

where Ṁ_acc is the accretion rate onto the protostar and M_⋆ is the mass of the protostar.

In order to relate the density at the launching point of a streamline $n_{\text{H}}^0$ to the mass-loss rate of the wind Ṁ_w, we assume that between R_in and R_out, the wind has a constant local mass-loss rate. This gives a density structure at the base of the wind of $$n_{\text{H}}{}^0=\frac{0.5{\dot{M}}_{\text{w}}}{2\pi 1.4m_{\text{H}}{V}_{\text{j}}{R}_0\left(R_{\text{out}}-R_{\text{in}}\right)},$$(18)

where Ṁ_w is the (two-sided) mass-loss rate of the wind. It is important to note that for winds launched from a narrow region of the disk (R_out ≃ few R_in), the choice of the power-law index of the local mass-loss rate has a weak influence on $n_{\text{H}}{}^0$. In this work, we follow the modeling results of the SiO jet in HH212 (Tabone et al. 2017) and we fix R_in = 0.05 au and R_out = 0.3 au, leading to a density at the base of the streamlines of $$\begin{array}{ll}{n}_{\text{H}}{}^0=2.5\times {10}^{10}\hfill & \frac{{\dot{M}}_{\text{w}}}{{10}^{-6}{M}_{\odot } {\text{yr}}^{-1}}\left(\frac{0.25 \text{au}}{R_{\text{out}}-R_{\text{in}}}\right)\hfill \\ \hfill & {\left(\frac{V_j}{100 \text{km} {\text{s}}^{-1}}\right)}^{-1}\left(\frac{0.15 \text{au}}{R_0}\right){\text{cm}}^{-3}.\hfill \end{array}$$(19)

To further reduce the parameter space, we follow the universal correlation between accretion and ejection observed from Class 0 to Class II jets and set the (two-sided) wind mass-flux to Ṁ_w = 0.1 Ṁ_acc.

In the end, the parameter space in the present study is thus reduced to only three free parameters: Ṁ_acc, M_⋆, and Q∕Q_ref. To investigate how the chemical content of a dust poor, laminar jet evolves in time with the decline of the accretion rate and the increase in stellar mass, we computed six sets of parameters summarized in Table 3. First, four models representative of Class 0 winds were computed. At this stage, the young embedded source has not reached its final mass yet and we chose M_* = 0.1 M_⊙. The accretion rate is varied from Ṁ_acc = 5 × 10⁻⁵ M_⊙ yr⁻¹ to 5 × 10⁻⁶ M_⊙ yr⁻¹ to model sources with various accretion luminosities. A model representative of a Class I source was also computed. At this stage, the protostar has accumulated most of its mass and we chose M_* = 0.5 M_⊙ with an accretion rate of Ṁ_acc = 10⁻⁶ M_⊙ yr⁻¹. Lastly, a model representative of a Class II source with an accretion rate of Ṁ_acc = 10⁻⁷ M_⊙ yr⁻¹ and M_* = 0.5 M_⊙ was computed. This accretion rate is representative for an actively accreting TTauri star. To study the impact of surviving dust, the six wind models were computed for the dust fraction Q∕Q_ref = 0, 10⁻³, 10⁻², and 10⁻¹.

The initial chemical abundances at the base of the streamline are computed as in Sect. 3, that is, assuming chemical equilibrium and no attenuation of the radiation field. Hence, they only depend on the adopted T_K = 1000 K, T_vis = 4000 K, and the initial ratios $n_{\text{H}}{}^0/G_0^0$ and $n_{\text{H}}{}^0/W^0$ given in Table 3 for all models.

4.2 Results: dust-free winds

In this section, we present the results on dust-free wind models (Q = 0). Chemical abundances and the local FUV radiation field of the selected wind models ⓐ, ⓑ, and ⓒ (see Table 3) are presented in Fig. 9. These specific models allow one to highlight the influence of the wind parameters, namely the density of the wind (e.g., mass-loss rate) and the radiation field. Models ⓐ and ⓑ have a different n_H and G₀, but they share the same n_H∕G₀ ratio. Models ⓑ and ⓒ have a similar FUV radiation field but a different density (Table 3). Asymptotic abundances for the full set of dust-free wind models are also given in Fig. 10 (solid line).

4.2.1 FUV field

One of the main difference between single point models and wind models is the inclusion of the attenuation of the radiation field along streamlines. In dust-free winds, only gas-phase species can shield the radiation field and decrease photodissociation rates of molecular species.

Figure 9 (right panels) shows that the attenuation of the radiation field results in sharp absorption patterns that are characterized by thresholds below which the radiation field is heavily extincted. Those thresholds correspond to ionization thresholds of the most abundant atomic species. This specific attenuation pattern has already been pointed out by Glassgold et al. (1989, 1991) in the context of stellar winds, although without focusing on the sharpness of the attenuation patterns that are at the root of the chemical richness of dust-free jets. The specific shielding mechanism that causes these unique attenuation patterns can be understood by the inspection of the model ⓑ (Fig. 9b). In the absence of attenuation, carbon is expected to be ionized. However, at z ≃ 2 au, the abundance of C⁺ drops by several orders of magnitude (right panel) and neutral carbon, which is not shown here, becomes the main carbon carrier. Figure 9b, right panel, shows that above this transition, photons below the photoionization threshold of carbon (λ ≤ 1100 Å) are attenuated by more than 8 orders of magnitude. In this region of the spectrum, the opacity of the gas is dominated by carbon. The steep decrease of C⁺ is thus due to the attenuation by carbon itself, triggering a C⁺/C transition when the gas becomes optically thick to ionizing photons. Because of the increase of C at the C⁺ /C transition, the local opacity of the gas increases even more. This leads to an attenuation of the FUV field below λ = 1100 Å that is much stiffer than attenuation by dust. This process, which is called continuum self-shielding and included in most of the dusty PDR models (Röllig et al. 2007), turns out to be of paramount importance in dust-free and dust-poor winds.

The resulting attenuation of the FUV radiation field is very sensitive to the wind model. For model ⓐ, with a higher mass-loss rate, the radiation field at z = 1000 au is strongly attenuated down to the ionization wavelength threshold of sulfur, which has a relatively large value (λ = 1600 Å); whereas for model ⓒ, with a lower mass-loss rate, the radiation field is barely extincted across the FUV spectrum. As any self-shielding process, it depends on the column density of the neutral atom X. In the inner ionized and unattenuated part of the jet, this column density increases with z as $$N_{\text{X}}(z)=n_{\text{X}}^0{z}_0{\displaystyle {\int }_0^{z/z_0}\frac{{(1+u)}^2+{(\frac{z_0}{R_0}u)}^2}{{(1+u)}^4}}\text{d}u,$$(20)

where $n_{\text{X}}^0$ is the density of X at the launching point and is proportional to ${(n_{\text{H}}{}^0)}^2/G_0^0$. The self shielding occurs if N_X(z) becomes typically larger than ${\sigma }_{\text{X}}^{-1}$, where σ_X is the FUV absorption cross section⁴ of atom X. Because the integral term on the right hand side converges, Eq. (20) reveals a threshold effect: If ${(n_{\text{H}}{}^0)}^2/G_0^0$ is too low, N_X(z) is found to never rise above the critical value required to trigger the self-shielding, regardless of z.

Models ⓐ and ⓑ share the same unshielded $n_{\text{H}}{}^0/G_0^0$ ratio and consequently exhibit similar atomic abundances at the base of the wind. However, as model ⓐ is denser, its ${(n_{\text{H}}{}^0)}^2/G_0^0$ ratio is larger. This increases the total column density of S and Si, which exhibit self-shielding transitions. This results in attenuation of the radiation field at much longer wavelengths. In contrast, as model ⓒ has both a lower $n_{\text{H}}{}^0/G_0^0$ ratio and a lower density, the column density of carbon is not sufficient in attenuating the radiation field and the shape of the radiation field remains unaltered.

In other words, in the absence of dust, the radiation field is attenuated by a continuum self-shielding of atoms. This process is very efficient at attenuating the radiation field below the photoionization thresholds of the atomic species. The wavelength below whichthe radiation field is attenuated depends on the species that are self-shielded. Carbon can attenuate the radiation field at a short wavelength, whereas S and Si attenuate the radiation field at a longer wavelength. Self-shielding by a specific species is very sensitive to the density of the wind. The result is that Class II and I models are not dense enough to be shielded by any atom, whereas Class 0 models are shielded by carbon for Ṁ_w ≥ 5 × 10⁻⁷M_⊙ yr⁻¹ and by carbon, silicon, and sulfur for Ṁ_w ≥ 2 × 10⁻⁶M_⊙ yr⁻¹.

Fig. 9 Computed chemical abundances and local FUV radiation field for dust-free isothermal wind models for TK = 1000 K. Left panels: chemical abundances relative to total H nuclei. Right panels: local mean intensity of the FUV radiation field at various positions along the wind (position indicated on the curves). Photoionization thresholds of C, S, Si, Mg, and Fe are also indicated by gray, yellow, orange, black, and blue dashed lines, respectively. Row a: Class 0 model with Ṁw = 2 × 10−6 M⊙ yr−1 and M* = 0.1 M⊙. Row b: Class 0 model with a lower mass-loss rate Ṁw = 10−6 M⊙ yr−1 and same mass. Row c: Class I model, with lower accretion rate Ṁw = 10−7   M⊙ yr−1 but higher mass M* = 0.5 M⊙.

Fig. 10 Abundances at z > 1000 au for a streamline anchored at 0.15 au in the disk as a function of the mass-loss rate for various dust fractions. For Ṁw ≥ 5 × 10−7 M⊙ yr−1, the mass of the central object is 0.1 M⊙ (Class 0 model) and 0.5 M⊙ for lower mass-loss rates (Class I and II models). Dust-free models are plotted in solid lines, and dust-poor models with Q∕Qref = 10−3, 10−2, and 0.1 are plotted in dashed, dashed-dotted, and dotted lines as indicated in each panel, respectively. Horizontal black dashed lines indicate the elemental abundance of carbon (panels on CO and C), silicon (panel on SiO), and oxygen (panel on H2O).

4.2.2 H₂

Figures 9 (left panels) and 10 show that dust-free wind models are poor in H₂ and consequently they are mostly atomic. H₂ is also smoothly increasing with an increasing mass-loss rate. For Class 0 wind models, typical values of 10⁻³ are found above z = 10 au. This global trend with a mass-loss rate is due to both an increase in the formation efficiency and a decrease in destruction efficiency.

Regarding formation pathways, the formation of H₂ is dominated by H⁻ for all models, except at the base of the wind where H₂ is formed by CH⁺ in Class I and II models. As seen in Sect. 3, the efficiency of the H⁻ route depends on the ratio n_H∕W that quantifies the ability of H⁻ to form H₂ instead of being photodestroyed by the visible field. All models share the same visible radiation field so that the ratio n_H ∕W only depends on the density of the wind (i.e., mass-loss rate). Consequently, the efficiency of the formation of H₂ by H⁻ increases progressively with the mass-loss rate.

Along with the increase of the formation efficiency, the destruction efficiency decreases with the mass-loss rate. For Class II models, destruction of H₂ are dominated by photodissociation. Class I models exhibit a sufficiently high column density of H₂ to insure anefficient self-shielding, quenching photodissociation route. Alternatively, destruction by C⁺ according to C⁺ + H₂ → CH⁺ → C⁺ + H takes over from photodissociation with a smaller efficiency. For even larger mass-loss rates, the self-shielding of atomic carbon leads to a drop of C⁺, quenching the former destruction route. Alternatively, H₂ is destroyed by atomic oxygen according to O + H₂ → OH → O + H.

In other words, dust-free wind models are mostly atomic. The H₂ abundance increases with the mass-loss rate due to an increase in the efficiency of the formation route by H⁻ and a decrease of the destruction route with the progressive self-shielding of H₂ and the recombination of C⁺.

4.2.3 Other molecules

Figure 10 shows that along with the increase of H₂ with the accretion rate, molecular abundances of interest increase with the mass-loss rate. As shown is Sect. 3, H₂ abundance regulates the formation of OH, CO, SiO, and H₂O. On the other hand, continuum self-shielding of atomic species is a crucial process for the survival of these molecules: it reduces photodissociation rates by attenuating the FUV radiation field and quenches destruction routes by C⁺. However, the precise impact of the shielding by atoms depends on the photodissociation threshold of each molecular species relative to the threshold below which the radiation field is attenuated.

OH is the first neutral species to be formed after H₂, and it is an important intermediate for the synthesis of other molecules. Figure 10 shows that the abundance of OH increases smoothly with the mass-loss rate following the increase of the abundance of H₂. Its abundance remains low, reaching 6 × 10⁻⁷ for model ⓐ. The wavelength dissociation threshold of OH is larger than the ionization thresholds of C, S, Si, Fe, and Mg, and thus OH is not efficiently shielded by those species. Consequently, OH is destroyed by photodissociation at all Ṁ_w and the increase of OH with Ṁ_w is mostly driven by the smooth increase of H₂ and the increase of formation rates with density. Interestingly, at the highest mass-loss rate, the abundance of OH is also limited by the reverse reaction OH + H → O + H₂.

CO is found to exhibit a much steeper increase with Ṁ_w from Class II-I models to Class 0 models (Fig. 10). For Class 0 models, the CO abundance is high at ≥ 10⁻⁵. Interestingly, the abundance ratio CO/H₂ is ≃ 2−5 × 10⁻² in Class 0 models. This value is much larger than the canonical value derived in dusty molecular gas (≃ 10⁻⁴) and constitutes one of the most striking characteristics of dust-free jets. The global behavior of CO is related to the shielding mechanism of the wind by carbon. In contrast with OH, the CO wavelength dissociation threshold lays below the wavelength ionization threshold of carbon. CO can thus be efficiently shielded by carbon when the column density of C is sufficient in triggering the self-shielding of C. This process is notable in both Class 0 models presented in Fig. 9 where a jump in CO by several orders of magnitude is seen across the C⁺ /C transition. Since the self-shielding of carbon only operates in Class 0, CO is only abundant in dust-free Class 0 models.

SiO exhibits an even steeper increase with the mass-loss rate and it is only abundant for the highest mass-loss rates (Ṁ_w ≥ 2 × 10⁻⁷M_⊙ yr⁻¹). Such a high sensitivity of SiO to the mass-loss rate is related to the self-shielding of silicon that ultimately controls the formation and destruction of SiO. Regarding destruction, SiO has an intermediate wavelength dissociation thresholds (λ_th = 1500 Å) that lies between C and Si ionization thresholds. As for OH, shielding by C does not significantly reduce photodissociation rates but in contrast to OH, shielding by S and Si is very efficient. As a consequence, the SiO photodissociation rate falls from model ⓑ to model ⓐ. Regarding formation routes, as shown in Sect. 3 and Fig. 7, SiO formation is more efficient when Si is in neutral form. Thus, the self-shielding of Si in model ⓐ activates a direct and efficient formation route of SiO by Si + OH → SiO + H. On the contrary, in models with a lower mass-loss rate, SiO synthesis proceeds via the much less efficient Si⁺ route.

The H₂O abundance is found to be low, reaching 10⁻⁶ for the highest mass-loss rate. As for OH, the H₂O wavelength dissociation threshold is longer than the ionization threshold of C, S, Si, and even Fe or Mg. H₂O is not efficiently shielded by the gas and its increase with the mass-loss rate is mostly due to the increase of the H₂ abundance. At the highest mass-loss rates (Ṁ_w ≥ 2 × 10⁻⁶ M_⊙ yr⁻¹), H₂O is also destroyed through the reverse reaction H₂O + H → OH + H₂.

4.3 Results: Dust-poor winds

The inclusion of a nonvanishing fraction of dust activates H₂ formation on dust (see Sect. 3), and it introduces a new source of opacity for the radiation field. Figure 10 summarizes the influence of the dust fraction Q∕Q_ref on the asymptotic molecular abundances for the full set of wind models.

The inclusion of a small fraction of dust Q∕Q_ref = 10⁻³ increases H₂ abundances by a factor ~2 (Fig. 10a). Over the explored range of wind mass-loss rates, this specific dust fraction is indeed close to the critical value above which H₂ formation on grains takes over gas phase formation (see Sect. 3). The chemistry of H₂ is consequently weakly affected. Similarly, the CO abundance is also weakly affected. This is because the opacity of the gas below λ ≤ 1100 Å is still dominated by carbon over most wind models. Consequently, dust does not affect CO photodissocitation rates and the CO abundance only changes by a factor of ≃3. On the contrary, this small amount of dust has a strong impact on the abundance of SiO and H₂O. Attenuation of the radiation field by this small amount of dust above λ ≥ 1100 Å shields SiO and H₂O, decreasing photodissociation rates. The attenuation by dust also triggers self-shielding of atomic species for a lower mass-loss rate (Ṁ_w = 10⁻⁶ M_⊙ yr⁻¹), leading to an even stronger attenuation of the radiation field at long wavelengths. Furthermore, regarding SiO, the Si⁺ /Si transition activates the very efficient SiO formation route through Si. As a consequence, the critical mass-loss rate Ṁ_w above which the wind is rich in SiO is lowered.

For larger Q∕Q_ref ratios, the formation of H₂ on dust takes over from gas-phase formation. Enhanced abundances of H₂ increase the formation rate of molecules, whereas attenuation of the radiation field by dust reduces destruction rates. Figure 10 shows that CO, SiO, and H₂O abundances increase by several orders of magnitude with increasing Q∕Q_ref above the critical value Q∕Q_ref = 10⁻³. The overall impact of dust on CO and SiO content is to lower the critical mass-loss rate above which the gas is rich in CO and SiO. For example, while dust-free wind models are rich in CO for Ṁ_w ≥ 5 × 10⁻⁷ M_⊙ yr⁻¹ and rich in SiO for Ṁ_w ≥ 2 × 10⁻⁶ M_⊙ yr⁻¹, dusty-wind models with Q∕Q_ref = 10⁻² are rich in CO for Ṁ_w ≥ 10⁻⁷ M_⊙ yr⁻¹ and rich in SiO for Ṁ_w ≥ 5 × 10⁻⁷ M_⊙ yr⁻¹. Interestingly, above those critical mass-loss rates, CO and SiO constitute the main carbon and silicon carriers.

In contrast, the H₂O abundance has a more complex dependency on Q∕Q_ref and Ṁ_w. H₂O only reaches large abundances at high mass-loss rates and for a rather large fraction of dust (Q∕Q_ref ≥ 10⁻²). When the wind is sufficiently shielded by dust, destruction via the reverse reaction H₂O + H → H₂ + OH limits the abundance of H₂O. Being rich in atomic hydrogen, dust-free and dust-poor winds are hostile to the formation and survival of H₂O.

4.4 Time dependent chemistry

For all models, chemistry is found to be out-of-equilibrium, leading to asymptotic abundances that are somewhat smaller than steady-state abundances. As the density and the radiation field drop with z due to geometrical dilution and attenuation, the ratio between the chemical and the dynamical timescales increases as V_j ∕(z n_H(z)) for two-body reactions and as V_j∕(z F(z)) for photoreactions, where F(z) is the local FUV photon flux. The flow thus necessarily undergoes transitions beyond which part or all of the chemistry is “frozen”. In the models presented here, we find that these transitions occur around z ≃ z₀. We note, however, that out-of-equilibrium effects do not explain the global trend with the mass-loss rate (i.e., wind density) and dust fraction, and they only reduce the overall asymptotic abundances by a factor less than 4.

4.5 Summary

In this section, the chemistry of dust-free and dust-poor winds are investigated by the use of parameterized wind models that include time-dependent chemistry and the attenuation of the radiation field. Our results for warm wind models (T_K = 1000 K) are summarized in Fig. 10. The overall molecular content of wind models increases with the mass-loss rate of the wind and with the dust fraction. Dust-free and dust-poor winds are atomic; however, the small fraction of H₂, which formed essentially via H⁻, regulates the synthesis of other molecules. The survival of these molecules is insured by the attenuation of the radiation field by atomic species, and by dust, if any. The attenuation of the radiation field by atomic species proceeds through self-shielding, a process thatdepends critically on the column density of the absorbing species. It results in the presence of density or equivalently, mass-loss rate thresholds above which specific molecules are very abundant. Those thresholds depend on both the specific species and the dust fraction.

5 Discussion

Our results on disk wind chemistry have been obtained from a simple parametric wind model. The main advantage of this wind model is due to a simple geometry that allows for a deep exploration of the parameter space and a detailed treatment of the radiative transfer. A number of limitations regarding the model presented here have to be taken into account before comparing our results to observations.

5.1 Attenuation in the wind

The radiative transfer is reduced to a 1D geometry by assuming that the radiation field is attenuated along the computed streamline. In a 2D geometry, the shielding of a streamline is provided by inner streamlines of the wind. This is especially true at the base of the wind, where photons coming from the accreting central object are impinging the streamline almost transversely. However, for most of the models presented here, the decline of the radiation field at the base of the wind z < R₀ is mostly caused by geometrical dilution, an effect that is properly taken into account by our 1D model. Species that contribute to the attenuation of the radiation field are mostly formed around z ≃ 25 R₀. At this distance, FUV photons propagate almost parallel to the streamlines and our 1D approximation is valid.

Our model also computes the attenuation of the radiation field by a limited number of species. As shown is Sect. 4, attenuation by atomic species is a key process for the survival of molecules. Because of their long-wavelength dissociation thresholds, molecules, such as OH or H₂O, cannot be efficiently shielded by the considered atoms, namely C, Si, S, Mg, and Fe. In this context, other elements that exhibit longer-wavelength dissociation thresholds, such as Al, Ca, and Li, could shield important molecules. However, our models show that in the absence of dust, the attenuation by Mg and Fe is negligible for Ṁ_w ≤ 10⁻⁶ M_⊙ yr⁻¹. Thus, rarer elements that also have similar photoionization cross-sections and similar recombination rates are expected to have a negligible contribution to the opacity of the gas.

5.2 Thermal structure and shocks

In this exploratory study, we only consider isothermal laminar disk wind models. The computation of thermal-balance is beyond the scope of this paper. The isothermal assumption is motivated by detailed thermal balance calculations in self-similar MHD disk winds, showing that ambipolar diffusion is a robust heating mechanism that is able to balance adiabatic and radiative cooling, and it naturally leads to a rather flat temperature profile (Safier 1993; Garcia et al. 2001). Depending on the disk wind mass-loss rate, the asymptotic temperature for R₀ ≃ 0.2 au varies from 500 to 4000 K (Panoglou et al. 2012; Yvart et al. 2016), which is in line with our assumed fiducial temperature of 1000 K. Similar calculations for a nonself-similar “X-wind” from the inner disk edge (Shang et al. 2002) also show quasi-isothermal behavior on inner streamlines; hence it appears to be a generic property of MHD disk winds that undergo large-scale collimation. In contrast, models of stellar winds heated by ambipolar diffusion exhibit a decline of temperature at large distance (Ruden et al. 1990) due to the lack of magnetic collimation and the much smaller acceleration scale of a few R_⋆ (see discussion in Sect. 4.4.3 of Garcia et al. 2001). We note that the shielding of the wind is not expected to be sensitive to the thermal structure of the jet since it relies on electron recombination and charge exchange rate coefficients that do not strongly depend on the temperature. The local radiation field computed in this work thus remains valid for warmer or cooler winds. On the other hand, it should be kept in mind that asymptotic molecular abundances would be significantly lower at wind temperatures below 1000 K (see Fig. 6a).

In reality, shocks are likely to also contribute to the heating and compression of the gas and impact molecule synthesis. Millimeter and submillimeter interferometric observations of Class 0 molecular jets show that high velocity components of molecular lines are spatially resolved as a series of knotty structures. The detection of proper motions of these knots (Lee et al. 2015) as well as detailed kinematical properties (Santiago-García et al. 2009; Tafalla et al. 2017) suggest that they are tracing internal shocks produced by time-variability in ejection velocity (Raga et al. 1990; Stone & Norman 1993). Numerical simulations that include dust-free chemistry of H₂, but neglect the UV radiation field, show that the H₂ abundance is indeed increased due to the increase of density in the post-shock gas (Raga et al. 2005).

The effect of time-variability in the flow is not included in our simplified laminar wind model. However, our results allow us to investigate the possible impact of internal shocks on molecule formation qualitatively. In Sects. 3 and 4, we show that high molecular abundances are reached in dust-free or dust-poor winds if they are warm with temperatures between ≃ 800 and ≃ 3000 K (see Fig. 6a), dense (see Fig. 6b), and well shielded by the gas (see Fig. 9). The impact of shocks is to locally increase the temperature and the density from their initial values. Both effects drive up molecular abundances compared to a laminar wind with the same equilibrium temperature and an average mass-flux. At a larger distance, all the successive internal shocks located between the protostar and z would also increase the gas attenuation of the radiation field, increasing the molecular richness of the wind even more. However, on a large scale, the interstellar radiation field or the radiation field generated locally by strong shocks (e.g., Tappe et al. 2012) is expected to take over from the protostellar radiation field. Hence, asymptotic chemical abundances presented in Fig. 10, which neglect the effect of shocks, should be considered as lower limits if the wind equilibrium temperature is T_K ≥ 800 K, and if the protostellar radiation field dominates over other sources of the FUV field. Otherwise, the peak abundances that are reached behind shocks could be lower than in Fig. 10 if molecule formation in the shocked gas is slower than the timescale to cool below ≃800 K.

Therefore, the precise quantitative impact of shocks remains to be constrained by detailed shock modeling, including self-consistent thermal balance. However, the overall decrease of molecular abundances with the decline of the mass-loss rate and dust fraction remains a robust prediction of our model.

5.3 Observational perspective

5.3.1 Mass-loss rates

Clarifying the contribution of high velocity jets to mass and energy extraction during protostar formation remains an important observational challenge. Mass, momentum, and energy fluxes of deeply embedded Class 0 jets are usually derived from spectrally and spatially resolved CO rotational lines, assuming a standard reference CO abundance of x(CO) = x_ref ≃ 10⁻⁴ (e.g., Podio et al. 2016). Therefore, the observable quantity is not the true wind mass flux Ṁ_w, but the “observed” mass-flux Ṁ_obs derived from CO observations assuming a standard CO abundance, namely $${\dot{M}}_{\text{obs}}=\frac{x(\text{CO})}{{10}^{-4}}{\dot{M}}_{\text{w}}.$$(21)

Our models show that if the jet is dust-free or dust-poor, this assumption on the CO abundance may lead to large errors on the derived energetic properties of jets. Indeed, the mass-flux is overestimated by a factor ~ 3.6 when all of the elemental carbon is in the form of CO (i.e., for the highest Ṁ_w or Q∕Q_ref ≃ 0.1), but underestimated by orders of magnitude at lower Ṁ_w or lower Q∕Q_ref. Interestingly, we find that when the CO abundance is x(CO) ≥ 10⁻⁶, neutral carbon and CO are the main carbon carriers (see Fig. 10). It suggests that submillimeter CI lines could be used together with CO to constrain the true CO abundance versus H and measure the true mass-flux of the jet.

Alternatively, our models suggest that ratios of molecular abundances could be used in conjunction with the observed Ṁ_obs to estimate thedust content of the wind and its true mass-flux. This is shown in the next section.

5.3.2 Dust content

Our model results have shown that for a given wind temperature T_K, the molecular content is very sensitive to the precise dust fraction, suggesting that chemistry could be used as a unique tool to constrain the presence of dust in protostellar jets, and as a consequence, constrain their launching radius. Here we illustrate this point using the grid of laminar and isothermal disk wind models presented in Sect. 4, recalling that our exact results depend on specific choices of fixed parameters, which are given in the footnote of Table 3.

We first display, in Fig. 11, the correlation between abundance ratios among the brightest molecular jet tracers, namely CO/H₂, SiO/CO, and H₂O/CO, obtained for 0 ≤ Q∕Q_ref ≤ 10⁻¹ and 5 × 10⁻⁷M_⊙ yr⁻¹ ≤ Ṁ_w ≤ 5 × 10⁻⁶M_⊙ yr⁻¹, that is, for objects with high mass-loss rates. For all dust poor models, the x(CO)∕x(H₂) ratio is found to be larger than 5 × 10⁻⁴, which is in contrast with the standard ratio of ≃10⁻⁴ derived in dusty molecular environments. Such a high ratio would thus directly point toward a small dust fraction ≪ 1. However, it may also be seen that all our models follow almost the same well-defined “path” in these ratio–ratio plots; in other words, there is a strong degeneracy between Ṁ_w (in effect $n_H^0$) and Q∕Q_ref in the sense that the same ratio-ratio pair may be obtained by increasing one parameter and decreasing the other. Such plots are therefore not able to determine the dust fraction Q∕Q_ref and true mass-flux precisely by themselves. However, they do provide a good validation test of the underlying hypotheses in our simple laminar wind models, as all observations should lie close to the predicted curves.

In Fig. 12, we present a more powerful observational diagnostic of the dust fraction, where the same ratios are plotted as a function of the “observable” (two-sided) mass-flux Ṁ_obs (see Eq. (21)). The plots start at Ṁ_obs = 10⁻⁷ M_⊙ yr⁻¹, which is theminimum detectable mass-flux in CO, assuming a threshold beam-averaged column density of N_CO = 3 × 10¹⁵ cm², a beam size of 3 × 10¹⁴ cm = 20 au, and a jet speed of 100 km s⁻¹. It may be seen that the degeneracy between the density and dust fraction is lifted in different ways for the different ratios.

The curve of x(CO)∕x(H₂) vs Ṁ_obs is found to provide a good diagnostic of the dust fraction, except for small Ṁ_obs near the detection limit or for very dust-poor winds with Q∕Q_ref ≤ 10⁻³ where curves come too close to discriminate. Measuring N(CO)∕N(H₂) in prototellar jets is challenging since pure rotational and ro-vibrational lines of H₂ and submillimeter lines of CO have very different upper energy levels, and so they do not necessarily trace the same region of the jet. Simultaneously probing pure rotational lines of H₂ and high-J lines of CO would help to circumvent this problem. In that perspective, future JWST observations, in synergy with Herschel data, would help to probe the dust content of high velocity jets.

In contrast to CO/H₂, the SiO/CO ratio is very sensitive to the dust fraction even for small Ṁ_obs and small Q∕Q_ref. SiO is routinely observed through its submillimeter rotational transitions along with low-lying rotational transitions of CO (Bachiller et al. 1991; Guilloteau et al. 1992; Tychoniec et al. 2019). These lines having similar upper energy levels, they can be used to derive the abundance ratio, though opacity effects can complicate the analysis (Cabrit et al. 2012). Because of the increase of $n_{\text{H}}{}^0$ with Ṁ_w (see Eq. (19)), our 1D laminar wind model predicted that SiO would only reach its maximum abundance above a critical mass-loss rate that depends on the precise dust fraction (Fig. 10). Figure 12 shows that even though the rise of SiO/CO with Ṁ_obs is not as steep as x(SiO) versus Ṁ_w in Fig. 10, it would be possible to constrain the dust fraction from these curves if SiO/CO < 0.1. However, when the maximum possible amount of gas-phase SiO is reached for a given dust fraction, the SiO/CO abundance ratio is found to be remarkably constant with a value of about 0.1. In that case, the observed Ṁ_obs only sets a lower limit to the dust fraction. In HH212, Cabrit et al. (2012) estimated a ratio SiO/CO > 0.04, while Lee et al. (2007) estimated Ṁ_obs ≃ 2 × 10⁻⁶ M_⊙ yr⁻¹. This suggests a lower limit of Q∕Q_ref ≥ 10⁻³. Joint observations of the SiO/CO and CO/H₂ abundance ratios together with Ṁ_obs could be used to derive the combination of the mass-loss rate and dust fraction in protostellar jets.

We recall that if SiO is formed in shocks and if the wind equilibrium temperature is T_K ≥ 800 K, the above diagnostic curves for laminar wind models overestimate the true dust fraction. Since SiO is shielded by S and Si, we also expect to find a correlation between high column density of atomic sulfur and silicon as well as high abundances of SiO. In that perspective, observations of the [SI] line at 25.2 μm by JWST in synergy with ALMA will be of interest.

The ratio H₂O/CO versus Ṁ_obs appears to be another good diagnostic for the dust content of the jets. Figure 11 shows that the H₂O/CO abundance ratio depends on both the mass-loss rate and the dust fraction, and we predict rather low H₂O abundances in dust-poor jets. On-source water spectra have been acquired toward a large sample of low-mass protostars over the course of the WISH and WILL Key Programs (van Dishoeck et al. 2011; Mottram et al. 2017). HIFI observations have unveiled high-velocity bullets toward 4 out of 29 Class 0 sources of the WISH sample (Kristensen et al. 2012). The H₂O/CO abundance ratio is estimated to be ≃1−10⁻¹ for mass-loss rates of about Ṁ_obs ≃ 10⁻⁵ M_⊙ yr⁻¹ (Kristensen et al. 2011; Hirano et al. 2010). Figure 12 shows that this is in line with our laminar wind models with Q∕Q_ref ≥ 3 × 10⁻³.

Finally, since dust-poor winds are shielded by S and Si below λ = 1500Å, we predict that H₂O is photodissociated by longer-wavelength photons. Interestingly, it corresponds to the photodissociation of water through the A state that produces OH in a vibrationally hot but rotationally cold state (van Harrevelt & van Hemert 2001). Thus, in the absence on an extra source of UV photons, OH would be expected to produce a little mid-IR line from rotationally excited levels but strong ro-vibrational lines in the near-IR. This prediction can also be tested by future JWST MIRI and NIRSPEC observations. Again, shock models should be developed to properly derive the dust content in these molecular bullets. A comparison with other oxygen bearing species, such as OH or atomic oxygen, would also allow to better contain dust-poor models.

Fig. 11 Asymptotic abundance ratios for a streamline launched at R0 = 0.15 au in the disk for dust-free models (cyan) and dust-poor models with Q∕Qref = 10−3 (blue), 10−2 (green), and 0.1 (red) as well as for different values of Ṁw = 5 × 10−7 (crosses), 10−6 (triangles), 2 × 10−6 (squares), and 5 × 10−6 (circles) M⊙ yr−1. All other parameters are kept constant to the values given in the footnote of Table 3. We note that the variations of the ratios with Ṁw is due to the variations of $n_{\text{H}}{}^0$, $n_{\text{H}}{}^0/G_0^0$, and $n_{\text{H}}{}^0/W^0$ with Ṁw as indicated in Eqs. (19), (17), and (16).

Fig. 12 Asymptotic abundance ratios for a streamline launched at R0 = 0.15 au in the disk as a function of the “observed” mass-flux Ṁobs (see Eq. (21)) for dust-free models (cyan) and dust-poor models with Q∕Qref = 10−3 (blue), 10−2 (green), and 0.1 (red) as well as for different values of Ṁw = 10−7 (straight crosses), 5 × 10−7 (crosses), 10−6 (triangles), 2 × 10−6 (squares), and 5 × 10−6 (circles) M⊙ yr−1. All other parameters are kept constant to the values given in the footnote of Table 3. In particular TK = 1000 K.

5.4 Comparison with stellar wind models

Our work demonstrates that disk winds launched within the silicate dust sublimation radius should be poor in H₂, with correspondingly high ratios of CO/H₂. This unique feature of dust-poor chemistry was already identified in the context of stellar winds by previous authors (Rawlings et al. 1988; Glassgold et al. 1989, 1991; Ruden et al. 1990). With the exception of Ruden et al. (1990), who had H₂ formation dominated by H${}_2^{+}$ in their network, most of the authors also found that photodetachement of H⁻ limits the gas-phase formation rate of H₂. However, our dust-free disk wind models predict H₂ abundances that are two orders of magnitude larger than those in dust-free stellar wind models with a similar mass-loss rate. For example, our Class 0 model (a) predicts an asymptotic H₂ abundance of ≃ 5 × 10⁻³, whereas the Case 1 model of Glassgold et al. (1991) predicts only ≃2 × 10⁻⁵. This difference is ultimately due to the geometry of the wind. For an equal mass-loss rate and constant V_j ≃ 100 km s⁻¹, the initial n_H ∕W ratio in our disk wind model is similar to isotropic stellar winds (≃10¹³ cm⁻³ for Ṁ_w ≃ 2−3 × 10⁻⁶ M_⊙ yr⁻¹; see Case 1 in Glassgold et al. 1991). However, the asymptotic n_H ∕W ratio in the disk wind is increased by a large factor ${(z_0/r_0)}^2$, whereas it remains constant in the stellar winds. This allows H⁻ to survive and form H₂ in larger abundances. On the other hand, a dust-free stellar wind that accelerates from the sonic point can reach an H₂ abundance closer to our model as density at the wind base is higher for the same mass flux, leading to more efficient H₂ synthesis (see Case 2 of Glassgold et al. 1991).

Regarding the formation and survival of OH, CO, SiO, and H₂O, most of the previous studies on dust-free stellar winds have discarded the impact of a UV excess, which turns to be of key importance in low-mass protostars with an accretion shock. A notable exception is Glassgold et al. (1991), who have shown that the inclusion of a strong FUV excess leads to a dramatic decrease of the abundances of H₂O, OH, and SiO, whereas the abundance of CO is weakly affected. Our modeling explored, for the first time, a level of FUV excess adjusted self-consistently to the wind mass-flux, assuming a canonical ratio of Ṁ_w to M_acc of 0.1. It shows that dust-free disk winds remain very rich in molecules above Ṁ_w ≃ 10⁻⁶ M_⊙ yr⁻¹ (for our adopted values of T_K, R_in, R_out, V_j, R_⋆, etc). Even though an effect of the temperature is not excluded, our disk wind model also exhibits a higher n_H ∕G₀ ratio and higher total column densities (for similar Ṁ_w) than stellar wind models, causing more efficient self-shielding by S and Si and favoring the survival of molecular species against FUV photodissociation.

6 Conclusion

In this work, the chemistry and resulting molecular composition of disk winds launched within the dust sublimation radius is explored by the use of single point and stationary wind models. Our results show that the formation of H₂ through H⁻ (electron catalysis) is a viable and dominant pathway to produce H₂ in dust-free disk winds. However, this route is not efficient enough to convert all atomic hydrogen into H₂. It results in dust-free winds that are atomic, which is in line with stellar wind models (Glassgold et al. 1991). The inclusion of surviving dust has a significant chemical impact on H₂ for Q∕Q_ref ≥ 10⁻³, where Q_ref ≡ 6 × 10⁻³ is a reference dust-to-gas mass ratio. Above this value, the H₂ abundance in dust-poor winds depends on the dust fraction.

Our results also demonstrate that despite the low abundance of H₂, dust-free and dust-poor disk winds can be rich in molecules, such as CO, SiO, or H₂O. Molecular richness of dust-free winds critically depends on the properties of the wind. High temperatures (T_K > 800 K) are required to activate efficient and direct formation routes of these molecules through OH. The overall molecular content of wind models increases with the mass-loss rate of the wind and with the dust fraction. The survival of these molecules is insured by the attenuation of the radiation field by atomic species, and by dust, if any. The attenuation of the radiation field by atomic species proceeds through continuum self-shielding, a process that critically depends on the column density of the absorbing species. It results in the presence of density or, equivalently, mass-loss rate thresholds above which specific molecules are very abundant. Those thresholds depend both on the specific species, and on the dust fraction. The order in which molecules are abundant as the density and/or the dust fraction increases is CO, SiO, and H₂O. CO in abundant when continuum self-shielding of atomic carbon is triggered, whereas SiO and H₂O are abundant when self-shielding of S and Si is activated.

Our results allow us to propose observational diagnostics to probe the presence of dust in warm molecular jets (T_K ≥ 800 K). As already suggested by Glassgold et al. (1991) and in the context of stellar winds, a CO/H₂ abundance ratio that is higher than that standard ISM value of 10⁻⁴ would be the primary piece of evidence for dust-poor winds. The SiO abundance appears to be a promising diagnostic to estimate the precise dust fraction in the wind: the SiO/CO abundance ratio exhibits a step increase from a critical Ṁ_w value that depends on the dust fraction. The H₂O/CO abundance ratio is also very sensitive to the dust fraction and to the density of the wind and could be used in combination with other oxygen bearing species.

The main limitation of our predictions resides in the assumption of a laminar and isothermal flow at T_K = 1000 K. Indeed, observations at high angular resolution show that molecules are present in shocks. Our models, together with Raga et al. (2005), suggest that internal shocks can locally enhance molecular abundances by increasing the temperature and the density of the gas. In the next paper of this series, shock models including the thermal balance will be investigated. Synthetic predictions of line intensities will allow one to identify the most promising atomic and molecular transitions to constrain the presence of dust in jets. In that perspective, the present modeling constitutes a necessary step toward shock modeling by computing preshock conditions (chemical abundances and local UV field) from which shock models can be run.

Our work shows that chemistry is a powerful tool to unravel the dust content of protostellar jets. Future JWST observations, in combination with ALMA and Herschel data, will be able to provide key information on the launching region of the jets. If jets are dust-free or dust-poor, the chemical inventory of jets will give access to elemental abundances of the inner regions of protostellar disks and uncover key physical processes at stake during this early phase of star formation.

Acknowledgements

This work is part of the research programme Dutch Astrochemistry Network II with project number 614.001.751,which is (partly) financed by the Dutch Research Council (NWO). This work was also supported by the French program Physique et Chimie du Milieu Interstellaire (PCMI) funded by the Conseil National de la Recherche Scientifique (CNRS) and Centre National d’Etudes Spatiales (CNES). BG and GPdF would like to acknowledge the support from the European Research Council, under the European Community’s Seventh framework Programme, through the Advanced Grant MIST (FP7/2017-2022, No 787813).

Appendix A Radiation field

Throughout this work, we model the radiation field by means of two components.

– The FUV component is considered to be equal to the FUV part of the standard interstellar radiation field (ISRF, Mathis et al. 1983) rescaled by the parameter G₀. Following the fitting formula adopted in the Meudon PDR code⁵ (Le Petit et al. 2006), the FUV part is modeled by a mean specific intensity of $$J_{\nu }^{\text{FUV}}=G_0\times 107.2\left(\text{tanh}(4.07 {10}^{-3}\lambda -4.6)+1\right) {\lambda }^{-2.89},$$(A.1)

with the wavelength λ in Å and $J_{\nu }^{\text{ISRF},\text{FUV}}$ in ${\text{ergs}}^{-1}{\text{cm}}^{-2}{\AA }^{-1}$. This formula fits the ISRF below λ ≃ 2000Å (see Fig. A.1b).

– The visible-NIR component is modeled by a black body emission at a temperature T_vis = 4000 K, diluted by afactor W: $$J_{\nu }^{\text{visible}}=W\times B_{\nu }(T_{\text{vis}}).$$(A.2)

In Sect. 3, W is used as a parametrization of the visible-NIR intensity. In Sect. 4, this parameter is the dilution factor of the photospheric emission and can thus be written as $$W(R)=\frac{1}{2}\left(1-{(1-R_{*}^2/R^2)}^{1/2}\right),$$(A.3)

where R_* is the stellar radius and R is the distance to the central object. For R ≫ R_*, the dilution factor yields $$W(R)=\frac{1}{4}{\left(\frac{R_{*}}{R}\right)}^2.$$(A.4)

Fig. A.1 Mean intensity of the adopted radiation field for two sets of parameters {G0, W}. The FUV and the visible components are in dotted blue and dotted red lines, respectively and the total spectrum is in a black solid line. (a) Mean intensity of the radiation field adopted in Sect. 3 corresponding to {G0 = 104, W = 5 × 10−3}. (b) Mean intensity with the same FUV component but a weaker visible one ({G0 = 104, W = 5 × 10−9}) overlaid on the Mathis ISRF scaled by G0 = 104.

Figure A.1a shows the radiation fields adopted in Sect. 3. Figure A.1b shows a radiation field with the same FUV component, but a weaker visible field {G₀ = 10⁴, W = 5 × 10⁻⁹} overlaid by the Mathis et al. (1983) ISRF scaled by G₀ = 10⁴. The parametrized radiation field reproduces the ISRF within a factor 2 from the visible (λ ≃ 4000 Å) to the NIR (λ ≃ 15 000 Å) well. Thus, our parametrization also gives a good proxy for the shape of the ISRF for W = 5 × 10⁻¹³.

Appendix B Updated chemical network

Here, we review chemical reactions and their adopted rate coefficients added to the chemical network of Godard et al. (2019). Reaction rate coefficients that are outdated or inaccurate in standard astrochemical databases (KIDA, UMIST, OSU) were fit from theoretical and experimental work in the form of a modified Arrhenius law.

Gas-phase formation routes of H₂ involving two-body reactions are initiated by radiative associations. As shown in Appendix C, the efficiency of these routes is directly proportional to the radiative association rate coefficient. It is thus crucial to get accurate rate coefficients, especially for the radiative association that is found to be dominant, namely between hydrogen and electron and between hydrogen and C⁺. Table B.1 compiles the radiative association rate coefficients adopted in this work together with references used to fit the rate coefficients. Original data and adopted fits are plotted in Fig. B.1. Other reactions that are related to H₂ formation and that were added to the original network are reported in Table B.2. Lastly, chemical reactions controlling the ionization state of magnesium added to the original chemical network are reported in Table B.3.

Fig. B.1 Adopted radiative association rate coefficients. Red lines and red dots are computed or measured rate coefficients taken from references indicated in Table B.1, and blue lines are the Arrhenius fits used in this work.

Appendix C Analytical approach of H₂ formation ingas-phase

In this appendix, we propose an analytic approach for the formation of H₂ in dust-free and dust-poor environments. Our goal is to derive an analytical expression of boundaries defining the dominant formation route in the {n_H, Q∕Q_ref} plane (see Fig. 3) and to derive the analytical expression of the abundance of H₂ plotted in Fig. 4. To do so, formation rates of H₂ by H⁻ and by any XH⁺ ion are first derived by taking the photodestruction of these intermediates into account. Then, a comparison between various routes, including a three-body reaction, is given.

C.1 Formation rate by gas phase catalysis

Formation of H₂ though H⁻ and any ion XH⁺ (with X = C, S, H...) is a two stage catalytic process. In this section, we derive the resulting formation rate of H₂, assuming a steady state for the intermediate species H⁻ and XH⁺ (Bodenstein approximation). Throughout this section, the gas is assumed to be neutral and atomic (x(H) ≃ 1), which is in line with our modeling results (Sects. 3 and 4).

Electron catalysis

H₂ can be formed through the intermediate anion H⁻ via a slow radiative attachment and a fast associative detachment: $$\begin{array}{l}\text{H}+{\text{e}}^{-}\to {\text{H}}^{-}+h\nu k_1^e(T) ,\\ {\text{H}}^{-}+\text{H}\to {\text{H}}_2+{\text{e}}^{-} k_2^e(T) .\end{array}$$

The efficiency of this route depends on the survival of the intermediate H⁻. In neutral medium, the main destruction route that can compete with the associative detachment (Eq. (C.2)) is the photodetachment of the fragile H⁻ anion by visible photons, $${\text{H}}^{-}+h\nu \to \text{H}+{\text{e}}^{-}.$$(C.3)

For a diluted black body at 4000 K, H⁻ decays preferentially through the associative detachment (Eq. (C.2)) to produce H₂ if $$n_{\text{H}}/W>n_{\text{H}}{}^{\text{crit}}\equiv \frac{k_{\varphi }^{0,e}}{k_2^e}=4.6\times {10}^{14} {\text{cm}}^{-3},$$(C.4)

where $k_{\varphi }^{0,e}$ is the photodetachment rate by an undiluted black body at 4000 K. Above this critical value, the efficiency of the electron catalysis is optimal. Defining $$\eta \equiv \frac{n_{\text{H}}}{W{n}_{\text{crit}}^W},$$(C.5)

a fraction $$\frac{\eta }{1+\eta }$$(C.6)

of the H⁻ formed by radiative association is converted in H₂. The associative detachment or photodetachment is very fast, and thus H⁻ reaches steady-state abundances much faster than H₂, in a typical time scale of ⁶ $$t_{{\text{H}}^{-}}\text{ lse}{(k_1^e{n}_{\text{H}})}^{-1}\simeq 8{\left(\frac{n_{\text{H}}}{{10}^8{\text{cm}}^{-3}}\right)}^{-1} \text{s}.$$(C.7)

For $t\gg t_{{\text{H}}^{-}}$, the resulting H₂ formation rate of the full catalytic process, which is valid even when H₂ abundance is out-of-equilibrium, is: $$\begin{array}{ll}{R}_{\text{elec}}^{H_2}=\hfill & 3.4\times {10}^{-19}\left(\frac{x_e}{4.8 {10}^{-4}}\right){\left(\frac{T_{\text{K}}}{1000K}\right)}^{0.9}\hfill \\ \hfill & \times \frac{\eta }{1+\eta }{\left(\frac{n_{\text{H}}}{1{\text{cm}}^{-3}}\right)}^2{\text{cm}}^{-3}{\text{s}}^{-1}.\hfill \end{array}$$(C.8)

The efficiency of this route is directly proportional to the electron fraction. For η ≫ 1, this rate does not depend on the rate of the fast associative attachment (Eq. (C.2)), but only on the rate of the limiting slow radiative association (Eq. (C.1)), leading the a rate of $$R_{\text{elec}}^{H_2}=3.4\times {10}^{-19}\left(\frac{x_e}{4.8 {10}^{-4}}\right){\left(\frac{T_{\text{K}}}{1000\text{K}}\right)}^{0.9}{\left(\frac{n_{\text{H}}}{1{\text{cm}}^{-3}}\right)}^2 {\text{cm}}^{-3}{\text{s}}^{-1}.$$(C.9)

If H₂ is predominantly formed through H⁻ and destroyed by photodissociation at a rate k_photo, the steady-state H₂ abundance is $$x({\text{H}}_2)=R_{elec}^{{\text{H}}_2}/(n_{\text{H}}{k}_{photo}).$$(C.10)

For an unsheilded Mathis radiation field rescaled by a factor G₀, the steady-state H₂ abundance yields⁷ $$x({\text{H}}_2)=5.6\times {10}^{-9}\frac{x_e}{4.8\times {10}^{-4}}{\left(\frac{T_{\text{K}}}{1000\text{K}}\right)}^{0.9}\frac{n_{\text{H}}}{G_0}\frac{\eta }{1+\eta },$$(C.11)

with η given in Eq. (C.5).

Ionic catalysis

H₂ formation can also be catalyzed by any ion noted here as X⁺, through the intermediate ion XH⁺. The process is made of a slow radiative association and a fast ion neutral reaction: $$\begin{array}{lll}\hfill & {\text{X}}^{+}+\text{H}\to {\text{XH}}^{+}+h\nu \hfill & k_1^X,\hfill \\ \hfill & X{\text{H}}^{+}+\text{H}\to {\text{H}}_2+{\text{X}}^{+}\hfill & k_2^X.\hfill \end{array}$$

This route is analogous to the former in which ions are the catalyst. As for the latter, the efficiency of the catalytic process depends on the survival of XH⁺. In protostellar winds, the main destruction routes that can compete with a reaction (Eq. (C.13)) are a dissociative recombination that reduces the efficiency of the catalysis by a factor less than 30% and by photodissociation. Unlike the fragile H⁻, ionic intermediates XH⁺ are photodestroyed by UV photons (λ < 400 nm). For an unshielded Mathis radiation field, reaction (Eq. (C.13)) dominates over photodissociation if $$n_{\text{H}}/G_0>n_{\text{H}}{}^{\text{crit},X}\equiv \frac{k_{\varphi }^{0,X}}{k_2^X} ,$$(C.14)

where $n_{\text{H}}{}^{\text{crit},X}$ is given in Table C.1 for the main ions and $k_{\varphi }^{0,X}$ is the photodissociation rate of XH⁺ by an unshielded Mathis radiation field. Below this critical value, photodissociation reduces the efficiency of the considered route. Defining $${\eta }_X\equiv \frac{{\eta }_X}{1+{\eta }_X},$$(C.15)

a fraction $${\eta }_X\equiv \frac{n_{\text{H}}}{G_0{n}_{\text{H}}{}^{\text{crit},X}}$$(C.16)

of the XH⁺ formed by radiative association is converted in H₂. As for H⁻, XH⁺ ions reach a steady-state in very short time scales and yield a total formation rate of H₂ by X⁺ that can be written as $$R_{X^{+}}=\alpha \left(\frac{x({\text{X}}^{+})}{x_X}\right){\left(\frac{T}{1000\text{K}}\right)}^{\beta }\frac{{\eta }_X}{1+{\eta }_X}{\left(\frac{n_{\text{H}}}{1{\text{cm}}^{-3}}\right)}^2 ,$$(C.17)

where x_X is a reference value for theabundance x(X⁺) of the catalyst X⁺ and α, β are given in Table C.1.

C.2 Comparison between routes

Dust-free

Formation rates of H₂ by H⁻ and XH⁺, which were deduced from parameters and are reported in Table C.1, allow one to directly compare the efficiency of the ionic and electron catalytic routes. At 1000 K, and for the ionic and electron abundances reported in Table C.1, the order in which species dominate is e⁻, C⁺, S⁺, Si⁺, and H⁺. We note that depending on the abundance of the catalyst and on the temperature, this order can vary. For example, formation by H⁺ takes over formation by C⁺ for x(H⁺) = 10⁻⁵ and T_K ≥ 10 000 K. However, in dust-free winds, the main physical parameter that changes this order is the photodestruction of the fragile intermediates anion H⁻ by the visible field. Indeed, as shown by Eq. (C.4), the photodetachment of H⁻ reduces the efficiency of this route for $$n_{\text{H}}/W<n_{\text{H}}{}^{\text{crit}}=4.6\times {10}^{14}{\text{cm}}^{-3}.$$(C.18)

This defines the boundary ② in Fig. 3. Though, at this critical n_H ∕W, formation though H⁻ still dominates over the CH⁺ route. Assuming that CH⁺ is not photodestroyed (n_H∕G₀ > 1.7 cm⁻³), formation by C⁺ takes over formation by e⁻ if $$R_{{\text{C}}^{+}}\ge R_{\text{elec}}.$$(C.19)

In combining Eqs. (C.17) and (C.8), this condition yields $$\begin{array}{ll}{n}_{\text{H}}/W\le \hfill & 6.7\times {10}^{12}\left(\frac{x({\text{C}}^{+})}{3.6\times {10}^{-4}}\right){\left(\frac{x_e}{4.8\times {10}^{-4}}\right)}^{-1}\hfill \\ \hfill & \times {\left(\frac{T_{\text{K}}}{1000\text{K}}\right)}^{-1.32}\text{}{\text{cm}}^{-3}.\hfill \end{array}$$(C.20)

This equation gives the boundary ③ of Fig. 3.

H₂ can also be formed by three-body reactions with a rate of $$R_{3B}=k^{3B}{n}_{\text{H}}^3.$$(C.21)

Assuming that electron catalysis is optimal and by using Eqs. (C.21) and (C.8), we find that a three-body reaction takes over formation through H⁻ if $$n_{\text{H}}\ge 1.9\times {10}^{13} \left(\frac{x_e}{4.8\times {10}^{-4}}\right){\left(\frac{T_{\text{K}}}{1000 \text{K}}\right)}^{1.24} {\text{cm}}^{-3}.$$(C.22)

This defines the boundary ①.

Dusty

The inclusion of dust does not change the efficiency of gas-phase formation routes but it adds a new formation route with a rate of $$R_{\text{dust}}=3.6\times {10}^{-17}\frac{Q}{Q_{\text{ref}}}\frac{S(T_{\text{K}})}{S(1000 \text{K})}\sqrt{\frac{T_{\text{K}}}{1000 \text{K}}}{n}_{\text{H}}^2 {\text{cm}}^{-3}{s}^{-1},$$(C.23)

where the formation rate of Hollenbach & McKee (1979), assuming a single grain size distribution of radius $a_{\text{g}}=\sqrt{\langle r_c^2\rangle }=20$ nm and a grain temperature of 15 K, is adopted. We note that S(T_K) is the sticking coefficient of H on grain: $$S(T_{\text{K}})=\frac{1}{1+0.4{\left(\frac{T_{\text{K}}}{100 \text{K}}+\frac{T_{\text{gr}}}{100 \text{K}}\right)}^{0.5}+0.2\frac{T_{\text{K}}}{100 \text{K}}+0.08{\left(\frac{T_{\text{K}}}{100 \text{K}}\right)}^2},$$(C.24)

where T_gr is the grain temperature. The sticking coefficient does not depend on T_gr as long as T_gr ≲ T_K and we adopt T_gr = 15 K for simplicity purposes.

When gas-phase formation is dominated by H⁻, the critical dust fraction ratio, above which formation on dust dominate, is obtained by combining Eqs. (C.23) and (C.8): $$\begin{array}{ll}Q/Q_{\text{ref}}\ge \hfill & 9.3\times {10}^{-3}\left(\frac{x_e}{4.8\times {10}^{-4}}\right){\left(\frac{T_{\text{K}}}{1000 \text{K}}\right)}^{0.4}\hfill \\ \hfill & \times {\left(\frac{S(T_{\text{K}})}{S(1000 \text{K})}\right)}^{-1}\frac{\eta }{1+\eta }\hfill \end{array}$$(C.25)

where η is defined in Eq. (C.5). Boundary ⑥ corresponds to the case η ≫ 1, and boundary ⑤ corresponds to the case η ≪ 1. When the gas-phase formation route is dominated by CH⁺, formation on dust takes over gas-phase formation for $$\begin{array}{ll}Q/Q_{\text{ref}}\ge 1.4\times {10}^{-4}\left(\frac{x({\text{C}}^{+})}{3.6\times {10}^{-4}}\right)\hfill & {\left(\frac{T_{\text{K}}}{1000 \text{K}}\right)}^{-0.92}{\left(\frac{S(T_{\text{K}})}{S(1000 \text{K})}\right)}^{-1}.\hfill \end{array}$$(C.26)

This relation defines the boundary ④.

References

All Tables

All Figures

Fig. 1 Schematic view of the 1D geometry used in this work. The code solves the chemical evolution of a slab of gas flowing at a constantvelocity Vgas and irradiated from the left (upstream) along a direction parallel to the flow. The proton density profile nH (z) is prescribedand the temperature TK is constant. The radiation field is the sum of a FUV part modeled by an ISRF scaled by G0 and a visible part modeled by a black-body radiation field at Tvis = 4000 K, which is diluted by a factor W (see Appendix A). The attenuation of the radiation field by gas-phase photoprocesses and by dust, if any, along z was consistently computed.

In the text

	Fig. 2 Dust-free formation routes for H2 where X stands for H, C, S, and Si. Blue, green, and red, arrows represent two-body reactions, the three-body reaction, and the photodestruction of key intermediates, respectively.
In the text

Fig. 3 Schematic view of the dominant H2 formation routes depending on the density, the dust fraction, and visible radiation field W summarizing our results presented in Sect. 3.1 and in Appendix C. The location of the boundaries are given for a temperature of TK = 1000 K, xe = 4.8 × 10−4, and x(C+) = 3.6 × 10−4. The schematic view remains valid from ≃100 K up to ≃ 5000 K. Dependencies on TK, xe, and x(C+) are given inAppendix C. We note that some limits depend on the visible flux W and others do not. Depending on the visible flux, boundaries ①, ②, and ③ can merge.

In the text

Fig. 4 Steady-state abundances relative to total H nuclei for H2 and chemical species involved in its formation for dust-free single point models with G0 = 104, W = 5 × 10−7, TK = 1000 K, and nH ranging from105 to 2 × 1012 cm−3. An analytical expression of the steady state abundance of H2, assuming destruction by photodissociation and formation by H− only, is alsoplotted in black dotted line (see Appendix C). Boundaries defined in Fig. 3 are also indicated on the upper axis. It is important to note that because of the decrease of the electron fraction following recombination at high density, a three-body reaction takes over from the formation via H− (boundary ①) at a lower density than indicated in Fig. 3.

In the text

Fig. 5 Steady-state abundances of H2 normalized to its value at Q∕Qref = 1 as a functionof Q∕Qref for two densities: nH =1010 cm−3 (solid red line) and nH = 105 cm−3 (solid blue line). Other parameters are constant and equal to their fiducial values (see Table 2). Gas-phase formation of H2 is dominatedby H− for nH = 1010  cm−3 (nH∕W = 2 × 1016 cm−3) and by CH+ for nH = 105  cm−3 (nH∕W = 2 × 1011 cm−3). Dashed linesindicate H2 abundance in the absence of dust. Abundances for each set of the model are normalized to their value at Q∕Qref = 1. We note that because of the low electron fraction at nH = 1010 cm−3 (xe ≃ 10−4), the formation on grains takes over from the formation via H− (boundary ⑥) at a lower dust fraction than indicated in Fig. 3.

In the text

	Fig. 6 Steady-state abundances relative to total H nuclei for relevant molecular species from single point models in the absence of dust. Panel a: abundances as a function of temperature for G0 = 104, W = 5 × 10−7, and nH = 109 cm−3, panel b: abundances as a function of nH (lower axis) and nH∕G0 (upper axis) for TK = 1000 K, G0 = 104, and W = 5 × 10−7.
In the text

Fig. 7 Dominant reactions controlling the abundance of CO, H2O, and SiO under warm (TK ≥ 800 K) and irradiated conditions. OH appears to be a key intermediate for the three species. The ionization state of carbon controls the destruction of SiO and H2O as well as the formation of CO. The ionization state of silicon controls the formation of SiO. For an unshielded ISRF FUV radiation field, carbon and silicon are ionized for nH ∕G0 ≲ 105 cm−3 and nH ∕G0 ≲ 3 × 106 cm−3, respectively.

In the text

Fig. 8 Wind model adopted in this work. Panel a: schematic view of the geometry of the model. Streamlines are assumed to be straight lines launched from the disk (see Fig. 1a). The wind velocity Vj is constant and equal to 100 km s−1. The wind launching region extends from Rin = 0.05 au out to Rout = 0.3 au. We focus on the chemical evolution of a representative streamline launched from 0.15 au in the disk and reduce the problem to 1D (see Sect. 4.1). Panel b: prescribed density nH and unshielded FUV flux G0 profile along the representative streamline launched from 0.15 au for wind solution ⓐ (see Table 3). For other models, nH and G0 are simply rescaled according to the Eqs. (19) and (17), respectively. Panel c: prescribed nH ∕G0 and nH ∕W ratio along the same streamline. We note that due to the collimation of the flow (z0 ≫ R0), these ratios increase with distance by a factor of ${(z_0/R_0)}^2=625$. In the absence of dust, the opacity of the gas in the visible is negligible so that the nH ∕W ratio is expected to be the true local ratio between the density and visible field.

In the text

Fig. 9 Computed chemical abundances and local FUV radiation field for dust-free isothermal wind models for TK = 1000 K. Left panels: chemical abundances relative to total H nuclei. Right panels: local mean intensity of the FUV radiation field at various positions along the wind (position indicated on the curves). Photoionization thresholds of C, S, Si, Mg, and Fe are also indicated by gray, yellow, orange, black, and blue dashed lines, respectively. Row a: Class 0 model with Ṁw = 2 × 10−6 M⊙ yr−1 and M* = 0.1 M⊙. Row b: Class 0 model with a lower mass-loss rate Ṁw = 10−6 M⊙ yr−1 and same mass. Row c: Class I model, with lower accretion rate Ṁw = 10−7   M⊙ yr−1 but higher mass M* = 0.5 M⊙.

In the text

Fig. 10 Abundances at z > 1000 au for a streamline anchored at 0.15 au in the disk as a function of the mass-loss rate for various dust fractions. For Ṁw ≥ 5 × 10−7 M⊙ yr−1, the mass of the central object is 0.1 M⊙ (Class 0 model) and 0.5 M⊙ for lower mass-loss rates (Class I and II models). Dust-free models are plotted in solid lines, and dust-poor models with Q∕Qref = 10−3, 10−2, and 0.1 are plotted in dashed, dashed-dotted, and dotted lines as indicated in each panel, respectively. Horizontal black dashed lines indicate the elemental abundance of carbon (panels on CO and C), silicon (panel on SiO), and oxygen (panel on H2O).

In the text

Fig. 11 Asymptotic abundance ratios for a streamline launched at R0 = 0.15 au in the disk for dust-free models (cyan) and dust-poor models with Q∕Qref = 10−3 (blue), 10−2 (green), and 0.1 (red) as well as for different values of Ṁw = 5 × 10−7 (crosses), 10−6 (triangles), 2 × 10−6 (squares), and 5 × 10−6 (circles) M⊙ yr−1. All other parameters are kept constant to the values given in the footnote of Table 3. We note that the variations of the ratios with Ṁw is due to the variations of $n_{\text{H}}{}^0$, $n_{\text{H}}{}^0/G_0^0$, and $n_{\text{H}}{}^0/W^0$ with Ṁw as indicated in Eqs. (19), (17), and (16).

In the text

Fig. 12 Asymptotic abundance ratios for a streamline launched at R0 = 0.15 au in the disk as a function of the “observed” mass-flux Ṁobs (see Eq. (21)) for dust-free models (cyan) and dust-poor models with Q∕Qref = 10−3 (blue), 10−2 (green), and 0.1 (red) as well as for different values of Ṁw = 10−7 (straight crosses), 5 × 10−7 (crosses), 10−6 (triangles), 2 × 10−6 (squares), and 5 × 10−6 (circles) M⊙ yr−1. All other parameters are kept constant to the values given in the footnote of Table 3. In particular TK = 1000 K.

In the text

Fig. A.1 Mean intensity of the adopted radiation field for two sets of parameters {G0, W}. The FUV and the visible components are in dotted blue and dotted red lines, respectively and the total spectrum is in a black solid line. (a) Mean intensity of the radiation field adopted in Sect. 3 corresponding to {G0 = 104, W = 5 × 10−3}. (b) Mean intensity with the same FUV component but a weaker visible one ({G0 = 104, W = 5 × 10−9}) overlaid on the Mathis ISRF scaled by G0 = 104.

In the text

	Fig. B.1 Adopted radiative association rate coefficients. Red lines and red dots are computed or measured rate coefficients taken from references indicated in Table B.1, and blue lines are the Arrhenius fits used in this work.
In the text