© A. Lehmann et al. 2022

1 Introduction

At first sight, the problem of galaxy formation might appear simple: a reservoir of gas in a potential well of dark matter must collapse under self-gravity to form the stars that make up the galaxy. This process is complicated by stellar and black-hole feedback mechanisms that individually impinge on the small-scale local environments, the dense molecular clouds in which stars ultimately form, and collectively produce a galaxy-scale feedback mechanism that affects the total reservoir of infalling material. Observational probes of these processes promise to uncover the energetics of the systems.

The complex galactic life cycle of baryonic matter always involves vast amounts of mechanical energy, often observed in the form of galaxy-scale high-velocity flows (V >500 km s⁻¹) that must result in shock waves. Such outflows have long been observed around starbursts and active galactic nuclei with atomic tracers (see the review by Veilleux et al. 2005). More recently, these kinematic signatures have been increasingly observed in molecular emission, including H₂ (Rupke & Veilleux 2013; Petric et al. 2018), CO (Gowardhan et al. 2018; Fluetsch et al. 2019), OH (González-Alfonso et al. 2017), HCN, HNC, and HCO⁺ (Aalto et al. 2012; Lindberg et al. 2016), and CH⁺ (Falgarone et al. 2017; Vidal-García et al. 2021), revealing the multi-phase nature of these flows. These molecular observations are challenging to explain given that molecules cannot survive in the hot gas (T >10⁶ K) generated by shock waves at these velocities. An accurate modelling of the molecular component is important because it carries a significant fraction of the outflow mass and momentum, directly concerns the star forming activity of galaxies, and can dominate the cooling budget over the X-ray cooling of the hot gas.

A recent paradigm for interpreting such puzzling molecular emission has emerged, wherein the mechanical energy of a large-scale high-velocity shock dissipates in myriad low- and intermediate-velocity shocks (V <100 km s⁻¹) (Guillard et al. 2009; Lesaffre et al. 2013; Appleton et al. 2017). These lower-velocity shocks are generated from the cascade of mechanical energy in the supersonic turbulence driven in the multi-phase plasma overrun by the large-scale shock. Little work has been done to use this framework to quantitatively connect observations to the energetics and physical conditions of the observed system. Accurate models of molecular shocks over a broad range of velocities and gas conditions are required to capture the diversity of dissipative structures within this framework. We present in this paper new shock models and make a first step in using grids of shock models to connect observed emission lines to the physical information of the system.

The first of the two main goals of this work is to model molecular shocks in a broad range of environments. The microphysics and computational methods required to accurately model atomic and molecular shocks have been rigorously studied for decades (see the introduction of Lehmann et al. 2020, hereafter Paper I). In Paper I we developed an iterative, post-processed radiative transfer algorithm to self-consistently compute the UV radiation from collisionally excited atomic H in weakly magnetised (pre-shock transverse magnetic field $B_0=0.1\text{ upmuG}{(n_{\text{H}}/{\text{cm}}^{-3})}^{1/2}$ for pre-shock density n_H), J-type molecular shocks at intermediate velocities (30 km s⁻¹  ≤ V_s ≤ 60 km s⁻¹) and a single pre-shock density, n_H = 10⁴ cm⁻³. We extend that work by computing a grid of shock models over a range of densities (10² cm⁻³  ≤ n_H ≤ 10⁶ cm⁻³) and broaden the range of velocities (5 km s⁻¹  ≤ V_s ≤ 80 km s⁻¹). We also compute two more grids at moderate magnetic field strengths ($B_0=1\text{ upmuG}{(n_{\text{H}}/{\text{cm}}^{-3})}^{1/2}$) and with or without an external UV radiation field. This work therefore presents molecular shock models with the broadest published range of physical conditions. The shock code we use and the grids of models are presented in Sect. 2. While we focus on an extragalactic application, shocks in this parameter space have applications to numerous local interstellar molecular environments that are the subjects of ongoing research, for example protostellar outflows (Dopita & Sutherland 2017; Gusdorf et al. 2017; Fang et al. 2018; van Dishoeck et al. 2021), infrared dark clouds (Pon et al. 2016), stellar winds (Tram et al. 2018), supernova remnants (Reach et al. 2019), or cloud-cloud collisions (Armijos-Abendaño et al. 2020; Enokiya et al. 2021).

Amongst the molecules observed in extragalactic environments, CH⁺ is particularly interesting. Its dominant formation pathway is highly endothermic, and yet it seems to take place in cold molecular environments. This fact, in combination with a destruction timescale too short to allow for transport far from the site of formation, makes the presence of CH⁺ in dense environments an indicator of bursts of turbulent dissipation. Recent observations of its rovibrational lines around high-redshift starburst galaxies (Falgarone et al. 2017; Vidal-García et al. 2021), in photodissociation regions in the Orion molecular cloud (Goicoechea et al. 2019), in the inner regions of NGC 1365 (Sandqvist et al. 2021), or in planetary nebulae (Neufeld et al. 2021) motivates a focused study on this molecule. In Sect. 3, we consider in detail the radiative transfer of CH⁺ rovibrational lines in different kinds of shocks. In Sect. 5, we also provide line intensities and column densities for species other than CH⁺. Many of these lines are observable with state-of-the-art observational facilities such as the James Webb Space Telescope (JWST) or the Atacama Large Millimeter/submillimeter Array (ALMA).

The second main goal is to develop a simple framework to use shock model results to interpret observations of environments containing an ensemble of shocks, as in the cascade of supersonic turbulence. This toy model is presented in Sect. 4, where we demonstrate how to use the results of the shock models to deduce the mechanical energy and physical conditions from an observed atomic or molecular line. As a template observation, we apply the ensemble model to the CH⁺ (1−0) broad (~ 1000 km s⁻¹) emission line from the high-z starburst galaxy the Cosmic Eyelash (Falgarone et al. 2017). Finally, in Sect. 5 we show how lines other than CH⁺, such as rovibrational lines of H₂ and CO andatomic fine-structure and metastable lines, can be used to further constrain our knowledge of the Eyelash galaxy and similar systems.

2 Shock model and grid

For this work, we use the Paris–Durham shock code¹, a publicly available state-of-the-art tool designed to study molecular chemistry in magnetised shocks. In this section we summarise some relevant theoretical aspects of magnetised shocks and introduce the grids of shocks that we analyse in this work.

2.1 Paris–Durham shock code

The Paris–Durham shock code solves the equations of magnetohydrodynamics (MHD) coupled with cooling functions and an extensive chemical reaction network appropriate to the molecular phases of the interstellar medium. The code simultaneously solves ordinary differential equations for the conservation of number density, mass density, momentum and energy for both the ion and neutral fluid, for each level of H₂, and each of the 141 species in the chemical network. To do this, it uses the thoroughly tested DVODE forward integration algorithm (Hindmarsh 1983) to give steady-state, plane-parallel shock models. The version of the shock code used in this work is that of Flower & Pineau des Forêts (2003), with updates described in Lesaffre et al. (2013), Godard et al. (2019), and Paper I.

This shock code allows the computation of stationary MHD shocks that, depending on the physical conditions, may be continuous C- and C^*-type or discontinuous J- and CJ-type shocks. It is important to capture the full diversity of shock types because differences in their thermal structure lead to wildly different observational predictions (see the review by Draine & McKee 1993). Representative temperatureand velocity profiles of all these types are shown in Figs. 1 and 2. The C-, C^* -, and CJ-typeshocks can only be modelled with a multi-fluid code, as their internal dynamical, thermal, and chemical structure requires the inclusion of ion-neutral drift. The basic division between these types of shocks arises from how the shock velocity compares tothe characteristic speeds in the medium and the interplay between heating and cooling. The three characteristic MHD speeds, the fast, intermediate and slow modes, allow for fast-, intermediate-, and slow-mode shocks (see e.g. Kennel et al. 1989; Lehmann & Wardle 2016), but we focus only on fast-mode shocks in this work. In order to have a fast-mode MHD shock, the fluid velocity must be greater than the total fluid fast magnetosonic speed ($c_{\text{f}}=\sqrt{c_{\text{s}}^2+B^2/4\pi {\rho }_{\text{n}}}$, for magnetic field strength B and where the neutral mass density is much greater than the charged fluid mass density ρ_n ≫ ρ_c in a weakly ionised gas, and sound speed $c_{\text{s}}=\sqrt{k_{\text{B}}T/\mu m_{\text{H}}}$ for temperature T, Boltzmann constant k_B, mean particle mass μ and proton mass m_H). However, if the shock velocity is below the magnetosonic speed of the charged fluid ($c_{\text{ims}}=\sqrt{c_{\text{s}}^2+B^2/4\pi {\rho }_{\text{c}}}$), the charged fluid decouples from the neutrals and streams ahead to heat and decelerate the gas via collisions in a magnetic precursor (Mullan 1971; Draine 1980). In this case, if the cooling is efficient enough to keep the local sound speed below the local fluid velocity, then all fluid variables remain continuous throughout resulting in a C-type shock solution. If the cooling cannot keep up with the heating there will be a transition from the supersonic fluid to a subsonic fluid. In some conditions there exists a unique physical solution that crosses this sonic point smoothly giving a C^* -type shock. Otherwise a discontinuous transition occurs, mediated by viscosity over length scales of the order of the mean free path, resulting in CJ-type shock. These C^* and CJ shock solutions require atypical numerical methods due to large gradients near the sonic point (e.g. Chernoff 1987; Roberge & Draine 1990) and were recently implemented in the Paris–Durham code (Godard et al. 2019). Finally, if the shock velocity is greater than the ion magnetosonic speed, there can be no drift between the fluid types and the fluid transitions without a magnetic precursor discontinuously due to viscosity resulting in a J-type shock. Viscous heating in J-type shocks can generate temperatures large enough to produce UV photons that can heat and ionise the gas ahead of the shock in a radiative precursor. As the self-generated UV photons affect the shock structure, and the shock structure is required to compute the UV radiative transfer, iterative methods are required to self-consistently model these shocks paying particular attention to non-equilibrium effects as gas is advected through the shock generated UV field (e.g. Shull & McKee 1979; Sutherland & Dopita 2017). This method has been implemented in the Paris–Durham code in Paper I.

