Modeling C-shock Chemistry in Isolated Molecular Outflows

As extensively studied previously, the physical conditions within shocks and turbulent sources are found to have a significant impact on the subsequent chemistry (Gusdorf et al. 2015; Burkhardt et al. 2016; Codella et al. 2017; Lefloch et al. 2017). While several studies looked into these effects by coupling a chemical network to magnetohydrodynamics (MHD) simulations (Hincelin et al. 2013), the chemical networks used were often limited in scope in order to maintain a reasonable computation time. In order to leverage the advanced gas-grain chemical network of osu.kida within nautilus, the physical conditions are simplified to a single, 1D slab. To do this, we adopted a system of parametric equations developed by Jiménez-Serra et al. (2008) that have been fit to match the simulated physical conditions by the MHD shock structure calculated by Flower & Pineau des Forêts (2003) and Kaufman & Neufeld (1996). By using these parametric equations, we are able to accurately reproduce these physical conditions while at the same time treating the chemistry in more detail than in previous studies.

Here we consider a plane-parallel C-shock where accelerated material interacts with a cold, quiescent cloud. Because these equations were fit to data from the numerical simulations, the dependent variable is the spatial location within a snapshot of the shock, z. Since the pre-shocked material is assumed to be homogeneous, these physical locations can be translated to corresponding times it would take to propagate through the shock to achieve the spatial coordinate given by assuming a single velocity to describe the aggregate speed of the shock (i.e., the shock velocity, v_s). It should be noted that the time described here refers to the post-shock time (i.e., beginning in Figure 1, right panel). Across this shock, the magnetic fields within the source will accelerate the charged particles. Because the dust grains tend to be charged in these sources, this causes the predominantly neutral gas to temporarily decouple from the dust within the shock. These velocities of the gas and grains can be described in the frame of reference of the downstream material as

$$\begin{eqnarray}&&{v}_{\mathrm{gas}/\mathrm{grains}}=({v}_{s}-{v}_{0})\left(1-\displaystyle \frac{1}{\cosh \left[(z-{z}_{0})/{z}_{\mathrm{gas}/\mathrm{grains}}\right]}\right),\end{eqnarray} \tag{ 5 }$$

where v₀ is the post-shock velocity of both the gas and dust in the comoving shock frame of reference (i.e., the difference between the peak drift velocity and the overall shock velocity), v_s is the shock velocity, z₀ is the distance where the decoupling begins, and z_gas and z_grains are parameters that set the length scale over which neutral gas and grains are decelerated in the shock, respectively. The relative velocity between these two populations, referred to the drift velocity, is then simply

$$\begin{eqnarray}&&{v}_{d}\equiv | {v}_{\mathrm{gas}}-{v}_{\mathrm{grains}}| .\end{eqnarray} \tag{ 6 }$$

With these velocities, we can covert the spatial coordinate within the shock to the flow of time a gas particle experiences as the shock propagates through it as

$$\begin{eqnarray}&&{t}_{\mathrm{postshock}}=\int \displaystyle \frac{1}{{v}_{s}-{v}_{\mathrm{gas}}}\,{dz}.\end{eqnarray} \tag{ 7 }$$

From the subsequent time-dependent velocity profile, an initial gas density, n₀, and an initial gas temperature, T_gas,0, we can develop the temporal evolution of the density and gas temperature. By assuming conservation of mass, one can describe the density as a function of the velocity (and thus time) such that

$$\begin{eqnarray}&&{n}_{\mathrm{gas}}({t}_{\mathrm{postshock}})=\displaystyle \frac{{v}_{s}}{{v}_{s}-{v}_{\mathrm{gas}}}{n}_{0}.\end{eqnarray} \tag{ 8 }$$

Since the extinction, A_V, can be estimated to be proportional to the gas density (Bergin et al. 1998), it can be similarly scaled from an initial extinction value A_V,0.

For the purposes of matching the MHD simulations, the gas temperature is described as a "Planck-like" function such that

$$\begin{eqnarray}&&{T}_{\mathrm{gas}}={T}_{\mathrm{gas},0}+\displaystyle \frac{{\left({a}_{T}(z-{z}_{0})\right)}^{{b}_{T}}}{{e}^{(z-{z}_{0})/{z}_{T}}-1},\end{eqnarray} \tag{ 9 }$$

where a_T and z_T are fitting parameters to the location and value of the peak gas temperature and b_T is an integer fitting the temperature's sensitivity to shock location.

For the purposes of this paper, we seek to model the prototypical shocked outflow L1157. Therefore, while this model is capable of simulating a range of shock velocities and gas densities, we will focus on a detailed analysis of the physical conditions based on derived values from Burkhardt et al. (2016). Hence, we assume a v_s of 20 km s⁻¹ and an n₀ of 5 × 10⁴ cm⁻³. The effects of varying these parameters are discussed further in Section 4. Assuming standard dark cloud conditions in the pre-shock gas, our initial gas and dust temperatures were both 10 K (Hincelin et al. 2015). From Table 4 of Jiménez-Serra et al. (2008), these initial conditions correspond to fitted values of z_gas = 1.4 × 10¹⁶ cm, z_grains = 3.2 × 10¹⁵ cm, z_T = 5.0 × 10¹⁵ cm, and a_T = 2.9 × 10⁻¹⁶ K^1/6 cm⁻¹, assuming z₀ = 0 cm.

