Self-generated ultraviolet radiation in molecular shock waves I. Effects of Lyman $\alpha$, Lyman $\beta$, and two-photon continuum

Shocks are ubiquitous in the interstellar and intergalactic media, where their chemical and radiative signatures reveal the physical conditions in which they arise. Detailed astrochemical models of shocks at all velocities are necessary to understand the physics of many environments including protostellar outflows, supernova remnants, and galactic outflows. We present an accurate treatment of the self-generated UV radiation in intermediate velocity, stationary, weakly magnetised, J-type, molecular shocks. Shock solutions computed with the Paris-Durham shock code are post-processed using a multi-level accelerated $\Lambda$-iteration radiative transfer algorithm to compute Ly$\alpha$, Ly$\beta$, and 2-photon continuum emission. The subsequent impacts on the ionisation and dissociation of key atomic and molecular species as well as on the heating by the photoelectric effect take the wavelength dependent cross-sections and the fluid velocity profile into account. We analyse shock models with velocities $V=25-60$ km/s, propagating in dense ($n \geq 10^4$ ${\rm cm}^{-3}$), shielded gas. Self-absorption traps Ly$\alpha$ photons in a small region in the shock, though a large fraction escapes into the line wings. We find a critical velocity $V\sim 30$ km/s above which shocks produce a Ly$\alpha$ photon flux exceeding that of the standard ISRF. The escaping photons generate a warm slab (T~100 K) ahead of the shock as well as pre-ionise C and S. These shocks are traced by bright atomic fine structure (e.g. O and S) and metastable (e.g. O and C) lines, substantive molecular emission (e.g. H2, OH, and CO), enhanced column densities of several species (e.g. CH+ and HCO+), as well as a severe destruction of H2O. As much as 13-21% of the initial kinetic energy of the shock escapes in Ly$\alpha$ and Ly$\beta$ photons if the dust opacity in the radiative precursor allows it.