The shock code assumes that the magnetic field is transverse to the direction of shock propagation and so we cannot compute oblique shocks, where there is an angle between the magnetic field and shock plane. However, these transverse shock models accurately describe oblique shocks for J-type shocks with the large Alfvénic Mach numbers considered in this work (Hollenbach & McKee 1979). For the C-type shocks, oblique shocks can be approximated by purely transverse models with a modified total field strength (Wardle & Draine 1987).

Fig. 1 Shock profiles of ion (red solid) and neutral (blue solid) fluid velocities, local sound speed (black solid), and neutral fluid temperature (green dash-dotted). The shocks propagate from right to left, and the velocities are given in the reference frame in which the shock is stationary. These shocks were computed with magnetic parameter b = 1 and no external radiation field, G0 = 0. It should be noted that both the distance and temperature axes are different for the C- and J-type shocks, and that the logarithmic scaling can give the false impression that both types of shocks have similar sizes for a given pre-shock density. We have excluded profiles for shocks with velocities Vs = 40 and 60 km s−1 for brevity.

Fig. 2 Same as Fig. 1, but for shocks with external radiation field G0 = 103.

2.2 Grids of shock models

The interstellar and circumgalactic medium displays an extraordinary diversity of physical conditions. In order to give a robust interpretation of observations we therefore model shocks over a broad range of conditions, including a range of densities, velocities, magnetic field strength, and external UV radiation fields. We consider shock velocities in the range 5–80 km s⁻¹, with models every 5 km s⁻¹ up to 30 km s⁻¹ and then every 10 km s⁻¹ up to 80 km s⁻¹, with pre-shock densities ranging from 10² cm⁻³ every half dex to 10⁶ cm⁻³. In this section, we introduce three grids of shock models varying the magnetic parameter b defining the pre-shock transverse magnetic field strength $B_0=b{\left(n_{\text{H}}/{\text{cm}}^{-3}\right)}^{1/2}$, and the radiation parameter G₀, which is a scaling factor² applied to the standard UV radiation field of Mathis et al. (1983). The three grids³ have: (i) a moderate magnetic field strength (magnetic parameter b = 1) and no external radiation (radiation parameter G₀ = 0), with profiles shown in Fig. 1; (ii) a weak magnetic field strength (magnetic parameter b = 0.1) and no external radiation (radiation parameter G₀ = 0), with profilesshown in Appendix A; and (iii) a moderate magnetic field strength (magnetic parameter b = 1) and strong external radiation⁴ (radiation parameter G₀ = 10³), with profilesshown in Fig. 2.

In total, these grids comprise 297 shock models. The elemental abundances and chemical network for these models are the same as in Paper I. The initial abundances of all chemical species are then computed by evolving a parcel of fixed density gas until thermal and chemical equilibrium for each of the nine densities. For the models with external radiation, there is an additional computation where the parcel isslowly advected away from the source of radiation until the intervening material gives an attenuation of the UV field of A_V = 0.1 magnitudes.The shock computation therefore begins with an initial state that depends on the density and radiation field considered. An analysis of the radiation field dependence for the initial temperature, electronic fraction, molecular fraction, and grain depletion can be found in Fig. 2 of Godard et al. (2019).

Here, we focus on the comparison of grids with or without external radiation, and leave the weak magnetic field models to Appendix A. To prefigure the focus on CH⁺ in this work, we present here its column densities in the shocks of these grids. The key shock parameters are summarised in Table 1.

2.2.1 Shocks with no external UV irradiation: G₀ = 0

Figure 1 shows profiles of the neutral and ionised fluid velocities and neutral fluid temperature for the shocks with moderate magnetic field strength, b = 1, and no external radiation field, G₀ = 0. For this grid the ion magnetosonic speed c_ims ~ 24 km s⁻¹, and so there is a separation of shock types into C- (V_s ≤ 20 km s⁻¹) and J-type (V_s ≥ 25 km s⁻¹). Following Godard et al. (2019), both neutral and charged grains are assumed to be fully coupled with the charged fluid. The mass density of the charged fluid is dominated by dust grains rather than atomic or molecular ions. It follows that the ion magnetosonic speed is roughly constant and low compared to a dust-free gas (Guillet et al. 2007). The neutral and ion fluid temperatures and velocities are decoupled in the C-type shocks, whereas they act as a single fluid in J-type shocks. This decoupling has important chemical implications, for example the kinetic energy provided by the ion-neutral velocity difference boosts endothermic reactions. It should be noted that the logarithmic scale gives the impression that the J-type shocks are of comparable thickness to the C-type shocks. However, the length scales for the J-type shocks to cool down to 100 K range from 10¹⁴ –10¹⁶ cm, while for the C-type shocks they range from 10¹⁵ –10¹⁸ cm. The corresponding shock timescales are less dependent on shock type, with values ranging from 100 years at high density to 10⁵ years for low density shocks.

The top panel of Fig. 3 shows the CH⁺ column densities through each shock in this grid, for the region of the shock before it cools down to a threshold temperature of 100 K. The column densities range from 3 × 10³ to 5 × 10¹³ cm⁻². For the J-type shocks (V_s ≥ 25 km s⁻¹), the column density increases with either density or velocity. There is a clear change in behaviour in the C-type shocks (V_s ≤ 20 km s⁻¹), with the column density being weakly dependent to velocity but inversely proportional to density.

Fig. 3 Column density of CH+ for shocks with magnetic parameter b = 1 and no external radiation field G0 = 0 (top) or strong external radiation field G0 = 103 (bottom). The dashed lines mark the boundaries between C-, C*-, CJ-, and J-type shocks. The coloured circles show exact values for each model in the grids.

2.2.2 Shocks with external UV irradiation: G₀ = 10³

In the presence of an external UV field, we model the shocks as propagating perpendicularly to the source of UV. In practice this means that we take the external UV field to have been attenuated by a layer of gas with A_V = 0.1 magnitudes,and the field is held constant at each point in the shock. Figure 2 shows profiles of the neutral and ionised fluid velocities and neutral fluid temperature for the grid with b = 1, G₀ = 10³. The ion magnetosonic speed c_ims ~ 21–24 km s⁻¹ for these shocks, and so the shocks with velocity V_s ≥ 25 km s⁻¹ are all J-type as in the G₀ = 0 grid. Below this velocity we find the full variety of C-, C^*- and CJ-type shocks, in agreement with Godard et al. (2019). Despite this variety of shock types, the length scales to cool down to 100 K are all similar for a given density, ranging from 10¹⁴ –10¹⁷ cm over the whole grid with no strong dependence on shock type. The corresponding shock timescales range from 100 yr to 10⁵ yr.

The bottom panel of Fig. 3 shows the CH⁺ column densities through each shock in this grid, for the region of the shock before it cools down to a threshold temperature of 100 K. In shocks where the external radiation field is able to maintain a temperature above 100 K in the post-shock, we instead take the temperature threshold to be 10% larger than the final temperature. The column densities range from 6 × 10⁹ to 8 × 10¹³ cm⁻². The column density roughly increases with preshock density until n_H = 10⁴ cm⁻³ before dropping off at densities larger than 10⁵ cm⁻³. There is a surprising continuity across the border between C^* - and J-type shocks.

2.2.3 Comparison of shock grids

A consequence of the broad parameter space that we explore is the occurrence of the full diversity of MHD shock types across the grids of models. All MHD shock types, J-, C-, C^* -, and CJ-type, arise in the grid with moderate magnetic field strength and a strong external radiation field. The UV field increases the ionisation fraction and hence the ion-neutral collisional heating rate, while strongly photodissociating molecules such as H₂ and CO, reducing the cooling rate. Thus the cooling cannot keep the local sound speed below the fluid velocity in sections of the parameter space and a transition takes place across the sonic point giving C^* -, and CJ-type shock solutions. Taking away the radiation field removes this possibility, and the shock grid contains just C- and J-type shocks. Finally, taking a weak magnetic field (as in the grid in Appendix A) reduces the total fluid ion-magnetosonic speed to below 5 km s⁻¹, and so that shock grid contains just J-type shocks.