Using this formalism, the physical conditions within the shock modeled here can be seen in Figure 2. It should be noted that there is a distinct delay in the heating of the gas relative to the peak drift velocity and dust temperature (whose calculation will be described in Section 2.2.2) by roughly an order of magnitude in time. The impact of this delay will be discussed in Section 3.

Zoom In Zoom Out Reset image size

In order to accurately assess the significance of enhancement through desorption versus chemical processing in shocks, it is crucial to consider the effects of the dust temperature on the system. Therefore, we adopted the formalism from Miura et al. (2017) and Aota et al. (2015), which considers the effects from the adiabatic heating of high-velocity particle collisions and the conductive heating from the dust residing in hot gas.

To start, the prescription for gas-grain energy transfer in the free flow approximation (Probstein & Shorenstein 1968; Kitamura 1986) is used to describe the evolution of the dust temperature as a balance between the radiative cooling, collisional heating, and radiative heating.

The temperature of the dust, T_dust, can thus be described as

$$\begin{eqnarray}&&\displaystyle \frac{4}{3}\pi {a}_{\mathrm{dust}}^{2}{\rho }_{\mathrm{mat}}\displaystyle \frac{\,{{dT}}_{\mathrm{dust}}}{\,{dt}}=4\pi {a}_{\mathrm{dust}}^{2}\left({\rm{\Gamma }}-{\epsilon }_{\mathrm{em}}{\sigma }_{\mathrm{SB}}{T}_{\mathrm{dust}}^{4}\right),\end{eqnarray} \tag{ 10 }$$

where a_dust is the average dust radius, ρ_mat is the mass density within a dust grain, _em is the emission coefficient for the dust given as

$$\begin{eqnarray}&&{\epsilon }_{{em}}={{Aa}}_{\mathrm{dust}}{T}_{\mathrm{dust}}^{2}\end{eqnarray} \tag{ 11 }$$

(where A is the dust emissivity given by Draine & Lee 1984 as A = 0.122 cm⁻¹ K⁻²), σ_SB is the Stefan–Boltzmann constant, and Γ is the function describing the rate at which energy is transferred from the gas to the dust per unit area (i.e., energy flux between gas and dust). Here Γ is given by

$$\begin{eqnarray}&&{\rm{\Gamma }}={\rho }_{\mathrm{gas}}{v}_{\mathrm{drift}}\left({T}_{\mathrm{recovery}}-{T}_{\mathrm{dust}}\right){C}_{{\rm{H}}},\end{eqnarray} \tag{ 12 }$$

where T_recovery is the recovery temperature, or adiabatic equilibrium temperature, defined in Probstein & Shorenstein (1968) and Kitamura (1986) as

$$\begin{eqnarray}&&{T}_{\mathrm{recovery}}=\displaystyle \frac{2}{\gamma +1}\left(\gamma +(\gamma -1){s}_{a}^{2}-\displaystyle \frac{\gamma -1}{1+2{s}_{a}^{2}+\tfrac{{s}_{a}}{\sqrt{\pi }{e}^{{s}_{a}^{2}}\mathrm{erf}({s}_{a})}}\right),\end{eqnarray} \tag{ 13 }$$

and C_H is the heat transfer function defined by

$$\begin{eqnarray}&&{C}_{{\rm{H}}}=\displaystyle \frac{\gamma +1}{\gamma -1}\displaystyle \frac{{k}_{{\rm{B}}}}{8\,\mu {m}_{{\rm{H}}}{s}_{a}^{2}}\left[\displaystyle \frac{{s}_{a}}{\sqrt{\pi }{e}^{{s}_{a}^{2}}}+\left(\displaystyle \frac{1}{2}+{s}_{a}^{2}\right)\mathrm{erf}({s}_{a})\right].\end{eqnarray} \tag{ 14 }$$

In these equations, γ is the specific heat ratio, k_B is the Boltzmann constant, μ is the mean molecular weight of a gas particle, m_H is the mass of a hydrogen atom, and s_a is a ratio of the drift velocity and the thermal gas velocity defined by

$$\begin{eqnarray}&&{s}_{a}=\displaystyle \frac{{v}_{d}}{\sqrt{2{k}_{{\rm{B}}}{T}_{\mathrm{gas}}/\mu {m}_{H}}}.\end{eqnarray} \tag{ 15 }$$

Equation (10) can be solved for the dust temperature at which the heating of the dust by the gas is in equilibrium with the radiative cooling,

$$\begin{eqnarray}&&{\rm{\Gamma }}-4{\epsilon }_{\mathrm{em}}{\sigma }_{\mathrm{SB}}{T}_{\mathrm{dust}}^{4}=0,\end{eqnarray} \tag{ 16 }$$

which can be further rewritten as

