Simulations of the symbiotic recurrent nova V407 Cyg. I. Accretion and shock evolutions
aa r X i v : . [ a s t r o - ph . H E ] M a r SIMULATIONS OF THE SYMBIOTIC RECURRENT NOVAV407 CYG. I. ACCRETION AND SHOCK EVOLUTIONS
Kuo-Chuan Pan ( 潘 國 全 ) , Paul M. Ricker , and Ronald E. Taam , Physik Department, Universit¨at Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland;[email protected] Department of Astronomy, University of Illinois at Urbana − Champaign, 1002 West GreenStreet, Urbana, IL 61801, USA; [email protected] Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road,Evanston, IL 60208, USA; [email protected] Academia Sinica Institute of Astronomy and Astrophysics, P.O. Box 23-141, Taipei 10617,Taiwan; [email protected] ; accepted 2 –
ABSTRACT
The shock interaction and evolution of nova ejecta with a wind from a redgiant star in a symbiotic binary system are investigated via three-dimensionalhydrodynamics simulations. We specifically model the March 2010 outburst ofthe symbiotic recurrent nova V407 Cygni from the quiescent phase to its eruptionphase. The circumstellar density enhancement due to wind-white dwarf inter-action is studied in detail. It is found that the density-enhancement efficiencydepends on the ratio of the orbital speed to the red giant wind speed. Unlikeanother recurrent nova, RS Ophiuchi, we do not observe a strong disk-like den-sity enhancement, but instead observe an aspherical density distribution with ∼
20% higher density in the equatorial plane than at the poles. To model the2010 outburst, we consider several physical parameters, including the red giantmass loss rate, nova eruption energy, and ejecta mass. A detailed study of theshock interaction and evolution reveals that the interaction of shocks with thered giant wind generates strong Rayleigh-Taylor instabilities. In addition, thepresence of the companion and circumstellar density enhancement greatly alterthe shock evolution during the nova phase. The ejecta speed after sweeping outmost of the circumstellar medium decreases to ∼ −
300 km s − , dependingon model, which is consistent with the observed extended redward emission in[N II] lines in April 2011. Subject headings: hydrodynamics — methods: numerical, — novae: general —stars:individual: V407 Cygni
1. INTRODUCTION
The symbiotic recurrent nova (SyRN) V407 Cyg consists of a white dwarf (WD)and a Mira-type red giant (RG, Munari et al. 1990). In March 2010, the gamma-rayemission from its outburst was detected by
Fermi -LAT (Large Area Telescope, Atwood etal. 2009), showing the first evidence for a gamma-ray nova (Abdo et al. 2010). Anothergamma-ray SyRN, V745 Sco, was also detected by
Fermi -LAT in February 2014 (Cheung etal. 2014). Since the RGs in symbiotic binaries could produce a high-density circumstellarmedium (CSM) through their stellar wind, gamma-ray emission in SyRNe could originatein high-energy particle acceleration that happens when blast waves pass through thehigh-density CSM (Abdo et al. 2010; Martin & Dubus 2013).Recently, classical novae (CNe) have also been classified as a new type of gamma-raysource (Ackermann et al. 2014). However, the lack of a high-density CSM in CNe cannotexplain the origin of gamma-ray emission as in SyRNe, requiring a new mechanism forgamma-ray CNe. Chomiuk et al. (2014) studied the high-resolution radio maps of thegamma-ray CN V959 Mon, proposing that the fast binary motion in CNe could shape thenova ejecta by transferring its angular momentum. Therefore, high-energy particles couldbe accelerated at the interface between equatorial and polar regions.Since particle acceleration in CNe and SyRNe is highly associated with shockpropagation in their nova eruptions, these new discoveries motivate us to study the shockevolution in nova systems. Additionally, SyRNe have long been considered as Type Iasupernova (SN Ia) progenitors (Hachisu & Kato 2001; Kato & Hachisu 2012), although thenature of the progenitor systems of SNe Ia remains uncertain (Maoz et al. 2014). However,evidence of interaction between SN Ia ejecta with a CSM has been found in some SNe Iain late-type spiral galaxies (Patat et al. 2007; Dilday et al. 2012; Silverman et al. 2013),suggesting that at least some SNe Ia could originate in symbiotic binaries. 4 –The gamma-ray emission from the 2010 outburst of V407 Cyg has been studied byMartin & Dubus (2013) through a semi-analytical study of diffusive shock accelerationand non-thermal emission in V407 Cyg. They found that the observed gamma-ray lightcurve could be fitted by requiring a circumbinary density enhancement (CDE) around theWD. Multi-dimensional hydrodynamics simulations by Orlando & Drake (2012) furthersupport the need for a CDE around the WD. They simulated a blast wave passing througha dense CSM in V407 Cyg including a sophisticated treatment of thermal conduction andradiative cooling (Orlando & Drake 2012), finding that the observed X-ray light curvescould be reproduced if a CDE exists around the WD. However, they do not include theasymmetric effect of binary motion, and the distribution of CDE is artificially imposed intheir simulations without modeling the WD-wind interaction in the quiescent phase.On the other hand, Walder et al. (2008) performed three-dimensional hydrodynamicssimulations of another SyRN, RS Oph, from the quiescent phase to a nova eruption. Intheir simulations, they found a disk-like CDE concentrated in the orbital plane, and theirshock evolution is consistent with the observations of the 2006 nova outburst (Sokoloski etal. 2006; Bode et al. 2006). Although some similarity between RS Oph and V407 Cyg hasbeen found in their eruptions’ spectra (Munari et al. 2011; Shore et al. 2011), the physicalparameters of the binary systems are quite different (Drozd et al. 2013), especially for theirorbital periods (separations). The orbital period in RS Oph is much shorter ( P orb ∼ .
25 yr)than in V407 Cyg ( P orb ∼
43 yr). Therefore, the amount and distribution of densityenhancement could be different in these two systems.In this paper, we investigate the evolution of the SyRN V407 Cyg from the quiescentphase to a nova eruption. In the next section, we review and summarize some observationalfacts about the binary system of V407 Cyg and its previous outbursts. The numericalmethods and the initial setup are presented in Section 3. Our simulation results with 5 –different red giant wind loss rates and nova eruption parameters are reported and discussedin Section 4. In the final section, we summarize our results and conclude.
2. OBSERVATIONS OF V407 CYG
As mentioned in the above section, the SyRN V407 Cyg consists of a massive WD anda Mira-type M6 III RG (Munari et al. 1990). The Mira-type RG has a pulsation period of763 day (Kolotilov et al. 2003). From the dust obscuration and the long-term optical lightcurve, an orbital period of ∼
43 yr was inferred by Munari et al. (1990), corresponding toan orbital separation of 16 AU for a 1 M ⊙ RG and a 1 . M ⊙ WD. The distance of V407 Cygis still unclear. Some distance estimates using the period-luminosity relation give a rangefrom 1 . . & . H α emission line-width, an ejecta mass M ej ∼ − M ⊙ , eruption energy E ej ∼ ergs, and ejecta speed v ej ∼ , ±
345 km s − could be estimated (Abdo et al. 2010; Schwarz et al. 2011). An ejecta mass of M ej ∼ − M ⊙ could also reproduce the observed velocity evolution, but may underestimate the RG massloss rate (Chomiuk et al. 2012). 6 –In general, the observed mass loss rate ˙ M wind of a Mira giant spans a wide range,varying from ∼ − M ⊙ yr − to ∼ − M ⊙ yr − with an average h ˙ M wind i ∼ × − M ⊙ yr − (Knapp & Morris 1985). The radio light curve from the 2010 outburst of V407 Cyg suggestsa mass loss rate ˙ M wind ∼ × − M ⊙ yr − (for v wind = 10 km s − and d = 3 kpc; Chomiuket al. 2012), but this mass loss rate is one order of magnitude higher than the estimatedvalue from the X-ray light curve ( ˙ M wind ∼ − − − M ⊙ yr − , Nelson et al. 2012; Martin& Dubus 2013). To obtain better agreement between the X-ray and radio estimates, abinary separation wider than 16 AU may be required (Chomiuk et al. 2012).The wind speed of a giant star correlates with the star’s luminosity (De Beck et al.2010; Politano & Taam 2011). Taking account of the observed luminosity L ∼ L ⊙ (Munari et al. 2011) and the uncertainty of the distance of V407 Cyg, a RG wind speedof 10 −
20 km s − can be estimated. In addition, a wind speed of ∼
10 km s − could bealso estimated from the optical P Cygni line profiles (Abdo et al. 2010). However, a higherwind speed of 30 km s − is estimated from the equivalent width of the Na I D lines in thequiescent phase (Tatarnikova et al. 2003) and in the eruption phase (Shore et al. 2011;Nelson et al. 2012).
3. NUMERICAL METHODS3.1. Simulation Code
The simulation code used in this study is FLASH version 4 (Fryxell et al. 2000; Dubeyet al. 2008; Lee 2013). FLASH is a grid-based, parallel, multidimensional hydrodynamicscode based on block-structured adaptive mesh refinement (AMR). The split piecewiseparabolic method (PPM) solver (Colella & Woodward 1984) is used to solve for the motion http://flash.uchicago.edu To model V407 Cyg and its 2010 eruption, we consider a three-dimensional (3D)simulation box with a size of 400 AU ×
400 AU ×
400 AU. Initially, a 1 M ⊙ RG is located atthe center (origin) of the simulation box, and a 1 . M ⊙ WD is located on the positive x − axiswith a binary separation of 16 AU, corresponding to an orbital period of about 43 yr. Theorbital plane is set on the x − y plane and the direction of angular momentum is set to be thepositive z − axis. Particle clouds representing the WD and the RG have a radius of 10 cm,which is slightly larger than three smallest zone spacings (∆ x min = 2 . × cm). We use9 levels of refinement based on the magnitudes of the second derivatives of gas density andpressure. Each AMR block contains 8 zones in the 3D simulation box, corresponding to aneffective uniform resolution of 2048 . To save computation time, we reduce the maximumAMR level based on the distance to the center of mass, giving an effective resolution of∆ x i ∼ . × d at the i th level of AMR, where d is the distance to the center of mass ofthe WD-RG binary. The highest resolution region contains the central 65 AU and thefirst forced AMR decrement occurs at d ∼
32 AU. The second AMR decrement occurs at d ∼
63 AU, the third at d ∼
127 AU, and so on. Outflow boundary conditions for fluidsand isolated boundary conditions for the Poisson solver are used. 8 –We assume the RG has a radius of 2 . v wind = 20 km s − in the co-rotating frame, and a constant mass loss rate ˙ M wind . Therefore, the initial gasdensity ρ is distributed using Equation 1: ρ ( r ) = ˙ M wind πr v wind = 1 . × − ˙ M wind − M ⊙ yr − ! (cid:16) v wind
20 km s − (cid:17) − (cid:16) r (cid:17) − g cm − , (1)where r is the distance to the RG. We assume a gamma-law equation of state with γ = 1 . T wind = 7 ,
000 K to mimic the thermodynamic behavior of gasbased on photoionization models of symbiotic binary systems (Nussbaumer & Vogel 1987;Nussbaumer & Walder 1993).We perform quiescent-phase simulations up to ∼ . − . . r < M ej , eruption energy E ej , and averaged ejectaspeed v ej . To avoid grid effects, the eruption covers a small spherical region with a radiusequal to fifteen smallest zone spacings (Pan et al. 2012). Within the erupting region, auniform density distribution and linear velocity profile in radius are used. We assume thekinetic energy is three-fourths of the total eruption energy (Dohm-Palmer & Jones 1996).The mass and kinetic energy of the CDE within the exploding region and the orbital speedof the WD are added to the eruption. Once the SyRN eruption is imposed, we reducethe resetting radius for the RG wind to the radius of the RG, and assume a completelyabsorbing surface on the RG. In reality, this assumption is not precisely correct, and areverse shock should be generated when the ejecta reach the surface of the RG. We ignorethese effects since the purpose of the paper is not to focus on the response of the RG, buton the global shock evolution. 9 –
4. RESULTS AND DISCUSSION
We perform runs with several different sets of values for the RG wind and SyRNeruption parameters, such as ˙ M wind , E ej , and M ej (see Table 1), and present our simulationresults in this section. In the first subsection, the evolution of wind focused by the WDduring the quiescent phase is described. The asymmetry due to WD-wind interaction isinvestigated as well. In Section 4.2, we study the shock evolution dependence on differentRG wind loss rates and different eruption energies and ejecta masses. We perform two quiescent-phase simulations by varying the mass loss rate from˙ M wind = 10 − M ⊙ yr − (model M6) to ˙ M wind = 10 − M ⊙ yr − (model M7; see Table: 1) tostudy the effects on the CDE. The quiescent phase evolution of both models is qualitativelythe same, but the gas density distributions scale with their mass loss rates. In thissubsection, we take model M6 as the reference simulation and describe the quiescent-phaseevolution in detail.Figure 1 demonstrates a typical evolution of the gas density in the orbital plane duringthe quiescent phase (model M6). Since the initial wind density distribution (see Equation 1)neglects the binary interaction, a turbulent and hot spiral arm quickly forms when theRG wind interacts with the WD. This turbulent spiral arm takes about an orbital period( ∼
40 yr) to reach a quasi-steady large-scale state in the co-rotating frame. We also showthe corresponding gas temperature distribution in Figure 2.The WD gravitationally deflects matter in its vicinity. Figure 3 shows the mass withincontrol surfaces for different radii for our models M6 and M7. It is clear that the enclosedmass increases for the first few years and then reaches stability after about 10 yr. Since the 10 –RG wind speed is higher than the orbital speed of the WD, the WD experiences a strongwind which limits the amount of mass enclosed within the control surface. Therefore, weobserve no significant wind-capture disk around the WD as reported by Walder et al.(2008), since the wind speed relative to the orbital speed is much higher than that inRS Oph. However, we do observe a small, elongated, disk-shaped structure around the WDwith a size about 5 AU (See Figure 4). This result is consistent with the fast simulation case( v wind = 60 km s − ) of RS Oph in Walder et al. (2010), since the orbital speed of RS Oph( v orb ∼
15 km s − ) is much faster than that in V407 Cyg ( v orb ∼ − ), giving a similarorbital to wind speed ratio of v orb /v wind ∼ /
4. We note, however, that the enclosed mass isnot gravitationally bound to the WD. As the WD ( R W D ∼ × cm) is not resolved inour simulations, we cannot determine the mass accreted by the WD.For r < − ∼ − P ∼ H p ∼ c s / Ω) in V407 Cyg, where Ω is the angularvelocity. Figures 5 and 6 show the gas column density and gas density-weighted temperaturedistribution as viewed from within the orbital plane. The density enhancement is onlyabout ∼ − ∼ ,
000 K. The opening angle ofthe spiral arm is about ∼ −
50 degrees. As described above, when the RG wind collideswith the WD, some gas periodically puffs out from the WD and forms the ring-like density 11 –distribution in Figure 5.Figure 7 and 8 show the angle-averaged density profiles of models M6 and M7 withinsurface of cones having different inclination angles. Each density profile is averaged from 180uniformly distributed radial rays in a given inclination angle and each radial ray is sampledfrom the surface of the RG to the boundary of the simulation box. From the upper panels,the averaged density profiles still follow the 1 /r distribution, but with about 20 − r < −
80 AU is highly concentrated in the region with θ < ◦ . Beyond this scale,the asymmetry is small and only affected at high inclination angles.We have also performed a low-resolution run by reducing the maximum AMR level to 8and enlarging the finest zone size to ∆ x min = 5 . × cm, corresponding to a factor of 2larger than the standard resolution. In this low-resolution run, the spiral arm qualitativelylooks the same as in the standard run, and the averaged density distribution looks the sameas well, but the turbulence is slightly suppressed due to a higher numerical viscosity. Theenclosed mass within the same control surface ( r < R c ) is about 20% higher than thestandard run.It should be noted that in our model we assume a constant RG wind speed during thewhole quiescent-phase simulation. However, in reality the RG wind is driven by radiationpressure on dust formed above the stellar surface during the Mira giant’s pulsation (Willson2000). The wind acceleration zone could extend to about 5 − R RG (Bowen 1988), whichis comparable to the binary separation in V407 Cyg. Therefore, the wind speed could beless than what we have assumed when it collides with the WD, implying that our estimateof the CDE could be low. 12 – After running for two orbital periods in the quiescent-phase simulations, we impose anova eruption on the WD. We perform two eruption simulations with two different eruptionenergies ( E = 1 . × erg and E = 1 . × erg) and ejecta masses ( M ej = 10 − M ⊙ and M ej = 10 − M ⊙ ) for each quiescent-phase simulation. Since we assume that theaveraged ejecta speed is h v ej i = 3 ,
000 km/sec and the kinetic energy ( E kin = M ej h v ej i ) isthree-fourths of the total eruption energy at the onset of eruption (Dohm-Palmer & Jones1996), a higher eruption energy implies a larger ejecta mass in that eruption, and vice versa.In total, we have four different sets of 3D eruption simulations as described in Table 1. InSection 4.2.1, we adopt model M6E44 as the reference simulation to describe the overallevolution in the eruption phase. We also discuss the shock evolution in all our models inSection 4.2.2. Figure 9 shows a typical gas density distribution for V407 Cyg in the orbital planeafter a nova eruption. In this case (model M6E44 in Table 1), the mass loss rate of the RGwind is M wind = 10 − M ⊙ yr − , the ejecta mass is M ej = 10 − M ⊙ , and the eruption energy is E = 1 . × erg. At about one week after the eruption, the forward shock sweeps out theRG wind and the spiral arm (label A), and forms a reverse shock propagating backward tothe WD (label B). The two shocks are separated by a contact discontinuity. At ∼
15 days,the forward shock reaches the RG. The propagation of the forward shock and the reverseshock can be seen in Figure 10. Note that although we assume the RG completely absorbsthe forward shock when it reaches the surface of the RG, it does not affect the propagationof the reverse shock from the eruption. 13 –Subsequently, Rayleigh-Taylor instabilities quickly develop (label C in Figure 9).After the forward shock passing through the RG, the RG wind and the CDE are shockedand heated (label D). These hot and dense plasma could be the source of X-ray emission(Orlando & Drake 2012). The self-interaction of Rayleigh-Taylor bubbles and the interactionof the forward shock with the spiral arm produce additional reverse shocks as well. Whenthese reverse shocks collide with each other, filament-like shocks form at about 70 daysafter the eruption (label E).After 70 days, the forward shock has slowed down because of the positive densitygradient toward the RG, forming a heart-shaped structure. Once the forward shock haspassed the RG, a high-density tail at the back side of the RG is formed (label F). Within theforward shock, several filaments are formed and destroyed repeatedly due to the interactionof the reverse shocks (label G). The turbulent spiral arm is completely destroyed by theforward shock. The corresponding gas temperature distribution is shown in Figure 11.It is clear that the temperature distribution is highly asymmetric due to the RG windinteraction and Rayleigh-Taylor instabilities. High Mach number shocks are associatedwith the forward shock, starting from M &
10 in the first week and then decreasing withtime. Once the Rayleigh-Taylor instabilities have developed, the Mach number decreases to M ∼ t = 137 day) a nova eruption is shown in Figure 12. The orange color shows the locationof the forward shock. Several filaments and Rayleigh-Taylor bubbles can also be seen.In Figure 13 and 14, we show the gas column density and gas density-weightedtemperature distribution in the r − z plane (edge-on view), to demonstrate the asphericalevolution of the shock radius due to the CDE. As we described in the previous subsection,the CDE makes the density ∼
20% asymmetric at different inclination angles. Therefore, 14 –although the ejecta speed is much higher than the orbital speed and the gravitationalbinding energy of the CDE is much less than the eruption energy, the ejecta speed greatlydecreases in the equatorial direction due to the CDE and the RG wind. The ratio of theshock radius in the poleward and equatorial directions is about 1 . − .
7, depending on thelocation of the RG.Orlando & Drake (2012) have performed simulations of the 2010 eruption of V407 Cygwith an artificial CDE. They found that the observed X-ray emission could be reproducedif there is a CDE around the WD, but is not required if the outburst energy and ejectamass are near the upper end of the range for these characteristics for classical novae,which would be extreme for SyRN. Comparing the results of their best-fitting simulationmodel with CDE, E44.3-NW7-CDE6.3-L40 with our results for model M7E44, we find thattheir forward shock expands more spherically and about twice as fast as that found here.Orlando & Drake (2012) found that most X-ray emission originates from the shocked CSM.The observed X-ray emission reaches a maximum at around 40 days after the eruption anddeclines after 50 days (Shore et al. 2011). This corresponds to the time when the RG windand the CDE are heated by the shock in our simulations (see label D in Figure 9). A lowershock expansion velocity in our simulations in comparison with Orlando & Drake (2012)may affect the shape of simulated X-ray lightcurve. However, we note that the initialaveraged ejecta velocity h v ej i ∼ ,
000 km s − in Orlando & Drake (2012) is probably toohigh based on recent observations (Abdo et al. 2010).Since the nova eruption completely destroys the CDE and creates a cavity thatcontains the binary system at the end of the simulation, the subsequent circumstellardensity distribution will need several years to decades to reestablish the original windprofile (Equation 1). The last recorded nova eruption of V407 Cyg prior to the 2010eruption was in 1998. Assuming that there were no undetected eruptions between the 1998 15 –and 2010 eruptions and the RG wind has a constant wind speed of 20 km s − , the newlydeveloped wind profile could only reach to about ∼
50 AU in distance. Therefore, ourinitial wind profile (Equation 1) may not be valid at large distances, which may affect ourshock evolution at late times ( t > r <
150 AU, butthe Rayleigh-Taylor instabilities are less obvious due to a higher numerical viscosity. Thefilaments produced by the interaction of reverse shocks are slightly less clear than in thestandard run, but overall the evolution does not change significantly.
Understanding the shock evolution in nova eruptions is crucial, since the high-energycharged particles required for the nonthermal gamma-ray and radio emission are likelyproduced in the shock front. Overall, the simulation results for the shock positions andevolution during the eruption phase runs are qualitatively similar, but quantitatively verydifferent. Therefore, in this subsection we describe the shock evolution of our four eruptionsimulations in detail.Figure 15 to Figure 18 show the shock evolution in all the eruption simulations. Atthe very beginning of a nova eruption ( t < ∼ t ∼
200 days(Figure 16). Similarly, with the same eruption energy, the shock would propagate muchfaster if the surrounding CDE were less dense (Figure 17).Nelson et al. (2012) and Martin & Dubus (2013) have described a simple semi-analyticalmodel for the shock evolution by comparing the circumstellar mass distribution with theejecta mass plus swept mass using momentum conservation. In their model, the densitydistribution is described by a 1D wind profile (Equation 1) plus a CDE. The CDE isassumed to be a gaussian distribution which can be described by ρ CDE ( r, θ ) = ρ , CDE exp − (cid:18) r sin θb (cid:19) − (cid:18) r cos θl (cid:19) ! , (2)where ρ , CDE , l , and b are density and length scale parameters. We perform similarcalculations by comparing the shock location of our 3D hydrodynamic simulations with thesemi-analytical model described by Martin & Dubus (2013).To find the most suitable sets of parameter values, we apply two constraints. First,we assume the ρ , CDE in model M6 is one order of magnitude higher than that in M7,i.e. ρ , CDE (M6) /ρ , CDE (M7) = 10, based on the analysis from Figure 7 and 8. Second,we assume the length-scale parameters l and b are equal, since we do not observesignificant asymmetry in the local density enhancement around the WD. By applyingthese two constraints, a best-fit set of parameter values is: for M7, l = b = 20 AU and ρ , CDE = 8 × − g cm − ; and for M6, l = b = 20 AU and ρ , CDE = 8 × − g cm − . Thebest-fit shock evolution is shown using the red dotted lines in Figures 15 - 18.The best fit model in Orlando & Drake (2012), which can reproduce the X-raylightcurve, corresponds to their model E44.3-NW7-CDE6.3-L40. It is characterized by anejecta mass of 2 × − M ⊙ , eruption energy E = 2 × erg, ρ , CDE ∼ × − g cm − ,and b = l = 40 AU. The CDE in this model is comparable to, in order of magnitude, with 17 –our model M7, but is less concentrated. However, as described above, the ejecta speed inOrlando & Drake (2012) is much faster than adopted in our work.Figure 19 shows the averaged forward shock evolution in the orbital plane and polewardregions of all considered models in Table 1. The shock radius evolution can be characterizedas r sh ∝ t α , where α is an evolution stage-dependent constant. Depending on the model, α is around 0 . − . α decreasea little for 10 < t <
100 days. After about 100 days, the forward shock is far beyond theRG, and α increases again. This behavior is consistent with the result in Walder et al.(2008), but for a different time scale due to differences in binary parameters. Furthermore,if we increase the RG mass loss rate (from model M7E44 to M6E44 or model M7E43 toM6E43) or decrease the eruption energy (model M6E43 to M6E44 or model M7E43 toM7E44), α becomes lower.The angle-averaged shock speed decreases greatly due to the sweeping out of the RGwind. Starting from a value of v ej ∼ ,
000 km s − at the onset of the nova eruption, theaverage shock speed in the radial direction drops to v ej ∼ −
300 km s − after about ayear. The evolution of shock speed can be characterized as another power-law relation withan index of − ∼ − (Figure 19). The speed evolution in Figure 19 is roughly consistentwith the spherically symmetric simulation of Moore & Bildsten (2012). In their simulation,a better treatment of the radiative cooling during shock expansion is implemented, but onlyin 1D. Our shock speed is also comparable with the observed ejecta speed (star symbols inFigure 19) of ∼ − on day +2 . ∼
200 km s − on day +196 from the FWHMof the broad component of Hα (Munari et al. 2011). By fitting the simulated ejecta speed 18 –with observations, the best-matched models are M6E44 and M7E43. Furthermore, the2011 April observation also shows that the maximum velocity is .
300 km s − (Shore et al.2011), which is also consistent with our simulations. The decrement of nova ejecta speedwith the existence of a CDE in our simulations is also consistent with the observed ejectavelocity decrement in RS Oph.
5. CONCLUSIONS
We have investigated the symbiotic recurrent nova V407 Cyg from the quiescent phaseto a nova eruption via three-dimensional hydrodynamical simulations. In quiescent-phasesimulations, we examined two different mass loss rates of the Mira-type RG. We studiedthe CDE produced due to the wind-WD interaction in V407 Cyg. It is found that theinduced spiral accretion wake in V407 Cyg is more turbulent and the CDE is less efficientthan found in RS Oph (Walder et al. 2008) due to a relatively high wind-to-orbital speedratio. In addition, we observe no significant wind-captured disk around the WD as reportedby Walder et al. (2008) in RS Oph, though the angle-averaged radial density profile in theequatorial plane is about 20% higher than that in the poleward direction.The shock evolution in the eruption phase also is studied. It is found that the forwardshock location is highly dependent on the inclination angle, nova eruption energy, andcircumstellar density distribution. The forward and reverse shock interactions with theRG wind and CDE are also crucial factors in the overall evolution. In addition, theshock radius evolution can be characterized by a power law, r sh ∝ t α , where α is about0 . − .
0, depending on model and evolution stage. Our model M6E44 ( M ej = 10 − M ⊙ and E ej = 1 . × erg) and model M7E43 ( M ej = 10 − M ⊙ and E ej = 1 . × erg)show good agreement with the observed ejecta speed from Munari et al. (2011). The shockevolution in the presence of CDE in our simulations is similar to what has been observed in 19 –V407 Cyg and RS Oph.Further work in this series of papers will include investigation of thermal and non-thermal emission in V407 Cyg and other symbiotic nova systems at different wavelengths,in particular in gamma-ray, X-ray, and radio. Similar analysis can be applied to other RNeand CNe as well.KCP would like to thank continuous support from F.-K. Thielemann. This workwas supported by the European Research Council (ERC) grant FISH, by PASC HighPerformance Computing Grant DIAPHANE, and by the Theoretical Institute for AdvancedResearch in Astrophysics (TIARA) in the Academia Sinica Institute of Astronomy andAstrophysics (ASIAA). PMR and RET acknowledge support by the NSF Division ofAstronomical Sciences under AAG 14-13367. FLASH was developed largely by theDOE-supported ASC/Alliances Center for Astrophysical Thermonuclear Flashes at theUniversity of Chicago. Analysis and visualization of simulation data were completed usingthe analysis toolkit yt (Turk et al. 2011). 20 – REFERENCES
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Science, 329, 817Ackermann, M., Ajello, M., Albert, A., et al. 2014, Science, 345, 554Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071Bode, M. F., O’Brien, T. J., Osborne, J. P., et al. 2006, ApJ, 652, 629Bowen, G. H. 1988, ApJ, 329, 299Chomiuk, L., Krauss, M. I., Rupen, M. P., et al. 2012, ApJ, 761, 173Chomiuk, L., Linford, J. D., Yang, J., et al. 2014, Nature, 514, 339Cheung, C. C., Jean, P., & Shore, S. N. 2014, The Astronomer’s Telegram, 5879,Colella, P., & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174De Beck, E., Decin, L., de Koter, A., et al. 2010, A&A, 523, AA18Dilday, B., Howell, D. A., Cenko, S. B., et al. 2012, Science, 337, 942Dohm-Palmer, R. C., & Jones, T. W. 1996, ApJ, 471, 279Drozd, K., Swierczynski, E., Ragan, E., & Buil, C. 2013, Advances in Astronomy and SpacePhysics, 3, A TEX macros v5.2. 23 –Table 1. Simulation Parameters
Abbreviation ˙ M a wind A b v c wind T d wind E e ej M f ej ( M ⊙ yr − ) (AU) (km sec − ) (K) (erg) ( M ⊙ )Quiescent phaseM6 † −
16 20 7,000 — —M7 10 −
16 20 7,000 — —Eruption phaseM6E43 10 −
16 20 7,000 1 . × − M6E44 † −
16 20 7,000 1 . × − M7E43 10 −
16 20 7,000 1 . × − M7E44 10 −
16 20 7,000 1 . × − † Reference simulation a RG mass loss rate b Binary separation c RG wind speed d Effective wind temperature e Eruption energy f Ejecta mass
24 –Fig. 1.— Gas density distribution in the orbital plane for model M6 at different labeledtimes. The color scale indicates the logarithm of the gas density in g cm − . 25 –Fig. 2.— Gas temperature distribution in the orbital plane for model M6 at different labeledtimes. The color scale indicates the logarithm of the gas temperature in K . 26 – -2 -1 Time (yr)10 -11 -10 -9 -8 -7 -6 M a ss ( M ⊙ ) c c c c c Model M6Model M7
Fig. 3.— Mass in the vicinity of the WD companion as determined for spherical controlsurfaces of different radii ( R c = 10 cm) as a function of time during the quiescent phase.Different line styles indicate different mass loss rates of the RG wind (see Table 1). 27 –Fig. 4.— 3D volume rendering of gas density around the WD for model M6. The whiteregion represents the disc-like accretion flow around the WD and the yellow region on theleft indicates the colliding shock front from the RG wind. The RG cannot be seen in thisfigure, but the bright region at the bottom shows the location of the RG. 28 –Fig. 5.— Gas density projection in the r − z plane perpendicular to the orbital plane formodel M6 at different labeled times, where the r -axis represents the axis passing throughthe RG and the WD. Note that the origin of the r − axis in this plot is at the center of mass,which is different from the simulation coordinate origin. In this figure, the RG is located onthe left. The color scale indicates the logarithm of the gas density in g cm − . 29 –Fig. 6.— Similar to Figure 5 but for gas temperature. 30 – -20 -19 -18 -17 -16 -15 -14 D e n s i t y ( g c m − ) Time = 82.5 yr θ = 90 ◦ θ = 60 ◦ θ = 45 ◦ θ = 30 ◦ θ = 0 ◦ Distance (AU)0.51.01.52.0 R a t i o Fig. 7.— Top: Averaged density profile of model M6 with different inclination angles as afunction of radius at 82 . -20 -19 -18 -17 -16 -15 -14 D e n s i t y ( g c m − ) Time = 82.5 yr θ = 90 ◦ θ = 60 ◦ θ = 45 ◦ θ = 30 ◦ θ = 0 ◦ Distance (AU)0.51.01.52.0 R a t i o Fig. 8.— Similar to Figure 7 but for model M7. 32 –Fig. 9.— Similar to Figure 1 but for the eruption phase of model M6E44. The labeled timesindicate the time after a nova eruption. Labels are important features that are described inSection 4.2. 33 – -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 D e n s i t y ( g / c m ) t = 6.9 day t = 15.5 day t = 26.7 day t = 333.9 day t = 398.7 dayDensityPressure -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 P r e ss u r e ( d y n e / c m ) -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 Fig. 10.— Density (solid lines) and pressure (dashed lines) profiles along the line extendingfrom the WD to the RG for model M6E44. Different colors represent different times afterthe RN eruption. The red shaded region shows the location of the RG. 34 –Fig. 11.— Similar to Figure 2 but for the eruption phase of model M6E44. 35 –Fig. 12.— 3D volume rendering of gas density for models M6 and M6E44 at different simula-tion phases. Left: In the quiescent phase (model M6), the labeled time indicates the simula-tion time since the start of the simulation. Right: In the eruption phase (model M6E44), thelabeled time indicates the simulation time after nova eruption. Orange colors in the rightpanel show the location of the nova shock. Red and yellow colors represent high-densityregions; green and deep blue colors denote low-density regions. (movies of simulations areavailable online). 36 –Fig. 13.— Similar to Figure 5 but for the eruption phase of model M6E44. The labeledtimes indicate the time after a nova eruption. 37 –Fig. 14.— Similar to Figure 13 but gas temperature. 38 – −150 −100 −50 0 50 100 150X (AU)−150−100−50050100150 Y ( A U ) This WorkMD13 ithout CDEMD13 ith CDE −150 −100 −50 0 50 100 150R (AU)−150−100−50050100150 Z ( A U ) This WorkMD13 ithout CDEMD13 ith CDE
Fig. 15.— Shock radius evolution. Left: in the orbital plane. Right: in the r − z planeperpendicular to the orbital plane. Blue contours indicate the shock radius evolution ofmodel M6E44 at different times. Black contours show a comparison of shock radius evolutionwith a simple analytical model (denoted by the label MD13) by Martin & Dubus (2013).Red dots represent the same simple analytical model, but with the CDE. 39 – −150 −100 −50 0 50 100 150X (AU)−150−100−50050100150 Y ( A U ) This WorkMD13 ithout CDEMD13 ith CDE −150 −100 −50 0 50 100 150R (AU)−150−100−50050100150 Z ( A U ) This WorkMD13 ithout CDEMD13 ith CDE
Fig. 16.— Similar to Figure 15 but for model M6E43. −200 −150 −100 −50 0 50 100 150 200X (AU) 200 150 100 50050100150200 Y ( A U ) This WorkMD13 without CDEMD13 with CDE −200 −150 −100 −50 0 50 100 150 200R (AU) 200 150 100 50050100150200 Z ( A U ) This WorkMD13 without CDEMD13 with CDE
Fig. 17.— Similar to Figure 15 but for model M7E44. 40 – −150 −100 −50 0 50 100 150X (AU)−150−100−50050100150 Y ( A U ) This WorkMD13 witho t CDEMD13 with CDE −150 −100 −50 0 50 100 150R (AU)−150−100−50050100150 Z ( A U ) This WorkMD13 witho t CDEMD13 with CDE