The peak temperature for the J shocks is independent of pre-shock density and temperature. It is accurately estimated using the Rankine–Hugoniot jump conditions, assuming an adiabatic index of 5∕3 and shock velocity much larger than the Alfvén velocity ($v_{\text{A}}=B/\text{}\text{}\sqrt{4\pi \rho }$ for fluid mass density ρ), as $$T_{\text{max}}~53\text{K}\left(\frac{\mu }{\text{2}.\text{33}\text{a}.\text{m}.\text{u}.}\right){\left(\frac{{\text{V}}_{\text{s}}}{\text{km} {\text{s}}^{\text{-1}}}\right)}^{\text{2}},$$(1)

where μ is the mean particle mass in a.m.u. In a fully molecular medium, μ = 2.33 a.m.u. and we recover the peak temperatures in the J-type shocks in all the shocks with no external irradiation. In the shocks with external irradiation, the picture is more complex. At densities lower than 10^3.5 cm⁻³ the external UV field fully photodissociates H₂ and the atomic gas has mean mass μ ~ 1.3 a.m.u. At densities larger than 10⁵ cm⁻³ the H₂ is fully self-shielded and μ ~ 2.33 a.m.u. However, at high velocities (V_s ≥ 50 km s⁻¹) the combined external and self-generated UV field is strong enough in the radiative precursor to fully dissociate H₂ and the pre-shock fluid enters the shock front as an atomic medium with mean particle mass μ ~ 1.3 a.m.u. Hence, the peak temperatures in these externally irradiated J-type shocks can be up to ~2 times lower than the same shocks with no external irradiation. In addition, H₂ is a key coolant in a warm, fully molecular gas. Hence the pre-shock gas in the radiative precursors with UV fields that fully dissociate H₂ are heated to 10⁴ K, two orders of magnitude larger than any pre-shock temperatures in shocks with no external irradiation.

The chemistry controlling the production of CH⁺ in shocks has been explained in detail in Godard et al. (2019). In short, CH⁺ formation requires the presence of H₂ and is linked to the available C⁺, while its destruction is dominated by H or H₂. The C⁺ abundance is enhanced by UV photoionisation of C. Hence the shocks with the external radiation field produce larger column densities of CH⁺ than those shocks with no external radiation field at all densities and velocities. This is especially true at low velocities (e.g. V_S ≤ 25 km s⁻¹), where the column densities of CH⁺ are boosted by up to eight orders of magnitude with an external radiation field.

3 Radiative transfer of CH⁺

The formation and destruction pathways of CH⁺ are described in detail in Godard et al. (2019). We note here that the shock code calculates the formation rate taking into account the excited rovibrational states of H₂. In this section we investigate the line emission from CH⁺ generated in these shocks, which we calculate using an updated version of the large velocity gradient (LVG) code of Gusdorf et al. (2008). We first summarise the details of the excitation and radiative transfer calculation. Then we show how the character of the emission is affected by, for example, the shock type, pre-shock density, shock velocities, and presence or absence of an external radiation field.

3.1 Excitation processes

The full details of the included processes, extrapolations, and testing of the excitation and radiative transfer of CH⁺ are given inNeufeld et al. (2021), but we give a summary here for completeness. Updated Einstein A coefficients are taken from Changala et al. (2021). We include the excitation to rovibrational states by non-reactive collisions with e⁻ (Hamilton et al. 2016), H (Faure et al. 2017), and He (Hammami et al. 2009). Rates for collisions with H₂ are deduced from the CH⁺-He data within the framework of the rigid rotor approximation as done by Godard & Cernicharo (2013). CH⁺ can also form in an excited state in the reaction $${\text{C}}^{+}+{\text{H}}_2\to {\text{CH}}^{+}(v,j)+\text{H},$$(2)

which is the dominant formation pathway of CH⁺. Since CH⁺ is short lived it is important to include this formation pumping in the excitation model. Furthermore, since the destruction in the reverse reaction occurs on timescales longer than its radiative decay, chemical de-pumping is not necessary to include. Radiative pumping in the infrared rovibrational lines is modelled within the LVG approximation. Radiative pumping in the optical electronic lines followed by a fluorescent cascade is modelled taking into account self-shielding processes using the FGK approximation (Federman et al. 1979).

With all these excitation processes in place, the radiative transfer is computed from the outputs of the Paris–Durham shock code, for example, the local temperature, ion fluid velocity, sound speed, CH⁺ density, and densities of collisional partners. This calculation is an update on that in Godard et al. (2019), with the main differences being the CH⁺– e⁻ collisional rates that are orders of magnitude different at low temperatures and the inclusion of chemical pumping and the impinging isotropic radiation field. In the end the only free parameters are the number of CH⁺ rovibrational levels to include in the computation and the microturbulent velocity dispersion, σ_turb, which we set to 100 levels and 2.5 km s⁻¹, respectively,for the rest of this work.

Fig. 4 Specific intensities of CH+ pure rotational lines v = 0 − 0 J = n → n − 1 up to n = 10 emitted parallel to the propagation direction of a (top) C-type shock with velocity Vs = 20 km s−1 or (bottom) J-type shock with velocity Vs = 50 km s−1, both with pre-shock density 104 cm−3, magnetic parameter b = 1, and no external radiation field G0 = 0.

3.2 CH⁺ emission line profiles

In Fig. 4, we show the specific intensity line profiles for the pure rotational lines of CH⁺ J =n → n − 1 up to n = 10 emitted parallel to the shock propagation direction for a 20 km/s C-type shock (upper panel) and a 50 km s⁻¹ J-type shock (lower panel) with no external radiation field. The Einstein A coefficients, central line frequency and upper energy level of these lines are given in Table 2. The emission from the C-type shock gives a skewed line profile, showing emission at velocities between the shock velocity and 0 km s⁻¹. This reflects the character of C-type shocks in which the gas is heated and decelerated smoothly. The line shape near the peak is given by the microturbulent velocity dispersion, because the peak CH⁺ production and excitation occurs where the thermal velocity dispersion is less than σ_turb = 2.5 km s⁻¹. In contrast, deceleration in the J-type shock occurs in the high temperature viscous jump that strongly dissociates all molecules. The CH⁺ emission therefore uniquely arises from the region where the shock has cooled, the velocity is ~ 0 km s⁻¹, and molecules have reformed. This region has cooled down to a temperature ~ 10³ K and so the thermal velocity dispersion is also smaller than σ_turb.

Fig. 5 Line intensities of CH+ J = n → n − 1 up to n = 10 emitted parallel to the propagation direction of shocks with velocity Vs = 30 km s−1, densities nH = 102–106 cm−3, and radiation parameter G0 = 0 (top) or G0 = 103 (bottom). We note that each figure has a different y-axis scale.

3.3 CH⁺ integrated line intensities

Figure 5 shows the integrated line intensities of the pure rotational lines of CH⁺ up to J = 10 → 9, emitted parallel to the shock propagation direction, for shocks with fixed velocity V_s = 30 km s⁻¹ over the density range, comparing the shocks with no external radiation field (upper panel) to those with an external radiation field G₀ = 10³ (lower panel). Externally irradiated shocks produce CH⁺ much more efficiently than non-irradiated shocks, and so their line intensities at a given density are always orders of magnitude stronger than shocks with no external UV field. The intensities for a given shock are remarkably similar over the rotational ladder, meaning an observation of multiple CH⁺ pure rotational lines does not distinguish different pre-shock densities.

Figure 6 shows the integrated line intensities of the pure rotational lines of CH⁺ up to J = 10 → 9 for shocks with fixed density n_H = 10⁴ cm⁻³ over the velocity range, comparing the shocks with no external radiation field (upper panel) to those with an external radiation field G₀ = 10³ (lower panel). In the upper panel there is strong CH⁺ emission for shocks with velocity V_s ≥ 40 km s⁻¹. At these velocities the shocks generate a strong Lyα flux, emphasising the need for a UV radiation field to initiate the CH⁺ chemistry by ionising C. With such a UV field provided externally, the low-velocity shocks in the lower panel produce just as much CH⁺ emission at shock velocity V_s = 25 km s⁻¹ as the non-irradiated shocks with shock velocity V_s = 60 km s⁻¹. With the exception of the lowest-velocity shocks that give negligible CH⁺ emission, the intensities for a given shock are again remarkably similar. Hence the CH⁺ rotational ladder does not strongly trace the shock velocity.

The 1–0 line of CH⁺ is the most commonly observed line, and so in Fig. 7 we show the line intensity over the entire grid of shocks with no external radiation field (upper panel) and for externally irradiated shocks (lower panel). This line is generally brighter in externally irradiated shocks over the whole parameter space, and is thus a good tracer for the presence of a strong UV field. In particular at low velocities (V_s ≤ 30 km s⁻¹) the intensity is up to seven orders of magnitude stronger in externally irradiated shocks than shocks with no external UV field.

In Fig. 8, we show the fraction of the shock kinetic flux reprocessed into all CH⁺ rovibrational lines for shocks with no external radiation field (upper panel) and for externally irradiated shocks (lower panel). These panels show that CH⁺ can be a significant pathway of the reprocessing of mechanical energy in shocks. In shocks with an external radiation field, up to 1% of the shock kinetic flux can escape in the rovibrational lines of CH⁺.

Fig. 6 Line intensities of CH+ J = n → n − 1 up to n = 10 emitted parallel to the propagation direction of shocks with velocities Vs = 5–80 km s−1, density nH = 104 cm−3, and radiation parameter G0 = 0 (top) or G0 = 103 (bottom).

4 Ensembles of shocks