$$\begin{eqnarray}&&{C}_{1}{v}_{d}^{3}{\alpha }_{1}({s}_{a}){\alpha }_{2}({s}_{a})-{C}_{2}{T}_{\mathrm{dust}}{v}_{d}{\alpha }_{2}({s}_{a})-{C}_{3}{T}_{\mathrm{dust}}^{6}=0\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{dust}}^{6}+\displaystyle \frac{{C}_{2}}{{C}_{3}}{v}_{d}{\alpha }_{2}{T}_{\mathrm{dust}}-\displaystyle \frac{{C}_{1}}{{C}_{3}}{v}_{d}^{3}{\alpha }_{1}{\alpha }_{2}=0,\end{eqnarray} \tag{ 18 }$$

where the C values are constants with respect to s_a defined by

$$\begin{eqnarray}&&{C}_{1}=\displaystyle \frac{1}{8}{\rho }_{\mathrm{gas},0}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{C}_{2}=\displaystyle \frac{1}{8}\displaystyle \frac{\gamma +1}{\gamma -1}\displaystyle \frac{{k}_{{\rm{B}}}}{\mu {m}_{{\rm{H}}}}{\rho }_{\mathrm{gas},0}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{C}_{3}=A{\sigma }_{\mathrm{SB}}{a}_{\mathrm{dust}}\end{eqnarray} \tag{ 21 }$$

and α values are functions of s_a defined by

$$\begin{eqnarray}&&{\alpha }_{1}({s}_{a})=1+\displaystyle \frac{\gamma }{\gamma -1}{s}_{a}^{-2}-\displaystyle \frac{1}{1+2{s}_{a}^{-2}+2\tfrac{{s}_{a}}{\sqrt{\pi }{e}^{{s}_{a}^{2}}\mathrm{erf}({s}_{a})}}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\alpha }_{2}({s}_{a})=\left(1+\displaystyle \frac{1}{2{s}_{a}^{2}}\right)\mathrm{erf}({s}_{a})+\displaystyle \frac{1}{\sqrt{\pi }{s}_{a}{e}^{{s}_{a}^{2}}}.\end{eqnarray} \tag{ 23 }$$

With this reduction of the equation, Equation (18) is solved for T_dust using the minpack's hybrd and hybrj routines with the modified Powell method (Moré et al. 1980). For the purposes of guiding the algorithm, initial guesses were computed by

$$\begin{eqnarray}&&{T}_{\mathrm{dust},\mathrm{guess}}={T}_{\mathrm{dust},0}+\displaystyle \frac{{v}_{d}}{{v}_{s}}\left(10\ {\rm{K}}\right).\end{eqnarray} \tag{ 24 }$$

Applying the temporal evolution of the gas temperature, density, and drift velocity determined in the previous section, we are able to compute the dust temperature across the evolution of the shock. In the pre-shock gas, the gas and dust are assumed to be in thermal equilibrium, or T_dust = T_gas.

In addition, we deviate from the formalism used in Jiménez-Serra et al. (2008), by setting an additional constraint of the velocity when used to calculate the rate of dust heating. Specifically, we assume that as the gas and grains begin to recouple, the drift velocity will eventually fall below the random thermal velocities of the gas. In the hot gas regime, this thermal velocity may be nontrivial. As such, the velocity of collisions may be more accurately described for these purposes as

$$\begin{eqnarray}&&v=\max \left[{v}_{d},\sqrt{\displaystyle \frac{2{k}_{{\rm{B}}}{T}_{\mathrm{gas}}}{\mu {m}_{{\rm{H}}}}}\right].\end{eqnarray} \tag{ 25 }$$

As can be seen in Figure 2, the grain temperature is found to initially track with the velocity structure within the shock, peaking at the same corresponding time as the peak velocity owing to adiabatic heating. However, as the gas continues to heat up, the grain temperature undergoes a secondary heating regime owing to the conductive heating from the hot gas. As such, while the sputtering rate (see Section 2.2.3) will fall off rather rapidly, the ability for species to redeposit onto the grains is delayed until the gas can cool sufficiently.

One major physical process that is common in astrophysical shocks is the nonthermal desorption of molecules from the ice into the gas phase through the collision of moderately fast gas-phase particles with the dust grains (e.g., sputtering). This process is a key mechanism for lifting complex molecules into the gas phase nondestructively. For typical C-shock velocities (e.g., 10–40 km s⁻¹), the sputtering is significant enough to desorb most species from the ice mantle without destroying them or having significant grain core erosion, both of which are common in high-velocity J-shocks (Requena-Torres et al. 2006).

In order to include these processes, the sputtering yields must be computed for each species on the ice, which is dependent on the conditions of the shock, the nature of the incident particle, and the binding properties of the target molecule. Utilizing the procedure discussed in Caselli et al. (1997) and Jiménez-Serra et al. (2008), with the initial formalism introduced in Draine & Salpeter (1979), we can compute the rate at which some ice species, J, is sputtered off the surface by some incident particle, P, such that