Introduction
Shocks are the fingerprints of the dynamical state and evolution of the interstellar medium. From protostellar outflows and supernova explosions to galactic outflows, massive amounts of mechanical energy are injected at supersonic velocities into the entire range of galactic scales and interstellar phases (e.g. Raymond et al. 2020). Such flows inevitably form shocks which immediately alter the state of the gas via compression and viscous heating (Draine & McKee 1993), but they also give rise to a turbulent cascade and the generation of lower velocity shocks (e.g. Elmegreen & Scalo 2004). By reprocessing the initial kinetic energy into atomic and molecular lines, shocks produce invaluable tracers which carry information on the physical conditions of the shocked gas and on the source itself, including its energy budget, lifetime, and mass ejection rate (e.g. Reach et al. 2005;McDonald et al. 2012;Nisini et al. 2015). The comparison of detailed models with observations, therefore, allows one to ad-dress complex and unsolved issues such as the role of stellar and active galactic nuclei (AGN) feedback on the cycle of matter (Ostriker & Shetty 2011) and galaxy evolution (Schaye et al. 2015;Hopkins et al. 2018;Richings & Faucher-Giguère 2018) as well as the transfer of mass, momentum, and energy from the large scales to the viscous scale in turbulent multiphase media.
New examples of the ubiquity and modus operandi of interstellar shocks have emerged from recent observations of extragalactic environments. In the Stefan's Quintet collision, for instance, a large scale ∼ 1000 km/s shock produces X-ray emission from the T > 10 6 K shock-heated gas. Warm molecular hydrogen, which should be destroyed at these temperatures, dominates the cooling in this region (Guillard et al. 2009), revealing the multiphase nature of the environment along with the necessary transfer of kinetic energy to lower velocity shocks. Interestingly, CO and [CII] lines have linewidths as large as 1000 km/s, bearing the signature of the turbulent cascade driven in the Article number, page 1 of 24 arXiv:2010.01042v1 [astro-ph.GA] 2 Oct 2020 post-shock gas by the large-scale shock (Appleton et al. 2017). Similarly, shocked, warm molecular hydrogen emission has also been observed in a sample of 22 radio galaxies in which star formation is quenched (Lanz et al. 2016). The authors suggest that shocks driven by the radio jets inject turbulence into the interstellar medium, which powers the luminous H 2 line emission. Further opening occurred with the discovery of ubiquitous, powerful molecular outflows in local and high-redshift starburst galaxies (Veilleux et al. 2020;Hodge & da Cunha 2020). Molecular lines in galactic outflows are found with velocity dispersions in excess of 1000 km/s, which can be seen as the signature of powerful turbulence generated by galactic winds (Falgarone et al. 2017). All of these examples reinforce the need for sophisticated and publicly accessible models capable of following the complex chemical and thermal evolution of shocks, down to the viscous scale, in a variety of physical conditions and for different input of mechanical energy.
Detailed shock models have long been developed for a variety of applications in the interstellar medium. Many works have been dedicated to modelling the heating, ionisation and radiative signatures of high velocity shocks (typically V s > 100 km/s) propagating into the atomic environments (with ambient gas densities n < 1 cm −3 ) found in the intercloud medium (Pikel'Ner 1957), supernovae remnants (Cox 1972;Dopita 1976;Shull & McKee 1979) and Herbig-Haro objects (Raymond 1979). Particular attention has been given to accurately treating the radiative precursor, the region ahead of the shock affected by photons generated by the shock itself. Photoionisation in this region can alter the state of the gas before it enters the shock front and thus determines the shock initial conditions. Most recently, state of the art models draw from increasingly large atomic databases to accurately predict the spectra of radiative atomic shocks (Hartigan & Wright 2015;. There have also been extensive studies of shocks focusing on molecular environments. Early work by Field et al. (1968) and Aannestad (1973) considered molecular chemistry in low velocity shocks (V s < 20 km/s) in dense environments (n = 10-100 cm −3 ), with a focus on the survival of molecules through the shock-heated gas. These shocks are not hot enough to produce significant UV photons, and hence no treatment of a radiative precursor is required. In a series of works, Hollenbach & Mc-Kee (1979, 1980 treated the reformation of molecules in the cooling flow behind a fully dissociated shocked gas. By using the UV field calculated by a specialised atomic shock code, they were able to make an extensive study of shocks over both low and high velocities (V s = 25 − 150 km/s) and a broad range of densities (n = 10-10 9 cm −3 ). This work thoroughly detailed the key physical processes-such as photoionisation and dissociation, formation of H 2 on dust, cooling from both atoms and molecules-which determine the structure and radiative characteristics of shocks in these environments. Further advancements were made by Neufeld & Dalgarno (1989) focusing on an accurate treatment of the UV radiative transfer, particularly Lyα and two-photon continuum, and its impact on molecular chemistry in intermediate velocity shocks (V s = 60-100 km/s) in dense media (n = 10 4 -10 6 cm −3 ).
In the present paper, we revisit these pioneering and sophisticated developments by implementing the physics and chemistry of self-irradiated shocks in the Paris-Durham code 1 , a public and versatile state-of-the-art tool initially designed for the study of molecular shocks. Building on the recent updates of Lesaffre et al. (2013) and Godard et al. (2019) who focused on low velocity shocks (V s 25 km s −1 ) irradiated by an external radiation field, we explore here an intermediate velocity range (25 km s −1 V s 60 km s −1 ) where shocks are hot enough to generate significant UV radiation, yet cool enough to prevent the production of multi-ionised species. We present tracers for this velocity range, including atomic line emission, the Lyα and Lyβ counterparts, as well as tracers resulting from the complex molecular chemistry in the cooling flow. Infrared atomic and molecular lines are particularly timely given the wealth of data coming from the Atacama Large Millimeter/submillimeter Array (ALMA) and in anticipation of the James Webb Space Telescope (JWST).
The Paris-Durham shock code is introduced in section 2 where we describe the updates implemented for this work. As a first step to compute the self-generated UV field we only treat the emission from atomic hydrogen. The corresponding threelevel atomic model is outlined in section 3. In section 4 we describe the post-processed radiative transfer treatment to compute the shock generated radiation field. The effects of this radiation field on shocks at different velocities are shown in section 5. We finally conclude by summarising our results in section 6.

Shock model
In order to self-consistently compute a shock and the UV radiation field it generates, we iteratively solve these two components separately. We first generate a shock solution with the Paris-Durham shock code, which we describe in this section, then post-process this output to compute the radiation field. We then run the shock code again in the presence of this radiation field and repeat until convergence.

J-type shocks
Interstellar shocks come in a variety of flavours, depending on the shock speed in relation to the signal speeds in the neutral and ionised components of the fluid (Draine & McKee 1993). We focus solely on single-fluid J-type (discontinuous) shocks in this work, because C-type (continuous) shocks do not reach the high temperatures required to generate significant UV radiation. In addition, for the conditions we consider (Table 2), in particular for weak magnetic fields, C-type shock velocities are bounded to velocities V s < 25 km/s (see figure D1 in Godard et al. (2019)). Figure 1 shows the typical evolution of the temperature and atomic and molecular hydrogen abundances in a J-type molecular shock. Viscous heating mediates an initial adiabatic jump in temperature up to a first plateau, whose temperature can be estimated using the Rankine-Hugoniot jump conditions, assuming an adiabatic index of 5/3, shock velocity much larger than Alfven velocity (v A = B/ 4πρ for magnetic field B and fluid mass density ρ) and fully molecular medium (mean molecular weight of 2.33 a.m.u.), as where V s is the shock speed. At velocities greater than ∼30 km/s interstellar shocks can reach temperatures sufficient to collisionally dissociate H 2 , with dissociation temperature of 5.2×10 4 K, and to collisionally excite atomic hydrogen, with electronic levels starting at ∼1.2×10 5 K. Atomic H is therefore produced in the hottest parts of the shock. Strong cooling due to Lyα excitations causes the temperature to drop from the peak to settle Article number, page 2 of 24 Lehmann, Godard, Pineau des Forêts & Falgarone: Self-generated UV in molecular shocks at T 10 4 K. This plateau is maintained until H 2 reforms on grains deeper in the shock. In the transition between the adiabatic and second temperature plateau the excitation of H is strong enough to produce significant amounts Lyα photons. Flower & Pineau des Forêts (2010) estimated that a 30 km/s shock with preshock density n H = 2×10 5 cm −3 would produce a Lyα photon flux three orders of magnitude more than that of the standard UV interstellar radiation field (ISRF, Mathis et al. (1983)). Such an effect calls for a detailed treatment of the UV radiative transfer in molecular shocks, and hence as a first step in this direction we consider the addition of just one source of UV photons, the hydrogen atom. After summarising key aspects of the Paris-Durham shock code and outlining the updates used in this work we describe the 5-level model of hydrogen in section 3.

Paris-Durham shock code
The Paris-Durham shock code gives the steady-state solution of the plane-parallel magnetohydrodynamics equations coupled with cooling functions and an extensive chemical reaction network appropriate to the molecular phases of the interstellar medium. The version used in this work is that of , with updates described in Lesaffre et al. (2013), Godard et al. (2019) and here.
The code solves a set of coupled, first-order, ordinary differential equations-given by Flower (2010)-using the DVODE forward integration algorithm (Hindmarsh 1983) from some initial conditions. The dynamical variables, that is to say the density and velocity, are often parameters with initial values that we choose to explore, but the initial temperature and chemical abundances are calculated by solving the chemical and thermal equations in a uniform slab for 10 7 years so that the material entering the shock is roughly in chemical and thermal equilibrium. We adopt the elemental abundances of Flower & Pineau des Forêts Table 1. Fractional gas-phase elemental abundances used in this work. We adopt the elemental abundances of  and put all of the species depleted on grain mantles into the gas phase. Numbers in parentheses denote powers of 10.  (2003) and put all of the species depleted on grain mantles into the gas phase (see Table 1). In order to integrate through the discontinuity in J-type shocks, the code employs an artificial viscosity treatment (Richtmyer & Morton 1957). This treatment has been verified to uphold the Rankine-Hugoniot relations, producing a smooth adiabatic jump in which momentum and energy are conserved.
In J-type shocks with molecular initial conditions the dominant source of atomic H in the hottest parts of the shock is the collisional dissociation of H 2 . The treatment of this process in the code has been detailed in Flower et al. (1996) and Le Bourlot et al. (2002), but we summarise it here in order to stress its importance for the accuracy of Lyα production in molecular shocks. The code takes into account the reduced energy required to dissociate H 2 when excited in its rovibrational states. The dominant collision partners are H, H 2 , H + , and e − . Excitations to the dissociative triplet state b 3 Σ + u by collisions with electrons has a rate coefficient given by where E T = 116300 K is the excitation energy of the triplet state and E (v, J) and f H 2 (v, J) are the energy and fractional population of H 2 in the rovibrational state v, J, respectively. For collisions with H, the dissociation rate coefficient is given by where E D = 52000 K (4.48 eV) is the dissociation energy of H 2 . Rate coefficients for collisions with H 2 are unknown for astrophysical conditions, but we use a rate 8 times lower than for H collisions, as indicated by shock tube experiments (Jacobs et al. 1967;Breshears & Bird 1973). The time-dependent populations of the rovibrational levels are solved in parallel with the dynamical and chemical variables as described in Le Bourlot et al. (2002) and , allowing for collisional and radiative transitions between the levels. In Appendix A we check the effect of the number of levels treated has on the shock structure. We found that treating 150 levels is both computationally feasible and accurately models the H 2 cooling and dissociation.

Updates
As we post-process the UV radiative transfer (section 4), the main update to the shock code is that it reads a given UV field at any point during the shock integration. Here we describe the modifications to account for this local field in addition to other updates on previous versions.

Photoreactions
Our treatment of photo-reactions (dissociation and ionisation) is slightly modified from Godard et al. (2019), because we don't have to account for the attenuation of the radiation as this is done in the post-processed radiative transfer. Given the local radiation field energy density-calculated in the post-processed radiative transfer described in section 4-defined as where c is the speed of light and J ν is the angle averaged specific intensity at frequency ν, the photo-reaction rate is given by where h is Planck's constant and σ ν is the frequency dependent cross-section. We use the cross sections of Tabone et al. (2020) for the most important photoreactions, the photo-ionisation of C, S, and Si, and both photo-ionisation and -dissociation of CH, OH, H 2 O, O 2 , and C 2 . For all other photoreactions, we use the form where α is the rate in the unattenuated standard ISRF, and G eff is the effective radiation parameter integrated over the whole frequency range (corresponding to wavelengths between 911-2400 Å), defined as G eff = dν cu ν /hν 1.55 × 10 8 ph/s/cm 2 .
We note that G eff is simply the flux of photons normalised to the Mathis et al. (1983) ISRF with modifications described in Tabone et al. (2020).

Atomic H cooling
The energy lost due to the collisional excitation of the electronic levels of H is given by where the C i j and ∆E i j are, respectively, the collision rates and energy differences between levels i and j, given in the next section, and the summation is taken over all pairs of levels. During the first shock iteration we compute the level populations n i by solving for statistical equilibrium treating only collisions and spontaneous emission. With no absorption or stimulated emission, this is an optically thin treatment. In subsequent iterations of the shock solution, we skip this computation and instead use the populations found during the post-processed radiative transfer described in section 4. In this way the shock solution has a self-consistent treatment of cooling due to H.  Three level atomic hydrogen model and its radiative transitions. Statistical weights are g 1 = 2, g 2a = 6, g 2b = 2, g 3a = 12 and g 3b = 6.

Chemical network
Starting with the reaction network of Godard et al. (2019), we remove all adsorption and desorption reactions as we are interested in shocks fast enough to produce a UV field that removes all mantles via photodesorption. We then extend the network to include collisional ionisation and dissociation rates from Hollenbach & McKee (1989). These include ionisation of He, all metals and select molecules-OH, H 2 O, O 2 , CH, and CO-by collisions with H, H + , He, and electrons. We also include collisional dissociation by electrons for H 2 O, OH, and CO. These updates result in a chemical reaction network including 1256 reactions for 141 species.

Collisional rates
In order to study the UV field generated inside a molecular shock, we consider a three-level atomic hydrogen model including the Lyα (n = 2 → 1), Lyβ (n = 3 → 1), and Hα (n = 3 → 2) line transitions. The second level is divided into 2 sub-levels in order to model the two-photon emission from the 2s metastable state. The third level is also divided into two, grouped by orbitals with radiative transitions to either the 2p or 2s orbitals. This division of sub-levels is schematically shown in Fig. 2 and the values of the transition wavelengths, Einstein A coefficients, and energy differences ∆E for the 6 radiative transitions use atomic data from the NIST database (Kramida et al. 2019) and are listed in Table 2.
Collisions with electrons dominate most of the rates so we consider mostly H-e − collisions. The collisional de-excitation rates coefficients are given by where Υ i j is the effective collision strength. We fit Υ i j data of Anderson et al. (2002) (given in Appendix B) as a temperature dependent power-law so that the de-excitation rate can be written Article number, page 4 of 24 Lehmann, Godard, Pineau des Forêts & Falgarone: Self-generated UV in molecular shocks Collisional rate coefficients for the 2s to 2p orbitals are stronger for H-H + collisions than H-e − collisions, so we also include these rates as given by Osterbrock & Ferland (2006) fit by the same power law equation (10). The collision rate parametres α, β, and energy differences ∆E for the 10 collisional transitions are listed in Table 2.

Non local thermodynamic equilibrium parameter
A first characterisation of a radiative transfer problem is given by the non local thermodynamic equilibrium (LTE) parameter for a line with collisional de-excitation rate C i j = k i j n C where n C is the number density of the collisional partner. characterises the competition between the re-emission of an absorbed photon versus its destruction due to collisional de-excitation, and so approximates a photon destruction probability. If = 1, collisions dominate so that LTE holds and the radiation is thermal. As → 0 the radiative transfer becomes increasingly difficult as the effects of scattering become important. In the fiducial shock we consider in section 5.1 the Lyα emission is generated in gas with temperature T ∼ 10 4 K and electron density n e ∼ 10 3 cm −3 , giving a non-LTE parameter ∼ 10 −14 . This parameter appears in the solution of the statistical equilibrium equations and at such a low value can cause important rounding errors with double precision variables. We thus go to quad precision to avoid this problem. In the next section we turn to the Accelerated Lambda Iteration method to overcome numerical difficulties caused by extreme optical depths.

Post-processed radiative transfer
In this section we describe the post-processing of the Paris-Durham outputs to obtain the UV field at every point inside the shock as well as extending into the unshocked region ahead of the shock. This region ahead of the shock influenced by these UV photons is called the radiative precursor.
To briefly summarise, we first use the Accelerated Lambda Iteration (ALI) algorithm rendered necessary by the extreme optical depths in the Lyα transition. We then use the ALI calculated excited level populations to solve the radiative transfer for the two-photon continuum emission from the 2s metastable state. Finally the radiation field is extended into the radiative precursor. We give a summary of the ALI algorithm here, but for a detailed review see Hubeny (1992) and references therein.

Accelerated Lambda Iteration
We seek to calculate the angle-averaged UV intensity at each position, z, in the shock in order to compute the UV energy density (equation 4) that is used to calculate photoelectric heating, photo-ionisation, and dissociation rates in the Paris-Durham shock code. For a plane-parallel semi-infinite slab, the specific intensity on a ray, I µν , with angle µ = cos θ with respect to the slab normal satisfies the Radiative Transfer Equation where the line emission coefficient for the transition i → j and line opacity where n i and n j are the densities of the upper and lower levels of atomic hydrogen respectively, and B i j and B ji are the Einstein B coefficients. The Gaussian line profile is given by where the thermal velocity and frequency shift from the Doppler-shifted line centre, ν i j , where v z is the flow velocity in the shock propagation direction. The level populations are needed to compute the emission coefficients, and are obtained by solving the Statistical Equilibrium (SE) equationṡ where theJ i j are the mean intensities, J ν , averaged over the line profile for that transition. In order to directly solve this at each position in the shock one would need to invert an enormous matrix with dimensions determined by the discretisation choices for the number of grid positions, frequencies, and angles. As this is computationally unrealistic, the usual strategy is to iterate between solving the RTE and SE equations. Our algorithm is summarised as follows: 1. We initialise the populations by solving the SE equations with theJ i j set to zero. The populations are used to compute the source function where we have assumed complete frequency redistribution, that is that the emission and absorption profiles are equal. 2. The source function is used to solve the RTE along some ray, giving the specific intensity at each point in the shock for every desired frequency. The solution can be written in terms of an integral operator, Λ µν , acting on the source function On a discrete grid this can be written emphasising that the intensity at any one point is coupled to the source function at all points, and that this coupling is encoded in Λ µν . 3. The intensity is used to compute the profile-integrated angleaveraged intensity for each transition We compute the intensity over two rays with angles θ = 0 and 180 • , which has two advantages. Firstly it is coherent with the energetics of the plane-parallel shock in which the energy flux changes only along this direction. Secondly it avoids the non-convergence of the integral caused by emission adding up over infinite length rays near θ = 90 • . 4. The SE equations are solved to update the level populations.
Without modification if we return to the first step and iterate to convergence we have the classical Lambda Iteration method, because the mean intensity in the SE equations is replaced by combining Eqs. (21) and (23) where Λ i j -the profile-integrated angle-averaged Λ µν -is the Λ-operator for this transition. This iteration scheme is known to be pseudo-convergent in systems with extreme optical depths. This effect is shown in Appendix C where we also verify the method against analytic solutions. The ALI algorithm overcomes such pseudo-convergence by splitting the Λ i j operator where Λ * i j is the approximate Λ-operator. The iterative scheme is chosen such that at iteration n Olson et al. (1986) demonstrated mathematically that a nearly optimal choice of Λ * i j is the diagonal of Λ i j . By inspecting Eq. (22), the diagonal encodes the local contribution of the source function to the intensity at a given point. The splitting in Eq. (27) can therefore be interpreted as solving the intensity exactly for the propagation within adjacent grid points while using the long range contribution from the previous iteration. With this choice of Λ * i j = diag Λ ij the radiative terms in the SE equations become Solving the SE equations with this modification, iteratively with Eq. (23) is the ALI algorithm. 5. Before going back to step 1 we further increase the rate of convergence of the iterative scheme by way of the Ng acceleration algorithm (Ng 1974) as formulated by Olson et al. (1986). This algorithm extrapolates the excited populations based on the previous three iterations. We only take this step if the column densities of all levels vary monotonically over those previous iterations.
The approximate Λ-operator, Λ * i j , as the diagonal of Λ i j emerges from the calculation of the intensities. In step 2, to solve the RTE along a given ray we use the method of short characteristics formulated by Paletou & Léger (2007) in which the intensity at position z is computed from the intensity at the upstream grid position, z u , and the source function as where the angle and frequency dependent upstream incremental optical depth is The integral in Eq. (30) is solved by assuming the source function is quadratic in the interval, interpolating on the grid points upstream (z u ) and downstream (z d ) of the point in question Article number, page 6 of 24 Lehmann, Godard, Pineau des Forêts & Falgarone: Self-generated UV in molecular shocks where the coefficients where and The diagonal of the Λ-operator is then given by Ψ o integrated over all angles and frequencies weighted by the line profile. In Appendix D we verify that the algorithm reproduces previous work on multi-level radiative transfer in idealised slabs.

Two-photon emission
The 2s metastable state of hydrogen cannot decay by a single photon process, and instead the 2s-1s transition proceeds by the emission of two photons. During the ALI computation this transition is treated as a resonance line with an Einstein coefficient A 2ph = 8.2 s −1 in order to determine the population of the 2s orbital. This population is then used to compute the continuum emission by solving the RTE including only dust opacity and an emission coefficient where ψ ν is the functional form of the spectrum used in Shull & McKee (1979).

Radiative precursor and postshock extension
After computing the radiation field within the shock with the ALI method, we extend the radiation into the preshock by solving the RTE with the intensity at the shock front as a boundary condition. The ALI algorithm is not necessary in this region because of the low temperatures. In the first iteration, this escaping intensity impinges on a homogeneous slab entering the shock with flow velocity equal to the shock velocity. In subsequent iterations we use values of temperature and chemical composition computed with the Paris-Durham shock code. To decide where to start the Paris-Durham integration the preshock is extended over a size, L pre , such that the dust attenuates the field that escapes the shock, G eff,0 , down to a negligible level compared to the ISRF. To do this we solve where κ D is the dust opacity at the Lyα central wavelength. We similarly extend the radiation further into the postshock by solving the RTE with the ALI output intensity as a boundary condition. This boundary is chosen at a location in the shock where the line radiative transfer has forced the Lyα energy into the wings but before dust attenuation takes effect.
In solving the RTE for the pre and postshock regions we include the opacity due to H, H 2 , and dust. While many H 2 lines of the Lyman (B 1 Σ + u -X 1 Σ + g ) and Werner (C 1 Π u -X 1 Σ + g ) band systems overlap with the Lyα or Lyβ emission, in practice we found only the v = 6 − 0 P(1) line to be of importance. In the first iteration, the H 2 level populations are assumed thermal, and then we use the output of the shock code for self-consistent populations in subsequent iterations.

Results
To understand the impact of the new treatment of the selfgenerated UV field we first consider a single fiducial case of a typical shock. We highlight the impact of this treatment by comparing to a shock model run without self-generated UV. After considering the effects on one shock, we compute shocks at a range of velocities to analyse trends with velocity and give observable predictions for molecular shocks at intermediate velocities.

Fiducial case
As a fiducial case we consider a 40 km/s shock propagating into gas with total hydrogen density n H = 10 4 cm −3 and magnetic field strength 10 µG. In order to emphasise the adiabatic plateau, we adopt here an artificial viscous length of 10 9 cm, 2 orders of magnitude smaller than the real viscous length deduced from the H 2 -H 2 collision cross section. We note, however that the results are weakly dependent on this choice as long as the chosen viscous length is smaller than the real viscous length. In Table 3 we list the shock parameters. The shock acts to reprocess energy through heat into various atomic and molecular lines, dust emission, and to overcome energy barriers in chemical reactions. Some of this energy is also converted into magnetic energy. For the fiducial shock, all the components of this reprocessed energy are shown in Fig. 3. Excitation of molecules and atoms results in emission which escapes the shock and can be used to probe the shock properties. Excitation of H-though it is the largest component reprocessing almost half of the total energy fluxresults in the emission of Lyα, Lyβ, and two-photon continuum UV radiation that is eventually absorbed by dust. Roughly half of this emission escapes ahead of the shock and is absorbed over a large distance set by the lengthscale of dust attenuation. The other half is reprocessed, in the cooling flow of the shock, into thermal energy via the photoelectric effect. Interestingly, the local cooling rate induced by collisional excitations of H and the Article number, page 7 of 24  (Table 3). This shows the energy lost due to excitation of atomic H, other atoms, molecules, and other processes as a percentage of the total energy flux. H 2 chemistry involves cooling due to collisional dissociation and heating due to reformation. Other chemistry is mostly cooling due to collisional ionisation. subsequent interactions of Lyα, Lyβ, and two-photon continuum emission with the surrounding gas are not very different from the results obtained with an optically thin treatment of Lyα and Lyβ. The convoluted radiative transfer presented in the previous section is important only to determine the exact line profile and the spatial asymmetry of the H emission, that is the exact fractions of the radiative energy that travel ahead of the shock and in the postshock. Profiles of temperature, density, and radiation field for this shock are shown in Fig. 4 with (solid) and without (thin dashed) the self-generated UV treatment included. The computation typically converges by 3 iterations of shock and radiative transfer cycles. With no radiation field, the shock propagates into cold gas at ∼10 K. After the adiabatic jump, seen in the middle panels of Fig. 4, dissociation of H 2 mostly due to collisions with electrons produces atomic H with abundance ∼ 1. Cooling due to H excitation by electron collisions-as discussed in section 2.3.2-then determines the transition down to a plateau at T ∼ 10 4 K. The cooling is peaked in this transition as the temperature quickly drops too low for significant excitation of H, at which point O becomes the dominant coolant. This plateau is maintained until H 2 reforms on dust grains and efficiently cools the gas along with other molecules-such as OH and CO-whose production follows the presence of H 2 .
With the self-generated UV treatment included, the photons that escape the shock front form a radiative precursor, heating the gas ahead of the shock to ∼100 K over a distance of ∼10 17 cm (top left panel of Fig. 4), which is the length scale over which dust absorption fully attenuates the field. Photoionisation of C and S produces an increase in electron abundance entering the shock ∼ 3 orders of magnitude more than the initial abundance. An increase in H + also generates more electrons, resulting in stronger dissociation of H 2 in the adiabatic plateau. Hence the cooling due to H excitation and the transition to the second tem-Article number, page 8 of 24 Lehmann, Godard, Pineau des Forêts & Falgarone: Self-generated UV in molecular shocks  (Table 3). For the radiation field, we show the total emission computed with equation 7 (violet) and the contributions to this total by Lyα (red), Lyβ (cyan), and two-photon continuum (blue). The left column shows quantities in the radiative precursor, with distance increasing towards the left (i.e. distance from the shock front), the middle column shows postshock quantities with distance increasing towards the right in log scale while the right column is the same as the middle but in linear scale. Thick lines show results including the UV treatment, while thin dashed lines have no UV included.
perature plateau occurs earlier than without the UV field. The position of this transition converges by the third iteration.
Strong self-absorption traps most of the Lyα photons near a peak of emission, seen in the lower centre panel of Fig. 4, between z = 10 11 -10 14 cm where the UV flux is an order of magnitude more intense (G eff ∼ 450) than that which escapes (G eff ∼ 55). Emerging from the shock front, the Lyα, Lyβ, and two-photon fluxes are roughly 39, 1, and 14 times the ISRF. The flux escaping the trapped region is not symmetric. For Lyα, the rightward escaping photon flux is ∼14% more than leftward. The rightward escaping photons are attenuated by dust in the postshock, generating an extended tail of warm molecular gas due to Article number, page 9 of 24 A&A proofs: manuscript no. 3141593  (Table 3) including (solid) and without (dashed) the self-generated UV. The left column shows quantities in the radiative precursor, with distance increasing towards the left (i.e. distance from the shock front), the middle column shows postshock quantities with distance increasing towards the right in log scale while the right column is the same as the middle but in linear scale.
photoelectric heating. This tail remains above 10 K for a factor of ∼3 longer than the first iteration.  Figure 6 shows a few representative spectra at the shock front (z = 0), in the precursor (z pre = 10 16 cm), at the peak of Lyα (z ∼ 10 12 cm) and Lyβ (z ∼ 2 × 10 11 cm) emission, and deep into the postshock (z = 10 14 cm). The peaks of the wings in the precursor are shifted further away from line centre than the wings in the deep postshock after the peak of emission. This is because the H in the fluid before the peak is at higher temperatures than after the peak, and so it has a wider absorption profile. The precursor spectrum at Lyα shows deep absorption due to the cold H in this region at negative of the shock velocity. The left shows profiles in the radiative precursor, with distance increasing towards the left (i.e. distance from the shock front), the middle shows postshock profiles with distance increasing towards the right in log scale while the right is the same as the middle but in linear scale. Peak temperatures increase with increasing shock velocity (see Eq. 1).
This is also seen in the Lyβ profile, as well as a deep absorption line due to H 2 v = 6 − 0 P(1) Lyman band absorption, whose line centre Doppler-shifted by fluid velocity is shown in the Lyβ panels of Fig. 6 by the dashed grey vertical lines. The spectrum at the peak of Lyα emission shows the typical flat-top profile due to saturation effects of optically thick emission. The same saturation effects make the Lyβ spectrum start to show the flat-top profile. These peaks of emission emerge from the hot regions of the shock at T∼ 50000 K.
Profiles of selected carbon-, oxygen-, sulphur-, and siliconbearing species are shown in Fig. 7. The shock-generated UV strongly changes the preshock abundances and ionisation fraction of C, S, and Si. Without the UV field these 3 species all enter the shock mostly neutral, whereas they all become mostly ionised in the preshock. Oxygen chemistry is also strongly affected, with the photodissociation cross-sections of oxygenbearing molecules O 2 , H 2 O, and OH overlapping Lyα, Lyβ, and/or the two-photon continuum. On the other hand, the photodissociation of CO proceeds via the indirect predissociation mechanism, requiring line absorption at specific wavelengths that do not overlap the H emission (van Dishoeck & Black 1988). Hence as the gas begins to be shocked most of the oxygen is contained in atomic O and CO. H 2 is another molecule that survives the strong UV radiation, however in this case it is because there is not enough time spent in the radiative precursor for strong photodissociation to take place. When the fluid enters the hot shock front the temperatures and densities of the adiabatic jump and transition to the second plateau are not so different with or without the UV field, so the chemistry produces very similar abundances by z ∼10 12 cm. However, the UV persists with fluxes stronger than the typical ISRF, that is G eff > 1, until z ∼10 15 cm strongly affecting the chemistry in the region z 10 13 cm. For example, the atomic ions C + , S + , and Si + all remain more than 2 orders of magnitude higher than when UV photo-ionisation is not taken into account. Oxygen remains mostly atomic whereas without the UV treatment there is a very strong formation of H 2 O in the this region. Finally, the peak of CH + production at z ∼10 14 cm gives an abundance 3-4 orders of magnitude larger than without the UV treatment, before it settles at close to the same level after the UV has been attenuated. We address in the following sections how these differences in local abundances af- fect integrated quantities such as line emission and column densities.

Critical velocity
As a first step in determining the velocity at which the selfgenerated UV becomes important to treat we computed a small grid of shock models with velocities V s = 25-60 km/s at every 5 km/s, propagating into gas with total hydrogen density n H = 10 4 cm −3 . All other parameters are the same as the fiducial shock, listed in Table 3. We leave the exploration of the dependence on density and other parameters to a forthcoming work. The self-consistent converged temperature profiles are shown in Fig. 8. The peak temperature of the adiabatic jump rises with velocity according to Eq. (1). At V s = 25 km/s the peak temperature is not strong enough to dissociate H 2 to remove it as a significant coolant, and so the temperature drops to the postshock equilibrium without any plateaus. As the shock velocity increases, the dissociation occurs faster due higher tem-peratures. Thus atomic H is produced earlier with increasing velocity and so the transition to the second plateau occurs earlier due to the cooling by collisional excitations of H with electrons. Stronger H cooling leads to larger fluid densities in the second plateau. In turn, this leads to larger dust densities and a faster H 2 formation rate, so that H 2 becomes the dominant coolant and the second plateau ends earlier with increasing velocity. The attenuation due to enhanced dust densities balances the stronger UV radiation fields here to give a nearly constant distance at which the shocks cool down to, say, 10 K.
The temperature profiles of the radiative precursors-left panel of Fig. 8-show that heating due to the self-generated UV field starts to increase the preshock temperature above the initial value at a velocity between 30 and 35 km/s. We ran a shock at 32 km/s to sample this transition with more detail. The UV strength parameter G eff emerging from the shock front (z=0) is shown in Fig. 9. This shows a dramatic transition in the shock velocity range V s = 30-32 km/s at which the emergent UV photon flux is stronger than the ISRF, that is G eff > 1. The adiabatic jump at this velocity heats the gas to high enough temperatures to simultaneously remove H 2 as a strong coolant via collisional dissociation and excite H enough to produce strong Lyα fluxes. At V s = 35 km/s the Lyα photon flux alone is ∼ 10 times the ISRF. This presents the rather unique situation of having a slab of fully molecular material up against a boundary radiation field with G eff =10-400 for shocks with velocities between 35-60 km/s. This situation is generally not realised in photondominated regions with equivalent UV flux, because the broad spectrum of the ISRF is more effective at photodissociating H 2 than the shock UV field concentrated around Lyα and Lyβ wavelengths. In such slabs molecule formation thus takes place once the UV field has been strongly attenuated.

Shock length-and time-scales
In Figs. 10 and 11 we show the shock sizes and timescales at which the shock has cooled down to 200 K or 10 K. For the size scale, we also plot the size of radiative precursor that has been heated above 10 K. The sizes of both radiative precursor and postshock region are somewhat constant for shock velocities V s ≥ 35, varying by less than a factor of 2 depending on the size criteria. This is because these sizes are essentially determined by the dust attenuation of the UV field.

Line emission
The shock observables that we consider here, that is selected atomic emission lines, rovibrational lines of H 2 , and column densities of various species, are all quantities obtained by integrating local quantities through the shock. The choice of where to end that integration therefore determines these quantities and the interpration of observations. For this work we integrate until the gas falls down to 10 K, so that we integrate through shock-heated gas. There is also material ahead of the shock heated above 10 K in the radiative precursor. However, the precursor slabs computed here assume a uniform medium over large distances, which is not necessarily relevant in real astrophysical systems. Because it is unknown how realistic it is to include emission generated in the precursor, so we simply consider two cases: including this material or not.
In Fig. 12 we show the intensity of selected infrared fine structure lines generated in the shock heated gas. Most of these lines show a dramatic increase in intensity over the velocity  range 30-60 km/s compared to lower velocity shocks. For example, the C + (158µm) and Si + (34.8µm) lines shows increases of more than 2 orders of magnitude once the UV becomes important. The C + line also shows an increasing trend with velocity in this range, and has a strong contribution from the gas heated in the radiative precursor. The other fine structure lines show a remarkably constant intensity in this range. This figure  Intensities of selected metastable atomic lines are shown in Fig. 13. More optical lines are accounted for in the model, but we show just the strongest lines for each species. As with the fine structure lines there are dramatic increases in these intensities over lower velocity shocks. In addition, there are stronger trends with velocity, for example the N + (6583 Å) line varies by more than 3 orders of magnitude over the velocity range. A general result is that the intermediate velocity molecular shocks produce strong emission in the metastable lines O(6300 Å), C(9850 Å), and S + (6731 Å).
In Fig. 14 we plot the intensities of the pure rotational lines S(0) up to S(4) of H 2 as well as the v = 1 − 0 S(1) generated Article number, page 13 of 24 A&A proofs: manuscript no. 3141593  in the shock-heated gas. Unlike the atomic lines, these lines do not show a significant increase compared their intensity from the 25 km/s shock. They are also remarkably constant over the velocity range, except for the contribution to the S(0) by the radiative precursor. Hence combined observations of atomic lines and H 2 lines would be necessary to probe the different shock velocities in systems with an ensemble of shocks, such as in the turbulent cascade in the wakes galactic outflows or supernova shocks. These H 2 lines lie in the observational bands of JWST and so could be used to interpret planned observations of such astrophysical systems.
The UV emission from atomic H-Lyα, Lyβ or two-photon continuum-are also possible observables. In Table 4 we list their fluxes that escape into the radiative precursor as well as the shock kinetic flux 0.5ρV 3 s . For shocks with velocities V s ≥ 40 km/s, 13-21% of the kinetic energy entering the shock comes back ahead of the shock in Lyα, Lyβ or two-photon emission. This UV emission can be absorbed by dust in the material ahead of the shock, which is often an unknown quantity in astrophysical systems. Combined with the previously discussed shock tracers, Table 4 could then be used to give a prediction of the maximum contribution to the Lyα emission from intermediate velocity shocks.

Column densities
In addition to atomic and H 2 line emission there are shock tracers in the molecular chemistry. Column densities of selected species are shown in Fig. 15. We show the column densities for postshock regions integrated until the gas cools to 10 K (solid lines) and also with the material in the radiative precursor added to shocked gas (dotted). Due to the extended warm tail produced by the UV there is an increase in the total H column density-by a factor ∼3-4 in the postshock or an order of magnitude when including the preshock-for shocks with velocities V s > 30 km/s. Above this velocity, however, there is no significant correlation with velocity. This is due to the roughly constant shock-size, discussed in the previous section, over the velocity range once Lyα production becomes significant.
Some species show enhanced column densities when the UV field becomes strong at velocities V s ≥ 32 km/s. For example, column densities of CH + and HCO + are increased by more than 2 orders of magnitude by the end of the velocity range compared to their values at V s = 30 km/s. This is due to the larger abundances of C + maintained by photoionisation deep in the postshock. The updated UV treatment may therefore allow us to consider shock models to interpret extragalactic ALMA observations of CH + emission (Falgarone et al. 2017) as well as protostellar jet observations of HCO + (e.g. Tafalla et al. 2010). In the other direction, photodissociation of H 2 O strongly reduces its column densities compared to molecular shocks at lower velocities. Observations involving oxygen chemistry have long stimulated the development of shock models. For instance, recent Herschel observations of protostellar environments reveal H 2 O abundances far too low to be explained with simple chemical models (e.g. Goicoechea et al. 2012;Kristensen et al. 2013;Karska et al. 2014Karska et al. , 2018. The precise origins of the emission of O 2 (e.g. Yıldız et al. 2013;Chen et al. 2014;Melnick & Kaufman 2015) and O (Kristensen et al. 2017) have met with similar struggles. Mixtures of high velocity dissociative shocks, photodissociation regions, low velocity C-type shocks, and shocks irradiated by UV from either external sources or the shocks themselves have been invoked to explain these observations. The present work broadens the range of parameters of the Paris-Durham shock code to be used to study these problems, and shows how self-irradiated shocks could be a viable solution to explain low abundances of H 2 O and O 2 . The column densities of all these species (except CH + ) in the heated precursor are larger than in the postshock because the precursors are orders of magnitude larger than the postshock, as seen in Fig. 10. H 2 O, OH, CO, and total H have more than a factor of 3 larger column in the radiative precursor than in the shocked gas. In summary, molecular shocks at intermediate velocities (V s =32-60 km/s) have large column densities of CH + , HCO + , and CO while also having low column densities of H 2 O.
In Appendix F we compare the column densities of selected species in V s = 60 km/s shocks to the results of the shock models of Neufeld & Dalgarno (1989). The good agreement between the two works is a striking result given the decades of updates to chemical reaction networks and rates, computational methods, and the inclusion of the magnetic field in our work.

Conclusion
We have implemented a treatment of the UV radiative transfer including the atomic hydrogen lines Lyα, Lyβ, and the twophoton continuum in order to solve for self-irradiated molecular shocks at intermediate velocities using the Paris-Durham public shock code. The main results are summarised as follows: -A detailed treatment of radiative transfer is necessary to accurately compute the line profiles and escape of Lyα and Lyβ. However, to understand the energetic impacts of H excitation in these shocks it is sufficient to model the radiation with an optically thin treatment. -For preshock density n H = 10 4 cm −3 , shocks with velocity V s > 30 km/s produce a UV radiation field that escapes into the preshock gas with a Lyα photon flux stronger than the standard ISRF. -For shock velocities between 35-60 km/s the escaping UV photons give a radiation field parameter G eff ∼ 10-400. In other words, the line and continuum emission of H alone carries as much as 400 times more photons than those contained in the entire ISRF. These photons escaping ahead of the shock heat the preshock gas in a radiative precursor over ∼ 10 17 cm to ∼ 100 K. -The maximum contribution to Lyα emission for shocks with velocities 40-60 km/s is as large as 13-21% of the shock kinetic energy flux. -After HI emission, molecular emission is the second most important coolant in intermediate velocity shocks. In addition, these shocks have large column densities of CH + , HCO + , and CO while also having low column densities of H 2 O. These shocks may be useful in explaining low H 2 O abundances found in Herschel observations of protostellar environments.
-Compared to lower velocity shocks, atomic fine structure and metastable emission lines are boosted by many orders of magnitude. For example, the O(63.2 µm) and S(25.2 µm) fine-structure lines and O(6300 Å) and C(9850 Å) metastable lines are particularly bright. Most fine structure lines are remarkably constant at intermediate shock velocities, but ratios of metastable lines could be used as a probe of shock velocity. -We predict intensities of ro-vibrational lines of H 2 . Many of these lines will be observable by JWST and hence this treatment will be critical for the interpretation of observations of phenomena with shocks at these velocities.
The present work shows that it is indeed important to properly model the UV radiation generated by the shock heated gas in order to extend models of molecular shocks to higher velocities. A forthcoming work will present the application of the present models for the interpretation of observations.

Appendix A: Number of H 2 levels
To estimate the impact of the number of H 2 levels on the shock structure we run 4 shock models with shock velocity V s = 40 km/s and preshock density n H = 10 4 cm −3 , varying the number of H 2 levels treated: N lev =50, 100, 150 or 200. Otherwise, we use the same input parameters as the fiducial model (see Table 3). In Fig. A.1 we show the temperature and hydrogen abundance profiles for these shocks, which are reasonably converged when 150 levels are included. When treating 50 levels the dissociation of H 2 is underpredicted by ∼ 4 orders of magnitude, leading to the reformation of H 2 to occur earlier and reduce the size of the shock by a factor of ∼ 5. Using 150 levels of H 2 , Fig. A.2 shows abundances profiles of H 2 and H over the grid of shock models described in section 5.2. At V s =25 km/s, H 2 is strongly dissociatied in the adiabatic jump and H is the dominant species, but there is a transition at V s =30 km/s above which shocks dissociate H 2 to a similar abundance x(H 2 ) ∼10 −7 .

Appendix B: Atomic hydrogen parameters
The three level model of hydrogen described in section 3 is built by combining the data of the sub-levels with energies and statistical weights taken from the NIST database (Kramida et al. 2019) and listed in Table B  The energies of the combined levels are the average energies of their sub-levels weighted by the statistical weights Similarly the Einstein coefficients of the combined levels, listed in Table 2, are the weighted average Einstein coefficients where the sum is over all the transitions from the sub-levels of the upper combined level to any of the sub-levels of the lower combined level. The sub-level transition data are taken from the NIST database (Kramida et al. 2019)-with the exception of the two-photon Einstein coefficient taken from Nussbaumer & Schmutz (1984)-and listed in Table B.2. The effective collision strengths for de-excitation due to collisions with electrons are taken from Anderson et al. (2002), listed in Table B.3. The effective collision strengths, Υ A , have been modified in order to be applied to the fine structure of H. For all transitions n l j → nl ĵ Υ =Υ A n l → nl (2 j + 1) (2s + 1)(2l + 1) (2 j + 1) (2s + 1)(2l + 1) .
Hence the sum of the collision strengths for sublevel transitions gives the collision strength for the combined level transitions. We fit this combined collision strength with a power-law which allows us to express the rate coefficients as power-laws, Eq. 10, with the fit coefficients listed in Table 2.

Appendix C: Two-stream approximation test
Our ALI transfer code is tested with the two-stream approximation test where the RTE of a homogeneous slab is solved along two rays with µ = ±1/ √ 3. In addition, we impose a boundary condition that deep in the slab (when τ → ∞) the source function equals the Planck function. This setup has analytic solution and is the non-LTE parameter For the homogeneous slab we use the Lyα radiative parameters listed in Table 2, and densities and temperatures typical for shock conditions: T = 10 4 K, n e = 10 3 cm −3 . This gives a non-LTE parameter ∼ 10 −14 .
The result without using the ALI modification (that is, with Λ * = 0) is shown in Fig C.1. After 50 iterations the solution is still more than 6 orders of magnitude below the analytic solution. Using Λ * as defined by Eq. (34)

Appendix D: Three-level hydrogen test
To further test the implementation of our ALI transfer code, we compare with the three-level benchmark problem of Avrett (1968). In this problem the radiative transfer is solved for a plane-parallel semi-infinite atmosphere with constant collisional de-excitation rates C i j = 10 5 s −1 for all three transitions and temperature 5000 K. The statistical weights are g 1 = 2, g 2 = 8, and g 3 = 18. The transition frequencies are ν 21 = 2.47 × 10 15 Hz and ν 31 = 2.93 × 10 15 Hz. The Einstein coefficients are A 21 = 4.68×10 8 s −1 , A 31 = 5.54×10 7 s −1 , and A 32 = 4.39×10 7 s −1 . The line profiles are Gaussian (eq. 16), we use an eight-ray Gaussian quadrature scheme and the boundary condition lim τ→∞ S i j = B i j .

Appendix E: Tables of results
Here we tabulate results for shock models in the velocity range V s = 25-60 km/s described in section 5.2. In table E.1 we show the intensities escaping the shock front, which includes material  (1968). Dotted lines are solutions tabulated in Avrett (1968). in the shock down to 10 K as discussed in section 5.2.3, of hydrogen and atomic fine structure and metastable lines. We also tabulate a subset of H 2 rovibrational lines given as output by the code. Additional lines are available upon request. A selection of these lines are shown in figures 12-14. In table E.2 we give column densities of a variety of species. Some of these column densities are shown in figure 15.

Appendix F: Comparison with Neufeld & Dalgarno (1989)
In   Note. Numbers in parentheses denote powers of ten.

Velocities (km/s) Species
Note. Numbers in parentheses denote powers of ten. Column densities (cm −2 ) of warm gas, T > 200 K, predicted with the Paris-Durham shock code for shocks with velocity V s = 60 km/s and ambient densities n H = 10 4 , 10 5 , and 10 6 cm −3 (L20), and comparison with the column densities predicted by Neufeld & Dalgarno (1989) (ND89). Total column densities for ND89 are estimated from their figures 14 and 15.
Article number, page 23 of 24  . This shows the energy lost due to excitation of atomic H, other atoms, molecules, and other processes as a percentage of the total energy flux. H 2 chemistry involves cooling due to collisional dissociation and heating due to reformation. Other chemistry is mostly cooling due to collisional ionisation.
Article number, page 24 of 24