Observations of extragalactic molecular line emission (e.g. Guillard et al. 2009; Falgarone et al. 2017) with extremely large linewidths (~ 1000 km s⁻¹), to first approximation, imply temperatures at which such molecules cannot exist. These observations have inspired an explanatory framework in which energy injected at large scales enters a turbulent cascade that dissipates in structures at smaller scales and much lower velocities. Here we make a first step in using this framework to make quantitative predictions of the mechanical energy and physical conditions of an observed system, assuming these dissipative structures are MHD shocks.

We consider a toy model in which the emission is due to shocks with velocities much lower than the observed linewidth. In particular, the velocities are low enough to produce molecular emission lines, which are then spread over the observed velocity range by bulk flows. We stress that the emission model we develop can be used with any line whose radiative transfer can be calculated in these shocks. As an illustrative example we apply the model to a particular observation, the CH⁺(1−0) line from the high-redshift starburst galaxy The Cosmic Eyelash. We show how to use the ensemble model with one shock output, giving a prediction of the total rate of kinetic energy dissipated in shocks and the minimum number of shocks needed in the observational solid angle. Finally, we show how a grid of shock models can then be used to constrain physical conditions and mechanical energy of an observed system.

Fig. 7 CH+(1−0) line intensityfor shocks with magnetic parameter b = 1 and radiation parameter G0 = 0 (top) or G0 = 103 (bottom).

4.1 Toy model

All the symbols used in this section are given in Table 3. A cartoon of our toy model is shown in Fig. 9, in which $\mathcal{N}$ shock surfaces are captured within an observational solid angle Ω. We consider each of the shocks in the ensemble to have the same shock velocity V_S, pre-shock density, and to subtend the same solid angle, Ω₁, on the planeof the sky. This model is a highly simplified representation of a large-scale injection of energy (for example, the infall/outflow of material onto galaxies like the Eyelash, or the galaxy collision in Stephan’s Quintet) exciting a dominant shock velocity in clumps dispersed according to a Gaussian distribution in velocity-space with linewidth equal to the observed linewidth.

Each shock is randomly placed in the observational solid angle. The specific intensity of any line emitted by one shock (e.g. see Fig. 4) is centred at a velocity V randomly drawn from a Gaussian distribution. For simplicity we replace the shock emission profile with a top-hat profile with linewidth Δν, line intensity I (e.g. see Fig. 7), and optical depth τ₁. The amplitude of the top-hat specific intensity of one shock is then I_ν,1 = I∕Δν. Δν is derived from the approximate velocity width ΔV of a shock line profile, found in Sect. 3.2 to be dominated by the microturbulent velocity dispersion, and set to ΔV = 2σ_turb = 5 km s⁻¹. Two shocks placed at velocities within this linewidth will both contribute to the intensity at V. In this bin, if there is no spatial overlap of shock surfaces then we model the ensemble flux density with $$\begin{array}{ll}{F}_{\nu }(V)\hfill & =\stackrel{\text{unattenuated}\text{background}}{\overbrace{I_{\nu ,\text{c}}\left(\Omega -N(V){\Omega }_1\right)}}+\stackrel{\text{attenuated}\text{background}}{\overbrace{I_{\nu ,\text{c}}\exp \left(-{\tau }_1\right)N(V){\Omega }_1}}+\stackrel{\text{shock}\text{emission}}{\overbrace{I_{\nu ,1}N(V){\Omega }_1}},\hfill \end{array}$$(3)

where N(V) is the number of shocks in the velocity bin of size ΔV and I_ν,c is the background continuum specific intensity. The first term gives the background emission that passes through unattenuated by shock slabs at velocity V, the second term gives the background emission attenuated by a single shock layer scaled by the number of shocks in the bin, and the third term gives the emission generated in the shocks themselves. On the other hand, if there are many shocks overlapping in both velocity and space, the ensemble emission is given instead by $$F_{\nu }(V)=\stackrel{\text{attenuated}\text{background}}{\overbrace{I_{\nu ,\text{c}}\exp \left(-N(V){\tau }_1\right)\Omega }}+\stackrel{\text{shock}\text{emission}}{\overbrace{I_{\nu ,1}{\Omega }_1{\sum }_{n=0}^{N(V)-1}\exp \left(-n{\tau }_1\right)}},$$(4)

where now there is no unattenuated background. The first term gives the background attenuated by as many shock layers as there are in the velocity bin and the second term gives the emission generated in the shocks attenuated by the other shock layers between them and the observer in the same velocity bin. Equation (4) is given here for completeness, but we find for our application in the next section that we can always use Eq. (3). This must be verified for each application, and can be decided with the following considerations.

The single shock peak emission is shifted to velocities randomly drawn according to a Gaussian distribution $$f(V,\sigma )=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{V^2}{2{\sigma }^2}\right).$$(5)

After constructing many ensembles, the number of shocks in the velocity bin around V will follow a Poisson distribution with the average given by $$\overline{N}\left(V\right)=\mathcal{N}\Delta Vf(V,\sigma )$$(6)

and the standard deviation by $${\sigma }_{\Delta V}\left(V\right)=\sqrt{\overline{N}\left(V\right)}.$$(7)

Equation (3) can only be used if the shocks within a velocity bin are not overlapping in space. This can be verified by making sure the typical number of shocks overlapping in velocity and space, $${\overline{N}}_{\text{S}}\left(V\right)=\overline{N}\left(V\right)\frac{{\Omega }_1}{\Omega }<1;$$(8)

otherwise, we use Eq. (4). It is worth noting that if the cumulative optical depth through the shock layers, ${\overline{N}}_{\text{S}}\left(V\right){\tau }_1$, is small, then Eq. (4) reduces to Eq. (3). That is, both the background and shock emission are unattenuated by the shock layers in the ensemble. These conditions need to be checked when computing the shock ensemble emission to reproduce a particular observation.

Fig. 8 Total CH+ flux as a fraction of shock kinetic flux for shocks with magnetic parameter b = 1 and radiation parameter G0 = 0 (top) or G0 = 103 (bottom).

Fig. 9 Cartoon of ensemble of blue- and redshifted shocks within an observational solid angle (black and white circle). The left and right panels represent two scenarios with the same total shocked surface area but a different surface area per shock.

Fig. 10 CH+ (1−0) flux density for the Eyelash galaxy from ensembles of 1000 (left) or 10 000 (right) shocks each with velocity Vs = 50 km s−1 and pre-shock density nH = 104 cm−3. The peak of the average ensemble (black) has been fitted to 5.2 × 10−26 erg s−1 cm−2 Hz−1, and the particular ensembles (blue) have a resolution ΔW = 50 km s−1.

4.2 Example application: the Eyelash Starburst galaxy

In this section, we show how the ensemble toy model can be used to extract information from an extragalactic observation of extremely broad CH⁺(1-0) emission line towards the Eyelash galaxy. To apply the ensemble model, we need the peak value of the background subtracted flux density integrated over the galaxy $F_{\nu ,\text{p}}^{\text{obs}}$, the observational solid angle of the galaxy Ω_obs, the background continuum flux density $F_{\nu ,\text{c}}^{\text{obs}}$, and the line velocity dispersion σ_obs. For further analysis we also use the observational velocity resolution ΔW, gravitational lensing magnification factor μ, and signal-to-noise ratio (S/N), of the observed line. These values are all summarised in Table 4.

We illustrate the application by using the emission from a single shock with velocity V_s = 50 km s⁻¹ propagating into gas with density n_H = 10⁴ cm⁻³, taken from the grid of models with magnetic parameter b = 1 and no external radiation field. The top-hat shock specific intensity of the CH⁺(1−0) line I_ν,1 ~ 7.2 × 10⁻¹⁴ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹ and optical depth τ₁ ~ 0.5. The emission from successive shock layers overlaps in velocities of the order of the width of their emission line profiles. In Sect. 3.2, we found that the line broadening is dominated by the microturbulent velocity dispersion, so we chose an overlap velocity bin ΔV = 5 km s⁻¹ for all shocks.

With the observational parameters chosen and shock radiative transfer computed, we can then constrain the total shock surface area. Setting $N(V)=\overline{N}(V)$ in Eq. (3), the background subtracted ensemble flux density at the peak (V = 0) averaged over many instances is then $$F_{\nu }(0)-F_{\nu ,\text{c}}^{\text{obs}}=\left[\frac{F_{\nu ,\text{c}}^{\text{obs}}}{{\Omega }_{\text{obs}}}\left(\exp \left(-{\tau }_1\right)-1\right)+I_{\nu ,1}\right]{\Omega }_1\mathcal{N}\Delta Vf(0,{\sigma }_{\text{obs}}).$$(9)

Equalising the right hand side of this expression to the observed value $F_{\nu ,\text{p}}^{\text{obs}}$ gives the total shock solid angle $${\Omega }_1\mathcal{N}=\frac{F_{\nu ,\text{p}}^{\text{obs}}{\sigma }_{\text{obs}}\sqrt{2\pi }}{\Delta V}{\left[\frac{F_{\nu ,\text{c}}^{\text{obs}}}{{\Omega }_{\text{obs}}}\left(\exp \left(-{\tau }_1\right)-1\right)+I_{\nu ,1}\right]}^{-1}.$$(10)