$$\begin{eqnarray}&&\left[\displaystyle \frac{\,{{dn}}_{{\rm{J}}}}{\,{dt}}\right]=\displaystyle \frac{\pi {a}_{\mathrm{dust}}^{2}{n}_{P}}{2{s}_{P}}\sqrt{\displaystyle \frac{8{k}_{{\rm{B}}}{T}_{\mathrm{gas}}}{\pi {m}_{P}}}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\int }_{{x}_{\mathrm{th}}}^{\infty }{x}^{2}\left({e}^{-{\left(x-{s}_{P}\right)}^{2}}-{e}^{-{\left(x+{s}_{P}\right)}^{2}}\right)\left\langle Y({E}_{P})\right\rangle \,{dx},\end{eqnarray} \tag{ 27 }$$

where a_dust is the average dust grain radius, n_P is the gas-phase density of the projectile, k_B is the Boltzmann constant, T_gas is the gas temperature, m_P is the mass of the projectile. Further, s_P is the ratio between the kinetic energy of a particle traveling at the current drift velocity and the thermal energy of the gas, or

$$\begin{eqnarray}&&{s}_{P}^{2}\equiv \displaystyle \frac{{m}_{P}{v}_{d}^{2}}{2{k}_{{\rm{B}}}{T}_{\mathrm{gas}}}.\end{eqnarray} \tag{ 28 }$$

Similarly, x is the ratio between the kinetic energy of a projectile at some specific velocity, v_x, and the thermal energy of the gas. This is related to the impact energy of the projectile, E_P, as

$$\begin{eqnarray}&&{E}_{P}={x}^{2}{k}_{{\rm{B}}}{T}_{\mathrm{gas}}.\end{eqnarray} \tag{ 29 }$$

The integral, therefore, is equivalent to the total contribution from projectiles at all velocities corresponding to an x value higher than x_th, which is defined as

$$\begin{eqnarray}&&{x}_{\mathrm{th}}=\sqrt{\displaystyle \frac{{\epsilon }_{0}{U}_{0}}{\eta {k}_{{\rm{B}}}{T}_{\mathrm{gas}}}},\end{eqnarray} \tag{ 30 }$$

where U₀ is the binding energy of the species to be sputtered off and ₀ is defined as

$$\begin{eqnarray}&&{\epsilon }_{0}=\max [1,4\eta ].\end{eqnarray} \tag{ 31 }$$

Here η is a ratio describing the relative masses of the projectile (m_P) and the target species (m_T)

$$\begin{eqnarray}&&\eta =4\xi \displaystyle \frac{{m}_{P}{m}_{T}}{{\left({m}_{P}+{m}_{T}\right)}^{-2}},\end{eqnarray} \tag{ 32 }$$

where ξ is a material-specific efficiency factor assumed to be 0.8 for ices (Draine & Salpeter 1979). The angle-averaged value for the sputtering yield, Y, which is approximately equal to double the perpendicular yield value, is described in Draine & Salpeter (1979) such that

$$\begin{eqnarray}&&{\left\langle Y\right\rangle }_{\theta }\approx 2Y(\theta =0)=K\displaystyle \frac{{\left(\epsilon -{\epsilon }_{0}\right)}^{2}}{1+{\left(\tfrac{\epsilon }{30}\right)}^{4/3}}\ \mathrm{for}\ \epsilon \gt {\epsilon }_{0},\end{eqnarray} \tag{ 33 }$$

where K ≈ 8.3 × 10⁻⁴ and is defined as

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{E}_{P}}{{U}_{0}}\eta .\end{eqnarray} \tag{ 34 }$$

For the purposes of calculating these sputtering yields, the binding energies, U₀, need to be known for each species; these are already given within nautilus as part of the KIDA-2014 network. It should be noted that values for many of these species are estimated and have not been explicitly experimentally measured.

For each grain species in the network, a sputtering reaction was added for the following gas-phase projectiles: H₂, He, C, O, Si, Fe, and CO. To compute n_P, the gas-phase abundance within the model at each time step for each projectile was multiplied by the overall gas density (i.e., ${n}_{{{\rm{H}}}_{2}}$). In the nautilus three-phase model, sputtering can only occur for species on the surface of the ice (i.e., the top few monolayers), similar to the vast majority of other desorption mechanisms (Ruaud et al. 2016). The integral over x was computed with the dqags routine in the quadpack package.⁹ Computing this indefinite integral numerically was found to be highly sensitive to the sampling rate of the integrating routine, as much of the range in which the integral is evaluated over is trivially small. As such, to increase the efficiency of the computation of the integral, we converted the integral to be definite by approximating the upper and lower bounds such that the following applies:

• 1.  
  If x_th > s_P+20, then the vast majority of the projectiles at this drift velocity are traveling at too low of a velocity to effectively sputter. In this case, we assume that the integral will evaluate to zero.
• 2.  
  If x_th < s_P − 10, then the vast majority of the projectiles at this drift velocity are potentially good projectiles for efficient sputtering. In this case, we reset the lower bound of the integral to s_P − 10 and the upper bound to s_P + 20.
