Internal Shocks from Variable Outflows in Classical Novae

TL;DR

This paper uses hydrodynamical simulations to study internal shocks in classical novae outflows, explaining multi-wavelength emissions and complex light curve features through variable outflow velocities and shock interactions.

Contribution

It introduces a model of non-monotonic, fluctuating outflow velocities causing internal shocks, which better explains observed nova emissions and light curve complexities.

Findings

Simulations reproduce correlated gamma-ray and optical flares.

Delayed X-ray and radio emissions are linked to shock evolution.

Multiple internal shocks explain complex spectral line behaviors.

Abstract

We present one-dimensional hydrodynamical simulations including radiative losses, of internal shocks in the outflows from classical novae, to explore the role of shocks in powering multi-wavelength emission from radio to gamma-ray wavelengths. Observations support a picture in which the initial phases of some novae generate a slow, equatorially-focused outflow (directly from the outer Lagrange point, or from a circumbinary disk), which then transitions to, or is overtaken by, a faster more isotropic outflow from the white dwarf which collides and shocks the slower flow, powering gamma-ray and optical emission through reprocessing by the ejecta. However, the common occurence of multiple peaks in nova light curves suggests that the outflow's acceleration need not be monotonic, but instead can involve successive transitions between "fast" and "slow" modes. Such a time-fluctuating outflow velocity naturally can reproduce several observed properties of nova, such as correlated gamma-ray and optical flares, expansion of the photosphere coincident with (though lagging slightly) the peak flare luminosity, and complex time-evolution of spectral lines (including accelerating, decelerating, and merging velocity components). While the shocks are still deeply embedded during the gamma-ray emission, the onset of $\sim$ keV X-ray and $\sim 10$ GHz radio synchrotron emission is typically delayed until the forward shock of the outermost monolithic shell (created by merger of multiple internal shock-generated shells) reaches a sufficiently low column through the dense external medium generated by the earliest phase of the outburst.

Equations

$$\dot{M}={\dot{M}}_{\mathrm{f}}=10^{-5}{M}_{\odot }\mathrm{wee}{\mathrm{k}}^{-1},v_{\mathrm{w}}=v_{\mathrm{f}}=1400\mathrm{km}{\mathrm{s}}^{-1}$$

$$\dot{M}={\dot{M}}_{\mathrm{s}}=5\cdot 10^{-6}{M}_{\odot }\mathrm{wee}{\mathrm{k}}^{-1},v_{\mathrm{w}}=v_{\mathrm{s}}=200\mathrm{km}{\mathrm{s}}^{-1}.$$

$$\begin{array}{rl}\delta t& \lesssim t_{\mathrm{thin}}\approx \frac{\dot{M}{\kappa }_{\mathrm{opt}}}{4\pi v_{\mathrm{f}}{v}_{\mathrm{s}}}\\ \approx & 1\mathrm{day}\left(\frac{\dot{M}}{10^{-4}{M}_{\odot }\mathrm{w}{\mathrm{k}}^{-1}}\right)\left(\frac{1500\mathrm{km}{\mathrm{s}}^{-1}}{v_{\mathrm{f}}}\right)\left(\frac{500\mathrm{km}{\mathrm{s}}^{-1}}{v_{\mathrm{s}}}\right),\end{array}$$

$$T_{\mathrm{eff}}={\left(\frac{L_{\mathrm{tot}}}{4\pi \sigma R_{\mathrm{ph}}^2}\right)}^{1/4}.$$

$$T_{\mathrm{B},\nu }^{\mathrm{obs}}(t)=T_{\mathrm{B},\nu }(t)\exp [-{\tau }_{\mathrm{r}}(t)],$$