We note that if a shock is a net absorber of the background continuum then ${\Omega }_1\mathcal{N}$ will be negative. In this case, the observation cannot be reproduced no matter the value of $\mathcal{N}$ and we can rule out the ensemble with that particular shock model. For this shock the total shock solid angle, ${\Omega }_1\mathcal{N}=2.1\times {10}^{-10}$ sr, can then be used in Eq. (8), giving a maximum average overlap of shocks in space and velocity ${\overline{N}}_S(V=0)~2.3\times {10}^{-2}$. Hence Eq. (3) can be used at all velocities over the entire emission line.

The total shock surface gives the total rate of energy dissipated by the shocks in this ensemble $${\mathcal{L}}_K=\frac{1}{2}{\rho }_0{V}_{\text{s}}^3\mathcal{A},$$(11)

where ρ₀ is the pre-shock mass density, and $\mathcal{A}$ is the total shock surface area. The area is $$\mathcal{A}={\Omega }_1\mathcal{N}{D}_{\text{A}}^2,$$(12)

where the angular diameter distance D_A = 1734 Mpc is calculated assuming the Planck cosmological parameters (Planck Collaboration VI 2020) at redshift z = 2.3. For this, 50 km s⁻¹ shock and distance to this galaxy, the area $\mathcal{A}~640$ kpc² and the dissipation rate ${\mathcal{L}}_K~8.9\times {10}^{45}$ erg s⁻¹ ~ 2.3 × 10¹² L_⊙.

After using the average ensemble flux density we can use the constrained ${\Omega }_1\mathcal{N}$ to build a particular ensemble by replacing N in Eq. (3) with a different integer drawn from the Poisson distribution for each velocity bin. Two particular ensembles are shown in Fig. 10 to illustrate the effect of different total number of shocks. The ensemble shock emission has been further binned to ΔW = 50 km s⁻¹ to match the observational resolution. The obvious difference between the scenarios is the relative amplitude of the shot noise, given by $$\frac{{\sigma }_{\Delta W}(V)}{\overline{N}(V)}=\frac{1}{\sqrt{\overline{N}(V)}},$$(13)

where $\overline{N}(V)$ is computed with the observational resolution ΔW. Since this noise is a real feature of the discrete velocity sampling of shocks, it could in principle be measured if the noise in the line centre signal were observed to be larger than in the continuum far enough from the line. In Falgarone et al. (2017), after removing Gaussian fits to the emission and absorption lines, the amplitude of the residuals was found to be constant over the band, and therefore the noise was dominated by instrumental noise. This gives an upper limit on the shot noise and hence constrains the lower limit on the number of shocks in the galaxy, $\mathcal{N}$. The minimum number of shocks is then found by comparing the relative shot noise, Eq. (13), to the signal-to-noise ratio, leading to the inequality $$\mathcal{N}>\frac{{\sigma }_{\text{obs}}\sqrt{2\pi }}{\Delta W}{(\text{S/N})}^2.$$(14)

We note that this is independent of the radiative transfer calculation for any particular shock model and is just a result of discrete sampling from a Gaussian distribution. For this observation, the signal-to-noise ratio ~ 6, giving a minimum number of shocks of ${\mathcal{N}}_{\text{min}}~1000$. This value can be combined with the fit for the total shock solid angle (Eq. (10)) to give a maximum surface area of each shock in the beam, equivalent to a circle of diameter d_max ~ 900 pc. It is important that this number be much larger than the thickness of the shock layer L_sh ~ 250 au so that we indeed have a sheet geometry that can be accurately modelled with a one-dimensional shock code. For any particular observed system if further arguments can be made to constrain this maximum size of shock sheets, for example with regards to inhomogeneities of the medium, then the minimum number of shocks can be further constrained. Reversing the argument, if we demand shock surfaces to have a minimum equivalent diameter of d_min = 100 L_sh ~ 10⁻¹ pc, then we constrain the upper limit of the number of shocks at ${\mathcal{N}}_{\text{max}}~{10}^{11}$.

Fig. 11 Average number of shocks in the velocity bin ±2.5 km s−1 around V = 0 overlapping in space according to Eq. (8) (top), and peak optical depth τ1 for each shock for the grid with G0 = 0 and b = 1, applied to the Eyelash galaxy observation (bottom).

4.3 Application with a grid of shock models

We now repeat the treatment of the previous section over the wide range of shock models discussed in Sect. 2.2 on the same observational source, the Eyelash starburst galaxy. Here we compare the grids of shocks with magnetic parameter b = 1 and external radiation parameters G₀ = 0 or G₀ = 10³.

First, we consider the choice between Eqs. (3) and (4). As previously, we assume Eq. (3) in order to derive the total shock solid angle using Eq. (10). In Figs. 11 and 12, we show the maximum number of shock surfaces that overlap in velocity and space (top panels) and the peak optical depth through each shock (bottom panels). We note that the blank space in the top panels represents the shock models seen in absorption and which therefore cannot reproduce the observation. The top panels show that Eq. (3) is immediately verified (${\overline{N}}_{\text{S}}<1$) for the shocks at large velocities and densities for the G₀ = 0 grid, and for almost all the shocks in the externally irradiated grid. For the models requiring many overlapping shocks, if we consider also the peak optical depth there is no model in the two grids with a cumulative optical depth ${\overline{N}}_{\text{S}}{\tau }_1>1$. Hence, we consider Eq. (3) accurate.

In Fig. 13, we show the total kinetic energy dissipation rate required to reproduce the observation with the ensemble model for shocks with external radiation parameters G₀ = 0 (top) or G₀ = 10³ (bottom). We find that no matter the shock type, the presence of a UV radiation field is key to providing the C⁺, via photoionisation of C, required for the shock heated gas to produce significant levels of CH⁺. With no external radiation field only the higher-velocity shocks (V_s > 30 km s⁻¹) can produce enough UV, in Lyα and Lyβ lines, to sufficiently photoionise C (see Paper I). In contrast, the strong external radiation field allows lower-velocity shocks to efficiently produce CH⁺ emission over a wide range of densities. These figures can be understood as representing the efficiency of the production of CH⁺ emission with respect to the available shock mechanical energy. In the parameter domain, there are two distinct regions in which externally irradiated shocks are one to two orders of magnitude more efficient at producing CH⁺(1−0) emission than the high-velocity shocks with no external radiation (region A in Fig. 13): the C^* - and CJ-type shocks with velocities V_s =5–15 km s⁻¹ and proton densities n_H = 10⁴–10⁵ cm⁻³, and J-type shocks with velocities V_s = 25–40 km s⁻¹ and proton densities n_H = 10⁴–10⁵ cm⁻³ (regions B and C in Fig. 13, respectively). Both these regions require a mechanical energy injection rate ${\mathcal{L}}_K~{10}^{11}$ L_⊙.

The bolometric luminosity of the Eyelash galaxy ${\mathcal{L}}_{\text{bol}}=2.3\times {10}^{12}$ L_⊙ (Ivison et al. 2010). This value gives a strong upper limit of the mechanical energy budget. The lowest density medium possibly responsible for shocked CH⁺ emission is n_H ~ 10^3.5 cm⁻³. If the shocks occur in an environment with no external radiation field, only high-velocity J-type shocks (V_s > 50 km s⁻¹) can explain the CH⁺ emission (region A). It should be noted, however, that these shocks would require a mechanical energy injection rate comparable to the bolometric luminosity. If a strong UV field is present, then the full range of velocities could reproduce the line. The lowest-velocity shocks (V_s < 40 km s⁻¹, regions B and C) would require the lowest injection rate of energy, with an injection rate ~ 20 times smaller than ${\mathcal{L}}_{\text{bol}}$. Such low velocities could be expected theoretically, as simulations of MHD turbulence in circumgalactic and interstellar medium conditions indicate that the dissipation of turbulence by shocks is dominated by the lowest velocities with a sharp drop-off towards higher velocities (Lehmann et al. 2016; Park & Ryu 2019).

While we have gained a lot of information from the analysis of just one line, the emission from these shocks is not limited to CH⁺. They produce other tracers, such as strong emission from rovibrational lines of H₂ that could be observed with the JWST. Therefore, observations of these additional tracers can be used to disentangle these two irradiation scenarios and give a definite value of the mechanical energy injection rate. In the next section we provide predictions of such tracers.

Fig. 12 Average number of shocks in the velocity bin ±2.5 km s−1 overlapping in space according to Eq. (8) (top), and peak optical depth τ1 for each shock for the grid with G0 = 103 and b = 1, applied to the Eyelash galaxy observation (bottom).

Fig. 13 Total shock dissipation rate in the Eyelash galaxy for ensembles of shocks with no external radiation field G0 = 0 (top) and forshocks with a strong external radiation field G0 = 103 (bottom). We note that one solar luminosity L⊙ = 3.8 × 1033 erg s−1. The violet labels ‘A’, ‘B’, and ‘C’ mark the regions of most efficient models discussed in Sect. 5.2.

5 Predictions

In this section we consider shock tracers other than CH⁺, for example Lyα and H₂. First, we look at emission coming from a single shock layer. Then we consider these other tracers in the framework of our shock ensemble toy model applied to the Eyelash galaxy.

5.1 Lyα and H₂ from a single shock