• 3.  
  Otherwise, x_th is reasonably close to s_P and will not cause any significant issues integrating as a lower bound. In this case, we only reset the upper bound to s_P + 20 for computational efficiency instead of evaluating an indefinite integral.

For the first ∼100 yr following the shock, the collisions between the gas and grains begin to increase in frequency and impact energy, but the drift velocity has yet to reach its peak value. While the velocities are too low to induce sputtering, it is sufficient to begin heating the dust grain (ΔT_dust ≲ 15 K). This can have a number of effects, the foremost being the increased mobility of species to diffuse within the ice. Similar to what is seen in chemical models of hot cores, this warm-up phase can, in turn, rapidly accelerate the formation of complex molecules. This phase is, however, much shorter than the normal timescales for thermal desorption (Garrod 2013), which is found to be mostly inefficient in this regime. In this regime, we classify species into one of three types:

• 1.  
  Those for which gas and ice abundance profiles are constant over time, implying that the dominant ice and gas-phase reactions were already efficient prior to this increase in dust temperature. These species tend to be relatively simple (e.g., CO, H₂O, NH₃, CH₃OH) and thus would have relatively low binding energies and could already easily diffuse.
• 2.  
  Those for which the ice abundance rapidly increases in this short time period and the gas-phase abundance becomes enhanced owing to nonthermal desorption from exothermic reactions (Garrod et al. 2007). These species (e.g., CH₃CHO, OCN, CN, SO₂) tend to have fairly low ice abundances compared to what might be expected through observational constraints, which is likely because these species were unable to efficiently diffuse in the dark cloud conditions. The primary formation routes for the molecules were most efficient in the gas phase and thus were far too abundant for the ice enhancements to alter it significantly. For some species, like OH, the molecule is enhanced on the ice through the dissociation of species, such as HOCO, which also forms CO that reacts with N to form OCN. In the event that the ice chemistry for these species is found to be more efficient (e.g., comparable to or more efficient than the gas-phase chemistry), then it possible that these reactions may induce additional gas-phase enhancements in the peak sputtering regime.
• 3.  
  Finally, some species also appear to have an enhancement in the ice due to the heating of the dust. However, it only starts occurring once the sputtering is also efficient. This behavior manifests itself as a small increase in the ice surface abundance prior to the rapid depletion in the next regime. These species tend to be the more complex species in the network, such as HCOOCH₃, NH₂CHO, and HCOOH, since these require additionally high dust temperatures to effectively diffuse, or were enhanced by the initial production of the ice species in the previous group.

While a large number of species display some sort of ice enhancement in this short time period, the subsequent chemistry in the more abundant gas phase dominates any additional production during this regime. This implies that in the astrophysical shocks the secondary chemical effects prior to the peak sputtering are relatively minimal. And to this extent, low-velocity shocks can act as a useful tool for constraining the underlying ice chemistry without having to consider additional physical processes. It should be noted that an expansion of the limited radical–radical reaction network used could impact the results described here.

Here the drift velocity, and thus the sputtering rate, reaches its peak value. At this point, the species in the ice surface and mantle are rapidly lifted into the gas phase, leaving only negligible abundances on the grain (i.e., less than one complete monolayer). Indeed, this process is seen to be much more efficient than both thermal and reactive desorption within shocks. Again, here, we describe the major groups of species in this regime:

• 1.  
  For highly abundant gas-phase species, such as OH, CS, and CN, the ice abundance was only a small fraction of the total reservoir of those molecules owing to gas-phase formation routes being much more efficient. As such, the liberation of the ice population did not appreciably enhance the gas-phase abundance of these molecules.
• 2.  
  For other species, the new gas-phase abundances are nearly equivalent to the pre-shock ice abundance with a near-instantaneous transition. While some of these species may have initially formed in the gas phase and froze out, the primary reservoir of these species, such as H₂O, NH₃, and HCOOH, prior to the shock was in these ices. Hence, the enhancements in the gas-phase abundances are exclusively due to sputtering and not post-shock chemistry. Thus, these species are excellent molecules for constraining their initial ice abundances.
• 3.  
  For many species in the model, such as CH₃OH, H₂CO, and HNCO, while the post-sputtering gas-phase abundances are equivalent to pre-sputtering ice abundances, the change is not as rapid as for the previous group and the gas-phase enhancement of the molecules begins prior to the majority of the nonthermal desorption brought on by the sputtering. This can be due to either the species efficiently sputtering off before the peak sputtering velocities or additional induced reactions by the increased mobility on the grain and additional enhancements of related species. Observationally, this would mean that the shock-enhanced abundances would correspond to times earlier than when it is seen for the previously mentioned class of molecules.
• 4.  
  Finally, some molecules have gas-phase significant enhancements on top of the sputtering that produce abundances far larger than the pre-sputtering peak solid-phase abundances. Some of these species tend to be the same species that, during the pre-sputtering peak dust heating, would have ice enhancements without the corresponding gas-phase bump (e.g., OCN and CH₃CHO) or were secondary products from the enhancements of other species (e.g., H₂CCO). Many of these species are enhanced by the initial liberations into the gas-phase of species typically locked in the ice or the slightly heated gas temperatures.

