TY - JOUR
AU - Metzger, Brian, D
AB - ABSTRACT Radiative shocks, behind which gas cools faster than the dynamical time, play a key role in many astrophysical transients, including classical novae and young supernovae interacting with circumstellar material. The dense layer behind high Mach number |$\mathcal {M} \gg 1$| radiative shocks is susceptible to thin-shell instabilities, creating a ‘corrugated’ shock interface. We present 2D and 3D hydrodynamical simulations of optically thin radiative shocks to study their thermal radiation and acceleration of non-thermal relativistic ions. We employ a moving-mesh code and a specialized numerical technique to eliminate artificial heat conduction across grid cells. The fraction of the shock’s luminosity Ltot radiated at X-ray temperatures |$kT_{\rm sh} \approx (3/16)\mu m_\text{p} v_{\rm sh}^{2}$| expected from a 1D analysis is suppressed by a factor |$L(\gt T_{\rm sh}/3)/L_{\rm tot} \approx 4.5/\mathcal {M}^{4/3}$| for |$\mathcal {M} \approx 4\text{--}36$|⁠. This suppression results in part from weak shocks driven into underpressured cold filaments by hot shocked gas, which sap thermal energy from the latter faster than it is radiated. Combining particle-in-cell simulation results for diffusive shock acceleration with the inclination angle distribution across the shock (relative to an upstream magnetic field in the shock plane – the expected geometry for transient outflows), we predict the efficiency and energy spectrum of ion acceleration. Though negligible acceleration is predicted for adiabatic shocks, the corrugated shock front enables local regions to satisfy the quasi-parallel magnetic field geometry required for efficient acceleration, resulting in an average acceleration efficiency of εnth ∼ 0.005–0.02 for |$\mathcal {M} \approx 12\text{--}36$|⁠, in agreement with modelling of the gamma-ray nova ASASSN-16ma. stars: novae, stars: supernovae: general, Shock waves, radiation: dynamics, X-rays: bursts, instabilities 1 INTRODUCTION Non-relativistic shock waves power thermal and non-thermal radiation from a wide range of astrophysical transient events. The latter include massive star eruptions or young supernovae (SNe) in which the stellar ejecta collides with a dense pre-explosion outflow from the progenitor star (e.g. SNe IIn, Smith & McCray 2007; Chevalier & Irwin 2011; Smith 2016) or a relic protostellar disc (Metzger 2010); the merger of binary stars, or ‘common envelope’ events (e.g. Ivanova et al. 2013; MacLeod, Ostriker & Stone 2018), following which the dynamically ejected envelope collides with mass loss from earlier stages of the merger process (Metzger & Pejcha 2017; Pejcha et al. 2017); and internal collisions within the outflows from hot white dwarfs during classical nova eruptions (Mukai & Ishida 2001; Sokoloski et al. 2006; Martin & Dubus 2013; Metzger et al. 2014; Martin et al. 2017). In classical novae, relativistic particles accelerated by such shocks power GeV gamma-ray emission, as discovered by Fermi-LAT in temporal coincidence with the nova optical outburst (Ackermann et al. 2014; Cheung et al. 2016; Franckowiak et al. 2018) . A common feature of all of these systems is that the density of the shocked gas can be sufficiently high for the radiative cooling time to be much shorter than the radial expansion time; such shocks, in which radiative losses play a crucial role in their stability and dynamics, are referred to as ‘radiative’ (e.g. Raymond 1979; Drake 2005). As their name implies, the kinetic luminosity dissipated at radiative shocks is radiated with nearly unity efficiency, making them extremely effective at powering luminous transient emission. This is particularly true when the shocks occur in an expanding outflow that has already lost most of its initial internal energy to adiabatic expansion, a condition which is closely related to the shocks being strong (Mach number |$\mathcal {M} \gg 1$|⁠). The shock-powered thermal energy typically emerges from this environment at visual wavelengths, because the column of gas ahead of the shocks effectively absorbs ultraviolet (UV) and soft X-ray photons emitted by the hot post-shock medium and re-radiates this energy at lower frequencies, where the opacity is much lower. Nevertheless, thermal X-rays can sometimes be directly observed at sufficiently high photon energies ≳ 1–10 keV, above the cut-off set by photoelectric absorption (e.g. Nymark, Fransson & Kozma 2006; Metzger et al.2014). Though emphasized less frequently, the high gas densities near radiative shocks can also result in efficient non-thermal emission, by providing a dense field of target photons and ambient ions for relativistic particles accelerated near the shock front to interact with via inverse Compton and hadronic (pion creation) processes, respectively (e.g. Metzger et al. 2015; Martin et al. 2017). Rapid cooling of the thermal plasma also boosts the energy available to non-thermal particles in the post-shock cooling regions through adiabatic compression (Vurm & Metzger 2018). The dense material surrounding radiative shocks can therefore provide a ‘calorimeter’ for measuring both the total shock power (from the thermal visual or X-ray wavelength radiation) and the energy in accelerated relativistic particles (from non-thermal X-rays or gamma-rays). In this way, simultaneous broad-band monitoring of transients thus enables one to constrain or measure the fraction of the shock power which is placed into a power-law distribution of non-thermal particles, i.e. on the acceleration efficiency εnt (Metzger et al. 2015). The classical nova ASASSN-16ma displayed temporally correlated optical and gamma-ray light-curve behaviour, a clear signature that the optical radiation was the product of reprocessed radiative shock emission rather than being directly powered by the white dwarf (Li et al. 2017). Several considerations favour a hadronic origin for the gamma-ray emission from classical novae (Metzger et al. 2015, 2016; Martin et al. 2017). The observed ratio of the optical and gamma-ray flux in ASASSN-16ma implied an ion acceleration efficiency of εnt ≈ 0.005 for this event (Li et al. 2017). Depending on ones point of view, the acceleration efficiency of ∼1 per cent measured for ASASSN-16ma (Li et al. 2017) and other classical novae (Metzger et al. 2015) is either surprisingly high or surprisingly low. Particle-in-cell (PIC) plasma simulations of high Mach number non-relativistic shocks find values of εnt ≈ 0.05–0.2 when the magnetic field upstream of the shock is oriented quasi-parallel to the shock normal. However, no measurable acceleration is seen (εnt ≈ 0) when the magnetic field is in the plane of the shock (e.g. Caprioli & Spitkovsky 2014a). The latter case is the expected one in the radially expanding outflows from a transient explosion, because the magnetic field being advected outwards from the central ejecta or pre-explosion stellar wind will be dominated by its toroidal component on large scales due to flux conservation (e.g. as in the Parker spiral characterizing the solar wind) and thus would be oriented perpendicular to the velocity of the radially expanding shock front. Most models of radiative shocks and their emission are based on 1D steady-state hydrodynamical solutions for the flow properties (e.g. the MAPPINGS catalogue; Dopita & Sutherland 1996; Allen et al. 2008) . However, radiative shocks are prone to various instabilities, which imprint a rich dynamical and geometrical structure that could play an important role in shaping both their thermal and non-thermal emission. Chevalier & Imamura (1982) performed a linear stability analysis of radiative shocks with power-law cooling function (Λ∝Tα), finding them to be subject to an oscillatory instability for cooling exponents α ≲ 0.4. Gaetz, Edgar & Chevalier (1988) explored thermal instability with realistic interstellar medium cooling. Innes (1992) developed 1D models of the thermal instability, finding secondary shocks driven into underpressured, cool gas and finding that magnetic pressure cushions the instability. The ‘pressure-driven thin-shell overstability’ (Vishniac 1983; Bertschinger 1986; Mac Low & Norman 1993) occurs when the dense shell of gas behind a radiative shock, which is bounded on its other side by a high pressure, is decelerating, causing growing overstable non-radial oscillations. Of particular relevance to this work, the dense cool layer of gas behind a radiative shock is also subject to a ‘non-linear thin-shell instability’ (NTSI; Vishniac 1994). This occurs because lateral perturbations of the interaction front cause material to be diverted from convex to concave regions, in such a way that the direct compression from the oppositely directed flows creates multiple elongated regions, imparting the shock front with a corrugated geometrical structure. Numerical simulations of radiative shocks and the NTSI find a shock interaction front characterized by highly complex regions of dense, cooled gas (Stevens, Blondin & Pollock 1992; Strickland & Blondin 1995; Walder & Folini 2000; Pittard 2009; Parkin et al. 2011; McLeod & Whitworth 2013; Kee, Owocki & ud-Doula 2014). One consequence of the ‘corrugated’ shock front created by the NTSI is a reduction in the thermal X-ray luminosity as compared to the naive expectation from a 1D laminar compression analysis (e.g. Kee et al. 2014). While a head-on strong shock of velocity |$v$|sh heats the colliding gas to a temperature |$kT_{\rm sh} \approx (3/16)\mu m_\text{p} v_{\rm sh}^{2}$|⁠, where μ is the mean molecular weight, this is reduced by a factor ∝cos 2α for shocks of obliquity α > 0. Furthermore, as we will show here, the hot post-shock gas undergoes additional losses by performing PdV work on (driving weak shocks into) underpressurized cold regions, transferring energy to the cool gas before it can be directly radiated by the hot gas and thus further reducing the average temperature of the post-shock emission. Accurate predictions for the X-ray efficiency of radiative shocks are of crucial importance, as they affect what inferences are made about the power of the shocks from X-ray observations. Only a few works have attempted to quantify the X-ray suppression of radiative shocks and its dependence on the upstream parameters. Kee et al. (2014) performed 2D simulations of colliding flows and dual radiative shocks, with Mach numbers |$\mathcal {M} \approx 42$|⁠, finding suppression of the X-ray luminosity by a factor of ≈1/50 relative to the expectation for an otherwise identical 1D flow. Kee et al. (2014) make the ansatz that the reduction in the X-ray luminosity scales with the Mach number as |$\propto \mathcal {M}^{-1}$|⁠, but to our knowledge they provide no motivation for this specific dependency. In addition, few of the previous numerical studies demonstrate convergence of their results with grid resolution. Parkin & Pittard (2010) emphasize that ‘numerical conduction’, and other effects associated with limited grid resolution, can artificially lower the temperature of the shock-heated regions and thus reduce their X-ray emission relative to the physical case. A fundamental numerical challenge is to resolve the small cooling length behind the shock, which for high Mach number |$\mathcal {M} \gtrsim 10\text{--}30$| shocks of interest can be orders of magnitude smaller than the system size. A related challenge is to avoid the artificial exchange of thermal energy between adjacent hot and cold grid cells with densities that differ by factors up to |${\sim } \mathcal {M}^{2} \gtrsim 10^{2}\text{--}10^{3}$|⁠. Beyond its effect on the temperature of the thermal radiation, the irregular geometry of the radiative shock front could have an equally pronounced effect on their ability to accelerate relativistic ions and thus to power non-thermal X-ray and gamma-ray emission. The lower Mach numbers of the oblique shock fronts created by the NTSI could in principle act to steepen the power-law index of the energy distribution of accelerated relativistic particles, dN/dE∝E−β, relative to the expectation of β ≃ 2 in the |$\mathcal {M} \gg 1$| limit (e.g. Drury 1983). Furthermore, as discussed above, the local inclination angle of the shock normal relative to the direction of the upstream magnetic field (itself perpendicular to the inflow velocity) plays a crucial role in whether ion acceleration can occur at all (e.g. Caprioli & Spitkovsky 2014a). Here, we present high-resolution numerical simulations of radiative shocks using the moving-mesh hydrodynamical code rich (Yalinewich, Steinberg & Sari 2015) . We introduce for the first time a specialized ‘Lagrangian’ numerical technique, which eliminates artificial diffusion of thermal energy across cell boundaries, allowing us to obtain greater convergence on this problem than possible in previous work. Our main goal is to quantify the effects of thermal instabilities and the NTSI in radiative shocks on the temperature distribution of the thermal X-ray emission, and on enhancement of non-thermal particle acceleration due to local shock obliquity, as a function of the Mach number. This paper is organized as follows. In Section 2, we overview the physical properties of radiative shocks and define the domain of applicability of our simulations in the parameter space of the properties of the shock and upstream medium. In Section 3, we outline our numerical code and methods, including the ‘Lagrangian’ scheme to reduce artificial diffusivity, for which we provide numerical verification tests in Appendix A. In Section 4, we describe our results and their application to shock-powered transients, such as ASASSN-16ma. In Section 5, we summarize our conclusions and their implications. 2 OVERVIEW OF RADIATIVE SHOCKS Consider two oppositely directed locally planar flows of equal density ρ which collide head-on, creating dual shocks of velocity |$v$|sh. Initially (before cooling instabilities grow and distort the shock front), the shocks are effectively 1D and heat the gas to a temperature given by the usual jump condition \begin{eqnarray*} T_{\rm sh} \approx \frac{3}{16}\frac{\bar{\mu }m_\text{p}}{k_{\text{b}}} v_{\rm sh}^{2} \approx 1.4\times 10^{7}\, {\rm K}\left(\frac{v_{\rm sh}}{10^{3}\, {\rm km\, s^{-1}}}\right)^{2}, \end{eqnarray*} (1) where |$\bar{\mu } = 0.62$| for solar composition and |$\bar{\mu } = 0.74$| for that typical of classical novae (Vlasov, Vurm & Metzger 2016). We assume that gas cools behind the shock at the optically thin rate, \begin{eqnarray*} \dot{q} = \Lambda (T)n_\text{e} n_\text{p}, \end{eqnarray*} (2) where ne and np are the number densities of free electrons and protons, respectively, and Λ(T) is the cooling function (Fig. 1). For a strong shock (⁠|$\mathcal {M} \gg 1$|⁠) and adiabatic index γ = 5/3, the gas immediately behind the shock is compressed to a density ρsh = 4ρ and cools radiatively on the time-scale \begin{eqnarray*} t_{\rm cool} = \frac{\bar{\mu }}{\mu _\text{p}\mu _\text{e}}\frac{(3/2)k_{\text{b}}T_{\rm sh}}{n_{\rm sh}\Lambda }, \end{eqnarray*} (3) where nsh = ρsh/mp, μe = 2mp/(1 + X) ≃ 1.16, and μp = 1/X ≃ 1.39 for hydrogen mass fraction X = 0.72. The immediate post-shock gas, slowed to a velocity |$v$|sh/4, cools over a characteristic length-scale given by \begin{eqnarray*} \mathcal {L}_{\rm cool} = v_{\rm sh}t_{\rm cool}/4. \end{eqnarray*} (4) If the interacting flows possess a vertical or radial extent ∼R, then define a ‘cooling efficiency’ according to \begin{eqnarray*} \eta & \equiv & \mathcal {L}_{\rm cool}/R \nonumber \\ &\simeq & 10^{-3}\left(\frac{n_{\rm sh}}{10^{10}{\rm cm^{-3}}}\right)^{-1}\left(\frac{v_{\rm sh}}{10^{3}\, {\rm km\, s^{-1}}}\right)^{4}\left(\frac{t_{\rm exp}}{1\, {\rm wk}}\right)^{-1}, \end{eqnarray*} (5) where we have approximated Λ ≈ 1.3 × 10−22(T/107K)−1 erg cm3 s−1 for nova metallicity gas in the temperature range |$10^{5}\, {\rm K} \lesssim T \lesssim 10^{7}$| K. We normalize the system size to a value R = |$v$|shtexp, for a time-scale texp = 1 week, and particle density nsh = 1010 cm−3, values characteristic of Type IIn SNe and classical novae. Figure 1. Open in new tabDownload slide Cooling function used in our simulations for solar composition (solid red line) and a composition appropriate to the ejecta of a classical novae (solid blue line), as calculated assuming collisional ionization equilibrium, using cloudy version C17.00 (Ferland et al. 2017). For nova composition, we follow Vlasov et al. (2016) and adopt solar composition with enhanced abundances as follows: [He/H] = 0.08, [N/H] = 1.7, [O/H] = 1.3, [Ne/H] = 1.9, [Mg/H] = 0.7, and [Fe/H] = 0.7. Figure 1. Open in new tabDownload slide Cooling function used in our simulations for solar composition (solid red line) and a composition appropriate to the ejecta of a classical novae (solid blue line), as calculated assuming collisional ionization equilibrium, using cloudy version C17.00 (Ferland et al. 2017). For nova composition, we follow Vlasov et al. (2016) and adopt solar composition with enhanced abundances as follows: [He/H] = 0.08, [N/H] = 1.7, [O/H] = 1.3, [Ne/H] = 1.9, [Mg/H] = 0.7, and [Fe/H] = 0.7. A shock is defined as ‘radiative’ when η ≪ 1 (tcool ≪ texp), since in this limit, the shocked gas has sufficient time to radiate its thermal energy before adiabatic radial or lateral expansion reconverts the thermal energy back into bulk motion. In this limit, most of the upstream kinetic power dissipated by the shock, \begin{eqnarray*} L_{\rm tot} \sim \bar{\mu } m_\text{p} n_{\rm sh} v_{\rm sh}^{3}R^{2}, \end{eqnarray*} (6) ultimately emerges as radiation at some waveband, where the pre-factor in this expression depends on the precise geometry of the system. In the opposite, adiabatic limit η ≫ 1 (tcool ≫ texp), the radiative luminosity is reduced from Ltot by a factor of 1/(1 + 5η/2) for cases when free–free emission dominates the cooling (e.g. Metzger et al. 2014). If gas behind the shock cools and compresses at roughly constant pressure (as turns out to be a crude but useful approximation), then its density increases by a factor of \begin{eqnarray*} \frac{\rho _{\rm c}}{\rho } \sim \mathcal {M}^{2} \approx \frac{T_{\rm sh}}{T_{\rm min}} \sim 10^{3}\left(\frac{T_{\rm min}}{10^{4}\, {\rm K}}\right)\left(\frac{v_{\rm sh}}{10^{3}\, {\rm km\, s^{-1}}}\right)^{2}, \end{eqnarray*} (7) where |$\mathcal {M} = v_{\rm sh}/c_{\rm s}$| is the Mach number of the shock defined relative to the sound speed cs of the cooled gas. The latter is generally regulated to a value ≈(kbTmin/μmp)1/2 ≈ 10–20 km s−1, corresponding to the temperature floor Tmin ∼ few × 104 K set by photoionization balance and the minimum of the cooling curve (Fig. 1). This dense gas, which remains at least partially self-shielded from ionizing radiation, provides an ideal environment for dust nucleation upon further cooling and expansion of the system (e.g. Derdzinski, Metzger & Lazzati 2017). Note that equation (7) neglects pressure support from magnetic fields or non-thermal ions, which may prevent the gas from compressing fully (though even non-thermal particles may cool rapidly, and magnetic fields self-generated at the shock could have a small coherence length and may decay downstream of the shock). In a 1D model, the thickness of the shell, |$w$|s, grows linearly with time t as it collects mass, approximately as |$w_\text{s} \approx c_{\rm \text{s}}t/\mathcal {M}$|⁠, and thus reaches a radial thickness |${\sim } R/\mathcal {M}^{2}$| in a time texp, where R ∼ |$v$|shtexp is the approximate radial size of the flow if it is expanding at a velocity ∼|$v$|sh. Fig. 2 shows the parameter space of upstream density and shock velocity separating adiabatic from radiative shocks for different ejecta sizes R = |$v$|shtexp. Our implicit assumption that the shocked gas cools at the optically thin rate is a good approximation as long as the timescale for photons to diffuse away from the shock tdiff ∼ τ(R/c) is much shorter than the expansion time texp ∼ R/|$v$|sh, where τ ∼ κρR is the optical depth and κ is the opacity. For tdiff ≫ texp, the energy density of radiation generated near the shock affects its structure (e.g. the shock jump is mediated by photon scattering rather than collisionless plasma processes; e.g. Katz, Budnik & Waxman 2010), a case not studied here. The condition τ = c/|$v$|sh which separates ‘effectively transparent’ from radiation-mediated shocks are shown by dashed lines in Fig. 2. The triangular wedge between the solid and dashed lines, as characterizes many shock-powered astrophysical transients, represents the region studied in this paper. Figure 2. Open in new tabDownload slide The regime of shocks studied in this paper in the space of upstream gas density |$n = \rho /(\bar{\mu }m_\text{p})$| and shock velocity |$v$|sh. Solid lines show the conditions separating adiabatic from radiative shocks (tcool = texp, where tcool is the optically thin cooling time-scale of the post-shock gas and texp = R/|$v$|sh is the characteristic expansion or evolution timescale) for different values of the shock radius R = 1013–1015 cm. Dashed lines show the condition, τ < c/|$v$|sh, for the optically thin cooling assumption to be valid and for the shock to be mediated by collisionless plasma processes, where τ = ρκR is the optical depth to the shock and κ = κes ≃ 0.4 cm2 g−1 is the opacity, taken to be electron scattering of fully ionized hydrogen (this is just a convenient approximation, as the real opacity will likely be the result of Doppler-broadened absorption lines and will be wavelength-dependent). Figure 2. Open in new tabDownload slide The regime of shocks studied in this paper in the space of upstream gas density |$n = \rho /(\bar{\mu }m_\text{p})$| and shock velocity |$v$|sh. Solid lines show the conditions separating adiabatic from radiative shocks (tcool = texp, where tcool is the optically thin cooling time-scale of the post-shock gas and texp = R/|$v$|sh is the characteristic expansion or evolution timescale) for different values of the shock radius R = 1013–1015 cm. Dashed lines show the condition, τ < c/|$v$|sh, for the optically thin cooling assumption to be valid and for the shock to be mediated by collisionless plasma processes, where τ = ρκR is the optical depth to the shock and κ = κes ≃ 0.4 cm2 g−1 is the opacity, taken to be electron scattering of fully ionized hydrogen (this is just a convenient approximation, as the real opacity will likely be the result of Doppler-broadened absorption lines and will be wavelength-dependent). In this paper, we neglect thermal conduction, which introduces a cooling term to the hot post-shock gas of the form |$\nabla \cdot \boldsymbol {F}_{\text{c}}$| in addition to radiative losses, where |$\boldsymbol {F}_{\text{c}}$| is the conductive flux. Neglecting magnetic fields, the latter can be written as \begin{eqnarray*} \boldsymbol {F}_\text{c} = \nabla (\kappa T), \end{eqnarray*} (8) where κ = κ0T5/2 is the conductivity, |$\kappa _0 \approx (2\times 10^{-5}/{\rm ln}\, \Lambda _{\rm c}$|⁠) erg K−7/2 s−1 cm−1 is a constant, and |${\rm ln}\, \Lambda _{\rm c} \approx 10$| is the Coulomb logarithm. Under the assumption that temperature gradients are comparable to the cooling length, i.e. |$\nabla \sim 1/\mathcal {L}_{\rm cool}$|⁠, we can approximate |$\nabla \cdot \boldsymbol {F}_\text{c} \sim \kappa _0 T^{7/2}/(\mathcal {L}_{\rm cool})^{2}$|⁠, where |$\mathcal {L}_{\rm cool} \approx (k_\text{b} T v)/(n\Lambda) \sim (k_\text{b} T^{2} v_{\rm sh})/(n\Lambda T_{\rm sh})$|⁠. Here, we have generalized equation (4) to account for the expectation, for a 1D flow, that the velocity behind the shock |$v$| = |$v$|sh(T/Tsh) decreases due to mass flux conservation nv = constant as gas compresses at roughly constant pressure nT = constant. The ratio of the conductive to radiative cooling rates is therefore approximately given by \begin{eqnarray*} \frac{\nabla \cdot \boldsymbol {F}_{\text{c}}}{n^{2} \Lambda (T)} &\sim & \frac{\kappa _0 T_{\rm sh}^{2}\Lambda }{k_\text{b}^{2} T^{1/2} v_{\rm sh}^{2}} \nonumber \\ &\approx & 0.05 \left(\frac{T}{T_{\rm sh}}\right)^{-3/2}\left(\frac{v_{\rm sh}}{10^{3}\, {\rm km\, s^{-1}}}\right)^{-1}, \end{eqnarray*} (9) where we have again used the nova composition approximation for Λ. Conduction thus provides a moderate perturbation to radiative cooling in the immediate post-shock gas where T ≈ Tsh, but could become more important at lower temperatures T ≪ Tsh. However, also note that even a weak magnetic field greatly suppresses thermal conduction perpendicular to the local field direction (e.g. Narayan & Medvedev 2001). 3 NUMERICAL METHOD AND SIMULATION RUNS We simulate the evolution of radiative shocks using the moving-mesh hydrodynamical code rich (Yalinewich et al. 2015). The moving-mesh technique is ideal for this problem because of the large compression factors involved. rich solves the compressible Euler equations on a moving Voronoi mesh via a second-order Gudonov scheme. The boundary conditions are an inflowing gas with a velocity of |$v$|sh and an adiabatic index of γ = 5/3. The density and temperature of the upstream flow are in all cases set to ρ = 2.6 × 10−12 g cm−3 and T = Tmin = 104 K, respectively. Our suite of numerical runs are summarized in Table 1, labelled according to Mach number and numerical technique (see below). We cover a range of shock velocities (⁠|$v$|sh ≈ 50–500 km s−1), Mach numbers (⁠|$\mathcal {M} \sim 4\text{--}36$|⁠), cooling functions (both solar and classical nova ejecta abundances), and grid resolutions. Note that once the shocked gas undergoes the cooling instability, the velocity of the shocks decreases by a factor of 3/4, since the unstable shock front expands at a much lower velocity |$v_{\rm sh}/\mathcal {M} \ll v_{\rm sh}$| towards theupstream direction once it has lost pressure support rather than the value |$v$|sh/4 for an adiabatic shock. Since the shock starts out adiabatic at early times t < tcool, we therefore define Tsh to be a factor 16/9 higher than the one defined in equation (1) for all of our runs (except the lowest Mach number runs, which are stable). Also note that the grid resolutions in Table 1 are defined with respect to the initial state of the upstream gas, prior to when the gas has cooled and compressed. Although our grid cells compress with the flow, we typically achieve a resolution in the compressed regions that is a factor of ∼1/|$\mathcal {M}$| smaller than the initial upstream one. Most of our simulations are performed in two spatial dimensions, motivated by the largely 2D nature of the NTSI. We also perform a single 3D simulation, which, as we shall describe, undergoes a qualitatively similar evolution compared to the otherwise equivalent 2D case. Table 1. Simulation runs. Model |$v$|sh |$\mathcal {M}^{a}$| Coolingb Cell sizec Num.d (km s−1) (⁠|$\mathcal {L}_{\rm cool}$|⁠) M36_L 500 36 Nova 400 L M36Sol_L – – Solar – – †M363D_L – – Nova 300 – M36LowRes_L – – – 280 – M28_L 387 28 – 400 – M22_L 275 22 – – – M12_L 162 12 – – – M4_L 50 4 – 30 – M36_SL 500 36 – 400 SL M28_SL 387 28 – – – M22_SL 275 22 – – – M12_SL 162 12 – – – M4_SL 50 4 – 30 – Model |$v$|sh |$\mathcal {M}^{a}$| Coolingb Cell sizec Num.d (km s−1) (⁠|$\mathcal {L}_{\rm cool}$|⁠) M36_L 500 36 Nova 400 L M36Sol_L – – Solar – – †M363D_L – – Nova 300 – M36LowRes_L – – – 280 – M28_L 387 28 – 400 – M22_L 275 22 – – – M12_L 162 12 – – – M4_L 50 4 – 30 – M36_SL 500 36 – 400 SL M28_SL 387 28 – – – M22_SL 275 22 – – – M12_SL 162 12 – – – M4_SL 50 4 – 30 – Notes: aMach number. bComposition assumed in cooling function (Fig. 1). cInitial cell size, in units of the initial cooling length |$\mathcal {L}_{\rm cool}$|⁠. dNumerical technique: ‘Lagrangian’ (L) or ‘Semi-Lagrangian’ (SL); †3D simulation. Open in new tab Table 1. Simulation runs. Model |$v$|sh |$\mathcal {M}^{a}$| Coolingb Cell sizec Num.d (km s−1) (⁠|$\mathcal {L}_{\rm cool}$|⁠) M36_L 500 36 Nova 400 L M36Sol_L – – Solar – – †M363D_L – – Nova 300 – M36LowRes_L – – – 280 – M28_L 387 28 – 400 – M22_L 275 22 – – – M12_L 162 12 – – – M4_L 50 4 – 30 – M36_SL 500 36 – 400 SL M28_SL 387 28 – – – M22_SL 275 22 – – – M12_SL 162 12 – – – M4_SL 50 4 – 30 – Model |$v$|sh |$\mathcal {M}^{a}$| Coolingb Cell sizec Num.d (km s−1) (⁠|$\mathcal {L}_{\rm cool}$|⁠) M36_L 500 36 Nova 400 L M36Sol_L – – Solar – – †M363D_L – – Nova 300 – M36LowRes_L – – – 280 – M28_L 387 28 – 400 – M22_L 275 22 – – – M12_L 162 12 – – – M4_L 50 4 – 30 – M36_SL 500 36 – 400 SL M28_SL 387 28 – – – M22_SL 275 22 – – – M12_SL 162 12 – – – M4_SL 50 4 – 30 – Notes: aMach number. bComposition assumed in cooling function (Fig. 1). cInitial cell size, in units of the initial cooling length |$\mathcal {L}_{\rm cool}$|⁠. dNumerical technique: ‘Lagrangian’ (L) or ‘Semi-Lagrangian’ (SL); †3D simulation. Open in new tab For reasons described above (Fig. 2), the shocked gas is assumed to cool at the optically thin rate (equation 2). While this assumption is not valid for the UV and X-ray emission directly emitted by the hot gas, absorption and reprocessing of the emission occurs very rapidly to visual frequencies, where the opacity and diffusion time is much lower tdiff ≪ t. The colliding medium is therefore effectively optically thin for purposes of the hydrodynamic evolution. Fig. 1 shows the cooling functions used, calculated from version C17.00 of cloudy (Ferland et al. 2017) using the nova and gass10 abundances. The cooling is implemented via the exact cooling method proposed in Townsend (2009). Previous works on colliding radiative shocks (e.g. Kee et al. 2014) find a large reduction in the X-ray luminosity arising from radiative shocks (Vishniac 1994). However, some of this reduction may not be physical due to artificial numerical diffusion of thermal energy from hot to cold cells (Parkin & Pittard 2010). There are two ways to overcome this issue. First, one may take the brute force approach of increasing the number of grid cells per cooling length; however, the required resolution depends on the Mach number and can be prohibitively high for |$\mathcal {M} \gg 1$|⁠. The second method, which we have decided to adopt, is to eliminate the advection term in the numerical scheme altogether. Though this leads to moderate numerical inaccuracies, we find that it preserves the main qualitative features of the shock evolution while allowing numerical convergence on the reduction of the X-ray emission. We refer to this new approach as the ‘Lagrangian’ scheme, in contrast to our simulations employing the standard formulation with the advection term included, which we refer to as ‘Semi-Lagrangian’. Appendix A provides a detailed description of the ‘Lagrangian’ scheme and its validation. To diagnose the structure of the shock and its role in heating the gas, we locate the external shock fronts by identifying cells that have a density within 3 per cent and pressure within 10 per cent of the original inflow parameters, as well as a large pressure gradient and a negative velocity divergence. Once those cells have been identified, we locate the cell with the highest temperature that is within three cells distance and along the pressure gradient direction, thus also defining the local shock direction. From the temperature ratio between the upstream and downstream cells, we calculate the Mach number of the shock, using the regular jump conditions across shocks. We have also tried using the procedure outlined by Schaal & Springel (2015) and Schaal et al. (2016) to detect shocks; however, this gave many potential spurious detections in the turbulent region between the outer forward and reverse shocks. 4 RESULTS Our results are summarized in Figs 3–17 and Table 2. We begin with a walk-through of the time evolution of the shock interaction and development of the NTSI for our fiducial run employing the ‘Lagrangian’ technique with |$\mathcal {M} = 36$| (Section 4.1). We then discuss the implications of our suite of results for the temperature distribution of the shock luminosity (Section 4.2) and the acceleration of non-thermal ions (Section 4.3). Figure 3. Open in new tabDownload slide Snapshots of the post-shock temperature, in units of log10(T/K), for our fiducial model |$\text{M}36\_\text{L}$| (‘Lagrangian’ scheme and Mach number |$\mathcal {M} = 36$|⁠), shown times t = 0.3tcool (top panel), 0.75tcool (middle panel), and 1.45tcool (bottom panel), where tcool is the initial cooling timescale of the shocked gas. The first snapshot shows the shock, while it is still adiabatic and the shock front is stable, while the final snapshot shows the fully developed NTSI. Note the change in spatial scale of the final snapshot relative to the first two. Figure 3. Open in new tabDownload slide Snapshots of the post-shock temperature, in units of log10(T/K), for our fiducial model |$\text{M}36\_\text{L}$| (‘Lagrangian’ scheme and Mach number |$\mathcal {M} = 36$|⁠), shown times t = 0.3tcool (top panel), 0.75tcool (middle panel), and 1.45tcool (bottom panel), where tcool is the initial cooling timescale of the shocked gas. The first snapshot shows the shock, while it is still adiabatic and the shock front is stable, while the final snapshot shows the fully developed NTSI. Note the change in spatial scale of the final snapshot relative to the first two. Table 2. Key results of simulations. Model |$T_{\rm sh}^{a}$| (K) 〈T〉b (K) |$L(\gt T_{\rm sh}/3)/L_{\rm tot}^{c}$| |$\langle \epsilon _{\rm nth} \rangle _{\theta _{\rm max} = 40^\circ }^{d}$| |$\langle \epsilon _{\rm nth} \rangle _{\theta _{\rm max} = 55^\circ }^{d}$| 〈q − 4〉e |$L_{\rm shock}/L_{\rm tot}^{f}$| M36_L 7.2 × 106 K 3.3 × 105 K 3.6 × 10−2 1.5 × 10−2 2.7 × 10−2 6.5 × 10−3 0.82 M363D_L – 1.2 × 105 K 9.5 × 10−3 10−2 2.2 × 10−2 4.3 × 10−3 0.63 M36LowRes_L – 4.4 × 105 K 5.6 × 10−2 8 × 10−3 1.5 × 10−2 7 × 10−3 0.7 M36Sol_L – 3.4 × 105 K 3.3 × 10−2 9.6 × 10−3 1.9 × 10−2 7.5 × 10−3 0.77 M28_L 4.3 × 106 K 3 × 105 K 5.9 × 10−2 8.8 × 10−3 1.8 × 10−2 10−2 0.86 M20_L 2.2 × 106 K 1.75 × 105 K 7.38 × 10−2 1.8 × 10−3 4.9 × 10−3 2.2 × 10−2 0.79 M12_L 7.5 × 105 K 1.1 × 105 K 1.6 × 10−1 4 × 10−3 9.2 × 10−3 8.5 × 10−2 0.72 M4_L 7.2 × 104 K 3.3 × 104 K 0.72 0 0 – 1 M36_SL 7.2 × 106 K 5 × 104 K 1.2 × 10−3 2.2 × 10−2 3.2 × 10−2 8.1 × 10−3 0.59 M28_SL 4.3 × 106 K 4.7 × 104 K 2.9 × 10−3 2.2 × 10−2 3.3 × 10−2 1.1 × 10−2 0.69 M20_SL 2.2 × 106 K 6.6 × 104 K 1 × 10−2 8.8 × 10−3 1.6 × 10−2 2.2 × 10−2 0.52 M12_SL 7.5 × 105 K 1.1 × 105 K 1.5 × 10−1 2.7 × 10−3 6.7 × 10−3 8.6 × 10−2 0.57 M4_SL 7.2 × 104 K 3.3 × 104 K 0.72 0 0 – 1 Model |$T_{\rm sh}^{a}$| (K) 〈T〉b (K) |$L(\gt T_{\rm sh}/3)/L_{\rm tot}^{c}$| |$\langle \epsilon _{\rm nth} \rangle _{\theta _{\rm max} = 40^\circ }^{d}$| |$\langle \epsilon _{\rm nth} \rangle _{\theta _{\rm max} = 55^\circ }^{d}$| 〈q − 4〉e |$L_{\rm shock}/L_{\rm tot}^{f}$| M36_L 7.2 × 106 K 3.3 × 105 K 3.6 × 10−2 1.5 × 10−2 2.7 × 10−2 6.5 × 10−3 0.82 M363D_L – 1.2 × 105 K 9.5 × 10−3 10−2 2.2 × 10−2 4.3 × 10−3 0.63 M36LowRes_L – 4.4 × 105 K 5.6 × 10−2 8 × 10−3 1.5 × 10−2 7 × 10−3 0.7 M36Sol_L – 3.4 × 105 K 3.3 × 10−2 9.6 × 10−3 1.9 × 10−2 7.5 × 10−3 0.77 M28_L 4.3 × 106 K 3 × 105 K 5.9 × 10−2 8.8 × 10−3 1.8 × 10−2 10−2 0.86 M20_L 2.2 × 106 K 1.75 × 105 K 7.38 × 10−2 1.8 × 10−3 4.9 × 10−3 2.2 × 10−2 0.79 M12_L 7.5 × 105 K 1.1 × 105 K 1.6 × 10−1 4 × 10−3 9.2 × 10−3 8.5 × 10−2 0.72 M4_L 7.2 × 104 K 3.3 × 104 K 0.72 0 0 – 1 M36_SL 7.2 × 106 K 5 × 104 K 1.2 × 10−3 2.2 × 10−2 3.2 × 10−2 8.1 × 10−3 0.59 M28_SL 4.3 × 106 K 4.7 × 104 K 2.9 × 10−3 2.2 × 10−2 3.3 × 10−2 1.1 × 10−2 0.69 M20_SL 2.2 × 106 K 6.6 × 104 K 1 × 10−2 8.8 × 10−3 1.6 × 10−2 2.2 × 10−2 0.52 M12_SL 7.5 × 105 K 1.1 × 105 K 1.5 × 10−1 2.7 × 10−3 6.7 × 10−3 8.6 × 10−2 0.57 M4_SL 7.2 × 104 K 3.3 × 104 K 0.72 0 0 – 1 Notes: a Post-shock temperature for the head-on collision of a 1D flow with the simulation parameters (equation 1), as approximately achieved at early times when the shock is still adiabatic. b Average temperature of the radiating gas, as weighted by the radiated luminosity, evaluated once the NTSI has reached saturation. cFraction of the shock’s luminosity Ltot emitted at temperatures >Tsh/3, evaluated once the NTSI has reached saturation. dAverage non-thermal ion acceleration efficiency across the shock front, as weighted by the local shock power, assuming the maximum angle between the shock inclination angle and the upstream magnetic field for particle acceleration is 40° or 55°, respectively (these correspond to the approximate range of cut-off angles found by Caprioli & Spitkovsky 2014a). This quantity is also averaged over time, after the NTSI has reached saturation. eAverage power-law index of the momentum distribution of accelerated non-thermal ions, as weighted by the local shock power and measured relative to the value q = 4 predicted by linear DSA in the strong shock limit |$\mathcal {M} \gg 1$| (equation 16). fFraction of the kinetic energy of the upstream flow which is dissipated into thermal energy at the outer shock fronts; the remaining fraction is instead dissipated through weaker shocks or turbulence in the post-shock region. Open in new tab Table 2. Key results of simulations. Model |$T_{\rm sh}^{a}$| (K) 〈T〉b (K) |$L(\gt T_{\rm sh}/3)/L_{\rm tot}^{c}$| |$\langle \epsilon _{\rm nth} \rangle _{\theta _{\rm max} = 40^\circ }^{d}$| |$\langle \epsilon _{\rm nth} \rangle _{\theta _{\rm max} = 55^\circ }^{d}$| 〈q − 4〉e |$L_{\rm shock}/L_{\rm tot}^{f}$| M36_L 7.2 × 106 K 3.3 × 105 K 3.6 × 10−2 1.5 × 10−2 2.7 × 10−2 6.5 × 10−3 0.82 M363D_L – 1.2 × 105 K 9.5 × 10−3 10−2 2.2 × 10−2 4.3 × 10−3 0.63 M36LowRes_L – 4.4 × 105 K 5.6 × 10−2 8 × 10−3 1.5 × 10−2 7 × 10−3 0.7 M36Sol_L – 3.4 × 105 K 3.3 × 10−2 9.6 × 10−3 1.9 × 10−2 7.5 × 10−3 0.77 M28_L 4.3 × 106 K 3 × 105 K 5.9 × 10−2 8.8 × 10−3 1.8 × 10−2 10−2 0.86 M20_L 2.2 × 106 K 1.75 × 105 K 7.38 × 10−2 1.8 × 10−3 4.9 × 10−3 2.2 × 10−2 0.79 M12_L 7.5 × 105 K 1.1 × 105 K 1.6 × 10−1 4 × 10−3 9.2 × 10−3 8.5 × 10−2 0.72 M4_L 7.2 × 104 K 3.3 × 104 K 0.72 0 0 – 1 M36_SL 7.2 × 106 K 5 × 104 K 1.2 × 10−3 2.2 × 10−2 3.2 × 10−2 8.1 × 10−3 0.59 M28_SL 4.3 × 106 K 4.7 × 104 K 2.9 × 10−3 2.2 × 10−2 3.3 × 10−2 1.1 × 10−2 0.69 M20_SL 2.2 × 106 K 6.6 × 104 K 1 × 10−2 8.8 × 10−3 1.6 × 10−2 2.2 × 10−2 0.52 M12_SL 7.5 × 105 K 1.1 × 105 K 1.5 × 10−1 2.7 × 10−3 6.7 × 10−3 8.6 × 10−2 0.57 M4_SL 7.2 × 104 K 3.3 × 104 K 0.72 0 0 – 1 Model |$T_{\rm sh}^{a}$| (K) 〈T〉b (K) |$L(\gt T_{\rm sh}/3)/L_{\rm tot}^{c}$| |$\langle \epsilon _{\rm nth} \rangle _{\theta _{\rm max} = 40^\circ }^{d}$| |$\langle \epsilon _{\rm nth} \rangle _{\theta _{\rm max} = 55^\circ }^{d}$| 〈q − 4〉e |$L_{\rm shock}/L_{\rm tot}^{f}$| M36_L 7.2 × 106 K 3.3 × 105 K 3.6 × 10−2 1.5 × 10−2 2.7 × 10−2 6.5 × 10−3 0.82 M363D_L – 1.2 × 105 K 9.5 × 10−3 10−2 2.2 × 10−2 4.3 × 10−3 0.63 M36LowRes_L – 4.4 × 105 K 5.6 × 10−2 8 × 10−3 1.5 × 10−2 7 × 10−3 0.7 M36Sol_L – 3.4 × 105 K 3.3 × 10−2 9.6 × 10−3 1.9 × 10−2 7.5 × 10−3 0.77 M28_L 4.3 × 106 K 3 × 105 K 5.9 × 10−2 8.8 × 10−3 1.8 × 10−2 10−2 0.86 M20_L 2.2 × 106 K 1.75 × 105 K 7.38 × 10−2 1.8 × 10−3 4.9 × 10−3 2.2 × 10−2 0.79 M12_L 7.5 × 105 K 1.1 × 105 K 1.6 × 10−1 4 × 10−3 9.2 × 10−3 8.5 × 10−2 0.72 M4_L 7.2 × 104 K 3.3 × 104 K 0.72 0 0 – 1 M36_SL 7.2 × 106 K 5 × 104 K 1.2 × 10−3 2.2 × 10−2 3.2 × 10−2 8.1 × 10−3 0.59 M28_SL 4.3 × 106 K 4.7 × 104 K 2.9 × 10−3 2.2 × 10−2 3.3 × 10−2 1.1 × 10−2 0.69 M20_SL 2.2 × 106 K 6.6 × 104 K 1 × 10−2 8.8 × 10−3 1.6 × 10−2 2.2 × 10−2 0.52 M12_SL 7.5 × 105 K 1.1 × 105 K 1.5 × 10−1 2.7 × 10−3 6.7 × 10−3 8.6 × 10−2 0.57 M4_SL 7.2 × 104 K 3.3 × 104 K 0.72 0 0 – 1 Notes: a Post-shock temperature for the head-on collision of a 1D flow with the simulation parameters (equation 1), as approximately achieved at early times when the shock is still adiabatic. b Average temperature of the radiating gas, as weighted by the radiated luminosity, evaluated once the NTSI has reached saturation. cFraction of the shock’s luminosity Ltot emitted at temperatures >Tsh/3, evaluated once the NTSI has reached saturation. dAverage non-thermal ion acceleration efficiency across the shock front, as weighted by the local shock power, assuming the maximum angle between the shock inclination angle and the upstream magnetic field for particle acceleration is 40° or 55°, respectively (these correspond to the approximate range of cut-off angles found by Caprioli & Spitkovsky 2014a). This quantity is also averaged over time, after the NTSI has reached saturation. eAverage power-law index of the momentum distribution of accelerated non-thermal ions, as weighted by the local shock power and measured relative to the value q = 4 predicted by linear DSA in the strong shock limit |$\mathcal {M} \gg 1$| (equation 16). fFraction of the kinetic energy of the upstream flow which is dissipated into thermal energy at the outer shock fronts; the remaining fraction is instead dissipated through weaker shocks or turbulence in the post-shock region. Open in new tab 4.1 Time evolution and saturated state The collision between opposite flows results in the creation of strong ‘forward’ and ‘reverse’ shocks (top panel of Fig. 3), behind which the temperature is close to the value for a 1D planar shock, Tsh ≃ 7 × 106 K (equation 1). After one cooling time has elapsed, the temperature at the centre of the accumulating shell plummets from radiative losses (middle panel of Fig. 3). This loss of pressure generates a rarefaction wave, which travels from the centre out to the shock fronts. Once the rarefaction reaches the shocks, their loss of pressure support results in them being pushed back by the ram pressure of the upstream gas, collapsing around the cold dense central shell. After the collapse, the upstream flows now collide with the dense shell, which becomes unstable to the NTSI, driving the shock interface to an irregularshape. Due to the large bending angles in the central shell, most of the shocks are now oblique and, partially for this reason, produce gas with a temperature lower than Tsh (bottom panel of Fig. 3). The NTSI quickly saturates, and the average properties of the shocked gas reach a steady state. A movie of the evolution just described is available in the online version of the manuscript. In 3D, the qualitative picture remains the same; Fig. 4 shows the temperature topology of the 3D shock front at saturation. Figure 4. Open in new tabDownload slide Contour of the external shocked gas (set by density slightly larger than the inflow density) colour coded by temperature for our 3D simulation. The NTSI corrugates the shock front and creates a ‘mountain range’ topology. Figure 4. Open in new tabDownload slide Contour of the external shocked gas (set by density slightly larger than the inflow density) colour coded by temperature for our 3D simulation. The NTSI corrugates the shock front and creates a ‘mountain range’ topology. The above evolution qualitatively describes all of our simulations with sufficiently high Mach numbers (e.g. |$\mathcal {M} = 36$| for the example shown in Fig. 3). For our lowest Mach number runs with |$\mathcal {M} = 4$|⁠, the shock front largely remains laminar, a result explained by the fact that the NTSI requires a minimum compression ratio to grow (Vishniac 1994). Our lowest Mach number simulations also result in shocked gas which is thermally stable, which may explain some differences in the behaviour. Some aspects of the development of the NTSI in our simulations can also be understood from linear theory. Vishniac (1994) shows that for an isothermal thin shell of width |$w_\text{s} \approx c_{\text{s}}t/\mathcal {M}$|⁠, perturbations with wavelengths λ in the range |$w$|s ≤ λ ≪ cst are unstable,1 and their amplitude, A, grows to saturation at a rate \begin{eqnarray*} t^{-1}_{\rm NTSI}\approx \frac{A^{1/2}c_\text{s}}{\lambda ^{3/2}}. \end{eqnarray*} (10) In order for the perturbations to grow, an initial amplitude A ≈ |$w$|s is required. At early times, the small wavelengths ∼|$w$|s grow the fastest, but as the width of the central shell increases, they become stable and are overtaken by the largest unstable wavelength λmax = cst. The latter, which defines the largest length-scale characterizing the ‘corrugated’ shape of the shock front, is a factor of |$\mathcal {M}$| times larger than the naive estimate of the shell thickness |$w_\text{s} \approx c_{\text{s}}t/\mathcal {M}$| for a 1D flow, but a factor |${\sim } 1/\mathcal {M}$| times smaller than the shock radius R ∼ |$v$|sht for an expanding outflow of velocity ∼|$v$|sh. This hierarchy of scales is schematically illustrated in Fig. 5. Figure 5. Open in new tabDownload slide Schematic illustration of the hierarchy of length-scales |$w$|s ≪ λmax ≪ R in radially expanding high Mach number |$\mathcal {M} \gg 1$| radiative shocks in an expanding medium with radial velocity comparable to the internal shock velocity |$v$|sh. The radial scale of the flow on timescales ∼t is R ∼ |$v$|sht. In the naive 1D picture (shown on the right), the shell of cold gas bounded between the shocks is thin, with a radial width |$w_\text{s} \approx c_{\rm s}t/\mathcal {M} \sim R/\mathcal {M}^{2}$|⁠. However, in multiple dimensions (shown on the left), the shell is unstable to the NTSI, which in its saturated state creates structure in the shock front up to a characteristic length-scale |$\lambda _{\rm max} = c_{\rm s}t \sim R/\mathcal {M}$|⁠, where cs is the sound speed of the cold gas. In both cases, the thickness of the layer of hot gas bounding the cold is approximately given by the cooling length |$\mathcal {L}_{\rm cool}$| (equation 4). Figure 5. Open in new tabDownload slide Schematic illustration of the hierarchy of length-scales |$w$|s ≪ λmax ≪ R in radially expanding high Mach number |$\mathcal {M} \gg 1$| radiative shocks in an expanding medium with radial velocity comparable to the internal shock velocity |$v$|sh. The radial scale of the flow on timescales ∼t is R ∼ |$v$|sht. In the naive 1D picture (shown on the right), the shell of cold gas bounded between the shocks is thin, with a radial width |$w_\text{s} \approx c_{\rm s}t/\mathcal {M} \sim R/\mathcal {M}^{2}$|⁠. However, in multiple dimensions (shown on the left), the shell is unstable to the NTSI, which in its saturated state creates structure in the shock front up to a characteristic length-scale |$\lambda _{\rm max} = c_{\rm s}t \sim R/\mathcal {M}$|⁠, where cs is the sound speed of the cold gas. In both cases, the thickness of the layer of hot gas bounding the cold is approximately given by the cooling length |$\mathcal {L}_{\rm cool}$| (equation 4). For high Mach number shocks, radiative cooling results in an enormous density contrast |$\rho _{\rm c}/\rho _{\rm h} \sim \mathcal {M}^{2}$| between the cold (∼104 K) and hot (∼106–107 K) phases of the post-shock gas (equation 7), which can give rise to spurious cooling due to advection terms in the Euler equation (‘numerical conduction’). To illustrate this effect, Figs 6 and 7 compare the temperature and density of the shocked gas after the NTSI has saturated between our ‘Lagrangian’ and ‘Semi-Lagrangian’ simulations for the fiducial |$\mathcal {M} = 36$| case. In both schemes, rapid cooling and thermal instability in the post-shock gas creates filamentary cold structures. However, in the ‘Semi-Lagrangian’ case artificial conduction of thermal energy from hot to cold cells (which radiate the received energy immediately because of the high cooling rate at low temperatures) effectively eliminates the hot phase between the shock fronts and thus artificially suppresses the temperature of the radiating gas (Section 4.2). In contrast, for our runs which employ the ‘Lagrangian’ technique, the post-shock region maintains separate cold and hot phases, resulting in a higher overall temperature distribution. Going to lower Mach numbers reduces the difference between the outcomes of the two schemes, to the limit where for |$\mathcal {M} \simeq 4,$| the two techniques give very similar results. Figure 6. Open in new tabDownload slide Demonstration of how numerical conduction artificially suppresses the hot phase of the post-shock gas. Here, we show the temperature of the gas between the shocks, in units of log10(T/K), in the saturated state of the NTSI. We compare our fiducial model |$\text{M}36\_\text{L}$| which adopts the ‘Lagrangian’ numerical scheme to minimize numerical conduction (top panel) to the ‘Semi-Lagrangian’ run |$\text{M}36\_\text{SL}$| (bottom panel). Figure 6. Open in new tabDownload slide Demonstration of how numerical conduction artificially suppresses the hot phase of the post-shock gas. Here, we show the temperature of the gas between the shocks, in units of log10(T/K), in the saturated state of the NTSI. We compare our fiducial model |$\text{M}36\_\text{L}$| which adopts the ‘Lagrangian’ numerical scheme to minimize numerical conduction (top panel) to the ‘Semi-Lagrangian’ run |$\text{M}36\_\text{SL}$| (bottom panel). Figure 7. Open in new tabDownload slide Similar to Fig. 6, but comparing the density profiles of the |$\text{M36}\_\text{L}$| and |$\text{M}36\_\text{SL}$| models in the saturated state, in units of log10 relative to the upstream density. Figure 7. Open in new tabDownload slide Similar to Fig. 6, but comparing the density profiles of the |$\text{M36}\_\text{L}$| and |$\text{M}36\_\text{SL}$| models in the saturated state, in units of log10 relative to the upstream density. 4.2 Temperature distribution of post-shock gas The geometrical structure of radiative shocks driven by thermal instabilities and the NTSI results in a reduction in the temperature of the radiating gas (e.g. Kee et al. 2014). Figs 8 and 9 show for each of our simulations the time evolution of the average temperature, 〈T〉, weighted by the radiative luminosity. We have normalized the latter to immediate post-shock temperature for a head-on adiabatic 1D flow, Tsh (equation 1). At early times (t ≪ tcool), when the shock is still in the adiabatic state, we have 〈T〉 ≈ Tsh, as expected. For high Mach numbers however, the temperature plummets to a minimum after the shock collapses due to runaway cooling at times t ∼ tcool. Following this minimum, the temperature rises again slightly, eventually reaching an asymptotic value 〈T〉 which is typically one to two orders of magnitude smaller than Tsh, depending on the Mach number. Figure 8. Open in new tabDownload slide Average temperature of the post-shock gas, weighted by the radiative luminosity and normalized to that for a head-one adiabatic shock Tsh (equation 1), as function of time in units of the initial cooling time, for the different simulation runs with Mach number |$\mathcal {M} = 36$| (Table 1). This figure illustrates that the standard ‘Semi-Lagrangian’ simulation technique (red curve) produces much lower average post-shock temperatures than our new ‘Lagrangian’ technique (green curve; see also Fig. 6), a result we attribute to the suppression of numerical conduction in the latter. Also note that the 3D ‘Lagrangian’ simulations (lavender curve) saturates to somewhat lower average temperatures than the otherwise equivalent2D run. Figure 8. Open in new tabDownload slide Average temperature of the post-shock gas, weighted by the radiative luminosity and normalized to that for a head-one adiabatic shock Tsh (equation 1), as function of time in units of the initial cooling time, for the different simulation runs with Mach number |$\mathcal {M} = 36$| (Table 1). This figure illustrates that the standard ‘Semi-Lagrangian’ simulation technique (red curve) produces much lower average post-shock temperatures than our new ‘Lagrangian’ technique (green curve; see also Fig. 6), a result we attribute to the suppression of numerical conduction in the latter. Also note that the 3D ‘Lagrangian’ simulations (lavender curve) saturates to somewhat lower average temperatures than the otherwise equivalent2D run. Figure 9. Open in new tabDownload slide Average temperature of the post-shock gas, weighted by the radiative luminosity and normalized to that for a head-one adiabatic shock Tsh (equation 1), as function of time in units of the initial cooling time, for the different Mach number runs (Table 1). Figure 9. Open in new tabDownload slide Average temperature of the post-shock gas, weighted by the radiative luminosity and normalized to that for a head-one adiabatic shock Tsh (equation 1), as function of time in units of the initial cooling time, for the different Mach number runs (Table 1). Fig. 8 shows the impact of numerical resolution and technique on the temperature in our results. At early times, before cooling instabilities develop and while the post-shock structure is still well resolved, our ‘Semi-Lagrangian’ simulations converge with resolution. However, in the late-time saturated state of the NTSI, the average temperature of the ‘Semi-Lagrangian’ simulations is much lower than the otherwise equivalent ‘Lagrangian’ case (red versus green lines in Fig. 8). The temperature furthermore increases with resolution, albeit slowly, and is therefore not converged for the ‘Semi-Lagrangian’ runs. By contrast, for the simulations which employ the ‘Lagrangian’ technique to prohibit advection between cells, the average temperature is almost an order of magnitude higher for high Mach numbers than even the highest resolution ‘Semi-Lagrangian’ simulations. This implies that artificial (numerical) conduction plays a crucial role in the post-shock temperature distribution and confirms that our ‘Semi-Lagrangian’ simulations, and likely those employed in some past numerical works, are under-resolved. Our ‘Lagrangian’ simulations also show evidence for convergence at late times, as shown by the similarity of the brown and green lines in Fig. 8. Fig. 10 shows the probability density function of the radiated luminosity as function of the gas temperature in the saturated state for different runs. The ‘Lagrangian’ simulations again show a higher emission temperature than the otherwise equivalent ‘Semi-Lagrangian’ case, whereas the lowest Mach run which is stable, emits all its luminosity at temperatures comparable to Tsh. Table 2 gives for each simulation the emission-weighted temperature and fraction of the radiated luminosity, L(> Tsh/3)/Ltot emitted by cells at temperatures greater than Tsh/3. Figure 10. Open in new tabDownload slide Probability density function of the shock luminosity as a function of temperature, shown for a selected sample of runs from our suite of simulations (Table 1). Low Mach number shocks are not susceptible to the NTSI and therefore dissipate most power at T ≈ Tsh, while for high Mach number, the shock front becomes highly oblique and releases most of its energy at substantially lower temperatures. Figure 10. Open in new tabDownload slide Probability density function of the shock luminosity as a function of temperature, shown for a selected sample of runs from our suite of simulations (Table 1). Low Mach number shocks are not susceptible to the NTSI and therefore dissipate most power at T ≈ Tsh, while for high Mach number, the shock front becomes highly oblique and releases most of its energy at substantially lower temperatures. Fig. 11 shows the X-ray reduction as a function of the Mach number, which for our Lagrangian runs is reasonably fit by the functional form \begin{eqnarray*} \frac{L(\gt T_{\rm sh}/3)}{L_{\rm tot}} \approx \frac{4.5}{\mathcal {M}^{4/3}}. \end{eqnarray*} (11) This result is close to the conjecture of Kee et al. (2014) that the fraction of the radiation at X-ray emitting temperatures ∼Tsh is reduced by a factor that scales as |${\sim } 1/\mathcal {M}$| compared to the naive expectation of a 1D flow. Figure 11. Open in new tabDownload slide Fraction of the total luminosity radiated behind the shock with X-ray temperature T ≳ Tsh/3 as a function of the Mach number for our suite of Lagrangian simulations (Table 1), where Tsh (equation 1) is the temperature for a 1D adiabatic shock. Figure 11. Open in new tabDownload slide Fraction of the total luminosity radiated behind the shock with X-ray temperature T ≳ Tsh/3 as a function of the Mach number for our suite of Lagrangian simulations (Table 1), where Tsh (equation 1) is the temperature for a 1D adiabatic shock. Fig. 14 shows a 2D histogram of the shock power for our fiducial |$\mathcal {M} = 36$| Lagrangian run as a function of Mach number and incident angle, θB, between the shock front and the assumed tangential magnetic field (the Y-axis). We find that our shock detection algorithm successfully detects 70–86 per cent of the incoming kinetic energy as being dissipated (converted into heat) at the external shock fronts, while the rest is dissipated in the turbulent region between the forward and reverse shock. In the next section, we describe how the inclination and Mach number distribution of the shocks plays an important role in the non-thermal ion acceleration. What physical effects cause the suppression in the X-ray luminosity of the shocks shown in Fig. 11? Part of the reduction is related to the lower post-shock temperature resulting from the oblique shocks Tsh∝cos 2α; however, the distribution of shock power with angle shown in Fig. 14 makes clear that this can hardly account for all of the suppression that we find (e.g. a factor of ∼30 for |$\mathcal {M} = 36$|⁠). Some insight may be provided by Fig. 12, which shows the distribution of pressure and temperature of the post-shock gas, weighted by its radiated luminosity, for our fiducial 2D and 3D |$\mathcal {M} = 36$| runs. The faint patch at high T and high P represents the hot volume-filling post-shock gas, which cools radiatively towards the low T, low P corner of the diagram. Importantly, however, much of the cool emitting gas is underpressurized relative to the hot surrounding medium. This pressure difference enables the hot gas to perform P ·dV work on the cool gas, e.g. through a series of weak shocks, robbing the hot gas of its thermal energy, which is instead radiated with high efficiency at lower temperature. To strengthen this conclusion, we show in Fig. 13 a zoomed-in region of the interface between the hot shocked gas and the cold dense gas. It is clear that the cells that are emitting most of the radiation are the cold, underpressured cells adjacent to the hot shocked gas. The underpressure we observe, which presumably results from runaway (acausal) cooling, appears to be a multidimensional effect, as it does not appear in 1D simulations we have performed for the same shock parameters. Figure 12. Open in new tabDownload slide Temperature and pressure distribution of the post-shock gas, weighted by its radiated luminosity, as measured at a snapshot in time following saturation of the NTSI. The top panel shows the 2D run |$\text{M}36\_\text{L}$|⁠, while the bottom panel shows the 3D case |$\text{M36}_{\rm 3D}\_\text{L}$|⁠. The hot gas has a higher pressure than much of the emitting cold gas, allowing the hot gas to transfer thermal energy to the cold gas via PdV work (e.g. weak shocks) before it is radiated away. Figure 12. Open in new tabDownload slide Temperature and pressure distribution of the post-shock gas, weighted by its radiated luminosity, as measured at a snapshot in time following saturation of the NTSI. The top panel shows the 2D run |$\text{M}36\_\text{L}$|⁠, while the bottom panel shows the 3D case |$\text{M36}_{\rm 3D}\_\text{L}$|⁠. The hot gas has a higher pressure than much of the emitting cold gas, allowing the hot gas to transfer thermal energy to the cold gas via PdV work (e.g. weak shocks) before it is radiated away. Figure 13. Open in new tabDownload slide Zoomed-in view of the interface between the outer (hot) shocked gas and the inner (cold) dense phase. The hot gas has a higher pressure and drives weak shocks into the cold dense phase, which is able to efficiently radiate the energy away before it can be radiated by the hot phase. This effect is partially responsible for reducing the average temperature of the post-shock gas as compared to the naive expectation for a planar 1D compression analysis. Figure 13. Open in new tabDownload slide Zoomed-in view of the interface between the outer (hot) shocked gas and the inner (cold) dense phase. The hot gas has a higher pressure and drives weak shocks into the cold dense phase, which is able to efficiently radiate the energy away before it can be radiated by the hot phase. This effect is partially responsible for reducing the average temperature of the post-shock gas as compared to the naive expectation for a planar 1D compression analysis. We now consider an order-of-magnitude analytic estimate of this effect. The rate of shock heating of the cold gas by the hot gas is approximately given by \begin{eqnarray*} L_{\rm c} \approx P_{\rm h}v_{\rm c}A_{\rm c} \sim kT_{\rm h}n_{\rm h}A_{\rm c}v_{\rm sh}(\rho _{\rm h}/\rho _{\rm c})^{1/2}, \end{eqnarray*} (12) where Ph is the pressure of the hot phase, Ac is the surface area of the cold gas, and |$v$|c is the velocity of the shock driven into the cold gas by the hot phase. The latter is estimated as |$v$|c = |$v$|sh(ρh/ρc)1/2 by equating the pressure of the hot gas Ph with the ram pressure |$\rho _{\rm c}v_{\rm c}^{2}/2$| and using |$P_{\rm h}/\rho _{\rm h} \approx v_{\rm sh}^{2}$|⁠. The hot gas radiates the energy at a rate \begin{eqnarray*} L_{\rm rad} \approx \frac{(3/2)kT_{\rm h}n_{\rm h}V_{\rm h}}{t_{\rm cool}} \sim kT_{\rm h}n_{\rm h}R^{2}v_{\rm sh}, \end{eqnarray*} (13) where Vh = R2Δ is the volume of the hot gas across the shock front of area R2, and |$\Delta \approx \mathcal {L}_{\rm cool} = \eta R$| is the thickness set by the cooling length (equation 4). Combining results, we conclude that \begin{eqnarray*} \frac{L_{\rm c}}{L_{\rm rad}} \approx \frac{A_{\rm c}}{R^{2}}\left(\frac{\rho _{\rm h}}{\rho _{\rm c}}\right)^{1/2}. \end{eqnarray*} (14) Thus, depending on to what extent the surface area of the (topologically complex) structure of cold filaments Ac exceeds the naive planar expectation ∼R2, and on the ratio of the densities at the interface between the two media (Fig. 13), we can have Lc ≳ Lrad, i.e. the shock work done by the hot phase on the cold gas can indeed be comparable to that directly radiated by the hot gas. Since the cold filamentary gas has a larger surface area in 3D than in 2D, this could contribute to the lower average temperature of the radiation in our 3D simulations. 4.3 Non-thermal ion acceleration Diffusive shock acceleration (DSA; e.g. Drury 1983; Blandford & Eichler 1987) is the process by which a small fraction of charged particles crossing a shock front are accelerated to much higher energies than the thermal population by means of repeated reflections across the shock interface (usually due to their interaction with magnetic turbulence, often generated by the current of the non-thermal particles themselves). When active, linear DSA theory predicts that particles of mass m and momentum p are accelerated to a spectrum f(p)∝p−q, corresponding to an energy spectrum \begin{eqnarray*} \frac{\text{d}N}{\text{d}E} = 4\pi p^{2}f(p)\frac{\text{d}p}{\text{d}E}. \end{eqnarray*} (15) In the non-relativistic regime E = p2/2m, such that dp/dE∝E−1/2 and dN/dE∝E(1 − q)/2; in the relativistic limit, E∝p and dN/dE∝E2 − q. The index q is related to the compression factor r according to \begin{eqnarray*} q = \frac{3r}{r-1};\, \, \, \, \, r = \frac{(\gamma +1)\mathcal {M}^{2}}{(\gamma -1)\mathcal {M}^{2} + 2}, \end{eqnarray*} (16) where γ is the adiabatic index. In the limit of a strong shock |$\mathcal {M} \gg 1$|⁠, we have r = q = 4 and the spectrum at high energies is a power-law dN/dE∝E−β with β = q − 2 ≈ 2. However, for lower |$\mathcal {M}$| as obtained in the oblique shocks induced by the NTSI, r ≲ 4(q ≳ 4), and we expect β ≳ 2. Values of q > 4 can also result from non-linear effects at the shock due to the feedback of cosmic rays on the compression factor (e.g. Malkov & Drury 2001), but we ignore this effect here since the acceleration efficiencies we find are quite modest. The power-law index β of the non-thermal particle energy distribution depends on the Mach number of the shock. However, the overall efficiency εnt with which the kinetic power of the shock is converted into non-thermal particle energy instead depends mainly on the orientation of the magnetic field ahead of the shock. Caprioli & Spitkovsky (2014a) show that εnt ≈ 0.05–0.2 for upstream magnetic fields quasi-parallel to the shock normal (depending weakly on the Mach number), but that εnt rapidly decreases for magnetic field inclinations θB ≳ 40°, approaching εnt ≈ 0 for θB ≳ θmax = 40°–55°. Because in the transient astrophysical systems of interest the upstream magnetic field is expected to be nearly perpendicular to the radially propagating shock fronts (i.e. α − 90° = θB), little or no particle acceleration is expected if the shock front remains smooth, in conflict with observations of GeV gamma-rays from classical novae. However, the corrugated shock structure created by the NTSI for high Mach number results in localized regions where θB ≲ θmax (Fig. 14), as illustrated schematically in Fig. 15, resulting in a smaller but non-zero average particle acceleration efficiency. The impact of fluctuations in the direction of the upstream magnetic field relative to the shock normal on the injection of ions into the DSA cycle has been explored in previous works (e.g. Guo & Giacalone 2010); however, the instabilities explored here provide a natural and self-consistent mechanism for creating the requisitegeometry. Figure 14. Open in new tabDownload slide Fraction of the upstream kinetic power which is dissipated across the radiative shock front as a function of the local Mach number and the angle θB between the shock front and the assumed toroidal magnetic field, shown for our fiducial 2D run |$\text{M}36\_\text{L}$| (top panel) and the equivalent 3D case |$\text{M}36_{\rm 3D}\_\text{L}$| (bottom panel). A vertical dashed line shows the critical angle, θmax ≈ 40°–55°, below which the shock accelerates relativistic ions with high efficiency (Caprioli & Spitkovsky 2014a), under the assumption that the upstream magnetic field is perpendicular to the upstreamvelocity. Figure 14. Open in new tabDownload slide Fraction of the upstream kinetic power which is dissipated across the radiative shock front as a function of the local Mach number and the angle θB between the shock front and the assumed toroidal magnetic field, shown for our fiducial 2D run |$\text{M}36\_\text{L}$| (top panel) and the equivalent 3D case |$\text{M}36_{\rm 3D}\_\text{L}$| (bottom panel). A vertical dashed line shows the critical angle, θmax ≈ 40°–55°, below which the shock accelerates relativistic ions with high efficiency (Caprioli & Spitkovsky 2014a), under the assumption that the upstream magnetic field is perpendicular to the upstreamvelocity. Figure 15. Open in new tabDownload slide Schematic diagram of the impact of the irregular geometry of the high Mach number radiative shock front on diffusive acceleration of relativistic ions. PIC simulations show that ion acceleration is limited to shocks for which the angle θB between the upstream magnetic field |$\boldsymbol {B}$| and the shock normal |$\hat{n}$| is less than a critical angle θmax ≈ 40°–55° (‘quasi-parallel’ shocks; e.g. Caprioli & Spitkovsky 2014a). In most astrophysical transients, due to magnetic flux conservation in the outflow, |$\boldsymbol {B}$| is perpendicular to the radial upstream velocity |$\boldsymbol {v}_{0}$|⁠, and therefore no acceleration can occur because for a planar shock |$\hat{n}$| is everywhere parallel to |$\boldsymbol {v}_{0}$|⁠. However, the corrugated shock front which results from the NTSI creates localized regions where θB < θmax is satisfied (Fig. 14), resulting in a low but non-zero net ion acceleration efficiencies of εnt ∼ 0.01 (Fig. 16). Figure 15. Open in new tabDownload slide Schematic diagram of the impact of the irregular geometry of the high Mach number radiative shock front on diffusive acceleration of relativistic ions. PIC simulations show that ion acceleration is limited to shocks for which the angle θB between the upstream magnetic field |$\boldsymbol {B}$| and the shock normal |$\hat{n}$| is less than a critical angle θmax ≈ 40°–55° (‘quasi-parallel’ shocks; e.g. Caprioli & Spitkovsky 2014a). In most astrophysical transients, due to magnetic flux conservation in the outflow, |$\boldsymbol {B}$| is perpendicular to the radial upstream velocity |$\boldsymbol {v}_{0}$|⁠, and therefore no acceleration can occur because for a planar shock |$\hat{n}$| is everywhere parallel to |$\boldsymbol {v}_{0}$|⁠. However, the corrugated shock front which results from the NTSI creates localized regions where θB < θmax is satisfied (Fig. 14), resulting in a low but non-zero net ion acceleration efficiencies of εnt ∼ 0.01 (Fig. 16). Fig. 16 shows our calculation of the non-thermal particle acceleration efficiency, εnt, as a function of Mach number. These are obtained by time-averaging the output for each of the numerical simulations and weighing the kinetic power of the shocked cells which satisfy the θB ≲ θmax condition of Caprioli & Spitkovsky (2014a). For concreteness, we approximate the results of the latter by assuming εnt = 0.1 for θB ≤ 55° and εnt = 0 for θB > 55°; however, our results are not sensitive to the precise transition angle (varying by a factor ≲ 2 if we instead take θmax = 40°; Table 2).Additionally, we assume that the magnetic field is not altered during the passage of the shock wave. We find that the acceleration efficiency εnt is roughly ≈0.01, depending only weakly on Mach number provided the latter is large enough for growth of the NTSI. This value agrees remarkably well with the ion acceleration efficiency inferred from modelling the optical and gamma-ray emission from ASASSN-16ma (Li et al. 2017). Figure 16. Open in new tabDownload slide Average non-thermal particle acceleration efficiency, εnt, of radiative shocks as a function of the Mach number. The upper border of the shaded region is given by a maximum angle of 55° between the shock normal and the upstream magnetic field and the lower border is given by 40°. Also shown for comparison is the derived efficiency from calorimetry applied to the gamma-ray emitting classical nova ASASSN-16ma (Li et al. 2017). Figure 16. Open in new tabDownload slide Average non-thermal particle acceleration efficiency, εnt, of radiative shocks as a function of the Mach number. The upper border of the shaded region is given by a maximum angle of 55° between the shock normal and the upstream magnetic field and the lower border is given by 40°. Also shown for comparison is the derived efficiency from calorimetry applied to the gamma-ray emitting classical nova ASASSN-16ma (Li et al. 2017). Fig. 17 shows our calculation of the average non-thermal spectrum for each of our ‘Lagrangian’ runs. These are obtained by combining the shock luminosity-weighted Mach number and inclination distribution (Fig. 14) with the results of equations (15) and (16) for |$q(\mathcal {M})$|⁠, again only counting regions of the shock which obey the θB ≲ θmax condition (we also give in Table 2 the shock luminosity-weighted value of the ion momentum index q). Figure 17. Open in new tabDownload slide Energy spectrum of non-thermal ions accelerated at the outer shock fronts, calculated by applying equation (16) locally to each patch of the shock interface, for our ‘Lagrangian’ runs with different Mach numbers. For comparison with a dashed line is the much steeper slope inferred by modelling the gamma-ray spectrum of ASASSN-16ma (Li et al. 2017). As we discuss further in the text, this discrepancy could be due to the preferential escape of high-energy cosmic rays (before producing gamma-rays), or due to an additional source of particle acceleration with a steeper index in the turbulent post-shock region. Figure 17. Open in new tabDownload slide Energy spectrum of non-thermal ions accelerated at the outer shock fronts, calculated by applying equation (16) locally to each patch of the shock interface, for our ‘Lagrangian’ runs with different Mach numbers. For comparison with a dashed line is the much steeper slope inferred by modelling the gamma-ray spectrum of ASASSN-16ma (Li et al. 2017). As we discuss further in the text, this discrepancy could be due to the preferential escape of high-energy cosmic rays (before producing gamma-rays), or due to an additional source of particle acceleration with a steeper index in the turbulent post-shock region. Despite the reduction in Mach number due to the irregular shape of the shock front, we find that the energy spectrum of relativistic particles for our high |$\mathcal {M}$| simulations is very close to the dN/dE∝E−β with β ≃ 2 expected in the |$\mathcal {M} \gg 1$| limit. This behavior is inconsistent with the steeper value β ≈ 2.7 of the ion spectrum inferred by modelling the gamma-ray emission of ASASSN-16ma (Li et al. 2017). Part of this discrepancy could be due to the energy-dependent escape of cosmic ray ions from the shock front, such that higher energy particles with large Larmor radii have a greater chance of escaping the surrounding medium without producing gamma-rays (see below). Alternatively, we may be missing an additional source of particle acceleration with a softer spectrum than the outer shock, such as that due to lower Mach number shocks or turbulence occurring downstream. If we compare the kinetic power dissipated by the external shocks as measured from the simulation to the total upstream kinetic energy flux, we find that only ≈60–80 per cent is accounted for, as tabulated in the final column of Table 2. This suggests that a sizable fraction of the free kinetic energy of the flow is dissipated in the chaotic downstream, through oblique shocks or turbulence. We speculate that this could contribute a softer non-thermal ion spectrum than that accelerated by the relatively strong externalshocks. We perform an additional check before proceeding. We have calculated the properties of shock-accelerated ions by summing the contributions separately from each position along the shock front, as if the acceleration process acts independently at each local position. This is a reasonable approximation if the ions ‘see’ an approximately infinite planar shock during the acceleration process. This in turn requires that the length-scale over which cosmic rays diffuse upstream in the process of being accelerated, λdiff ≈ D/|$v$|sh, be small compared to the length-scale over which the shock or post-shock turbulent region in changing, where here D = rgc/3 is the diffusion coefficient (in the limit of Bohm diffusion, Caprioli & Spitkovsky 2014b) and rg = E/eB is the ion gyro radius, where e is the electron charge. When thermal instabilities and the NTSI is active, the length-scale for changes in the shock front or post-shock structure range from the radiative cooling length |$\mathcal {L}_{\rm cool}$| (equation 4) to the coherence size of the largest corrugations |$\lambda _{\rm max} = v_{\rm sh}t_{\rm exp}/\mathcal {M}$| (Fig. 5). Considering the cooling length as the relevant one,we find \begin{eqnarray*} \frac{\lambda _{\rm diff}}{\mathcal {L}_{\rm cool}} &=& 0.3\left(\frac{\epsilon _{B}}{0.01}\right)^{-1/2}\left(\frac{E}{100\, {\rm GeV}}\right)\times \nonumber \\ &&\times \,\left(\frac{n}{10^{10}\, {\rm cm^{-3}}}\right)^{1/2}\left(\frac{v_{\rm sh}}{10^{3}\, {\rm km\, s^{-1}}}\right)^{-7}, \end{eqnarray*} (17) where we have estimated the magnetic field strength near the shock as |$B \simeq (6\pi \epsilon _{\text{B}}\mu m_\text{p} n v_{\rm sh}^{2})^{1/2}$| and εB ≪ 1 is the fraction of the kinetic power of the shock placed into magnetic field energy. Equation (17) shows that, because the energy of ions at the injection scale is very low E ∼ few × kTsh, the sub-structure of the shock introduced by thermal instabilities and the NTSI should not affect the injection of ions and thus the acceleration efficiency. However, higher energy particles E ≲ 100–1000 GeV which produce the gamma-rays detectable by Fermi-LAT can in principle have |$\lambda _{\rm diff} \gtrsim \mathcal {L}_{\rm cool}$| for shock velocities and densities relevant to classical novae and other transients. The preferential escape of higher energy particles during the acceleration process could therefore play a role in the steeper than expected ion spectrum inferred for classical novae such as ASSASN-16ma (Fig. 17). 5 DISCUSSION AND CONCLUSIONS We have presented the results of 2D and 3D moving-mesh hydrodynamical simulations of the radiative shocks created by dual equal-density colliding flows. Our goal was to assess the impact of the thin-shell instability on the temperature distribution of the thermal radiation and the efficiency of non-thermal ion acceleration relevant to the non-thermal X-ray and gamma-ray radiation. Our main findings are summarized as follows: Radiative shocks play an important role in a variety of astrophysical transients, such as SNe IIn , classical novae, and stellar mergers. When the photon diffusion time away from the shock is much less than the expansion time (Fig. 2), the shock dynamics are well approximated assuming optically thin radiative cooling. High Mach number radiative shocks are susceptible to thermal instabilities (Chevalier & Imamura 1982) and the NTSI (Vishniac 1994), the latter of which in its saturated non-linear state creates a highly corrugated shock front (a ‘mountain range’ in 3D) and a complex turbulent downstream characterized by cold filaments immersed in a hot volume-filling gas (Figs 3, 4, and 6). The radial extent of the cold filamentary region is enhanced as compared to the naive 1D laminar compression picture, but it remains a factor of at least |${\sim } 1/\mathcal {M}$| thinner than it would otherwise be for an adiabatic shock (Fig. 5). This difference between the naive 1D laminar flow expectation, and the true multidimensional post-shock structure, could affect estimates of the clumping of the cool gas and the (eventual) process of dust formation (Derdzinski et al. 2017). Application of our results to these issues will be the aim of future work. Properly capturing the multiphase structure of the post-shock region requires suppressing artificial numerical conduction from the hot to the cold regions (Parkin & Pittard 2010), which is challenging because of the large density contrast involved and the high cooling rate of the cold gas. We mitigate this effect in our simulations by employing a new ‘Lagrangian’ numerical technique that eliminates the advective term in the numerical scheme at the expense of only moderate loss of accuracy (Appendix A). The impact of this technique on preserving the high-temperature phase of the post-shock region is striking (Fig. 6). The average temperature of the radiated thermal emission (Figs 8 and 9), and the fraction of the shock’s thermal luminosity radiated at the maximum temperature ≈Tsh (Fig. 10) is strongly suppressed compared to the naive expectation from a 1D planar compression analysis. The degree of suppression is well approximated by the expression |$L(\gt T_{\rm sh}/3)/L_{\rm tot} \sim 4.5/\mathcal {M}^{4/3}$| for |$4 \lesssim \mathcal {M} \lesssim 36$| (Fig. 11). Our findings broadly agree with those of previous work (e.g. Kee et al. 2014), but we illustrate the Mach number and 2D versus 3D dependence of this effect here for the first time. We find that the temperature suppression only partially results from the lower Mach numbers of the oblique shocks created by the NTSI. A more important effect appears to be the presence of weak shocks behind the main outer shocks, which are being driven into underpressurized cold filaments (presumably created by thermal instability; Figs 12 and 13), robbing the hot medium of its thermal energy before it can be radiated. The suppression appears to be greater in 3D than in 2D, a result which may possibly be attributed to the larger surface area of the cold gas exposed to the hot medium in the 3D case (equation 14). The direction of the upstream magnetic field in the dynamically expanding outflows of astrophysical transients is expected to lie close to the plane of the shock, resulting in the nearly complete suppression of the DSA of relativistic ions (Caprioli & Spitkovsky 2014a). However, the NTSI results in local regions where the obliquity angle is sufficiently large θB ≳ 40°–55° to allow ion acceleration (Figs 14 and 15). Averaged over the shock front, the acceleration efficiency is typically εnt ∼ 1 per cent (Fig. 16), depending only weakly on the Mach number. Though well below the ∼10 per cent maximum predicted for quasi-parallel shocks, this acceleration efficiency agrees remarkably well with those inferred from multiwavelength observations of gamma-ray novae (Metzger et al. 2015; Li et al. 2017). On the other hand, the ion particle energy spectrum we calculate using standard DSA theory is too flat compared to observations, e.g. dN/dE∝E−β with β ≈ 2 (Fig. 17). The preferential escape of higher energy ion from the shock region (before they produce gamma-rays), or an additional component of ion acceleration with a softer spectrum from e.g. weaker shocks within the turbulent region between the outer shocks, could contribute to the steeper spectrum β ≈ 2.7 inferred from the gamma-ray data. Our results do not address the implications of radiative shocks for electron acceleration, for which current PIC simulations suggest may be relatively independent of the upstream magnetic field geometry (e.g. Park, Caprioli & Spitkovsky 2015). Indeed, efficient electron acceleration is inferred indirectly in classical novae from radio synchrotron emission at later times than the observed gamma-ray emission (e.g. Chomiuk et al. 2014; Weston et al. 2016), but when the forward shock may still be in the radiative regime (Metzger et al. 2014; Vlasov et al. 2016). Our results for suppressed X-ray emission have implications for shock-powered astrophysical transients. As one example, many classical novae are accompanied by thermal X-ray emission of luminosity LX ∼ 1033–1035 erg s−1 and temperatures ≳ 1 keV (e.g. Sokoloski et al. 2006; Schwarz et al. 2011; Chomiuk et al. 2014). The fact that these luminosities are much lower than the power of the shocks ≳ 1037 erg s−1 required to explain the GeV gamma-ray emission has generally led to the conclusion that the shocks are giving rise to these phenomenon are distinct (e.g. Vlasov et al. 2016). However, from equation (11), we conclude that for Mach numbers |$\mathcal {M} \sim 30\text{--}200$| of relevance, the X-ray emission will be suppressed by a factor of ∼20–200. If the observed X-rays arise from the same radiative shocks responsible for the gamma-rays, then the true shock power could be several orders of magnitude higher than would be naively guessed from the measured X-rayluminosity. Beyond their application to transients, radiative shocks can play a role in a variety of other astrophysical systems, including colliding wind binaries (e.g. Stevens et al. 1992; Lamberts, Fromang & Dubus 2011; Hendrix et al. 2016), old SN remnants (Raymond et al. 1997; Blondin et al. 1998), young stellar object outflows (e.g. Hansen et al. 2017), and cloud–cloud collisions in active galactic nucleus (Daltabuit & Cox 1972). Our results demonstrate that expectations for the radiation and geometric structure of these phenomenon could require substantial revision when generalizing from expectations based on 1D simulations. We conclude by addressing some limitations of our work. First, though we believe our novel ‘Lagrangian’ technique represents an important advance in eliminating artificial conduction in the post-shock region, we have not demonstrated complete numerical convergence, as resolving the dense cold filamentary structures remains a challenge. Furthermore, the actual physical effect of electron conduction between the hot and cold phases could be relevant in some cases (equation 9), particularly for lower velocity shocks (though magnetic fields could substantially suppressconduction). Another limitation of our work is that we do not include pressure support in the post-shock gas from magnetic fields or non-thermal particles (though the cooling time of even relativistic ions may be short compared to the dynamical time in many cases of relevance; Vurm & Metzger 2018), which could impose a maximum compression ratio lower than that introduced by our temperature floor. Magnetic pressure, by reducing the density contrast between the cooling layers and the inflowing medium (for fields parallel to the shock normal) or by resisting shear (if the field is in the same plane of the shock), can act to decrease the effectiveness of the NTSI (Ramachandran & Smith 2005; Heitsch et al. 2007). On the other hand, ambipolar diffusion or other forms of reconnection may act to eliminate the magnetic field, and its effects on shear may be less important if its coherence length-scale is much smaller than other macroscopic scales in the problem. ACKNOWLEDGEMENTS ES and BDM are supported through the NSF grant AST-1615084, NASA Fermi Guest Investigator Program grants NNX16AR73G and 80NSSC17K0501; and through Hubble Space Telescope Guest Investigator Program grant HST-AR-15041.001-A. We thank Damiano Caprioli, John Raymond, and Lorenzo Sironi for helpful comments and conversations. We acknowledge computing resources from Columbia University’s Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement Grant 1G20RR030893-01, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation Contract C090171, both awarded 2010April 15. Footnotes 1 The maximum unstable wavelength λmax = cst results because the instability cannot grow supersonically, i.e. |$t^{-1}_{\rm NTSI}\lt c_\text{s}/\lambda$|⁠. REFERENCES Ackermann M. et al. , 2014 , Science , 345 , 554 https://doi.org/10.1126/science.1253947 Allen M. G. , Groves B. A. , Dopita M. A. , Sutherland R. S. , Kewley L. J. , 2008 , ApJS , 178 , 20 https://doi.org/10.1086/589652 Bertschinger E. , 1986 , ApJ , 304 , 154 https://doi.org/10.1086/164151 Blandford R. , Eichler D. , 1987 , Phys. Rep. , 154 , 1 https://doi.org/10.1016/0370-1573(87)90134-7 Blondin J. M. , Wright E. B. , Borkowski K. J. , Reynolds S. P. , 1998 , ApJ , 500 , 342 https://doi.org/10.1086/305708 Caprioli D. , Spitkovsky A. , 2014a , ApJ , 783 , 91 https://doi.org/10.1088/0004-637X/783/2/91 Caprioli D. , Spitkovsky A. , 2014b , ApJ , 794 , 47 https://doi.org/10.1088/0004-637X/794/1/47 Cheung C. C. et al. , 2016 , ApJ , 826 , 142 https://doi.org/10.3847/0004-637X/826/2/142 Chevalier R. A. , Imamura J. N. , 1982 , ApJ , 261 , 543 https://doi.org/10.1086/160364 Chevalier R. A. , Irwin C. M. , 2011 , ApJ , 729 , L6 https://doi.org/10.1088/2041-8205/729/1/L6 Chomiuk L. et al. , 2014 , ApJ , 788 , 130 https://doi.org/10.1088/0004-637X/788/2/130 Daltabuit E. , Cox D. , 1972 , ApJ , 173 , L13 https://doi.org/10.1086/180908 Derdzinski A. M. , Metzger B. D. , Lazzati D. , 2017 , MNRAS , 469 , 1314 https://doi.org/10.1093/mnras/stx829 Dopita M. A. , Sutherland R. S. , 1996 , ApJS , 102 , 161 https://doi.org/10.1086/192255 Drake R. P. , 2005 , Ap&SS , 298 , 49 https://doi.org/10.1007/s10509-005-3911-7 Drury L. O. , 1983 , Rep. Prog. Phys. , 46 , 973 https://doi.org/10.1088/0034-4885/46/8/002 Ferland G. J. et al. , 2017 , Rev. Mex. Astron. Astrof. , 53 , 385 Franckowiak A. , Jean P. , Wood M. , Cheung C. C. , Buson S. , 2018 , A&A , 609 , A120 https://doi.org/10.1051/0004-6361/201731516 Gaetz T. J. , Edgar R. J. , Chevalier R. A. , 1988 , ApJ , 329 , 927 https://doi.org/10.1086/166437 Guo F. , Giacalone J. , 2010 , ApJ , 715 , 406 https://doi.org/10.1088/0004-637X/715/1/406 Hansen E. C. , Frank A. , Hartigan P. , Lebedev S. V. , 2017 , ApJ , 837 , 143 https://doi.org/10.3847/1538-4357/aa5ca8 Heitsch F. , Slyz A. D. , Devriendt J. E. G. , Hartmann L. W. , Burkert A. , 2007 , ApJ , 665 , 445 https://doi.org/10.1086/519513 Hendrix T. , Keppens R. , van Marle A. J. , Camps P. , Baes M. , Meliani Z. , 2016 , MNRAS , 460 , 3975 https://doi.org/10.1093/mnras/stw1289 Innes D. E. , 1992 , A&A , 256 , 660 Ivanova N. , Justham S. , Avendano Nandez J. L. , Lombardi J. C. , 2013 , Science , 339 , 433 https://doi.org/10.1126/science.1225540 Katz B. , Budnik R. , Waxman E. , 2010 , ApJ , 716 , 781 https://doi.org/10.1088/0004-637X/716/1/781 Kee N. D. , Owocki S. , ud-Doula A. , 2014 , MNRAS , 438 , 3557 https://doi.org/10.1093/mnras/stt2475 Lamberts A. , Fromang S. , Dubus G. , 2011 , MNRAS , 418 , 2618 https://doi.org/10.1111/j.1365-2966.2011.19653.x Li K.-L. et al. , 2017 , Nat. Astron. , 1 , 697 https://doi.org/10.1038/s41550-017-0222-1 Mac Low M.-M. , Norman M. L. , 1993 , ApJ , 407 , 207 https://doi.org/10.1086/172506 MacLeod M. , Ostriker E. C. , Stone J. M. , 2018 , preprint (arXiv:1803.03261) Malkov M. A. , Drury L. O. , 2001 , Rep. Prog. Phys. , 64 , 429 https://doi.org/10.1088/0034-4885/64/4/201 Martin P. , Dubus G. , 2013 , A&A , 551 , A37 https://doi.org/10.1051/0004-6361/201220289 Martin P. , Dubus G. , Jean P. , Tatischeff V. , Dosne C. , 2017 , A&A , 612 , A38 McLeod A. D. , Whitworth A. P. , 2013 , MNRAS , 431 , 710 https://doi.org/10.1093/mnras/stt203 Metzger B. D. , 2010 , MNRAS , 409 , 284 https://doi.org/10.1111/j.1365-2966.2010.17308.x Metzger B. D. , Pejcha O. , 2017 , MNRAS , 471 , 3200 https://doi.org/10.1093/mnras/stx1768 Metzger B. D. , Hascoët R. , Vurm I. , Beloborodov A. M. , Chomiuk L. , Sokoloski J. L. , Nelson T. , 2014 , MNRAS , 442 , 713 https://doi.org/10.1093/mnras/stu844 Metzger B. D. , Finzell T. , Vurm I. , Hascoët R. , Beloborodov A. M. , Chomiuk L. , 2015 , MNRAS , 450 , 2739 https://doi.org/10.1093/mnras/stv742 Metzger B. D. , Caprioli D. , Vurm I. , Beloborodov A. M. , Bartos I. , Vlasov A. , 2016 , MNRAS , 457 , 1786 https://doi.org/10.1093/mnras/stw123 Mukai K. , Ishida M. , 2001 , ApJ , 551 , 1024 https://doi.org/10.1086/320220 Narayan R. , Medvedev M. V. , 2001 , ApJ , 562 , L129 https://doi.org/10.1086/338325 Nymark T. K. , Fransson C. , Kozma C. , 2006 , A&A , 449 , 171 https://doi.org/10.1051/0004-6361:20054169 Park J. , Caprioli D. , Spitkovsky A. , 2015 , Phys. Rev. Lett. , 114 , 085003 https://doi.org/10.1103/PhysRevLett.114.085003 Parkin E. R. , Pittard J. M. , 2010 , MNRAS , 406 , 2373 https://doi.org/10.1111/j.1365-2966.2010.16888.x Parkin E. R. , Pittard J. M. , Corcoran M. F. , Hamaguchi K. , 2011 , ApJ , 726 , 105 https://doi.org/10.1088/0004-637X/726/2/105 Pejcha O. , Metzger B. D. , Tyles J. G. , Tomida K. , 2017 , ApJ , 850 , 59 https://doi.org/10.3847/1538-4357/aa95b9 Pittard J. M. , 2009 , MNRAS , 396 , 1743 https://doi.org/10.1111/j.1365-2966.2009.14857.x Ramachandran B. , Smith M. D. , 2005 , MNRAS , 362 , 1353 https://doi.org/10.1111/j.1365-2966.2005.09406.x Raymond J. C. , 1979 , ApJS , 39 , 1 https://doi.org/10.1086/190562 Raymond J. C. , Blair W. P. , Long K. S. , Vancura O. , Edgar R. J. , Morse J. , Hartigan P. , Sanders W. T. , 1997 , ApJ , 482 , 881 https://doi.org/10.1086/304183 Schaal K. , Springel V. , 2015 , MNRAS , 446 , 3992 https://doi.org/10.1093/mnras/stu2386 Schaal K. et al. , 2016 , MNRAS , 461 , 4441 https://doi.org/10.1093/mnras/stw1587 Schwarz G. J. et al. , 2011 , ApJS , 197 , 31 https://doi.org/10.1088/0067-0049/197/2/31 Smith N. , 2016 , preprint (arXiv:1612.02006) Smith N. , McCray R. , 2007 , ApJ , 671 , L17 https://doi.org/10.1086/524681 Sokoloski J. L. , Luna G. J. M. , Mukai K. , Kenyon S. J. , 2006 , Nature , 442 , 276 https://doi.org/10.1038/nature04893 Steinberg E. , Yalinewich A. , Sari R. , 2016 , MNRAS , 459 , 1596 https://doi.org/10.1093/mnras/stw783 Stevens I. R. , Blondin J. M. , Pollock A. M. T. , 1992 , ApJ , 386 , 265 https://doi.org/10.1086/171013 Strickland R. , Blondin J. M. , 1995 , ApJ , 449 , 727 https://doi.org/10.1086/176093 Townsend R. H. D. , 2009 , ApJS , 181 , 391 https://doi.org/10.1088/0067-0049/181/2/391 Vishniac E. T. , 1983 , ApJ , 274 , 152 https://doi.org/10.1086/161433 Vishniac E. T. , 1994 , ApJ , 428 , 186 https://doi.org/10.1086/174231 Vlasov A. , Vurm I. , Metzger B. D. , 2016 , MNRAS , 463 , 394 https://doi.org/10.1093/mnras/stw1949 Vurm I. , Metzger B. D. , 2018 , ApJ , 852 , 62 https://doi.org/10.3847/1538-4357/aa9c4a Walder R. , Folini D. , 2000 , Ap&SS , 274 , 343 https://doi.org/10.1023/A:1026597318472 Weston J. H. S. et al. , 2016 , MNRAS , 457 , 887 https://doi.org/10.1093/mnras/stv3019 Yalinewich A. , Steinberg E. , Sari R. , 2015 , ApJS , 216 , 35 https://doi.org/10.1088/0067-0049/216/2/35 APPENDIX A: LAGRANGIAN HYDRODYNAMICS ON VORONOI MOVING MESH A1 New ‘Lagrangian’ Scheme While Voronoi-based moving-mesh codes greatly reduce the advection terms in the Euler equations, they are still present to some extent, resulting in a finite amount of mass transfer between cells. We provide below an adaption of this scheme which, at the expense of some loss of accuracy, is fully Lagrangian. We start by calculating the wave speed of the discontinuity, and solve the Riemann problem at that reference frame. This naturally gives a flux that has no advection terms. However, since the interface between cells moves at a different velocity, an error occurs. In order to minimize this error, we fix the energy flux to compensate for the PdV work that was actually done. In practice, this means that the energy flux is set to be P*|$v$|intr instead of P*|$v$|*, where P* and |$v$|* are the pressure and velocity at the discontinuity and |$v$|intr is the velocity of the interface between cells. Since for smooth flows the difference between |$v$|intr and |$v$|* should be inversely proportional to the resolution, we expect a first-order convergence, while still retaining the conservative character of the scheme. There are two caveats that need to be addressed. First, a scenario in which ||$v$|cell| > ||$v$|intr| can occur when there is a shock wave, where |$v$|cell is the velocity of the downstream cell. Since this leads to an unphysical cooling of the shocked cell, we set the energy flux to be P*max(⁠|$v$|intr, |$v$|cell) whenever a shock wave occurs. Secondly, it is possible for the sign of |$v$|intr and |$v$|* to be opposite. This can give rise to negative thermal energy in the upstream cell. In this case, we set the energy flux to be min(P*, Pup)|$v$|intr, where Pup is the pressure of the upstream cell. Since this scheme is handled on an interface basis, in this paper we limit its use to be only on interfaces which have a temperature ratio larger than 3 between the left and the right states. However, for the examples shown below, we employ this scheme on all the interfaces. A2 Lagrangian mesh steering In order for the error from the Lagrangian scheme be minimized, the Riemann problems should be solved at interfaces moving as closely as possible to the contact speed. We propose the following scheme of setting the velocities of the mesh generating points Calculate the contact speed on the jth interface, |$\boldsymbol {w}_j$|⁠. For each Voronoi cell, we calculate its predicted new centre of mass \begin{eqnarray*} \boldsymbol {{\rm CM}}_{\rm new} &=& \left[V\cdot \boldsymbol {{\rm CM}}_{\rm old}+\sum A_j|\boldsymbol {w}_j|\text{d}t\left(\boldsymbol {F}_j + 0.5\boldsymbol {w}_j \text{d}t\right)\right]\nonumber \\ &&\left(V+\sum A_j{\rm sgn}(\boldsymbol {w}_j) |\boldsymbol {w}_j|\text{d}t\right), \end{eqnarray*} (A1) where V is the volume of the cell, |$\boldsymbol {F}_j$| is the centroid of the jth interface, Aj is the area of the jth interface, dt is the time-step, and the summation is done over all of the Voronoi cell’s interfaces. Set the velocity of the mesh generating point to be \begin{eqnarray*} \boldsymbol {v}_{\rm mesh} = \left(\boldsymbol {{\rm CM}}_{\rm new}-\boldsymbol {{\rm CM}}_{\rm old}\right)/\text{d}t. \end{eqnarray*} (A2) Additionally, for each cell, we calculate its deviation from Lagrangian motion, i.e. |$\sum A_j\left(\boldsymbol {w}_j-\boldsymbol {v}_{{\rm intr},j}\right)\text{d}t$|⁠, and at the following time-step add an additional velocity to compensate for this deviation that is limited in its magnitude to be a fraction of the mesh generating point’s velocity. A3 Convergence for the Yee vortex We benchmark the performance of our new scheme for smooth flows by running the Yee vortex test problem on a randomly distributed mesh. This problem describes the evolution of isentropic vortex that is pressure supported. The set-up is identical to that described in Steinberg, Yalinewich & Sari (2016) and is run until time t = 10. We calculate the L1 error norm of the density and the angular momentum, defined as |$L_1=\frac{1}{V}\int |\phi (x)-\phi (x)_Y|\text{d}V$|⁠, where ϕ(x)Y is the analytical time-independent variable and ϕ(x) is the computed variable from the simulation. Fig. A1 shows the L1 error norms, demonstrating convergence of both the density and angular momentum, although at a somewhat lower rate of −0.75 compared to the expected first-order convergence. This is likely due to the fact that we add an additional small non-Lagrangian velocity to the mesh generating points in order to keep the Voronoi cells ‘round’. Figure A1. Open in new tabDownload slide The L1 error norm for the density and angular momentum as a function of numerical resolution, for the implementation of our ‘Lagrangian’ scheme to the Yee vortex test problem. Figure A1. Open in new tabDownload slide The L1 error norm for the density and angular momentum as a function of numerical resolution, for the implementation of our ‘Lagrangian’ scheme to the Yee vortex test problem. A4 Kelvin–Helmholtz To test our new scheme’s performance for discontinuities, we redo the Kelvin–Helmholtz test problem as presented in Yalinewich et al. (2015). We run the problem until time t = 1.5 using both the ‘Lagrangian’ scheme as well as the normal ‘Semi-Lagrangian’ scheme, both with a resolution of 1282. Fig. A2 shows the density at the completion of the run. While the ‘Semi-Lagrangian’ method produces a smoother result than our new ‘Lagrangian’ method, the discontinuity is nevertheless perfectly preserved in the ‘Lagrangian’ method and the system maintains thesamedynamicalevolution. Figure A2. Open in new tabDownload slide Kelvin–Helmholtz test problem, comparing our ‘Semi-Lagrangian’ (top panel) and ‘Lagrangian’ (bottom panel) schemes. The colour level denotes density. Figure A2. Open in new tabDownload slide Kelvin–Helmholtz test problem, comparing our ‘Semi-Lagrangian’ (top panel) and ‘Lagrangian’ (bottom panel) schemes. The colour level denotes density. A5 Sedov blast wave The Sedov blast wave problem provides an ideal way to benchmark how our new scheme handles strong shocks. We set-up an initial grid with 2500 randomly distributed mesh generating points with an initial uniform density and a hotspot at the centre. Fig. A3 shows density snapshots following some period of evolution, comparing the Semi-Lagrangian and the Lagrangian schemes. While the solution for the Semi-Lagrangian scheme is smoother, the Lagrangian scheme does a good job of catching the overall morphology of the blast wave and the shock location. This is shown further in Fig. A4, which compares the radial density profiles between the schemes. The peak density achieved in Lagrangian case is about 25 per cent higher than in the Semi-Lagrangian case, but the overall shape is similar. Figure A3. Open in new tabDownload slide Sedov blastwave test problem, comparing our ‘Semi-Lagrangian’ (top panel) and ‘Lagrangian’ (bottom panel) schemes. The colour level denotes density. Figure A3. Open in new tabDownload slide Sedov blastwave test problem, comparing our ‘Semi-Lagrangian’ (top panel) and ‘Lagrangian’ (bottom panel) schemes. The colour level denotes density. Figure A4. Open in new tabDownload slide Snapshot of the radial density profiles of the Sedov blastwave problem, comparing the results of the Semi-Lagrangian and Lagrangian schemes. Figure A4. Open in new tabDownload slide Snapshot of the radial density profiles of the Sedov blastwave problem, comparing the results of the Semi-Lagrangian and Lagrangian schemes. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)
TI - The multidimensional structure of radiative shocks: suppressed thermal X-rays and relativistic ion acceleration
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty1641
DA - 2018-09-01
UR - https://www.deepdyve.com/lp/oxford-university-press/the-multidimensional-structure-of-radiative-shocks-suppressed-thermal-ufDD8ldDcf
SP - 687
VL - 479
IS - 1
DP - DeepDyve
ER -