In Paper I, we defined a photon flux normalised to the standard interstellar radiation field (ISRF) of Mathis et al. (1983) given by $$G_{\text{eff}}=\frac{{\displaystyle \int \text{d}}\nu c{u}_{\nu }/h\nu }{1.55\times {10}^8{\text{ph/s/cm}}^{\text{2}}},$$(15)

where u_ν is the energy density of the radiation field at frequency ν, c is the speed of light, and h is the Planck constant. For the shock models considered in that work (pre-shock proton density of 10⁴ cm⁻³), we found a critical shock velocity ~ 32 km s⁻¹ above which this photon flux integrated over just the Lyα line escaping into the precursor region is greater than the entire photon flux of the standard Mathis interstellar radiation field, that is, G_eff > 1. The top panel of Fig. 14 shows the escaping Lyα photon flux for shocks with b = 1 and no external radiation field. This confirms the previous result for the critical velocity for shocks with pre-shock density of 10⁴ cm⁻³. This critical velocity depends inversely on density, and once this threshold is crossed the Lyα flux shows a linear proportionality with density. However, no matter how large the density, the critical velocity is never below ~ 30 km s⁻¹, which corresponds to the minimum velocity at which the temperature in the viscous jump is large enough to generate enough atomic H from the collisional dissociation of H₂.

The bottom panel of Fig 14 shows the escaping Lyα photon flux relative to the kinetic flux of the shock $1/2{\rho }_0{V}_{\text{s}}^3$. This shows that once shocks become strong enough to destroy H₂ in the initial viscous jump (e.g. V_s > 30 km s⁻¹), atomic H becomes a strong coolant. At these intermediate velocities, up to 30% of the shock mechanical energy is reprocessed and escapes ahead of the shock in the form of Lyα. For a broad range of densities and velocities, this fraction is between 10% and 20%. We note that initially about twice these percentages are converted to Lyα in the shocks, but that the half of the photons that escape further into the shock are quickly reabsorbed by dust. The fractions at 80 km s⁻¹ are in line with the Lyα fluxes of the V_s = 100 km s⁻¹ shock models in the database of Alarie & Morisset (2019) using shocks computed with the MAPPINGS V code (Sutherland & Dopita 2017).

In the panels of Fig. 15, we show the line intensities of a selection of H₂ rovibrational lines for the grid of shocks with b = 1 and no external radiation field (upper panels) or external radiation field G₀ = 10³ (lower panels). These figures show that in dissociative shocks at high velocities, despite the strong destruction of H₂ in the initial temperature jump, there is strong rotational and rovibrational H₂ emission. This is due to the reformation of H₂ in the post-shock cooling layer, which pumps even the highest H₂ levels. This process only weakly depends on the gas temperature, and so the emission of H₂ in its rovibrational lines is weakly dependent on the velocity at high shock velocity. In the absence of an external UV field there is a sharp transition between low- and high-velocity shocks, near V_s ~ 20 km s⁻¹ depending on density, in the vibrational lines of H₂. Indeed, low-velocity shocks do not reach temperatures high enough to either dissociate H₂ or collisionally excite its high energy levels. With an external UV field, the vibrational lines of H₂ are excited even at low velocities. This is because the electronic levels of H₂ are pumpedby UV photons and induce a fluorescent cascade in the ground electronic state. Finally, these figures show perhaps surprisingly that strong emission from rovibrational lines of H₂ and strong Lyα emission can emerge from the same shocks over a broad range of parameters.

In Appendix B, we tabulate the integrated line intensities for a selection of H₂ and CH⁺ rovibrational lines, CO rotational lines, atomic H lines and two-photon continuum, and atomic fine structure and metastable lines. The tables are given for shock models in all three grids, at velocities at every 10 km s⁻¹ and integer logn_H. These lines are a subset of those output by the code.

Fig. 14 Lyα contribution to the UV flux escaping the shock front for shocks with magnetic parameter b = 1 and no external radiation field, G0 = 0. Top: flux relative to the Mathis ISRF, with flux equal to the standard ISRF at the 0 contour. Bottom: flux relative to the kinetic flux of the shock $1/2{\rho }_0{V}_{\text{s}}^3$ as a percentage.

5.2 Application to the Eyelash galaxy

For the shocks with no external irradiation that give ensemble total dissipation rates below the bolometric luminosity of the Eyelash galaxy (V_s ≥ 50 km s⁻¹ and n_H = 10⁴–10⁵ cm⁻³), ~15–25% of the shock kinetic energy is re-emitted in Lyα photons that escape into the pre-shock medium, as shown in Fig. 14. This means there is up to ~ 4 × 10¹¹ L_⊙ in Lyα photons permeating this region, bathing the surrounding dense, fully molecular medium. This leads to widespread heating of both dust and gas.

In Sects. 4.2 and 4.3, we pushed the limits of what information could be extracted using just the CH⁺(1-0) rotational line. We were able to constrain regions of allowable shock models within the grids with or without external irradiation. Here we consider how to further observationally constrain the mechanical energy and physical conditions of extragalactic systems. We consider 3 regions of parameters for the most efficient ensemble shock models, denoted ‘A’, ‘B’, and ‘C’ in Fig. 13. Region A contains high-velocity (V_s ≥ 50 km s⁻¹) J-type shocks, propagating in dense gas (n_H = 10⁴⁻⁵ cm⁻³) with no external irradiation driven by a mechanical energy injection rate ~ 10¹² L_⊙, region B contains low-velocity (V_s ~ 10 km s⁻¹) C^* -type shocks with an external UV field G₀ = 10³ driven by a mechanical energy injection rate ~ 10¹¹ L_⊙, and region C contains low-velocity (V_s = 25–40 km s⁻¹) J-type shocks with an external UV field G₀ = 10³ driven by a mechanical energy injection rate ~ 10¹¹ L_⊙.

For each shock, we use the outcome of Sect. 4.2 in which the ensemble CH⁺(1−0) emission is fit to the observation. We then predict the ensemble emission counterpart of rovibrational lines of CH⁺, rovibrational lines of H₂, pure rotational lines of CO, and atomic metastable and fine-structure lines. The intensities of these lines are multiplied by the total shock solid angle given by Eq. (10) to give the total line fluxes for the Eyelash galaxy. We note that the rovibrational lines of H₂ and atomic lines are computed by the shock code with an optically thin assumption, while the CO lines are computed with the LVG approximation like CH⁺.

In Fig. 16, we show a selection of total line fluxes for the efficient ensemble shock models. There are many differences between the models, suggesting that even just one more observation is able to provide strong constraints on galactic environments. Taking note that the observed wavelength, λ, is redshifted by $\lambda ={\lambda }_0\left(1+z\right)$ where λ₀ is the rest-frame wavelength and z ~ 2.3 is the cosmic redshift of the Eyelash galaxy, many of the lines from this system fall in observable wavebands of current and near-future state-of-the-art observatories, for example JWST and ALMA.

In molecular emission, H₂ rovibrational lines are good at distinguishing between shock types. In J-type shocks (regions A and C) these lines are typically an order of magnitude weaker than in the C^*-type shocks (region B). Stronger differences are found in CO lines, with two to three orders of magnitude stronger fluxes in the high-velocity ensemble in both low- and high-J rotational lines. This is due to the efficient photodissociation of CO by the broad-spectrum external radiation field in the externally irradiated models. By contrast the high-velocity model does not efficiently photodissociate CO, despite the similarly strong UV flux, due to a lack of overlap between CO absorption lines and the Lyα line. The high temperatures reached in the J-type shocks allow for stronger emission in atomic metastable lines. The low-velocity shocks only weakly excite these lines, and so they are typically around four orders of magnitude weaker than from the high-velocity ensemble. So, H₂ and atomic metastable lines probe shock type, CO lines probe the external radiation conditions, and H lines distinguish high and low-velocity shocks.

In the end, a multi-wavelength study is able to probe the physical conditions and energetics of galaxy-scale molecular outflows. An extremely broad CH⁺(1-0) emission line may uncover the presence of shocks in a turbulent cascade. In this framework, we are able to disentangle whether the ensemble of shocks is composed of high-velocity self-irradiated shocks or low-velocity externally irradiated shocks. Even just one extra line to complement the CH⁺(1−0) observation would strongly tighten constraints on the density and radiation conditions of the turbulent gas as well as give a minimum mechanical energy input from the outflow. These results open new perspectives for observational programmes using state-of-the-art facilities such as the JWST and ALMA.

Fig. 15 Line intensities of a selection of H2 rovibrational lines, for shocks with no external radiation field (top) and for radiation field G0 = 103 (bottom).

Fig. 16 Line fluxes of a selection of H2 and CH+ rovibrational lines, CO rotational lines, and atomic fine-structure and metastable lines emitted from ensembles of shocks around the Eyelash galaxy. Ensemble models are grouped by regions in Fig. 13, denoted ‘A’ (blue) for high-velocity J-type shocks with no external irradiation driven by a mechanical energy injection rate ${\mathcal{L}}_K~{10}^{12}$ L⊙, ‘B’ (purple) for low-velocity C*-type shocks with an external UV field G0 = 103 driven by a mechanical energy injection rate ${\mathcal{L}}_K~{10}^{11}$ L⊙, and ‘C’ (green) for low-velocity J-type shocks with an external UV field G0 = 103 driven by a mechanical energy injection rate ${\mathcal{L}}_K~{10}^{11}$ L⊙. Markers indicate fluxes averaged over all the models in the defined region.