The major relic of this phase is clearly the desorption of essentially the entire ice surface and mantle into the gas phase. While some species lift into the gas phase without additional simultaneous chemical processing, some species are significantly enhanced with the introduction of newly, abundant gas-phase species. From these findings, it is clear that this nonthermal desorption both is highly efficient and produces a unique chemical makeup of the gas in the shock relative to other forms of desorption. Furthermore, for a shock velocity of ∼20 km s⁻¹, this phase would hypothetically occur over a roughly 500 au physical scale. And for nearby molecular outflows such as L1157, ∼250 pc away (Looney et al. 2007), this interior structure within the bow shock would be potentially resolvable with an observation with better than 2'' resolution, which is easily achievable with modern interferometers.

Here, with the ice surface and mantle now residing in the gas phase, the gas temperature reaches its peak temperature. With the chemically enriched gas in a heated environment, new chemistry can proceed that would have been highly inefficient. Due to the extremely low abundances of solid-phase species, the overall profiles of the non-gas-phase abundances will not be discussed here in detail. In general, for a given species the gas-phase abundances will display one of the following trends:

• 1.  
  For some, no appreciable change in the gas-phase abundance occurs in the high-temperature regime. Many of the species that were significantly enhanced by the sputtering in the previous regime, such as CO, H₂O, and NH₃, maintain a stable gas-phase abundance such that there is no significant additional production. This is also true for species that experienced additional enhancements during the peak sputtering, such as SO and SO₂. In particular, HNCO does not see any significant post-shock gas-phase enhancements in this regime, which may be due to inefficiencies in the ice chemistry of the model and will be discussed in a later section.
• 2.  
  Other gas-phase abundances are enhanced during the high-temperature regime for species. For some, like NH₂CHO, CH₃CHO, and HCOOCH₃, this is due to an increase of its precursor species due to the sputtering in the previous phase and certain reactions with large barriers or strong temperature dependencies now being accessible. For others, this enhancement is due to the dissociation of other, more complex species, such as OH from the processing of recently enhanced H₂CO.
• 3.  
  The gas-phase abundances were also diminished for some species, including OH, CN, and OCN. For these, the post-shock chemistry that is enhancing other species is likely depleting these species that were predominantly in the ice prior to this.

In the shocked outflow L1157, we are able to sample different age shocks by comparing B1 and B2, which differ by about 2 × 10³ yr (Gueth et al. 1996). Thus, the comparison between these two shocks should probe the hot gas phase of our shock chemistry model. This will be discussed in detail in Section 4.1.

Here the drift velocity and temperatures have fallen back down to the initial pre-shock conditions. Also at this point, the density has finally reached its final value, which is roughly a factor of 5 larger than the pre-shock conditions. Also, the reservoir of molecules on the ice can return to roughly what it was in the initial cold core phase. While the abundances on the grain did increase in the hot gas regime, it is likely due to the reduction of drift velocity and thus the repressing of sputtering rather than adsorption. Here, we see the following:

• 1.  
  Many species have the ice abundances surpass the gas-phase abundances. In addition, over the timescale in which the ice abundance became larger, the gas-phase abundance varies. In many cases, the majority of the gas-phase population is redeposited onto the ice with little additional processing, such as with CH₃OH, NH₃, and H₂O, and occurs very close to 10⁴ yr after the shock. Other species, such as HNCO and NH₂CHO, undergo additional gas-phase enhancements, as their precursors have not yet undergone dissociation or adsorption. As such, the ice abundances of these species tend to not surpass the gas phase until closer to 5 × 10⁴ yr. For many of these species, the ice abundance is either comparable to (e.g., H₂O and NH₃) or greater than the values prior to the shock (e.g., NH₂CHO, CH₃CHO, and HCOOCH₃).
• 2.  
  Other species, such as OH, OCN, and SO₂, also have their solid-phase abundances significantly enhanced, but not enough to become close to the gas-phase abundances. Many of these species were predominantly gas-phase species in the pre-shock conditions. And, as such, it is to be expected that the gas-phase reactions would come to dominate again. Meanwhile, some molecules such as SO were found to have comparable ice and gas-phase abundances, again due to significant gas-phase production routes in addition to the redeposition.

Overall, essentially every species shows significant rebuilding of the ice mantle. However, the timescale in which species fall back onto the grain was found to vary from species to species, which could be possible to observationally probe for additional post-shock enhancement. While the species with earlier redeposition times and no additional gas-phase production in this regime were found to have comparable ice abundances prior to the peak sputtering, many of the other species had the resulting ice abundances being strongly enhanced (e.g., HCOOCH₃ and NH₂CHO). As will be discussed in Section 4.2, this shock enhancement of the ice may prove to be a useful probe for studying the physical and chemical history of these sources.

Beyond 10⁵ yr, the physical conditions do not change over time and represent the effects on the chemistry long after the shock has passed. Excluding the change in density and extinction, the physical conditions are similar to the pre-shock conditions. As such, we should be able to determine whether any changes from pre-shock molecular abundances are due to this change in physical conditions or the aftereffects of the shock.