$$T_{\mathrm{B},\nu }\simeq 4.2\cdot 10^5\mathrm{K}{\left(\frac{n}{10^7\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{\frac{1+p}{4}}{\left(\frac{\nu }{10\mathrm{GHz}}\right)}^{-\frac{p+3}{2}}$$

$$\times {\left(\frac{v_{\mathrm{sh}}}{10^3\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{\frac{11+p}{2}}{\left(\frac{{ϵ}_B}{10^{-2}}\right)}^{\frac{p+1}{4}}\left(\frac{{ϵ}_e}{10^{-2}}\right){\left(\frac{\Lambda (T_{\mathrm{sh}})}{10^{-22}\mathrm{erg}\mathrm{c}{\mathrm{m}}^3{\mathrm{s}}^{-1}}\right)}^{-1},$$

$${\tau }_{\mathrm{r}}={\int }_r^{\infty }{\alpha }_{\mathrm{r}}dr;{\alpha }_{\mathrm{r}}\approx 0.06T^{-3/2}{n}^2{\nu }^{-2}\mathrm{c}{\mathrm{m}}^{-1},$$

$$\nu L_{\nu }=4{\pi }^2{R}_{\mathrm{sh}}^2\frac{2{\nu }^{3}k{T}_{\mathrm{B},\nu }^{\mathrm{obs}}}{c^2},$$

$$r_{\mathrm{dust}}={\left(\frac{L}{4\pi \sigma T_{\mathrm{cond}}^4}\right)}^{1/2},$$

$${\rho }_{\mathrm{w}}$$

$$p_{\text{ratio}}$$

$$T_{\text{ratio}}$$

$${\rho }_{\text{ratio}}$$

$$Q_{\text{h,WD}}=L_{\text{WD}}/13.6\text{eV}{\mathrm{s}}^{-1},$$

$$Q_{\mathrm{h},\mathrm{net}}=\max \left(Q_{\text{h,WD}}-{\int }_0^{r_{\mathrm{f}}}{\alpha }_{\mathrm{B}}{\mathrm{n}}_{\mathrm{e}}{\mathrm{n}}_{\mathrm{i}}4\pi {\mathrm{r}}^2\mathrm{dr},0\right),$$

$$\nu L_{\nu }=4{\pi }^2{R}_{\mathrm{ph},\nu }^2\frac{2{\nu }^{3}k{T}^{\mathrm{ion}}}{c^2}+{\int }_{r(\tau =1)}^{\infty }{\epsilon }_{ff}4\pi r^{2}dr,$$

Full text

††pubyear: 2015††pagerange: Internal Shocks from Variable Outflows in Classical Novae–B

1 Introduction

A classical nova occurs when a shell of hydrogen-rich material on the surface of a white dwarf, which is accreting from a main sequence binary companion, undergoes runaway nuclear burning (Gallagher & Starrfield, 1978; Prialnik, 1986; Bode & Evans, 2008; Yaron et al., 2005; Starrfield et al., 2016). The resulting energy release leads to the ejection of a shell of mass $\sim 10^{-3}-10^{-7}M_{\odot}$ at high velocities $\sim 500-5000$ km s*-1* over days to months or longer. The ejecta contains the original accreted material and its burning products, as well as heavier nuclei dredged up from the underlying white dwarf (e.g. by Kelvin-Helmholtz instabilities at the interface between the inert core and the burning envelope; e.g. Casanova et al. 2018).

Traditionally, models of novae predicted that their light curves are powered exclusively by radiation which diffuses outwards directly from the hot white dwarf envelope through the expanding ejecta (Gallagher & Starrfield, 1978; Yaron et al., 2005). Early in the outburst, the ejecta are optically thick, and this luminosity emerges primarily at visual wavelengths. As the ejecta becomes increasingly transparent with time, the spectral peak shifts into the UV and, ultimately, the soft X-rays, once the white dwarf surface becomes visible (e.g. Schwarz et al. 2011; Henze et al. 2014), powered by stable burning of the residual fuel (e.g. Wolf et al. 2013).

However, a number of complications have long afflicted this simple picture. The optical light curves of many nova show multiple peaks (e.g. Strope et al. 2010) and complex spectral components, suggestive of a number of discrete mass ejection events (e.g. Cassatella et al. 2004; Csák et al. 2005; Hillman et al. 2014). The photosphere appears to temporarily expand during flares (e.g. Tanaka et al. 2011a; Aydi et al. 2019), behavior often accompanied by the appearance of new absorption line systems (e.g. Jack et al. 2017). Furthermore, X-ray observations reveal the presence of hot $\gtrsim 10^{7}$ K gas indicative of shocks (O’Brien et al., 1994; Mukai & Ishida, 2001; Nelson et al., 2019). However, the prevalance, and energetic importance, of these shocks were not fully appreciated until the discovery by Fermi/LAT of $\gtrsim 100$ MeV gamma-rays from Galactic classical novae (Ackermann et al. 2014; Cheung et al. 2016; Franckowiak et al. 2018).

The gamma-rays observed from classical novae are very likely produced by relativistic particles accelerated at shocks internal to the ejecta (e.g. Martin & Dubus 2013; Metzger et al. 2015; Martin et al. 2018). The high columns of gas ahead and behind the shocks make nova outflows excellent “calorimeters” for converting cosmic ray energy into gamma-rays (Metzger et al., 2015). For typical acceleration efficiencies of at most tens of percent of the shock power going into the energy of relativistic particles (e.g. Caprioli & Spitkovsky 2014), the high observed gamma-ray luminosities, $L_{\gamma}\sim 10^{35}-10^{36}$ erg s*-1*, require shocks with kinetic powers approaching the bolometric output of the nova, typically at or near the Eddington luminosity of the white dwarf, $L_{\rm Edd}\sim 10^{38}$ erg s*-1*.

The high densities in nova outflows also imply that the hot shocked gas will cool rapidly compared to the expansion time, i.e. the shocks are radiative rather than adiabatic. The layer of gas behind a radiative shock in this temperature range are often subject to thermal instabilities (Chevalier & Imamura 1982) as well as dynamical instabilities in the thin shell generated by the gas’ compression (Vishniac 1994). In conflict with naive expectations based on the 1D jump conditions (which would predict mainly hard $\gtrsim$ keV X-ray emission for $\sim 10^{3}$ km s*-1* shocks), recent multi-dimensional hydrodynamical simulations of radiative shocks show that these instabilities can reduce the temperature of the bulk of the luminosity to soft X-ray or UV wavelengths (Kee et al., 2014; Steinberg & Metzger, 2018). However, as we now discuss, such direct emission from shock is anyways not typically visible to an external observer until late times due to absorption by the overlying ejecta.

Early in the nova outburst, when the mass loss rate from the white dwarf is greatest, any UV and X-ray emission from the shocks will be absorbed by gas ahead of the shocks, reprocessing its energy into optical line or continuum emission. In fact, the required shock powers are so high that it led Metzger et al. (2014) to conclude that a significant fraction of the total optical luminosity of the nova may be powered (indirectly) by shocks (see also Metzger et al. 2015; Martin et al. 2018) rather than the white dwarf. This possibility has received strong support from the observed temporal correlations between the optical and gamma-ray light curves (Li et al. 2017; Aydi et al., in prep). High gas densities behind the shocks also imply that primary and secondary relativistic leptons (electrons or positrons) accelerated at the shock transfer their energy through Coulomb scattering to the thermal plasma faster than it can be radiated (e.g. via Inverse Compton scattering). This suppresses the non-thermal X-ray emission relative to what one would expect based on a naive downward extrapolation of the $\sim$ GeV photon spectrum measured by Fermi (Vurm & Metzger, 2018), consistent with deep upper limits from NuSTAR on a component of non-thermal X-rays contemporaneous with Fermi detections (Nelson et al., 2019).

Further evidence regarding the nature of the shocks in novae comes from late-time ($\gtrsim$ several month to year timescale) radio observations, both of thermal emission from the expanding photo-ionized ejecta (e.g. Ribeiro et al. 2014a; Cunningham et al. 2015) as well non-thermal synchrotron emission from shocks (e.g. Krauss et al. 2011; Weston et al. 2016). Resolved radio imaging supports a picture in which the initial stages of the nova eruption generate a slow outflow that is geometrically-concentrated in the binary equator, which then transitions to$-$or is overtaken by$-$a faster outflow or wind with more isotropic angular distribution (Sokoloski et al., 2008; Chomiuk et al., 2014; Linford et al., 2015). This outflow sequence is consistent with the bipolar ring- and jet-like structures revealed by optical imaging and spectroscopy of nova shells (e.g. Gill & O’Brien 2000; Harman & O’Brien 2003; Woudt et al. 2009; Shore 2013; Ribeiro et al. 2013; Shore et al. 2013). The “slow" equatorial outflow plausibly arises from an outflow or circumbinary disk which is fed by Roche lobe overflow at velocities $\lesssim$ few 100 km s*-1* comparable to the binary orbital speed. In contrast, the “fast" ejecta likely originates directly from the white dwarf surface, e.g. through a momentum-driven wind (Kato & Hachisu 1994).

Collisions between “fast” and “slow” outflow components thus provide a natural site for shocks and bipolar collimation in novae. However, the physical mechanism giving rise to the transitions between these markedly distinct outflow modes remains unclear. Most previous theoretical work has considered the evolution of the outflow velocity from slow to fast to be monotonic in time (e.g. Metzger et al. 2014; Martin et al. 2018). However, the presence of multiple flares in nova optical and gamma-ray light curves, if interpreted as distinct collision events, requires multiple transitions between fast and slow modes in a single event (e.g.  Kato & Hachisu 2009). Rather than a single shock, such a scenario would generate an entire ‘train’ of collisions. At late times, both single- and multiple-transition scenarios will generate a broadly similar global hourglass geometry of a fast bipolar wind confined by the “waist" of the slow equatorial shell of shock-compressed gas (Fig. 1). However, at early times (e.g. during the observed gamma-ray flares), these scenarios make distinct predictions for the physical locations of the shocks and their impact on the observed emission.

This work explores the role of internal shocks in nova outflows and their emission using global one-dimensional hydrodynamical simulations which include radiative losses. 1D calculations cannot account for latitudinal or azimuthal variations in the shocks, nor can they resolve important multi-dimensional features of the shock front. The latter can be approximately accounted for in a sub-grid sense using results gained from previous local multi-dimensional simulations of radiative shocks (Steinberg & Metzger, 2018). In this paper we simply point out where multi-dimensional results are likely to have the greatest affects on our conclusions and thus where some caution in the interpretation is warranted.

This paper is organized as followed. In §2 we describe details of our numerical model (see also Appendix A). In §3 we present results of the simulations. We compare the single- and multiple-transition scenarios, focusing on their abilities to reproduce the observed variability in nova optical light curves as well as the evolution of the photospheric radius and line velocity evolution. We also present $\gamma-$ray, X-ray, and radio light curves from our models, and briefly discuss implications for dust formation. In §4 we summarize our conclusions.

2 Description of Numerical Model

We perform one-dimensional hydrodynamical simulations of nova outflows using the moving-mesh hydrodynamical code RICH on a spherical grid (Yalinewich et al., 2015). Our aim is to capture the shock interaction between fast and slow outflow components, which we expect to occur primarily in the binary equatorial plane (Fig. 1). We are therefore not accounting for additional shocks or emission which could potentially take place in the polar regions due to velocity fluctuations in the fast outflow itself. The latter could readily be studied in an approximate manner with a similar 1D calculation to that presented here in future work.

As discussed in greater detail below, we specify the time-dependent velocity and mass-loss rate from the central spherical source (white dwarf atmosphere or circumbinary disk) via an inner boundary condition on the radial velocity and density at the radius $r=5\times 10^{11}$ cm. We employ a radial grid extending to an outer radius of $2.5\times 10^{15}$ cm with a high enough resolution to capture the large compression behind the shocks. We assume ideal gas pressure (adiabatic index $\gamma=5/3$) and include optically-thin radiative losses, using a cooling function calculated assuming collisional ionization equilibrium (CIE) as tabulated in version C17.01 of CLOUDY (Ferland et al., 1998, 2017). In our fiducial model we adopt the GASS10 solar abundances, but we also explore models using CLOUDY’s default nova abundances to account for the effects of possible enrichment of the ejecta from nuclear burning (Starrfield et al. 1972).

It is reasonable to treat the cooling as being effectively optically-thin for purposes of capturing the dynamical evolution. This is because the timescale for radial diffusion of the reprocessed optical radiation is much shorter than the dynamical or expansion timescale of the outflow, both on global scales and on smaller scales across the post-shock cooling layer (e.g. Metzger et al. 2014). Since our simulations do not consider heating of the gas by the outwards diffusion of thermal radiation from the white dwarf or shocks, we place a minimum temperature floor of $T_{\rm min}=10^{4}$ K similar to that maintained by the radiation through the layers of greatest interest. Our results are not sensitive to the precise value of the floor, as its main effect is limit the maximum compression ratio of the shocked gas. Additional details of our numerical set-up are described in Appendix A.

Mass loss is a key ingredient in novae, which plays an important role in setting the timescale of the eruption by removing mass from the white dwarf envelope (Livio et al., 1990; Kato & Hachisu, 1994). However, the exact mechanisms giving rise to the unbound outflows, as needed to produce variable internal shocks, are uncertain and not addressed in this work. Relatively slow outflows, of up to several hundreds of km s*-1*, may be launched from near the binary $L_{2}$ Lagrange point. These could be driven for instance by frictional drag (e.g. Livio et al. 1990; Shankar et al. 1991; Lloyd et al. 1997; however, see Kato & Hachisu 1991), non-axisymmetric gravitational torques from the binary (Pejcha et al., 2016b), or in the outflow of a circumbinary disk formed from matter ejected from $L_{2}$ that initially remains gravitationally bound (Pejcha et al., 2016a). These slow outflows, with velocities of order the binary escape speed, would plausibly be concentrated in the binary equatorial plane.

Outflows can also occur from the white dwarf directly, with little or no influence from the binary. These could be driven by radiation pressure acting on the iron opacity bump (e.g. Kato & Hachisu 1994), or by any large (near- or super-Eddington) localized energy deposition in the white dwarf envelope, e.g. from radioactive heating or damping of convectively-excited waves (e.g. Quataert et al. 2016). Outflows launched from deeper locations in the potential well of the white dwarf would be expected to possess higher velocities $\sim 2000-5000$ km s*-1* and a more spherical geometry than the slower outflows from $L_{2}$.

Lacking a first principles model, we are motivated on phenomenological grounds to consider two outflow modes characterized by distinct mass-loss rates, $\dot{M}$, and asymptotic velocities, $v_{\rm w}$. Somewhat arbitrarily (as the precise values will vary between novae), we assume a “fast” mode outflow with

[TABLE]

and a “slow” outflow mode with

[TABLE]

These are spherically-equivalent mass-loss rates. To obtain the true value, $\dot{M}$ must be multipled by the fraction of the total solid angle covered by each outflow component, $f_{\Omega}$. For instance, $f_{\Omega}\approx 1$ for the fast spherical outflow, while $f_{\Omega}\approx 0.3$ for the slow equatorially-focused outflow.111This is consistent with $\sim 30\%$ of nova showing evidence for dust creation along the line of site assuming the dust forms most effectively in the shock-compressed slow equatorial ejecta (Derdzinski et al. 2017; see $\S\ref{sec:dust}$).

The top panel of Figure 2 shows the chosen time evolution of $\dot{M}$ and $v_{\rm w}$ for our two fiducial models (“single transition” and “multiple transition”). In the single transition model (red lines in Fig. 2) the outflow velocity is increased from slow ($v_{\rm w}=v_{\rm s}$) to fast ($v_{\rm w}=v_{\rm f}$) only once at $t\approx 0$, while in the multiple transition model (blue lines in Fig. 2) the outflow velocity fluctuates between slow and fast modes several times. Local maxima (“flares”) are observed in classical nova light curves on timescales ranging from hours to weeks (e.g. Strope et al. 2010; Walter et al. 2012; Henze et al. 2018; Aydi et al. 2019). To generate collisions on approximately this timescale, we impose transitions between the slow to the fast mode on a characteristic time interval $\delta t\approx$ hours$-$days and with delays between the transitions of a day to a week. In our multiple transition model, we space the transitions across a range of $\delta t$, in order to explore its effect on the behavior of the photosphere during flares. In both models, the outflow velocity $v_{w}$ returns back to the slow mode at late times, which we further taper along with $\dot{M}$ to mimic the end of the outburst.

We assume that by the start of our simulations ($t=0$) the initial slow outflow has already reached a characteristic radius $r_{0}$, outside of which we truncate the density as $\rho\propto r^{-4}$ (Fig. 2, bottom panel). Early spectra of novae showing narrow absorption line features indicative of the presence of slowly-expanding material already present around the binary by the time of the nova outburst (e.g. Williams et al. 2008). The origin of this material is unclear: it may represent the earliest stages of mass ejection after the thermonuclear runaway, or, more speculatively, pre-existing material initially orbiting the binary. In our “Fiducial" model we take $r_{0}=6\times 10^{13}$ cm for the radial extent of the pre-existing medium, in which case the total mass of the ejecta in the fast and slow components (assuming $f_{\Omega}=1$ for both) are $\approx 1.5\times 10^{-5}M_{\odot}$ and $\approx 5.5\times 10^{-5}M_{\odot}$, consistent with the measured range in classical novae (e.g. Roy et al. 2012). We also consider a “Low Mass” model with $r_{0}=2\times 10^{13}$ cm and thus a lower mass in the slow component by a factor of a few compared to the Fiducial case. As we will show later, the X-ray and radio light curves of the shocks are sensitive to the properties of the initial external medium. Finally, we consider a model (“Bumpy") in which the initial slow outflow extends to $r_{0}=6\times 10^{13}$ cm but is radially inhomogeneous.

Nova outflows are intrinsically multi-dimensional, both on large spatial scales (e.g. equatorial-polar asymmetry) and on small scales (e.g. the cooling length of the post-shock gas) due to instabilities of the radiative shocks. Using zoomed-in high-resolution 2D and 3D simulations of head-on colliding dense flows, Steinberg & Metzger (2018) demonstrate that the luminosity-weighted temperature of the emission from the gas behind high Mach number radiative shocks is substantially lower than would be guessed from the adiabatic Rankine-Hugoniot jump conditions. This is due to a combination of obliquities generated at the shock front by the non-linear thin-shell instability (Vishniac, 1994) and the transfer of thermal energy from the hot immediate post-shock gas to cool dense clumps downstream (which radiates energy away more efficiently). The present 1D treatment of nova outflows will therefore underestimate the thickness of the cool shells of shock-compressed gas and lead to overestimates of the emission temperature (see $\S\ref{sec:Xrays}$). The former deficiency will not alter our qualitative conclusions, but can lead to artificially high luminosity when two such thin-shells collide ($\S\ref{sec:luminosity}$).

3 Internal Shocks in Novae

In this section we present the results of our simulations and our predictions for the optical, gamma-ray, X-ray, and radio emission in comparison to observed nova properties.

3.1 Single Versus Multiple Transition Scenarios

We begin with a general discussion of the dynamics of the outflow self-interaction and the resulting optical light curves. In the single transition scenario, the nova outflow undergoes a single change from a slow to fast outflow at $t\approx 0$ (Fig. 2). The fast flow quickly catches up to the slow flow, generating a single forward-reverse shock structure that propagates outwards through the slow outflow. The shocks sweep up mass from both flows into a shell that decelerates in a momentum-conserving fashion (Metzger et al., 2014). The shell is extremely thin due to radiative losses; following cooling, the gas density increases by a factor of $\propto\mathcal{M}^{2}\sim 10^{4}$ compared to its original unshocked value, where $\mathcal{M}=v_{\rm sh}/c_{\rm s}$ is the Mach number of the shocks, $v_{\rm sh}$ is the shock speed, and $c_{\rm s}\sim 10$ km s*-1* is the sound speed of the upstream.

By contrast, in the multiple transition model each new flip back to the fast outflow mode generates a new forward/reverse shock structure sandwiched around its own cool shell, which continues propagating out into the slow ejecta (see Figure 3 for snapshots of the radial profiles of the ejecta properties). These thin shells themselves collide at larger radii, creating a second (or third) generation of shocks, and ultimately merging pairwise into a single shell. In all models, the final state is thus similar: a monolithic dense shell bounded in front by a forward shock propagating down the density gradient of the outermost layers of the initial slow upstream medium. Depending on the composition and density profile of the latter (bottom panel of Fig. 2), the forward shock can remain radiative$-$with a radiative cooling time much shorter than that of radial expansion$-$for timescales of months or longer (Fig. 4).

In all our models we arbitrarily shut off the injection after two weeks by having both the density and the velocity of the ejecta exponentially decay, the exact method of the decay does not influence our main results. The time scale of two weeks was chosen to give a peak luminosity time (T2) of order three weeks. Changing the onset of the decay would simply affect the T2 timescale.

Figure 5 shows the radiated luminosity, which tracks the outwards movement of the shocks, as a function of the column density $\Sigma(r)=\int_{r}^{\infty}\rho dr^{\prime}$ to the outflow surface, for the two fiducial models. At early times, just after the onset of the fast flow, the shocks take place behind large columns $\Sigma\gtrsim 10^{22}-10^{25}$ cm2, which are high enough to absorb the shock’s UV/X-ray luminosity across a wide photon energy range $\sim 10$ eV $-1$ keV for a month or longer after the outburst. Following absorption and re-emission (which takes place effectively instantaneously), the shock’s luminosity has been “reprocessed” into optical continuum or line emission. The emission will emerge in the optical band because the ejecta opacity at optical wavelengths is much lower than in the UV/X-ray and the photon diffusion time is much shorter than the outflow expansion rate (such that losses due to adiabatic expansion can be neglected).

Under the reasonable assumption of instantaneous reprocessing, Figure 6 shows the optical light curve of the nova outburst (solid lines) in our different models. The luminosity shown includes both the shocks as well as a contribution $L_{\rm wd}=10^{38}$ erg s*-1* from the central white dwarf (also reprocessed by the ejecta at early times) and has been somewhat arbitrarily held constant in time. Although the single transition scenario (red lines) produces a relatively smooth light curve, the multitude of collisions in the multiple-transition scenario (blue lines) generate several high amplitude flares. That such a complex variable light curve can be generated, even from relatively simple assumptions about the outflow time-evolution and 1D geometry, illustrates that radiative shocks are a promising mechanism to generate flares in classical novae (Strope et al., 2010).

Particularly bright and narrow flares occur in the multiple-transition model when two of the shock-generated thin shells collide and merge. However, our calculations are likely overestimating the luminosity of such second-generation shocks for two reasons. First, our 1D model does not account for angular spread in the outflow properties, i.e. in a realistic situation the shell collision will not take place simultaneously at all polar and azimuthal angles. Second, our simulations are underestimating the thickness of the shells when they collide and thus the duration of the shock interaction. Our 1D simulations do not capture the non-linear thin-shell instability (Vishniac, 1994; Steinberg & Metzger, 2018), which would act to corrugate the shell and thus increase its effective radial thickness.

We define the optical photosphere as the location where the radial optical depth obeys $\tau_{\rm opt}=\Sigma\kappa_{\rm opt}=1$, where $\kappa_{\rm opt}=0.03$ cm2 g*-1* is an estimate of the effective opacity of the nova outflow to electron scattering and Doppler-broadened absorption lines (e.g. Pinto & Eastman 2000).222Free-free opacity generally contributes negligibly at optical/UV wavelengths because the shocked gas rapidly cools to $\lesssim 10^{4}$ K, resulting in a moderate column of ionized gas. The dashed lines in Figure 6 show just that portion of the shock luminosity which is reprocessed above the optical continuum photosphere, which we may therefore take as a rough proxy for the shock’s contribution to emission lines and/or optically-thin free-free emission. Figure 7 shows the radii of the outermost shock and the photosphere (top panel) as well as the blackbody temperature calculated from the radius and luminosity (bottom panel).

In the single transition scenario, the shock quickly moves ahead of the photosphere. Thus, after only a short delay, all of the shock-powered emission happens in optically-thin regions. In such a scenario, in addition to the simple light curve behavior, it would be challenging to explain the correlated behavior between the light curve and photosphere properties (see $\S\ref{sec:photo}$ below). By contrast, in the multiple-transition scenario shell collisions can take place at small radii, near or beneath the optical photosphere, even well after the onset of the outburst. Specifically, the bulk of the shock’s emission will take place below the photosphere if the timescale over which the outflow undergoes the transition from slow-to-fast $\delta t$ is sufficiently short,

[TABLE]

where $t_{\rm thin}$ is the time required for a thin shell of velocity $\sim v_{\rm f}$ takes to reach the photosphere of the slow upstream wind of mass-loss rate $\dot{M}$ and velocity $v_{\rm s}$. Depending on the precise values of $\dot{M}$ and $\delta t$, it is possible to get shocks that release most of their energy below, or above, the photosphere.

Green lines in Figure 6 show the predicted light curves of the “Bumpy" model in which (as in the single transition model) the outflow velocity evolution is monotonic, but rather than a steady-wind the initial slow outflow is radially inhomogeneous (bottom panel, Fig. 2). In this case, one obtains a multi-peaked light curve, as the single shock structure generated by the fast transition propagates into the external-most slow medium. However, the amplitude of the flares are modest (tens of percent) compared to the larger variations observed in some novae and those seen in the multiple-transition model. Furthermore, the shock quickly moves ahead of the photosphere, which (as in the single-transition model) precludes again correlated behavior in light curve and photosphere properties. Variations in the $\dot{M}$ of the initial slow outflow alone do not appear sufficient to explain flares similar to those observed in novae.

3.2 Photosphere Evolution During Flares

Time-resolved optical spectroscopy show that new absorption line systems appear in novae at times coinciding with light curve flares (e.g. Jack et al. 2017). Tanaka et al. (2011a, b) argue that such appearances/disappearances of absorption features are caused by the temporary expansion/shrinking of the optical photosphere. During maxima, the optical photosphere expands in size and the absorption components appear as the ejecta shell expands to become optically-thin. During mimina, the optical photosphere shrinks and the absorption component weakens relative to the emission lines. As we now discuss, such behavior can arise naturally from the multiple-transition internal shock scenario.

When shocks start below the photosphere, this causes an almost instantaneous rise in the optical luminosity, $L_{\rm tot}$. This is because the timescale for photon diffusion from the shock through the ejecta, $\tau_{\rm opt}(R_{\rm sh}/c)$, is short compared to the outflow expansion timescale, $R_{\rm sh}/v$, provided that $\tau_{\rm opt}\lesssim c/v\sim 10^{2}-10^{3}$ ($\Sigma\lesssim 10^{25}-10^{26}$ cm*-2*), where $c$ is the speed of light and $v$ is the radial velocity of the emitting gas. The photosphere radius $R_{\rm ph}$ is initially unchanged before the shock reaches it. Therefore, the new contribution to the luminosity from the shocks causes an immediate increase in the effective temperature of the continuum,

[TABLE]

After a delay $\sim t_{\rm thin}$ (eq. 3) the shock emerges from the photosphere of the slow outflow; the dense shell behind it carries the photosphere temporarily outwards with it, causing $R_{\rm ph}$ to rise and thus $T_{\rm eff}$ to decrease. A time lag $\Delta t$ therefore separates maxima in $T_{\rm eff}$ and $R_{\rm ph}$, which is approximately given by the timescale $t_{\rm thin}$ required for the swept-up shell to reach the photosphere (eq. 3).

Figure 9 provides a side-by-side comparison of the evolution of $L_{\rm tot}$ and $R_{\rm ph}$ for the multiple-transition model. Each maxima in $L_{\rm tot}$ and $R_{\rm ph}$ corresponds to the emergence of a given shock-bounded shell, for which we find delays in the range $\Delta t\approx 0.5-2.3$ days. Such anti-correlated $L_{\rm tot}/R_{\rm ph}$ behavior with a delay is seen in novae such as ASASSN-17pf (Aydi et al. 2019), though the lag may in some cases be too short to observe.

Figure 8 shows the time evolution of shock luminosity in the single and multiple-transition models, but now broken down by the radial velocity of the emitting gas. For shocks which take place near or above the continuum photosphere, these provide a rough proxy for the Doppler widths of features in the optical spectrum (e.g. emission lines or P-Cygni line profiles). However, note that in detail the observed line widths will differ from those predicted by Fig. 8, as our 1D models do not account for angular spread and line-of-sight projection effects in the outflow. We are justified in assuming that the optical emission occurs locally near the cooling gas in our simulation, because of the very short mean free path over which the UV and soft X-ray photons (which carry most of the shock’s power; Steinberg & Metzger 2018) are reprocessed.

The single-transition model (top panel of Fig. 8) predicts smooth time evolution of line features, covering a range of velocities 500$-$800 km s*-1* centered around that of the cool shell. The emission is dominated by the layers of cooling gas behind the forward and reverse shock, which are geometrically narrow but span a range of speeds as the emitting gas settles sub-sonically onto the central cool shell.

In contrast to the relatively simple evolution of the single-transition model, the multiple-transition model (bottom panel of Fig. 8) shows multiple complex velocity features, including both decelerating and accelerating components. The first decelerating feature in Fig. 8 (marked A) is emission from the forward/reverse shock structure generated by the first fast outflow; the fast outflow turns off after a few days, but the forward shock from this shell continues to dominate the luminosity as it sweeps up mass and decelerates. The shocks generated from the second (B) and third (C) transitions generate accelerating line features as the emitting gas is being pushed by their respective fast outflows. Eventually (D), the front end of the shell generated by transitions B catches up to the rear of the shell produced by the first transition (A), leading to a collision between the shell and a brief flare (which as we discussed earlier, may not be accurately captured by a 1D model due to the artificially thin shells). Around the same time, the fourth transition (E) occurs as well as a collision between the second and third shells (B/C).

Acceleration, deceleration and even the apparent “merger” of velocity components are common features of time-resolved spectra of classical novae (e.g. Walter et al. 2012; Aydi et al. 2019). Although detailed radiative transfer on our simulation data would be required to make more definitive statements, our findings suggest that internal radiative shocks, near or above the optical photosphere, provide a promising explanation for this observed behavior.

3.3 Gamma-Ray Emission

Figures 10 and 11 show the predicted X-ray and gamma-ray light curves from our models. The $\sim$ GeV gamma-rays detected from novae by Fermi LAT (Ackermann et al., 2014) are likely produced either by (a) relativistic ions accelerated at the shock, which collide in-elastically with ambient ions to generate neutral pions that quickly decay into gamma-rays; (b) relativistic electrons accelerated at the shock (or secondary electron/positron pairs produced from charged pion decay), which in-elastically scatter ambient nova optical photons up into the gamma-ray band via the Inverse Compton process, or which interact with ambient gas and emit gamma-rays through the (relativistic) bremsstrahlung process. In either of these (a) “hadronic” or (b) “leptonic” models, the timescale for relativistic particles to lose their energy and generate gamma-rays is generally short compared to the outflow expansion timescale over which adiabatic losses occur (Metzger et al., 2016; Martin et al., 2018; Vurm & Metzger, 2018).

The main source of opacity for gamma-rays in the $\sim 0.1-10$ GeV energy range is due to photon-matter pair creation, for which $\kappa_{\gamma}\approx 0.004$ cm2 g*-1* (e.g. Zdziarski & Svensson 1989). Figure 5 shows that, in both single- and multiple-transition models, the shocks spend most of the outburst duration above the gamma-ray photosphere. Nevertheless, the shocks can be deeply embedded at early times after the fast outflow is launched ($\tau_{\gamma}=\Sigma\kappa_{\gamma}>1$; $\Sigma\gtrsim 10^{26}$ cm2 g*-1*) and this could generate a modest delay in the onset of gamma-ray emission relative to the optical peak of a flare (Li et al., 2017).

Given these two facts – fast-cooling cosmic rays and gamma-ray-transparent ejecta – we expect the observed gamma-ray emission to approximately track the instantaneous power of the shock, much as with the reprocessed optical emission. Indeed, correlated gamma-ray and optical light curves are seen in both cases with time-resolved gamma-ray data: ASASSN-15ma (Li et al., 2017) and ASASSN-18fv (Aydi et al., in prep). Metzger et al. (2015) show how this behavior allows for an estimate of the efficiency of relativistic particle acceleration at the shocks, $\epsilon_{\gamma}\sim L_{\gamma}/L_{\rm sh}$, where $L_{\gamma}$ and $L_{\rm sh}$ are, respectively, the gamma-ray luminosity and the portion of the optical luminosity powered by shocks. Application of this method to gamma-ray detected novae results in typical values of $\epsilon_{\gamma}\approx 3-5\times 10^{-3}$ (Li et al. 2017), consistent with theoretical expectations for the efficiency of diffusive shock acceleration.333Once one takes into account the average inclination angle of the upstream magnetic field relative to the shock front, taking into account its highly-corrugated shape when the shock is radiative (Steinberg & Metzger, 2018). The gamma-ray light curves in Figs. 10, 11 were calculated by multiplying the kinetic power of the shocks from the simulation by a constant value $\epsilon_{\gamma}=3\times 10^{-3}$. Our fiducial models predict peak gamma-ray luminosities of $\sim 10^{35}-10^{36}$ erg s*-1*, which are consistent with the range observed by Fermi LAT (Ackermann et al., 2014; Linford et al., 2015; Cheung et al., 2016). Our results are also consistent with Martin et al. (2018), who show that the gamma-ray properties of LAT-detected novae can be understood as hadronic interactions following particle acceleration in internal shocks between faster and slower outflows (the analog of our single-transition model).

3.4 Thermal X-Ray Emission

In contrast to the non-thermal gamma-rays, the dominant source of X-ray emission from the shocks is thermal free-free emission. Non-thermal X-rays (e.g. from relativistic bremstrahlung or inverse Compton scattering emission) are generally suppressed because the electrons that emit in the hard X-ray band possess relatively low energies compared to those emitting in the gamma-ray band and are subject to severe losses to Coulomb scattering off the dense background of the nova ejecta (Vurm & Metzger, 2018). Also in contrast to gamma-rays, the X-ray emission is significantly attenuated by photo-electric absorption in the ejecta. The X-ray light curves shown in Fig. 10 and 11 are calculated assuming each cell $i$ at radius $r_{i}$ in the simulation emits free-free emission at its local temperature $T_{i}$, accounting for photoelectric attenuation according to $e^{-\tau_{\rm X,i}}$, where $\tau_{\rm X,i}=\int_{r_{i}}^{\infty}\kappa_{\rm X}\rho dr$ and $\kappa_{\rm X}$ is the opacity at the X-ray energy of interest.

The X-ray light curves in Figure 10 are calculated assuming solar composition for the ejecta opacity. In our multiple-transition model, the gamma-ray/optical flares are accompanied by hard X-ray $\gtrsim 5-10$ keV flares of luminosity $L_{\rm X}\sim 10^{33}-10^{34}$ erg s*-1*. At face value, these luminosities appear to be broadly consistent with the NuSTAR detection of thermal X-ray emission contemporaneous with the LAT gamma-ray detection of V5855 Sgr (Nelson et al. 2019). However, Nelson et al. (2019) infer very high temperatures for the un-absorbed X-ray emitting-plasma ($kT\sim 5-20$ keV) which require a higher shock velocity $\gtrsim 2000-4000$ km s*-1* than is obtained in our fiducial model. If we had instead adopted a much faster white dwarf outflow velocity $v_{\rm f}$ (for otherwise similar kinetic energy), our predicted NuSTAR band X-ray fluxes would be several orders of magnitude greater $\gtrsim 10^{36}-10^{37}$ erg s*-1* than those shown in Fig. 10, too high compared to observations.

As discussed by Nelson et al. (2019), the X-rays observed from V5855 Sgr are thus in fact sub-luminous relative to expectations given the high observed gamma-ray luminosity. To reconcile this fact, one is forced to conclude that either the X-rays component which dominates in NuSTAR originates from distinct shocks from those which are responsible for powering the gamma-rays, i.e. if the latter are too deeply buried for even hard X-rays to escape. Alternatively, the hard X-ray luminosity of the gamma-ray-producing shocks could be suppressed from the naive predictions of a 1D model. Supporting the latter possibility, Steinberg & Metzger (2018) found that the X-ray luminosity of high-Mach number $\mathcal{M}\gg 1$ radiative shocks is suppressed by a factor of $\sim 4.5/\mathcal{M}^{2}\sim 1/100$ for radiative shocks relative to the 1D predictions (see also Kee et al. 2014). As mentioned above, a non-thermal origin for the hard X-rays from V5855 (e.g. bremsstrahlung emission from the relativistic particles generating the gamma-rays) is disfavored, due to both the luminosity and spectral slope of the emission (Vurm & Metzger, 2018).

In contrast to $\gtrsim 10$ keV X-rays, lower-energy X-rays (e.g. in the $\lesssim 5$ keV Swift XRT bandpass) are suppressed for several months or longer after the eruption. The $\sim$ keV X-ray light curves in our models peak only once the forward shock of the outermost monolithic shell (created by the merger of previous shells) reaches a sufficiently low external column $\tau_{X}\lesssim 1$. Our models predict peak keV luminosities from the shocks which range from undetectably weak to $L_{\rm X}\sim 10^{34}-10^{35}$ erg s*-1*, depending sensitively on the mass, radial distribution, and composition of the slowly expanding external medium which exists prior to the onset of the first fast outflow.

Figure 11 shows how the predicted X-ray light curves depend on the properties of the pre-existing external slow medium, for an otherwise identical outflow behavior at $t>0$ (after the onset of the first fast outflow). Specifically, we consider how the results change if we adopt the standard CNO-enhanced nova abundance set from CLOUDY (Ferland et al., 2017) instead of solar abundances. We also consider results for the Low Mass model, in which the pre-existing slow external medium is three times less massive than in the Fiducial model due to its smaller initial radial extent $r_{0}=2\times 10^{13}$ cm (Fig. 2). The greater keV opacity of enhanced CNO abundances compared to the solar metallicity case reduces the column $\Sigma$ ahead of the shock at the time of the X-ray peak ($\tau_{X}=1$), thus resulting in a lower shock power $\propto\Sigma v_{\rm sh}^{3}$ and correspondingly lower X-ray luminosity. The effect of smaller $r_{0}$ and lower mass in the external medium is to increase the forward shock velocity $v_{\rm sh}$ (since the shell has less mass to sweep up) and shorten the timescale over to achieve $\tau_{X}=1$, both which increase the peak X-ray luminosity. The wide range of X-ray luminosities our models find, even for moderate variations in the nova outflow properties, appears to be at least qualitatively consistent with the diverse range of X-ray light curve behavior in novae (Mukai & Ishida, 2001; Orio et al., 2001b, a; Mukai et al., 2008; Schwarz et al., 2011).

The late-time onset of soft X-rays in classical novae is commonly interpreted as emission from the white dwarf surface (e.g. Schwarz et al. 2011). However, we note that in some cases such as our Fiducial model with a high external mass, the forward shock velocity is sufficiently low at late times for the X-ray emission to peak at low energies $\lesssim$ 0.3 keV (bottom panel of Fig. 10) in which case shock emergence could masquerade as the onset of the supersoft phase. On the other hand, differences in the predicted spectral energy distribution (blackbody emission from the white dwarf surface versus optically-thin emission from shocks at much larger radii) should be discernible with good spectral coverage.

Also note that our 1D models do not take into account the suppression of X-ray emission from radiative shocks found by Steinberg & Metzger (2018) and therefore the luminosities shown in Figs. 10, 11 may be over-estimates. This caveat becomes less severe with time, as the forward shock transitions from being radiative to adiabatic over a month or longer (depending on the velocity of the shocks and the composition of the upstream medium; Fig. 4).

3.5 Radio Emission

The radio emission from classical novae was long thought to originate exclusively from thermal free-free emission generated by $\sim 10^{4}$ K photo-ionized ejecta (e.g. Hjellming et al. 1979; Cunningham et al. 2015). However, a growing sample of novae show additional maxima in their early-time radio light curves when the ejecta radius is still small and thus the brightness temperature $\gg 10^{4}$ K, which indicate an additional contribution of non-thermal synchrotron emission (e.g. Mioduszewski & Rupen 2004; Weston et al. 2016; Vlasov et al. 2016).

We calculate the synchrotron emission from the forward shock444Radio emission from the reverse shock is strongly attenuated by free-free absorption in the dense cool shell and is therefore negligible in comparison to the forward shock emission in a 1D model. in our models via the following procedure. The observed brightness temperature at radio frequency $\nu$ is given by

[TABLE]

where $\tau_{\rm r}(t)$ is the free-free optical depth from the shock surface to infinity (see below) and $T_{\rm B,\nu}(t)$ is the un-attenuated brightness temperature of synchrotron emission from the post-shock gas. We approximate the latter by the expression (Vlasov et al. 2016, their eq. 40)

[TABLE]

where $\Lambda(T_{\rm sh})$ is the cooling function evaluated at the post-shock temperature $T_{\rm sh}$; $n\equiv\rho/m_{p}$ is the density of the gas just upstream of the shock; $v_{\rm sh}$ is the shock velocity; $\epsilon_{B}/\epsilon_{e}$ are the fraction of the shock power placed into magnetic fields and relativistic electrons, respectively, normalized to characteristic values; $p$ is the index of the energy distribution of the accelerated electrons (e.g. $dN_{e}/dE\propto E^{-p}$), which we take to be $p=2.5$.555Although synchrotron emission can also originate from secondary electron/positron pairs downstream of the shock by charged pion decay in the hadronic model for the gamma-rays, (Vlasov et al., 2016) show that the observed emission is likely dominated by primary electrons directly accelerated at the shock (i.e. the radio emission is “leptonic” even if the gamma-rays are “hadronic”). In deriving equation (6) we have assumed that the shock is radiative, such that the width of the synchrotron-emitting region behind the shock responsible for the synchrotron emission is controlled by the thickness of the post-shock cooling layer $\propto 1/\Lambda$. Since the cooling layer thickness scales as $t_{\rm cool}/t_{\rm exp}$, Figure 4 shows that assuming the shocks are radiative is good for the first $\approx 40$ days in our low-mass external medium model ($r_{0}=2\times 10^{13}$ cm) and substantially longer in the fiducial case ($r_{0}=6\times 10^{13}$ cm).

The optical depth due to free-free absorption through the ejecta is given by (Rybicki & Lightman 1979)

[TABLE]

where $T$, $n$, and $\nu$ are in cgs units and the pre-factor of $\alpha_{r}$ depends only weakly on the ejecta composition (Metzger et al., 2014). We calculate the ionization state of the ejecta ahead of the shock at each timestep following the procedure described in Appendix B. In practice, complete ionization is a reasonable approximation at late times when $\tau_{\rm r}\lesssim 1$ given the large Lyman continuum flux from the UV radiation of the radiative shock incident on the dilute upstream medium. Finally, we convert brightness temperature back to radio luminosity using

[TABLE]

where $R_{\rm sh}$ is the shock radius. Again, the luminosity we present assumes a spherical outflow and therefore should be multiplied by $f_{\Omega}\sim 0.3$ to account for the smaller angular extent of the equatorial shock region.

In addition to the non-thermal synchrotron emission, we calculate the thermal free-free radiation from the hot photo-ionized gas following the procedure described in Appendix B. Except during brief intervals when the radio photosphere is passing through a shock, the brightness temperature of the thermal emission is at the temperature floor $10^{4}$ K of the simulation. This is reasonable physically because the latter is close to the temperature expected for gas photo-ionized by UV/X-rays from the shocks or the central white dwarf (Cunningham et al., 2015). Also note that the light curve shape from the true bipolar ejecta morphology, including both the fast polar and slower equatorial components, will differ in detail from the predictions of a 1D model presented here (e.g. Ribeiro et al. 2014b).

Figure 12 shows the 10 GHz radio light curves (top panel) and synchrotron brightness temperatures (bottom panel) for the multiple-transition model. We show results separately for the Fiducial ($r_{0}=6\times 10^{13}$ cm) and Low Mass ($r_{0}=2\times 10^{13}$ cm) cases. In both models, the thermal emission (blue lines) peak after roughly one year, consistent with that of observed thermal emission in novae (e.g. Hjellming et al. 1979; Chomiuk et al. 2014). By contrast, the non-thermal synchrotron emission (red lines) varies by two orders of magnitude between the models. The synchrotron component furthermore only peaks above the thermal free-free emission in the Low External Mass model ($r_{0}=2\times 10^{13}$ cm). Figure 13 shows similar light curves, but in the case of a higher frequency $\nu=100$ GHz. In this case the synchrotron emission is overwhelmed by the thermal component at all times in both models.

As in the X-ray case, although the power of the shocks peaks at early times, the peak in the synchrotron light curve is delayed by free-free absorption through the slow external medium until $\tau_{\rm r}\lesssim 1$. The wide range in the peak luminosity of the shocks results because of the sensitive dependence of the intrinsic synchrotron luminosity on the shock velocity ($T_{\rm B}\propto v_{\rm sh}^{6.8}$ for $p\approx 2.5$; eq. 6). The value of $v_{\rm sh}$ at late times is in turn sensitive to the amount of slow ejecta into which the outermost shock is decelerating. The large diversity in the predicted peak luminosities of the non-thermal radio emission appears broadly compatible with the range $\nu L_{\nu}\lesssim 3\times 10^{28}$ erg s*-1* (Chomiuk et al., 2014) to $\sim 3\times 10^{30}$ erg s*-1* (Weston et al., 2016) from observations given realistic variations in the properties of the external medium and of the fast outflow from the white dwarf.

3.6 Dust Formation

Dust formation is commonly inferred to take place in the ejecta of novae (e.g. Gehrz et al. 1998). This is evidenced by the delayed onset of infrared emission and the occurrence in many cases of deep minima in the optical light curves caused by the sudden formation of dust along the line of site (e.g. Evans et al. 2005; Strope et al. 2010; Evans et al. 2017; Harvey et al. 2018; Finzell et al. 2018; Aydi et al. 2019). The expanding ejecta becomes cool enough for solid condensation to take place within weeks to months of the outburst, once the ejecta reaches radii larger than,

[TABLE]

where $T=T_{\rm cond}\approx 1500$ K is the approximate condensation temperature and $L$ is the central irradiating luminosity. Even once the condition $r>r_{\rm dust}$ is satisfied, however, it is not clear how molecules and grains are able to form given the potentially harsh environment created by ionizing radiation from the white dwarf (e.g. Pontefract & Rawlings 2004), particularly if the ejecta shell is part of a smooth temporally-extended (and thus radially thick) outflow.

Derdzinski et al. (2017) argue that the cool, dense gas generated by radiative shocks provide ideal conditions capable of dust nucleation and growth (see also Aydi et al. 2019). We expand upon this possibility explicitly here in Figure 14, which shows the mass-weighted distribution of gas density for matter passing through the condensation radius $r_{\rm dust}$ in our single and multiple-transition models. Our 1D simulations cannot capture the full multi-dimensional, multi-phase structure of the thermally-unstable post-shock gas. To generate Figure 14 we have therefore convolved the relatively narrow density distribution derived from our 1D models with a lognormal distribution (with a mean $\sim$ variance $\sim\mathcal{M}^{2}$) motivated by the results of high-resolution 2D simulations of high Mach number $\mathcal{M}\gg 1$ radiative shocks in Steinberg & Metzger (2018). These simulations capture the wide range of gas densities in the multi-phase region generated by thermal instabilities in the post-shock gas.

It is fruitful to compare the density distribution of shocked gas to the alternative simpler scenario of a steady outflow in which no shocks occur. For an outflow of constant mass-loss rate $\dot{M}$, the outflow density at $r=r_{\rm dust}$ is

[TABLE]

Comparing a characteristic value of $\rho_{\rm w}$ (shown as a vertical dashed line in Fig. 14) to that obtained in the case of radiative shocks for a flow with otherwise similar mean mass outflow rates, we see that the shock case spans $\sim 4-5$ orders of magnitude in density.

The “clumpy” ejecta structure generated in the thin shells of radiative internal shocks could provide a natural explanation for the small average filling factors inferred in nova ejecta of $\lesssim$ few percent (e.g. Shore et al. 2013), which otherwise would be challenging to produce in a steady-wind or the smooth homologously-expanding outflow from a singular mass ejection event. These high-density pockets also provide local environments for the nucleation and growth of dust grains as they are shielded from the otherwise detrimental affects of X-ray/UV radiation (Derdzinski et al., 2017).

4 Conclusions

We have studied internal shocks in classical nova outflows by means of 1D hydrodynamical simulations including radiative losses. Our primary aim is to model the shock interaction between slow and fast outflows from these systems and to assess the ability of shocks to simultaneously explain multi-wavelength observations. Figure 15 summarizes the broad range of emission components (both thermal and non-thermal) which contribute to the spectral energy distribution of novae, many of which are fundamentally shock-powered.

Our conclusions may be summarized as follows:

• •

The hourglass/embedded ring geometry of nova ejecta on large spatial scales is suggestive of earlier shock interaction between a fast outflow from the white dwarf with a slower equatorially-concentrated outflow shaped by the binary orbit (Fig. 1). However, this structure would be relatively insensitive to the detailed time evolution of the outflow at early times, in particular whether the transition from slow-to-fast outflow modes occurs only once (“single-transition” scenario) or several times in succession (“multiple-transition” scenario).

• •

The interaction between consecutive fast and slow outflow components is mediated by a forward-reverse shock structure, which sweeps up and compresses the shocked gas into a narrow dense shell. Each transition back to a fast flow generates its own dense shell. In most cases these shells merge pairwise into a single monolithic shell, which continues to propagate down the density gradient of the initial slow outflow until late times (Fig. 3).

• •

For the first $\sim$ month or longer after the outburst, the shocks take place deep under the UV and soft X-ray photosphere. Direct radiation from the shock is absorbed and reprocessed into optical emission, powering time-variable optical luminosity on top of the more smoothly-varying light curve powered by nuclear burning on the white dwarf surface (Fig. 6). In single-transition scenarios, the shock quickly moves ahead of the optical photosphere, such that after days or less its reprocessed luminosity is channeled into emission lines or optically-thin emission.

By contrast, the fresh shocks continually generated at small radii in the multiple transition scenario result in much of the shock emission taking place below the photosphere, where it can contribute to the optically-thick continuum. This behavior leads to correlated maxima in the optical light curve and photospheric radii, with a possible short time lag of days or less, consistent with observations of some flares in novae (e.g. Tanaka et al. 2011a; Aydi et al. 2019). For shock-powered optical emission above the photosphere, the predicted evolution of the spectral line velocities can be complex, with accelerating, decelerating, and even merging velocity components (Fig. 8). The emission lines originate from material behind the forward and reverse shocks with a velocity width intermediate between those of the slowest and fastest outflows. The slow unshocked medium may instead be that which is probed by transient narrow absorption line features (e.g. the THEA systems described by Williams et al. 2008).

• •

Thermal instabilities behind radiative shocks naturally produce a clumpy medium in the post-shock central shells, spanning several orders of magnitude in the gas density (Fig. 14). Mason et al. (2018) argue, from high-resolution UV and optical spectroscopy, that gamma-ray producing novae can be understood by a clumpy ejecta shell expelled during a brief, singular ejection episode, rather than necessitating a wind-type geometry. Our models predict that the majority (of at least the equatorial component) of the ejecta has passed through at least one radiative shock on its journey to large radii. Such clumpiness is therefore a natural prediction of the model. The same clumpy thin shells provide natural sites for dust nucleation, once the ejecta shells reach sufficiently large radii for solid condensation.

• •

As both gamma-ray and optical emission track the shock power with little or no delay, correlated optical and gamma-ray flares are a natural prediction of the Multiple Transition model in which a significant fraction of the optical luminosity originates from shocks. The observed range of luminosities are obtained if the acceleration efficiency of relativistic particles at the shock is $\sim 0.3\%$, similar to those expected due to ion acceleration at the corrugated interface of radiative shocks (Steinberg & Metzger, 2018) and thus further supporting a hadronic origin for the gamma-ray emission (e.g. Metzger et al. 2016; Li et al. 2017; Martin et al. 2018).

• •

The Multiple Transition scenario is clearly favored by the gamma-ray/optical variability of novae. However, this does not address the physical origin of the abrupt and frequent transitions in the nova outflow properties between “fast" and “slow" modes, over a wide range of timescales. The fast outflow likely originates directly from the white dwarf, while$-$given its low velocity and equatorially-focused geometry$-$the slow outflow probably originates from Roche lobe overflow (by frictional heating, non-axisymmetric gravitational torques from the binary, or in a wind from a circumbinary excretion disk; Pejcha et al. 2016b, a). Kato & Hachisu (2011) model the structure of the white dwarf burning layer accounting in an approximate matter for the effects of the binary companion. They find distinct optically-thick wind (“fast") and static envelope (“slow") phases, transitions between which could give rise to the observed variability (and, in our picture, internal shocks). However, to advance to the stage of being truly predictive (e.g. regarding the timescale between outflow mode switches), models of novae must account self-consistently for the multi-dimensional boundary conditions imposed by the Roche Lobe and identify the dominant physical mechanism driving the slow outflow. .

• •

Due to the relatively low photoelectric opacity at higher photon energies, flares in the hard $\gtrsim 10$ keV band may accompany the optical/gamma-ray flares (Figs. 10, 11). This is particularly true in cases in which external medium has a low metallicity and solar-like composition (e.g. ejecta dominated by unburnt material from the companion star). On the other hand, the hard X-ray emission could be substantially suppressed from the naive predictions of the 1D jump conditions if the shocks remain radiative to late times (Steinberg & Metzger, 2018).

• •

In contrast to the hard X-rays, the soft $\lesssim 1$ keV emission is typically delayed for weeks or longer by photoelectric absorption in the external slow outflow. This emission originates not from the deeply-embedded shell collisions themselves, but from the forward shock of the outermost shell, once it reaches a sufficiently low external column. The luminosity range of the $\sim$keV X-ray emission spans several orders of magnitude, being particularly sensitive to the mass, composition, and radial distribution of the initial slow outflow, or pre-existing medium surrounding the binary, at the time of the first transition to a fast outflow. In some cases the X-ray light curve may peak at sufficiently low energies $\lesssim 0.3$ keV that it could be confused for emission from the white dwarf surface in the absence of good spectral coverage.

• •

Like the soft X-rays, radio synchrotron emission from shock-accelerated electrons is not generally visible until late times $\gtrsim$ weeks-months (in this case due to free-free absorption) from the forward shock of the outermost shell (Fig. 12). If the forward shock is radiative, the radio emission is extremely sensitive to the velocity of the forward shock, which in turn depends on the momentum input of the fast outflow and the mass profile of the initial slow outflow or pre-existing external medium.

The keV X-ray and synchrotron radio emission together provide a potentially sensitive probe of the initial (potentially pre-thermonuclear runaway) medium surrounding the white dwarf, the origin of which is debated (e.g. Williams et al. 2008). Unfortunately detailed inferences may be plagued by several degeneries, such as the microphysical parameters $\epsilon_{e}$, $\epsilon_{B}$ in the radio case and the precise ejecta composition in the X-ray case. Very early time radio detections, especially at low frequencies where free-free absorption is severe, cannot be explained by the slow equatorial outflow and may require an additional emission site, such internal shocks in the polar regions due to variability in the fast outflow itself (Vlasov et al. 2016).

Acknowledgements

We thank Elias Aydi, Laura Chomiuk, and Indrek Vurm for helpful information and discussions. This work is supported in part by NSF through grant #AST-1615084 and by NASA through the Fermi Guest Investigator Program (grant # 80NSSC18K1708). BDM received partial support from the Simons Foundation through the Simons Fellows Program (grant # 606260).

Appendix A Numerical Method

We solve the Euler equations using a 1D second-order Lagrangian Finite-Volume scheme with spherical geometry. The code employs an exact Riemann solver with an ideal gas equation of state with a polytropic index $\gamma=5/3$. The code has been validated against a suite of smooth (e.g. linear waves) and non-smooth (e.g. Sedov blast wave) test problems.

During the run we attempt to enforce cells to have a size $dr\approx 3\cdot 10^{-4}r$. However, once shocked cells cool, their size decreases by a factor of $\sim\mathcal{M}^{2}$, which greatly reduces the time step. Whenever a cell is smaller than the desired threshold, we merge it with its smaller neighbor, provided that it is in a smooth region. If the cell is not in a smooth region, we remove it only if it satisfies that the ratio in the pressure, temperature and density between both of its neighbors satisfies:

[TABLE]

We include optically-thin radiative losses using the cooling function calculated from version C17.01 of CLOUDY (Ferland et al., 1998, 2017) using the GASS10 (solar) abundances and the exact cooling integration method as described in Townsend (2009).

Appendix B Radio Opacity

In this section we outline our algorithm to determine the radio opacity. We start by identifying the location of the outermost shock, $r_{\rm f}$. We assume that the WD has an ionizing photon luminosity of

[TABLE]

and we subtract from it the total rate of recombinations up to the forward shock

[TABLE]

where $\alpha_{\rm B}=2.6\cdot 10^{-13}\left(\frac{T}{10^{4}\,{\rm K}}\right)^{-0.8}\,{\rm cm^{3}\,s^{-1}}$ is the case B recombination rate (Osterbrock & Ferland, 2006). The final ionizing photon flux emerging from the forward shock, $Q_{\rm h}$ is the net transmitted ionizing photon flux, $Q_{\rm h,net}$, with the addition of the ionizing photon flux generated at the forward shock assuming a free-free spectrum. Following the procedure described in Metzger et al. (2014) (eq. 42 and 43) we estimate the ionization state of the hydrogen above the forward shock.

If the radio photosphere is above the forward shock (i.e. $r(\tau=1)>r_{\rm f}$), we estimate the free-free radio emission to be

[TABLE]

where $\varepsilon_{ff}$ is the free-free emissivity and $T^{\rm ion}=10^{4}$K.

If the photosphere is below the forward shock, we calculate the ionization state of the gas between the surface of the Strömgren sphere and the photosphere using the values for CIE based taken from Gnat & Sternberg (2007).