6 Conclusion

We have presented models of magnetised, molecular shocks with the broadest published range of physical parameters, with shock velocities V_s = 5–80 km s⁻¹, pre-shock proton densities n_H = 10²–10⁶ cm⁻³, weak to moderate magnetic field strengths B = 1–1000 μG, and with or without a strong external radiation field, G₀ = 0 or 10³. Post-processed molecular emission lines from these shocks were used in a simple model of an ensemble of identical low-velocity shocks transported by high-velocity bulk flows, representing the dissipative structures in a turbulent cascade driven by a large-scale injection of mechanical energy. As an illustrative test case, we successfully used the ensemble model to interpret the observation of broad (σ_obs > 500 km s⁻¹) CH⁺(1−0) emission from the high-z starburst galaxy the Cosmic Eyelash. With just this one emission line, we were able to constrain the medium responsible for a shock origin of this line and suggested key atomic and molecular lines that could further constrain these conditions. We summarise the key results as follows:

1. For the physical conditions explored, the full diversity of stationary MHD shock types arise (C-, C^*-, CJ- and J-type), emphasising the requirement for a multi-fluid treatment and molecular chemical network to model molecular shocks in these conditions;

2. A significant fraction of the shock mechanical energy is reprocessed in the excitation of H. Up to ~30% of the shock kinetic flux escapes as Lyα photons;

3. Strong emission of H₂ rovibrational lines is found over a broad range of parameters and is the brightest in the same density and velocity regimes that show the strongest Lyα emission. At low velocities (V_s ≤ 40 km s⁻¹), vibrational lines are strong tracers of the presence of an external UV field;

4. CH⁺ can be a significant pathway for reprocessing the mechanical energy, radiating up to 1% of the kinetic flux in shocks irradiated by a strong external UV radiation field. The UV field, whether self-generated or already present in the shocked environment, is a key component enhancing the production and excitation of CH⁺ in shocks. With respect to the required mechanical energy input, low-velocity C* and J-type shocks (V_s ≤ 40 km s⁻¹) with a strong external radiation field (G₀ = 10³) are at least 10 times more efficient than any of the shocks with no external radiation field;

5. We have demonstrated how to use these shock models to interpret extragalactic observations of molecular emission. Our model of emission from an ensemble of shocks is able to extract the mechanical energy and physical conditions of the turbulent cascade triggered by high-velocity flows around external galaxies;

6. In the Eyelash galaxy, we find that at least 10¹¹ L_⊙ of mechanical energy must be transferred from the large-scale, high-velocity flows to 10³–10¹¹ dissipative structures at lower velocity via a turbulent cascade. This minimum amount of mechanical energy is 20 times smaller than the total bolometric luminosity and is comparable to estimates of the galaxy-scale input of mechanical energy due to, for example, supernova-driven galactic outflows (Veilleux et al. 2020). Non-irradiated shocks are barely allowed, while irradiated shocks could be a viable scenario for explaining the CH⁺ emission observed in the Eyelash galaxy. The stringent constraint given by the bolometric luminosity shows, however, that this scenario should be extended to shocks irradiated by a stronger radiation field and improved by including the contribution of photodissociation regions to the CH⁺ emission. In any case, a strong UV radiation field is mandatory to explain the CH⁺ emission observed in the Eyelash galaxy. The presence of this UV field is supported by the analysis of the spectral energy distribution of the Eyelash (Ivison et al. 2010) that is well fitted by two dust components at temperatures of 30 and 60 K;

7. We have estimated the H₂, CH⁺, CO, and atomic line counterpart from an ensemble of low-velocity C^*-type (V_s ~ 10 km s⁻¹) or J-type (V_s ~ 25 km s⁻¹) externally irradiated shocks, or high-velocity (V_s = 50 km s⁻¹) J-type shocks with no external irradiation. CO lines are particularly good at distinguishing between the two irradiation scenarios as strong photodissociation by the external UV field leads to CO emission that is more than two orders of magnitude weaker across the whole rotational ladder for the low-velocity ensemble. In addition, the high-velocity ensemble produces shock-heated gas temperatures high enough to excite atomic lines so that Lyα, Lyβ, and all atomic metastable lines are at least four orders of magnitude brighter than from the low-velocity ensemble.

We have shown how accurate shock modelling over a wide range of physical conditions allows for the interpretation of molecular and atomic line emission from galaxies. The tool developed in this paper has been applied to observations from the Eyelash galaxy. However, this new method opens up a promising avenue for probing the physical conditions and energetics of galaxies through observations of molecular outflows.

Acknowledgments

We are very grateful to the referee for their thorough reading and their clever comments which triggered important discussions and helped us correct an error in the interpretation of observations from high-z galaxies. The research leading to these results has received fundings from the European Research Council, under the European Community’s Seventh framework Programme, through the Advanced Grant MIST (FP7/2017–2022, No. 787813). We would also like to acknowledge the support from the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.

Appendix A Grid of weakly magnetised shocks

Here we consider shocks with magnetic parameter b = 0.1 and no external radiation G₀ = 0.

In Fig. A.1 we show profiles of the neutral fluid temperatures and velocities, and local sound speed. The initial ion magnetosonic speed is less than 5 km s⁻¹ for these conditions, and so all the models are J-type shocks. For shock velocities V_S ≥ 25 km s⁻¹, these shocks are thinner by an order of magnitude than the equivalent shocks with stronger magnetic field, shown in Fig. 1. This decrease in size is due to the weaker magnetic pressure support in the post-shock. The length scales for the shocks to cool down to 100 K are in range of 10¹² -10¹⁵ cm, or on timescales of 1-10⁵ yrs.

In Figs. A.2, A.3, and A.4 we give the CH⁺ column density, integrated intensity of the 1-0 pure rotational line, and fraction of kinetic flux reprocessed in any rovibrational line, respectively. The column densities show the same general trends as the shocks with stronger magnetic field strength shown in the upper panel of Fig. 3. At high velocities (V_s > 40 km s⁻¹), the column densities are reduced by an order of magnitude compared to the stronger B field shocks, with a maximum column of 10¹³ cm⁻². This is simply due to the shorter length scales. At low velocities (V_s < 30 km s⁻¹) the J-type shocks are much less efficient at producing CH⁺ than the C-type shocks at the same velocities and densities in the b = 1 grid of shocks.

To apply this grid to the Eyelash galaxy described in Sect. 4.3, we show in Fig. A.5 the maximum number of shocks overlapping in velocity and space in the observed beam (upper panel), and the peak optical depth for the CH⁺(1-0) line for each shock in the grid. The cumulative optical depth ${\overline{N}}_S{\tau }_1<1$ for all ensemble models, so Eq. (3) is appropriate over the whole grid. In Fig. A.6 we show the mechanical energy injection rate required to drive the ensembles of shocks. This figure does not significantly differ from Fig. 13 for the b = 1 grid.

For shock tracers other than CH⁺, we show in Fig. A.7 the Lyα photon flux relative to the Mathis ISRF (upper panel) and energy flux relative to the shock kinetic flux (lower panel). In addition, in Fig. A.8 we show a selection of intensities of H₂ rovibrational lines.

Fig. A.1 Shock profiles of ion (red solid) and neutral (blue solid) fluid velocities, local sound speed (black solid), and neutral(green dash-dotted) temperature. The shocks propagate from right to left, and the velocities are given in the reference frame in which the shock is stationary. These shocks were computed with magnetic parameter b = 0.1 and no external radiation field, G0 = 0. We have excluded profiles for shocks with velocities Vs = 40 and 60 km s−1 for brevity.

Fig. A.2 Column density of CH+ for shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0.

Fig. A.3 Line intensities of CH+(1-0) for shocks with magnetic parameter b = 0.1 and no external radiation field G0 = 0.

Fig. A.4 Total CH+ flux as a fraction of shock kinetic flux for shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0.

Fig. A.5 Average number of shocks in the velocity bin ±2.5 km/s around V = 0 overlapping in space according to Eq. (8) (top), and peak optical depth τ1 for each shock for the grid with G0 = 0 and b = 0.1, applied to the Eyelash galaxy observation (bottom).

Fig. A.6 Total shock dissipation rate in the Eyelash galaxy for ensembles of shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0. We note that one solar luminosity L⊙ = 3.8 × 1033 erg s−1.

Fig. A.7 Lyα contribution to the UV flux escaping the shock front for shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0. (Top) Flux relative to the Mathis ISRF, with flux equal to the standard ISRF at the 0 contour. (Bottom) Flux relative to the kinetic flux of the shock $1/2\rho V_s^3$ as a percentage.

Fig. A.8 Intensities of H2 rovibrational lines for shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0.

Appendix B Intensities from single shocks

In Tables B.1-B.15 we tabulate the integrated line intensities for a selection of H₂ and CH⁺ rovibrational lines, CO pure rotational lines, atomic H lines and 2-photon continuum, and atomic fine-structure and metastable lines. The tables are given for shock models in all three grids, at velocities at every 10 km s⁻¹ and integer logn_H. We note that these tabulated lines are a subset of those output by the code that contains

• 200 rovibrational lines of H₂,

• 100 rovibrational lines of CH⁺,

• 14 pure rotational lines of CO,