• 1.  
  Many species, at this point, have reached a stable abundance for both the gas phase and ice mantle. This includes the abundant gas-phase species (e.g., CO and HCN) and the main constituents of the ice (e.g., H₂O and NH₃). As would be expected, the abundance trends discussed previously for these species remain comparable to the pre-shock values.
• 2.  
  For many of the most complex species, such as HCOOCH₃ and CH₃CHO, the ice abundances decrease owing to the lack of efficient formation mechanisms on cold grains, as was seen in the pre-shock conditions. This can also been seen by the increase in ice abundance of CH₃OH and HCOOH. Meanwhile, the gas-phase species also decrease owing to processing and additional redeposition, but they eventually stabilize at an abundance that can be maintained by the gas-phase production routes the species may have.
• 3.  
  Finally, for species such as OCN and many of the sulfur-bearing species, both the gas-phase and ice abundances drop like the previous group of molecules. But because these lighter species are less reliant on the accelerated diffusion on warm grains to form, they eventually stabilize to a final abundance like the first group.

Here we can see that for many of the species that had significant ice-phase enhancements in the first regime of the model there exists a temporary surge in ice chemistry that lasts for a few times 10⁴ yr, which is longer than many of the gas-phase enhancements commonly seen in these shocks. For a select number of species, their ice abundances are permanently enhanced, which appear to generally be the species that underwent significant post-shock chemistry, like CH₃CHO and NH₂CHO or species that were already highly abundant in the ice but increased in abundance at the end of the model, such as CH₃OH, CH₃OCH₃, and HCOOH. More generally, the entire post-shock regime may be characterized by having ices that have a higher ratio of complex organics relative to H₂O.