• Lyα, Lyβ, Hα, and 2-photon continuum (<2400 Å),

• 81 atomic infrared and optical lines including

  • C [230.41 μm], [370.37 μm], [609.75 μm], [4621.57 Å], [4627.35 Å], [8727.13 Å], [9808.32 Å], [9824.13 Å], and [9850.26 Å],

  • C⁺ [157.68 μm], [2323.50 Å], [2324.69 Å], [2325.40 Å], [2326.93 Å], and [2328.12 Å],

  • N [1.04006 μm], [1.04010 μm], [1.04100 μm], [1.04104 μm], [1.147 mm], [6527.23 Å], [6548.05 Å], and [6583.45 Å],

  • N⁺ [76.45 μm], [121.80 μm], [205.34 μm], [3062.83 Å], [3070.55 Å], [6527.23 Å], [6548.05 Å], and [6583.45 Å],

  • O [44.06 μm], [63.19 μm], [145.53 μm], [2958.37 Å], [2972.29 Å], [5577.34 Å], [6300.30 Å], [6363.78 Å], and [6391.73 Å],

  • O⁺ [0.4995 mm],[2470.219 Å], [2470.341 Å], [3726.032 Å], [3728.815 Å], [7318.92 Å], [7319.99 Å], [7329.67 Å], and [7330.73 Å],

  • S [1.08241 μm], [1.13089 μm], [1.15407 μm], [17.43 μm], [25.25 μm], [56.31 μm], [4507.31 Å], [4589.26 Å], and [7725.05 Å],

  • S⁺ [1.0290 μm], [1.0323 μm], [1.0339 μm], [1.0373 μm], [0.31456 mm], [4068.60 Å], [4076.35 Å], [6716.44 Å], and [6730.82 Å],

  • Si [1.5876 μm], [1.6073 μm], [1.6459 μm], [44.812 μm], [68.472 μm], [129.682 μm], [6526.78 Å] and [6589.61 Å],

  • Si⁺ [34.814 μm], [2328.52 Å], [2334.40 Å], [2334.60 Å], [2344.20 Å], and [2350.17 Å].

References

All Tables

All Figures

Fig. 1 Shock profiles of ion (red solid) and neutral (blue solid) fluid velocities, local sound speed (black solid), and neutral fluid temperature (green dash-dotted). The shocks propagate from right to left, and the velocities are given in the reference frame in which the shock is stationary. These shocks were computed with magnetic parameter b = 1 and no external radiation field, G0 = 0. It should be noted that both the distance and temperature axes are different for the C- and J-type shocks, and that the logarithmic scaling can give the false impression that both types of shocks have similar sizes for a given pre-shock density. We have excluded profiles for shocks with velocities Vs = 40 and 60 km s−1 for brevity.

In the text

	Fig. 2 Same as Fig. 1, but for shocks with external radiation field G0 = 103.
In the text

	Fig. 3 Column density of CH+ for shocks with magnetic parameter b = 1 and no external radiation field G0 = 0 (top) or strong external radiation field G0 = 103 (bottom). The dashed lines mark the boundaries between C-, C*-, CJ-, and J-type shocks. The coloured circles show exact values for each model in the grids.
In the text

	Fig. 4 Specific intensities of CH+ pure rotational lines v = 0 − 0 J = n → n − 1 up to n = 10 emitted parallel to the propagation direction of a (top) C-type shock with velocity Vs = 20 km s−1 or (bottom) J-type shock with velocity Vs = 50 km s−1, both with pre-shock density 104 cm−3, magnetic parameter b = 1, and no external radiation field G0 = 0.
In the text

	Fig. 5 Line intensities of CH+ J = n → n − 1 up to n = 10 emitted parallel to the propagation direction of shocks with velocity Vs = 30 km s−1, densities nH = 102–106 cm−3, and radiation parameter G0 = 0 (top) or G0 = 103 (bottom). We note that each figure has a different y-axis scale.
In the text

	Fig. 6 Line intensities of CH+ J = n → n − 1 up to n = 10 emitted parallel to the propagation direction of shocks with velocities Vs = 5–80 km s−1, density nH = 104 cm−3, and radiation parameter G0 = 0 (top) or G0 = 103 (bottom).
In the text

	Fig. 7 CH+(1−0) line intensityfor shocks with magnetic parameter b = 1 and radiation parameter G0 = 0 (top) or G0 = 103 (bottom).
In the text

	Fig. 8 Total CH+ flux as a fraction of shock kinetic flux for shocks with magnetic parameter b = 1 and radiation parameter G0 = 0 (top) or G0 = 103 (bottom).
In the text

	Fig. 9 Cartoon of ensemble of blue- and redshifted shocks within an observational solid angle (black and white circle). The left and right panels represent two scenarios with the same total shocked surface area but a different surface area per shock.
In the text

	Fig. 10 CH+ (1−0) flux density for the Eyelash galaxy from ensembles of 1000 (left) or 10 000 (right) shocks each with velocity Vs = 50 km s−1 and pre-shock density nH = 104 cm−3. The peak of the average ensemble (black) has been fitted to 5.2 × 10−26 erg s−1 cm−2 Hz−1, and the particular ensembles (blue) have a resolution ΔW = 50 km s−1.
In the text

	Fig. 11 Average number of shocks in the velocity bin ±2.5 km s−1 around V = 0 overlapping in space according to Eq. (8) (top), and peak optical depth τ1 for each shock for the grid with G0 = 0 and b = 1, applied to the Eyelash galaxy observation (bottom).
In the text

	Fig. 12 Average number of shocks in the velocity bin ±2.5 km s−1 overlapping in space according to Eq. (8) (top), and peak optical depth τ1 for each shock for the grid with G0 = 103 and b = 1, applied to the Eyelash galaxy observation (bottom).
In the text

	Fig. 13 Total shock dissipation rate in the Eyelash galaxy for ensembles of shocks with no external radiation field G0 = 0 (top) and forshocks with a strong external radiation field G0 = 103 (bottom). We note that one solar luminosity L⊙ = 3.8 × 1033 erg s−1. The violet labels ‘A’, ‘B’, and ‘C’ mark the regions of most efficient models discussed in Sect. 5.2.
In the text

	Fig. 14 Lyα contribution to the UV flux escaping the shock front for shocks with magnetic parameter b = 1 and no external radiation field, G0 = 0. Top: flux relative to the Mathis ISRF, with flux equal to the standard ISRF at the 0 contour. Bottom: flux relative to the kinetic flux of the shock $1/2{\rho }_0{V}_{\text{s}}^3$ as a percentage.
In the text

	Fig. 15 Line intensities of a selection of H2 rovibrational lines, for shocks with no external radiation field (top) and for radiation field G0 = 103 (bottom).
In the text

Fig. 16 Line fluxes of a selection of H2 and CH+ rovibrational lines, CO rotational lines, and atomic fine-structure and metastable lines emitted from ensembles of shocks around the Eyelash galaxy. Ensemble models are grouped by regions in Fig. 13, denoted ‘A’ (blue) for high-velocity J-type shocks with no external irradiation driven by a mechanical energy injection rate ${\mathcal{L}}_K~{10}^{12}$ L⊙, ‘B’ (purple) for low-velocity C*-type shocks with an external UV field G0 = 103 driven by a mechanical energy injection rate ${\mathcal{L}}_K~{10}^{11}$ L⊙, and ‘C’ (green) for low-velocity J-type shocks with an external UV field G0 = 103 driven by a mechanical energy injection rate ${\mathcal{L}}_K~{10}^{11}$ L⊙. Markers indicate fluxes averaged over all the models in the defined region.

In the text

Fig. A.1 Shock profiles of ion (red solid) and neutral (blue solid) fluid velocities, local sound speed (black solid), and neutral(green dash-dotted) temperature. The shocks propagate from right to left, and the velocities are given in the reference frame in which the shock is stationary. These shocks were computed with magnetic parameter b = 0.1 and no external radiation field, G0 = 0. We have excluded profiles for shocks with velocities Vs = 40 and 60 km s−1 for brevity.

In the text

	Fig. A.2 Column density of CH+ for shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0.
In the text

	Fig. A.3 Line intensities of CH+(1-0) for shocks with magnetic parameter b = 0.1 and no external radiation field G0 = 0.
In the text

	Fig. A.4 Total CH+ flux as a fraction of shock kinetic flux for shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0.
In the text

	Fig. A.5 Average number of shocks in the velocity bin ±2.5 km/s around V = 0 overlapping in space according to Eq. (8) (top), and peak optical depth τ1 for each shock for the grid with G0 = 0 and b = 0.1, applied to the Eyelash galaxy observation (bottom).
In the text

	Fig. A.6 Total shock dissipation rate in the Eyelash galaxy for ensembles of shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0. We note that one solar luminosity L⊙ = 3.8 × 1033 erg s−1.
In the text

	Fig. A.7 Lyα contribution to the UV flux escaping the shock front for shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0. (Top) Flux relative to the Mathis ISRF, with flux equal to the standard ISRF at the 0 contour. (Bottom) Flux relative to the kinetic flux of the shock $1/2\rho V_s^3$ as a percentage.
In the text

	Fig. A.8 Intensities of H2 rovibrational lines for shocks with magnetic parameter b = 0.1 and no external radiation field, G0 = 0.
In the text