In addition, we found that the evolution of certain clusters of species tended to proceed similarly within the shock, even if they were not directly chemically related, as seen for abundance ratios of a single species at different times within the model. First, any overall production or destruction of a molecule during the early onset of the shock can be described through the ratio of the pre-shock ice abundance, X_ice,preshock, and the gas-phase abundance when the gas temperature is at its peak value, ${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$. In the case where sputtering occurs, but no significant post-shock chemistry, the ratio of these abundances, ${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$/X_ice,preshock, will be of order unity. When post-shock chemistry creates further enhancements, such as with NH₂CHO, this ratio should be large.

Second, the timescale in which any gas-phase enhancements or depletion occur can be studied by the ratio between ${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$ and the gas-phase abundance at a time after the initial shock has passed. In this case, we adopt a comparison time ${t}_{\mathrm{postshock}}={10}^{5}$ yr, as it is a key inflection point for many of the species studied here. This ratio, ${X}_{\mathrm{gas},{10}^{5}\mathrm{yrs}}$/${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$, should be small (e.g., <1) for species where little additional chemistry occurred after they sputtered from the grain. For species that still undergo gas-phase chemistry or appear to continue to redeposit well into 10⁵ yr after the shock, this ratio can range from order unity to much greater than 1 (e.g., CN). These two ratios stratify the species studied here into general groups defined by their various means of enhancements, as seen in Figure 11.

Zoom In Zoom Out Reset image size

Many species follow the general bulk evolution of the ice and do not display any significant post-shock gas-phase chemistry, including CO, H₂O, NH₃, H₂CO, CH₃OH, HCN, HCOOH, and CH₃OCH₃, as seen in Figure 11 (blue) to have small values for both ratios. Other species, however, display additional enhancements in the gas-phase during the shock events, and many have similar profiles. HNCO, CH₃COCH₃, H₂S, H₂CCO, and SO are abundant in the ice prior to the shock but are strongly enhanced during both the sputtering and redeposition regimes. These species can be seen as a cluster of magenta points in Figure 11 to have similarly small values for ${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$/X_ice,preshock, but ${X}_{\mathrm{gas},{10}^{5}\mathrm{yrs}}$/${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$ are of order unity. Meanwhile, CS, NO, and SO₂, shown in green, are primarily gas-phase species that display reasonable enhancements during or just after the redeposition regime. The red-colored species are characterized by strong, consistent post-shock chemistry. CH₃CHO shows two discrete phases of enhancements, while HCOOCH₃ and NH₂CHO show strong, continual enhancements from the start of sputtering to redeposition. Both of these groups of species have large values for ${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$/X_ice,preshock but a wider range of ${X}_{\mathrm{gas},{10}^{5}\mathrm{yrs}}$/${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$ owing to the different timescales of depletion and enhancement. Gas-phase species that tend to be highly reactive, such as OH, CN, and OCN, undergo significant changes over time, especially characterized by a drop in abundance during the hot gas regime, resulting in large ${X}_{\mathrm{gas},{10}^{5}\mathrm{yrs}}$/${X}_{\mathrm{gas},{T}_{\mathrm{gas}\mathrm{peak}}}$ values, and are displayed as orange in Figure 11. As all these species will be enhanced at different discrete stages in the model, these enhancements should be theoretically reproducible by comparing similar-velocity shocks of different ages or the internal structure of single bow shocks.

Methanol has long been used as a classical indicator of "complex" chemistry, as the six-atom species is thought to be efficiently formed on the ice through successive hydrogenation of the highly abundant CO (Herbst & Millar 2008). As such, for the purposes of studying the chemistry in shocks, CH₃OH can be a useful species to study the lifting of the ice into the gas phase. Through deep CARMA observations of L1157, Burkhardt et al. (2016) found that the consistent abundance enhancement throughout the shock event of CH₃OH indicates that its primary source is from the liberation from grains during the shock.

In this study, we see that the strongest source of enhancement of gas-phase CH₃OH is, in fact, the sputtering event with no other significant gas-phase enhancements. On the order of 10³ yr after the shock, as probed by the two shock fronts in L1157, the abundance falls off rather smoothly in a relatively log-linear nature between the peak sputtering and the hot gas phase. Thus, on the timescales discussed, it would appear to be a proper tracer to compare the relative shock enhancements. From the analysis of the full network, other species that may also serve as good tracers of overall shock enhancements include both the major constituent species, such as CO, H₂O, and NH₃, and additional complex species like H₂CO, HCOOH, and CH₃OCH₃. These species display similar falloff profiles and are all species that should fairly exclusively be produced on the grain in terms of post-shock enhancements.

While we refrain from discussing species abundances in this model owing to the lack of strong constraints on the ice chemistry, it should be noted that the peak phase abundance of CH₃OH is in agreement with what is typically seen in astrophysical shocks (10⁻⁶; McGuire et al. 2015; Burkhardt et al. 2016). However, it has been proposed that dissociation during sputtering may be important and is not considered in this model at this time (Suutarinen et al. 2014).

In order to quantify this enhancement observationally, the measured abundances will need to be compared to another value. While CO is a common species to compare to historically, the majority of the pre-shock abundance of CO is also locked up in ice. As such, a comparison of CO to other species described here would probe the rate of desorption and redeposition more than any specific enhancements. As such, a more suitable species to test this would be one (1) that is primarily in the gas phase prior to the shock and (2) whose abundance is either fairly stable or could be described well by a physical model within the source. It may also be possible to compare the CH₃OH abundance to other positions within the source to determine the degree of enhancement, taking into account the change in physical conditions like density and temperature. For example, if a source contains multiple shock events with similar velocities or a single shock event where CH₃OH is sufficiently bright where either the gas is not shocked or the shock has long passed, the overall abundance profile of the molecule can be compared to modeled profiles.

Due to the abundance enhancements of HNCO in the L1157 bow shocks, especially comparing the B1 shock front to the older B2 shock front (Δt ∼ 2000 yr), it has been proposed that HNCO may be enhanced by post-shock gas-phase chemistry in addition to the initial release from sputtering (Rodríguez-Fernández et al. 2010; Burkhardt et al. 2016). Specifically, HNCO was seen to have a factor of two enhancement in the abundance over this timescale. In the model, however, at 2 × 10³ yr after the shock (i.e., during the hot gas regime), by comparing the abundances between the beginning of this phase (∼200 yr) and at the peak (∼2.1×10³ yr), one can see that this enhancement is only ∼20%. While we do not reproduce this peak, it is likely because the dominant reactions to produce HNCO in this phase require CH₂ or OCN, both of which were depleted prior to the peak of HNCO. In addition, OCN seems to be underproduced on the ice relative to estimated abundances from infrared absorption ice studies (Boogert et al. 2015). It should be noted that the calculated ice abundance of OCN is comparable to observed gas-phase abundances in dense cores (Marcelino et al. 2018). As such, while we do not reproduce this post-shock enhancement, it is possible that it may still be feasible if the pre-shock ice abundance of OCN was closer to observed values.

Despite this lack of overall abundance enhancement in the hot gas chemistry regime, the chemistry of HNCO is unique compared to the sputtering-dominated species, as discussed in Section 3.3.6 and seen in Figure 11. Instead of a strong enhancement, as seen by NH₂CHO, the gas-phase abundance of HNCO merely remains elevated relative to the pre-shock value for much longer than the sputtering-dominated species, such as CH₃OH. As a result, the HNCO/CH₃OH abundance ratio increases as the shock gets older until HNCO fully redeposits onto the ice, as was described by Burkhardt et al. (2016).

It has been proposed that the chemistry of sulfur-bearing molecules may be a useful indicator of shock evolution (Hartquist et al. 1980; Neufeld & Dalgarno 1989; Codella et al. 2005; Wakelam et al. 2005), including an increase of the SO₂/SO ratio over time and the general conversion of H₂S into SO and SO₂. Here, while both SO and SO₂ are enhanced at timescales >10³ yr, we find that SO₂ increases significantly more rapidly over this time, providing further evidence to support this earlier hypothesis. However, H₂S also appears to be enhanced in the later stages of the shock, which suggests that it is not a primary source of sulfur for SO and SO₂. As it is apparent that sulfur chemistry can be strongly tied to the evolution of a shock, a more in-depth study of high-temperature sulfur chemistry may be warranted.