Early Hydrodynamic Evolution of a Stellar Collision
aa r X i v : . [ a s t r o - ph . H E ] N ov D RAFT VERSION A UGUST
28, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
EARLY HYDRODYNAMIC EVOLUTION OF A STELLAR COLLISION D ORON K USHNIR AND B OAZ K ATZ
Draft version August 28, 2018
ABSTRACTThe early phase of the hydrodynamic evolution following collision of two stars is analyzed. Two strongshocks propagate at a constant velocity (which is a small fraction of the velocity of the approaching stars) fromthe contact surface toward the center of each star. The shocked region near the contact surface has a planarsymmetry and a uniform pressure. The density vanishes at the (Lagrangian) surface of contact and the speed ofsound diverges there. The temperature, however, reaches a finite value, since as the density vanishes, the finitepressure is radiation dominated. For Carbon-Oxygen white dwarfs collisions this temperature is too low for anyappreciable nuclear burning at early times. The divergence of the speed of sound limits numerical studies ofstellar collisions, as it makes convergence tests exceedingly expensive unless dedicated schemes are used. Weprovide a new one-dimensional Lagrangian numerical scheme to achieve this. Self-similar planar solutions arederived for zero-impact parameter collisions between two identical stars, under some simplifying assumptions.These solutions provide rough approximations that capture the main features of the flow and allow a generalstudy as well as a detailed numerical verification test problem. The self-similar solution in the upstream frameis the planar version of previous piston problems that were studied in cylindrical and spherical symmetries. Wefound it timely to present a global picture of self similar piston problems. In particular, we derive new resultsregarding the non trivial transition to accelerating shocks at sufficiently declining densities (not relevant forcollisions).
Subject headings: hydrodynamics — self-similar — shock waves — supernovae: individual (Ia) INTRODUCTION
It was recently argued (Katz & Dong 2012) that as manyas ∼ of all stars may collide with each other duringtheir lifetime, due to the dynamics of typical field triple sys-tems. Especially interesting is that the rate of collision be-tween white dwarfs (WDs) in such systems may be as highas the rate of type Ia supernovae (SNe Ia). Although colli-sions of WDs were earlier believed to have rates which areorders of magnitude smaller than the rate of SNe Ia, they mo-tivated three-dimensional hydrodynamic simulations of suchcollisions and the possible resulting thermonuclear explosion(Benz et al. 1989; Raskin et al. 2009; Rosswog et al. 2009;Lor´en-Aguilar et al. 2010; Raskin et al. 2010; Hawley et al.2012). While the amount of Ni (the decay of which pow-ers the observed light, Colgate & McKee 1969) synthesized inmost of these simulations was non-negligible, the results werecontradictory, with inconsistent amounts of Ni and differentignition sites of a detonation wave for the same initial condi-tions. This discrepancy was resolved by Kushnir et al. (2013),where high resolution two-dimensional (2D) simulations witha fully resolved ignition process were employed . Moreover,it was shown that there is a strong correlation between the Ni yield and the total mass of colliding Carbon-Oxygen(CO) WDs (insensitive to their mass ratio) which spans theobserved range of SNe Ia yields for the observed range of COWDs masses. In all collisions the nuclear detonation is dueto a well understood shock ignition, devoid of the commonlyintroduced free parameters such as the deflagration velocityor transition to detonation criteria. The detonation triggered Institute for Advanced Study, Einstein Drive, Princeton, New Jersey,08540, USA John N. Bahcall Fellow We note that the recent simulations preformed by Garc´ıa-Senz et al.(2013) are at significant lower resolutions than those performed byKushnir et al. (2013). by the collisions results in explosions which match key obser-vational properties of SNe Ia. We believe that this is the mainchannel for these explosions.In this paper we analyze the early hydrodynamic evolu-tion of zero-impact parameter collisions between two identi-cal stars. Our results are applicable to a wide variety of stellarcollisions, and in particular they clarify the early evolutionin the case of collisions between CO WDs. We restrict ouranalysis to early times in which the velocity of the approach-ing stars, ± v , is roughly constant. The velocity is increasingsince the stars accelerate towards each other in the gravita-tional field of each star, g ≃ GM ⋆ /R ⋆ , where M ⋆ is thestellar mass and R ⋆ is the stellar radius. In what follows, welimit our analysis to early times, for which t ≪ t ≡ v /g .The approaching velocity v is much larger than the speedof sound near the stellar edge and immediately after contacttwo strong shock waves form that propagate from the con-tact surface towards the center of each of the stars. As weshow, the shock velocity in the collision frame, ˙ z sh , is muchsmaller than v for stellar collisions with ˙ z sh . v / . Inthe context of collisions between CO WDs this property leadsto a lack of any appreciable nuclear burning at early times, t . . v /g . At later times, when a detonation is ignited(typically at t ∼ t , see e.g. Kushnir et al. 2013), there is al-ready a significant amount of shocked material, which allowsan efficient synthesis of Ni, consistent with SNe Ia obser-vations. A detailed study of these detonations at later timesis beyond the scope of this paper and will be described in asubsequent publication (Kushnir & Katz in preparation).The paper is organized as follows. In Section 2 we discussthe properties of the contact region in an example of colli-sion between CO WDs. We show that the contact region hasa planar symmetry, a uniform pressure, a diverging speed ofsound, and a finite temperature. These main features do notdepend on the assumptions that the stars have equal masses,that the impact parameter is zero or on the specific densityprofile. The planar symmetry allows us to verify that the re-sults of a realistic 2D simulation of a collision (with a limitedresolution near the contact surface) are correct, by comparingit to a high resolution one-dimensional (1D) simulation. A nu-merical difficulty arises due to the diverging speed of soundand is overcome using a new 1D numerical scheme. In Sec-tion 3, we consider the ideal case of zero impact collisions ofidentical stars with a (pre-collision) power-law density pro-file, ρ = Kr − ω , and an ideal equation of state (the some-what confusing notations, r for the distance from the con-tact surface and ω < for increasing profile, are used to beconsistent with previous literature on propagating shocks incylindrical and spherical symmetries, see below). The exactself-similar solution of this problem as well as a simple andaccurate analytical approximation are derived. We find ex-cellent agreement between the 1D numerical results and theself-similar solutions, which shows that the self-similar solu-tions are achieved for this flow and that our numerical schemesolves the flow equations accurately.The self-similar solution in the upstream frame is the planarversion of previous piston problems that were studied in cylin-drical and spherical symmetries (Sedov 1945; Taylor 1946;Sedov 1959, and references therein). We found it timely topresent a global picture of self similar piston problems. InSection A, the solutions for planar, cylindrical and sphericalsymmetries, for all values of the density power law index ω are presented and compared. In particular, we derive new re-sults regarding the non trivial transition to accelerating shocksat sufficiently declining densities (not relevant for collisions)and point out interesting similarities and differences with thestrong explosion problem. THE CONTACT REGION HAS A PLANAR SYMMETRY, AUNIFORM PRESSURE, A DIVERGING SPEED OF SOUND, AND AFINITE TEMPERATURE
In this section we discuss the properties of the contact re-gion. We begin with an example of collision between COWDs (Section 2.1). We show that the contact region has aplanar symmetry, a uniform pressure, a diverging speed ofsound, and a finite temperature. These main features are gen-eral (Section 2.2), and do not depend on the assumptions thatthe stars have equal masses, that the impact parameter is zeroor on the specific density profile. The planar symmetry allowsus to verify that the results of a realistic 2D simulation of acollision (with a limited resolution near the contact surface)are correct, by comparing it to a high resolution 1D simula-tion (Section 2.3). A numerical difficulty arises due to thediverging speed of sound and is overcome using a new 1Dnumerical scheme. Finally, we show in Section 2.4 that forall CO WD collisions the temperature at the contact surfaceis generally too small for any appreciable nuclear burning atearly times t . . v /g . Example: early stages of the collision of two CO WDs
The early evolution of the collision of two . M ⊙ COWDs approaching with a Keplerian velocity of v = 2 . × km s − and with a zero impact parameter is shown inFigures 1 (density map) and 2 (density, pressure, speed ofsound, and temperature profiles along the axis of symme-try, x = 0 ). The Figures correspond to t = 0 . s ≃ . v /g ( g ≃ . × cm s − ). This collision was cal-culated by Kushnir et al. (2013) using high resolution 2DFLASH4.0 simulations with nuclear burning (Eulerian, adap- tive mesh refinement, 19 isotope alpha-chain reaction net-work, Dubey et al. 2009; Timmes 1999). The system ofequations is closed with the Helmholtz equation of state(Timmes & Swesty 2000) and a multipole gravity solver. Ini-tially the CO WDs are at contact with free fall velocities. Thestructure of each CO WD is obtained from an isothermal stel-lar model at T = 10 K and with a uniform composition of
Carbon and
Oxygen by mass. Since the
Helmholtz equation of state assumes complete ionization, the initial pro-file is not reliable for very low densities ρ ∼ < g cm − . Forsimplicity, we assume that the density profile at the lowestdensities is a power law density profile ρ ∝ d − ω , where d isthe distance from the stellar edge and ω ≃ − . determinedby a fit to the profile from the stellar model in the vicinity of ρ = 10 g cm − . The exact profile at low densities has a smallinfluence on our results, as discussed below.This particular problem has a cylindrical symmetry and theposition is described by cylindrical coordinates x (radius withrespect to the axis of symmetry) and z (distance from theplane parallel to the surfaces at contact). The fact that thestars are identical implies a mirror symmetry ± z allowing usto focus on one of the stars. We note that the main featuresthat are described below are not restricted to this 2D scenario.As can be seen, shocks are propagating into each of theidentical stars (reaching z sh ≈ . × cm at t = 0 . s)as evident by the jump in density, pressure, speed of sound,and temperature. The velocity of the shocks ˙ z sh ≃ z sh /t ≃ km s − is much smaller than v . As we show in Section 3,generally v / ˙ z sh ∼ > for stellar collisions. Several interest-ing features are apparent in these figures which are generic toearly phases of collisions. The evolution has planar symmetry at early times, t ≪ ( R ⋆ /v )( v / ˙ z sh ) , in the vicinity of the contact region, z, x ≪ R ⋆ — This is evident in figure 1 and results from the factthat the (cylindrical) radius of the contact region grows as x cont ∝ ( R ⋆ v t ) / while the depth of the shocked regiongrows as z sh ∝ ˙ z sh t implying that at early times the shockedregion is a thin disk with diverging aspect ratio x cont /z sh ∝ [( R ⋆ / ˙ z sh t )( v / ˙ z sh )] / ≫ . A quantitative illustration ofthe planar symmetry is provided in Figure 2, where the hy-drodynamic profiles from a 1D planner numerical scheme,described in Sections 2.3 and B, are shown to agree (to betterthan ) with the profiles along the axis of symmetry. Theinitial CO WD density of the 1D model equals to the densityon the axis of symmetry in the 2D model, and the initial ve-locity is the free fall velocity. The gravitational field is mim-icked by an adjustable acceleration, which is constant in timeand space. We choose to apply the surface acceleration of g , for which the close agreement between the two codes isfound. However, since we are interested in the early evolu-tion t ≪ v /g , our results are not sensitive to the exact valueof the applied acceleration. It is clear from Figure 2 that itis difficult to infer the behavior of the flow variables near thecontact surface directly from the 2D simulation because of thelimited resolution. The pressure is roughly uniform between the two shocks — This isa common feature of colliding mediums and is due to the shortsound crossing time z sh /c s compared to the evolution time t which allows the pressure to be evenly distributed. This, inturn, is true due to the fact that the sound speed is generically http://cococubed.asu.edu/code pages/adiabatic white dwarf.shtml F IG . 1.— A density map at t = 0 . s from a 2D FLASH4.0 simulation with ≃ km resolution of the zero-impact-parameter collision of 0.64-0.64 M ⊙ CO WDs, previously moving on a Keplerian orbit ( ρ ≡ ρ/ g cm − ).The contact surface is at z = 0 . The positions of the shocks on the symmetryaxis ( x = 0 ) are z sh ≈ ± . × cm. The profiles of the flow variableson the symmetry axis (dashed lines) are given in Figure 2. faster than the velocity in which the shock moves with respectto the shocked fluid in the downstream ( z < z sh ). In thescenario considered here the speed of sound is even larger atsmaller z as shown in Figure 2 and explained below. The speed of sound is diverging towards the contact surface — Given that the pressure is nearly uniform, this is directly re-lated to the fact that the density is decreasing. The reason thatthe density is decreasing toward the contact surface is thatthese mass elements are near the surface of the star where the(pre-shocked) density is approaching zero. As we next show,while each element is adiabatically compressed by the flow,the compression is not sufficient to compensate for its initiallow value.
The temperature reaches a finite value towards the contact surface — The vanishing density and the finite pressure near the con-tact surface imply that sufficiently close to the surface, thepressure is dominated by radiation, and therefore the temper-ature is given by T ∝ p / . As we show in Section 2.4, thistemperature is generally too small for any appreciable nuclearburning at early times. −1 z [cm] t = 0 . T c s, p ρ
2D @ x = 0Planar 1D z pr F IG . 2.— The profiles of the flow variables (temperature, speed of sound,pressure, and density) as function of the distance from the contact surface, z from the 2D FLASH4.0 simulation of colliding CO WDs shown in Figure 1(black, the profiles are shown for x = 0 , which correspond to the dashedlines in Figure 1). The hydrodynamic profiles from a 1D planner numericalscheme (red, see text for details) are shown to agree (to better than ) withthe 2D profiles. The outer position of the uniform pressure region at that time, z pr , is marked with red circles (the pressure at z < z pr is fixed to be uniformin the 1D simulation). The profiles are normalized as T ≡ T/ K, c s, ≡ c s / cm s − p ≡ p/ erg cm − , and ρ ≡ ρ/ g cm − . The speed of sound generally diverges at the contactregion
Consider a mass element which is in the vicinity of the con-tact region. It is useful to work with the column density m = Z z ρdz (1)as a Lagrangian coordinate, where m = 0 at the contact sur-face. The pressure of the element at the current time t is equalto the value of the pressure p ( t ) throughout the shocked regionwhere it is approximately uniform. Since the shock velocityat early time in the upstream frame is approximately constant(it is roughly v because for stellar collision v / ˙ z sh ∼ > , aswe show in Section 3), the pressure that the mass elementhad immediately after it was shocked is proportional to itspre-shocked density p , sh ( m ) ∝ ρ ( t = 0 , m ) ≡ ρ ( m ) .The increase in pressure since that time is p ( t ) /p , sh ( m ) ∝ p ( t ) /ρ ( m ) . Assuming that the adiabatic compression ofthe mass element behind the shock can be approximately de-scribed with an effective adiabatic index γ , the density andthe speed of sound of the mass element at the current time aregiven by ρ ( t, m ) ≃ ρ ( m ) (cid:18) p ( t ) p , sh ( m ) (cid:19) /γ ∝ p ( t ) /γ ρ ( m ) ( γ − /γ ,c s ( t, m ) ∝ p ( t ) ρ ( t, m ) ∝ (cid:18) p ( t ) ρ ( m ) (cid:19) γ − γ . (2)Since ( γ − /γ > (for fully ionized plasma / < γ < / ), the temperature diverges at the contact surface, m → ,where the pre-shocked density was vanishing. We emphasizethat the arguments leading to Eqeuation (2) do not depend onthe assumptions that the stars have equal masses, that the im-pact parameter is zero or on the specific density profile. Notethat the diverging speed of sound makes the approximation ofa uniform pressure more accurate at the contact region. Sincethe pressure is radiation dominated in the vicinity of the con-tact region, the temperature there, T c ( t ) , can be estimated as a T c ( t ) ≃ p ( t ) ≃ γ + 1 ρ ( t ) v ⇒ (3) T c ( t ) ≃ (cid:18) ρ ( t )4 × g cm − (cid:19) / (cid:18) v × km s − (cid:19) / K , where a is the black body radiation constant and ρ ( t ) is thepre-shocked density of the element being shocked at the time t . In the simple case that the initial density is a power law ρ ( t = 0 , z ) = Kz − ω ∝ m ω/ ( ω − (see section 3 for a de-tailed discussion), the pressure grows with time as p ( t ) ∝ t − ω . The density and the speed of sound profile near the con-tact surface scale as ρ ( t, m ) ∝ t − ω/γ m ω ( γ − γ ( ω − ,c s ( t, m ) ∝ t − ω ( γ − γ m − ω ( γ − γ ( ω − . (4)The finite (nonzero) values of the density and the speed ofsound at z = 0 which were obtained in the 2D simulation pre-sented in Figure 2 are the result of the finite resolution of thatsimulation. For higher resolution, the speed of sound (den-sity) at z = 0 increases (decreases). The diverging speed ofsound implies that pressure can be quickly equilibrated in thecontact plain and material is ejected in a thin layer parallelto the contact surface. The study of this “jet” is beyond thescope of this paper. A 1D numerical scheme that allows accuratecalculations
The planar geometry which is valid at early times allow theflow to be accurately solved for using a very high resolution1D code. In particular, a fully Lagrangian code can be appliedwhich allows the possible ignition of detonation process to bestudied (Kushnir et al. 2013; Kushnir & Katz in preparation).However, the compression of the elements close to the con-tact surface, combined with the fact that the speed of soundis diverging, significantly limits the time step allowed by theCurrant condition. One way around this is to use Eulerianschemes for the flow near the contact surface, and restrict thesize of the cells there. The disadvantage of this method is thenumerical “smearing” of the flow variables near the contactsurface on the scale of the cell size ∆ z (which, for example,limits the maximal temperature there). Here we describe a 1DLagrangian scheme, which allows to overcome the Curranttime-step condition constrain without limiting the size of theinnermost cells.As discussed above, the mass elements near the contact sur-face equalize their pressure efficiently due to the high speedof sound. The scheme uses this fact by approximating thepressure of a chosen Lagrangian region, < m < m pr ,in the vicinity of the contact surface m = 0 , to be strictlyuniform. In the outer regions, m > m pr , the hydrodynamicequations are solved as usual. The evolution of the hydro-dynamic profile in the region < m < m pr is completely determined by the evolution of the value of the pressure inthe region p pr ( t ) . Indeed, the density of each mass elementdepends on p pr ( t ) , and the position of each element can beexpressed as z ( m ) = R dmρ − . The pressure p pr ( t ) is deter-mined by expressing the position of the edge of the uniformpressure region (which is determined similarly as in the outerregion) in terms of the profiles in the region as z ( m pr ) = Z m pr dmρ − . (5)A discretized numerical scheme which implements such aregion is described in Section B. In the numerical calcula-tions presented below, this scheme was used. Nuclear burn-ing can be added in a straight forward way and is described inKushnir & Katz (in preparation).The choice of the (Lagrangian) position of the boundary m pr is updated throughout the simulation to account for thegrowing region of nearly uniform pressure. In practice, a cellis added to the uniform pressure region if its Currant time-stepbecomes smaller than an appropriate threshold (as comparedto outer cells) and convergence is verified by changing thethreshold. The validity of this code is confirmed below whereit is compared to the exact solutions which are obtained in thecase of an ideal gas and power-law density profiles (Figure 3). Lack of significant nuclear burning
The highest temperature in the vicinity of the contact sur-face are still too low to allow a significant nuclear burning.In particular, the ignition of a detonation is postponed to laterstages of the evolution. This is a general feature of all COWDs collisions. In the extreme case of M ⊙ − M ⊙ colli-sion, where v ≃ . × km s − , g ≃ . × cm s − ,and t ≃ . s, the upstream density at . t is ≃ × g cm − , leading to a temperature of ≃ . × K (seeEquation (3)) in the vicinity of the contact surface, which istoo low for any significant nuclear burning. Note that this es-timate derived under the condition v / ˙ z sh ≫ , which weshow below to generally hold for stellar collisions. SOLUTION IN THE SIMPLIFIED CASE OF AN IDEAL GAS AND APOWER LAW DENSITY PROFILE
We showed in Section 2 that at early times, where the veloc-ity of the approaching stars is roughly constant, the problemhas planar geometry. At these times the pressure in the pre-shocked region is negligible, since close to the surface of thestar the speed of sound decreases significantly and is muchsmaller than v , which is larger or equal to the free fall ve-locity and comparable to the typical speed of sound in thestar. By adopting the following approximations, the problemis significantly simplified and allows an exact solution: • The equation of state is that of an ideal gas, with anadiabatic index γ . This is not strictly correct for theshocked region, as near the shock the equation of stateis close to that of an ideal gas ( γ = 5 / ), while nearthe contact surface the pressure is radiation dominated( γ = 4 / ). • The density distribution has a power law dependenceon the distance from the contact surface, ρ ( t = 0 , z ) = Kz − ω . (6)While this is not strictly correct near the surfaceof WDs, by choosing appropriate values for ω = − d log ρ/d log z , the obtained solutions approximatethe profiles at any given region. This approximationis also useful in collisions of other type of stars, withconcrete examples including a radiative envelope with ω = − or efficiently convective envelope (or degener-acy pressure) with ω = − / (Chandrasekhar 1939).The solution describes the position of the shock wave, z sh ( t ) and the hydrodynamical profiles between z = 0 and z = z sh ( t ) at any given time t after contact. In particular, we showthat v / ˙ z sh & for all γ and ω values that represent stellarcollisions.Since there are only three dimensional variables in theproblem, v , K , and t , it is reasonable to assume that anydimensional quantity is given by the appropriate combina-tion of these (up to a dimensionless multiplication factor)and the solution is self-similar. The profiles depend on thenon-dimensional parameters ω , γ and normalized location z/z sh ( t ) . The self-similarity allows the hydrodynamic partialequations to be reduced to an ordinary differential equationwhich can be easily (numerically) solved. The self-similaritysolution is derived in Section 3.1. The equations are solved inthe frame of the upstream fluid, where the problem is equiv-alent to a piston moving into stationary fluid, and the resultsare transformed back to the laboratory frame of the stellar col-lisions. By comparing the obtained solutions to the results ofdirect 1D simulations of the same problem, the validity ofthe exact solutions and of the 1D numerical scheme is vali-dated. In addition to the exact solution, a simplified analyticmodel is derived by approximating the pressure to be exactlyuniform throughout the shocked region and is shown to pro-vide an excellent approximation (Section 3.2). In Section A,the solutions of the piston problem for planar, cylindrical andspherical symmetries, for all values of the density power lawindex ω are presented and compared. Exact self similar solution
We solve the problem in the upstream frame, where it be-comes identical to a planar piston moving with velocity v into a stationary, cold fluid. The position of a fluid elementlocated at a distance z from the piston ( = contact surface)is given in this frame by r = z + v t where the initial po-sition of the piston is at r = 0 . The initial density profileis ρ ( t = 0 , r ) = Kr − ω . The position of the shock in thisframe is given by R sh ( t ) = z sh ( t ) + v t and its velocity is ˙ R sh ( t ) = ˙ z sh ( t ) + v . The local fluid velocity in this frameis u = u z + v . The notations r and R sh are used to beconsistent with the discussion in Section A which includescylindrical and spherical coordinate systems.Before solving the hydrodynamic equations, we note thatthe velocity of the shock is exactly constant. Most simply,by using dimensional analysis, R sh can only be constructedfrom v t (we discuss the validity of the dimensional analysisin Section A.2). Note that the shock cannot decelerate due tothe zero velocity boundary condition at the piston and it can-not accelerate in the increasing density profile due to the lim-ited energy budget at any given time (more details are givenin Section A.2).Using dimensional arguments, it is possible to show thatin the case where the flow is independent of any characteristiclength scale, the flow fields must be of a self-similar form (e.g.Zel’dovich & Raizer 1966; Waxman & Shvarts 2010) whichwe choose to be u = ˙ R sh ξU ( ξ ) , c s = ˙ R sh ξC ( ξ ) , ρ = KR − ω sh G ( ξ ) , (7) where u ( r, t ) , c s ( r, t ) , and ρ ( r, t ) are the fluid velocity, soundspeed, and density, respectively, and ξ ( r, t ) = r/R sh ( t ) = r/ ( ˙ R sh t ) (8)is the similarity parameter which is unity at the shock posi-tion. Note that the density scales like R − ω sh since (for strongshocks) the density just behind the shock wave is a constantfactor, ( γ + 1) / ( γ − , times the pre-shocked density justahead of the shock which is given by KR − ω sh . The value of ξ at the piston’s position, r p = v t , is given by ξ p = v / ˙ R sh . (9)Finally, the pressure is given by p = ρc /γ = KR − ω sh P ( ξ ) ˙ R /γ , where P ( ξ ) = G ( ξ ) ξ C ( ξ ) .Using Equations (7) and (8) , the hydrodynamic equations,Equations (A2) (with n = 1 representing the planar case),can be expressed as a single ordinary differential equation,Equation (A4), dUdC = ∆ ( U, C )∆ ( U, C ) , and one quadrature, Equation (A5), d ln ξdU = ∆( U, C )∆ ( U, C ) or d ln ξdC = ∆( U, C )∆ ( U, C ) , where ∆ , ∆ , and ∆ are given by ∆ = C − f , ∆ = U f − C (cid:18) U − ωγ (cid:19) , ∆ = C (cid:20) f − C + ( γ − ω γ C f (cid:21) , (10)where f = 1 − U. (11)The boundary condition at the piston position is u ( t, r p ) = v which can be expressed as U ( ξ p ) = 1 . (12)The Rankine–Hugoniot relations at the shock front determinethe boundary conditions for the self-similar solutions to be(e.g. Zel’dovich & Raizer 1966) U (1) = 2 γ + 1 , C (1) = p γ ( γ − γ + 1 , G (1) = γ + 1 γ − . (13)As illustrated here and in Section A (see also Guderley1942; Meyer-Ter-Vehn & Schalk 1982; Waxman & Shvarts1993), many of the properties of the self-similar flows maybe inferred by analyzing the contours in the ( U, C ) -plane de-termined by Equation (A4). Numerical integration of Equa-tions (A4) shows that for solutions starting at the strong shockpoint, Equations (13), C diverges close to the piston as U ap-proaches 1 (see Figure 6 for a representative case), implyingthat the speed of sound is diverging near the piston as ex-pected. We next analyze the behavior of the solution near thepiston point ( U, C ) = (1 , ∞ ) . Equation (A4) is given, toleading order in f , by d ln fd ln C = ν fC ⇒ f ∝ C ν , (14)where ν = 2 ( γ − ω ) ω ( γ − < . (15)The quadrature, Equation (A5), gives to leading order in ff = (cid:18) − ωγ (cid:19) ln (cid:18) ξξ p (cid:19) . (16)Using these results and Equation (A6) we find G ∝ f ω ( γ − / ( ω − γ ) , (17)and by using Equation (14), we find that P ( ξ p ) is finite(nonzero) implying that the pressure near the piston is notvanishing or diverging, as expected. We may now determinethe dependence of the density and the speed of sound near thepiston on mass and time. For a given mass element, we applyEquation (16) to its trajectory ξ (which coincides with a C characteristic), drdt = u ⇒ d ln ξ d ln R = − f ( ξ ) , (18)to get f ∝ ln (cid:18) ξ ξ p (cid:19) ∝ R − (1 − ω/γ )sh ∝ t − (1 − ω/γ ) . (19)At a given time, the mass scales as m ∝ ρξ . Using thisresult and Equation (19) with Equation (17) we find ρ ∝ t − ω/γ m ω ( γ − /γ ( ω − , and with Equation (14) we find c s ∝ t − ω ( γ − /γ m − ω ( γ − /γ ( ω − , the same dependence that wasderived in Equation (4). Note that since the pressure, density,and velocity are finite for the whole flow, the energy containedin the self-similar solution diverges as t − ω when t → ∞ ,in accordance with the work done by the piston on the gas.This proves the consistency of the self-similar solution, as de-scribed with more details in Section A.2.The solutions in the collision frame can be expressed by thesolutions in the upstream frame with the relation zz sh = r − v tR sh − v t = ξ − ξ p − ξ p . (20)Note that in the collision frame the self-similar solution de-scribes the whole space between the shock and the contactsurface ( ≤ z/z sh ≤ ). The density, speed of sound andpressure are presented for two cases ( γ = 5 / , ω = − / and γ = 4 / , ω = − ) in Figure 3. While the pressure inthe shocked region is roughly uniform, the speed of sound di-verges near the contact surface and the density vanishes there.The (normalized) shock velocity, ˙ z sh v = ˙ R sh − v v = ξ − p − , (21)is provided for some values of γ and ω in Table 1 and in Fig-ure 4. As can be seen, for cases relevant to stellar collisions( ω ≤ − / ), the shock velocity is a small fraction of v .This fraction increases with ω , up to a value of ( γ − / for ω = 0 (the solution of the corresponding Riemann problem).The behavior for ω > is discussed in Section A.6. An ana-lytic approximate expression for this velocity, Equation (24),is derived below by approximating the pressure to be exactlyuniform and is in excellent agreement with the exact result.Next we compare the self-similar solution to the results ofthe direct D numerical simulations described in section 2.3 −4 −2 z/z sh c s /c s,d ρ/ρ d p/p d γ = 5 / , ω = − / −4 −2 z/z sh c s /c s,d ρ/ρ d p/p d γ = 4 / , ω = − Self-similarNumerical z pr /z sh F IG . 3.— The self-similar profiles (speed of sound, pressure, and density)in the collision frame, normalized by their values immediately behind theshock front (denoted by a subscript d ), are presented in black solid lines fortwo cases: γ = 5 / , ω = − / (left) and γ = 4 / , ω = − (right).While the pressure in the shocked region is roughly uniform, the speed ofsound (density) diverges (vanishes) near the contact surface. The profiles ofthe flow variables, from a numerical calculation, at the time the shock arrivedto the edge of the computational grid, are presented in red dashed lines. Thenumerical profiles agree with the self-similar profiles to better than in therange z/z sh > − . The outer positions of the uniform pressure region atthat time, z pr /z sh , are marked with red circles. The transition between theuniform pressure region and the regular region is smooth, and the numericalsolution agrees with the self-similar solution in both regions. The excellentagreement between the results validates both the self-similar solution and thenumerical scheme. TABLE 1S OME VALUES OF ˙ z sh /v = ξ − p − ASFUNCTION OF ω AND γγ ω = − ω = − / / / for the same ideal gas equation of state and density profiles.The initial mesh consists of N cells with a uniform spacing, ∆ z . The initial pressure was chosen such that the outgoingshock wave is always strong. The shock trajectories for tworepresentative cases, calculated with N = 8000 , are shownin Figure 5. The agreement of the numerical trajectories withthe self-similar trajectories is better than . . The profiles ofthe flow variables, at the time the shock arrived to the edge ofthe computational grid, are shown in Figure 3. The numericalprofiles agree with the self-similar profiles to better than in the range ¯ ξ > − . The outer positions of the uniformpressure region at that time, z pr /z sh , are marked with circles.The transition between the uniform pressure region and theregular region is smooth, and the numerical solution agreeswith the self-similar solution in both regions. The excellentagreement between the results validates both the self-similarsolution and the numerical scheme. −10 −8 −6 −4 −2 0 110 −2 −1 ω ˙ z s h / v = ξ − p − γ = 5 / γ = 4 / p ( γ − / F IG . 4.— The shock velocities in the collision frame (normalized by v )as function of ω for γ = 4 / , / are shown in black. For cases relevant tostellar collisions ( ω ≤ − / ), the shock velocity is a small fraction of v .This fraction increases with ω , up to a value of ( γ − / for ω = 0 (thesolution of the corresponding Riemann problem). The behavior for ω > is discussed in Section A.6. The simple analytic solution given by Equa-tion (24), which is derived under the assumption that the pressure is uniform,is shown in red. The analytic expression provides an excellent approximationto the exact solution and is accurate to better than for / < γ < / and − < ω < − . t/ (∆ z/v ) z s h / ∆ z γ = 5 / , ω = − / γ = 4 / , ω = − F IG . 5.— The shock trajectories in the collision frame are presented for twocases: γ = 5 / , ω = − / (solid line) and γ = 4 / , ω = − (dashedline). The trajectories from a numerical calculation are presented in circles(the initial mesh consists of cells with a uniform spacing, ∆ z ). Theagreement of the numerical trajectories with the self-similar trajectories isbetter than . . . Simple analytic solution
We next provide a simple analytic solution by using the factthat the pressure is nearly uniform. Assuming that the pres-sure is exactly uniform between the contact surface and theshock, we can use the same arguments that lead to equation(2), to express the velocity of the fluid element immediatelyat the downstream of the shock. The pressure p pr ( t ) in theshocked region is given by its value in the immediate down-stream and grows with time as p pr ( t ) ∝ ρ ( t ) ∝ t − ω . Thedensity of a given mass element in the downstream grows withtime as it is adiabatically compressed, ρ ( t, m ) ∝ p pr ( t ) /γ ∝ t − ω/γ . The size of each element thus shrinks according to dx/dm = ρ − ∝ t ω/γ . This implies that the entire La-grangian region between and m scales in the same way z ( t, m ) ∝ t ω/γ and thus the fluid velocity of each element satisfies u z ( m ) = ∂z ( t, m ) ∂t = ωγ z ( t, m ) t . (22)Applying this to the element which is immediately down-stream of the shock we obtain the following equation: u z | sh = ωγ ˙ z sh . (23)By using the strong shock jump condition, ˙ z sh − u z | sh = ( γ − / ( γ + 1)( ˙ z sh + v ) , we can solve for ˙ z sh in terms of v , ˙ z sh v = γ ( γ − γ − ω ( γ + 1) . (24)As can be seen in Figure 4, equation (24) provides an ex-cellent approximation to the exact solution and is accurate tobetter than for / < γ < / and − < ω < − .This expression does not capture the growing shock velocity(normalized by v ) as ω approaches since the pressure issignificantly non-uniform in the downstream region at thesevalues. The hydrodynamic profiles under this approximationare exact power-laws with indexes as in Equation (4), and am-plitudes set by the shock jump conditions. SUMMARY AND DISCUSSION
The early phase of the hydrodynamic evolution followingthe collision of two stars is analyzed, focusing on the regionnear the contact region. It was shown in Section 2 that theshocked region has a planar symmetry, a uniform pressure,and a diverging (vanishing) speed of sound (density) whenapproaching the contact surface (Equation (2)). The tempera-ture reaches a finite value towards the contact surface (Equa-tion (3)), since the vanishing density and the finite pressurenear the contact surface imply that the pressure is dominatedby radiation. We showed in Section 2.4 that for all CO WDscollisions this temperature is generally too small for any ap-preciable nuclear burning at early times t . . v /g , beforethe velocity increases due to the gravitational acceleration. Inparticular, the ignition of a detonation is postponed to laterstages of the evolution. This is tightly related to the fact thatthe shock moves very slowly in the collision frame comparedto the fast approach speed v (see Figure 4 and Table 1).The divergence of the speed of sound has an important con-sequences for numerical studies of the stellar collisions. Thenumerical Currant condition will require a rapidly decreasingtime step for higher resolutions, making convergence tests ex-ceedingly expensive unless dedicated schemes are used. Wedescribed in Sections 2.3 and B a new 1D Lagrangian numer-ical scheme to achieve this.We provided self-similar planar exact solutions for the sim-plified case of a power-law density profile and an ideal equa-tion of state in Section 3.1. These solutions provide roughapproximations that capture the main features of the flow andallow a general study as well as a detailed numerical verifica-tion test problem. Finally, we derived an approximate analyticexpression for the shock velocity (Equation (24)) which is ac-curate to a few precent over a wide range of density profiles.We thank S. Dong, E. Waxman, and E. Livne for useful dis-cussions. D. K. gratefully acknowledges support from MartinA. and Helen Chooljian Founders’ Circle. FLASH was in partdeveloped by the DOE NNSA-ASC OASCR Flash Center atthe University of Chicago. Computations were performed atPICSciE and IAS clusters. APPENDIXA. SELF-SIMILAR PISTON DRIVEN FLOWS
The problem of a piston propagating into a medium with a power law density profile which is studied in Section 3.1 (in theupstream frame) is the planar version of previous piston problems that were studied in cylindrical and spherical symmetries(Sedov 1945; Taylor 1946; Sedov 1959, and references therein). We found it timely to present a global picture of such selfsimilar piston problems. In Section A.1 we write down the hydrodynamic equations of the flow for planar, cylindrical andspherical symmetries ( n = 1 , , , respectively) along with the resulting ODEs assuming self similarity. The solutions for allgeometries and all values of the density power law index ω are presented in the sections that follow. Particular emphasis is givento the non trivial transition to accelerating shocks at sufficiently declining densities, for which we derive new results and point outinteresting similarities and differences with the strong explosion problem. A detailed discussion of the solutions for the planarcase with growing density profiles ( ω < ), which are relevant for stellar collisions, is presented in Section 3.1. A.1. The equations describing self-similar flows
In the problems considered, a piston moves with a constant velocity v into a medium with an ideal gas equation of state withan adiabatic index γ , and an initial power-law density profile ρ ( t = 0 , r ) = Kr − ω , (A1)where r is the radial coordinate. The medium is assumed to have zero pressure initially and a strong shock propagates ahead ofthe piston. The equations describing the adiabatic 1D flow behind the shock are (e.g. Landau & Lifshitz 1987) ( ∂ t + u∂ r ) ln ρ + r − ( n − ∂ r ( r n − u ) = 0 , ( ∂ t + u∂ r ) u + ρ − ∂ r ( γ − ρc s ) = 0 , ( ∂ t + u∂ r )( c s ρ − γ ) = 0 , (A2)where n = 1 , , are for planar, cylindrical, and spherical symmetry, respectively.Self similar solutions are obtained by substituting Equations (7) and (8) in the hydrodynamic Equations (A2), with a shockvelocity scaling (the shock velocity is not constant in general), ˙ R sh = AR δ sh . (A3)The partial differential equations, Equations (A2), are replaced with a single ordinary differential equation (Zel’dovich & Raizer1966; Waxman & Shvarts 2010), dUdC = ∆ ( U, C )∆ ( U, C ) , (A4)and one quadrature d ln ξdU = ∆( U, C )∆ ( U, C ) or d ln ξdC = ∆( U, C )∆ ( U, C ) . (A5)The normalized density, G , is given implicitly by ( ξC ) − n − ω ) | − U | λ G ( γ − n − ω )+ λ ξ nλ = const (A6)with λ = ω ( γ −
1) + 2 δ. (A7)The functions ∆ , ∆ , and ∆ are ∆ = C − (1 − U ) , ∆ = U (1 − U )(1 − U − δ ) − C (cid:18) nU + 2 δ − ωγ (cid:19) , ∆ = C { (1 − U )(1 − U − δ ) − γ − U [( n − − U ) + δ ] − C + 2 δ + ω ( γ − γ C − U } . (A8) A.2. The self similarity assumption
There are two types of similarity solutions (see e.g. Zel’dovich & Raizer 1966). Following Kushnir & Waxman (2010), solu-tions of the first-type may be defined as solutions that are valid over the entire ( r, t ) -plane (or the part of which where the flowtakes place). Such solutions must satisfy the global conservations laws of mass, momentum, and energy, and hence the valuesof the similarity exponents of such solutions may be determined by dimensional considerations. Solutions of the second-typemay be defined as solutions which describe only part of the flow. Such solutions should be required to allow the existence of acharacteristic line, ξ c ( R sh ) , along which the self-similar solution is matched to another solution, and to comply with the globalconservation laws within the region of the ( r, t ) -plane described by the self-similar solution (note that it is commonly acceptedthat the similarity exponents of a second-type solution are determined by the requirement that the solution passes through asingular point of the hydrodynamic equations, but this condition is not general enough, see Kushnir & Waxman 2010).Before solving the hydrodynamic equations, we can use simple arguments to derive some properties of the self-similar solu-tions. To begin with, the shock must propagate with a constant velocity ( δ = 0 ) or accelerate ( δ > ), since if it is decelerating,the piston reaches it at some finite time. Next, let us assume that for the cases where the mass near the piston is finite ( ω < n ),the self-similar solution is valid everywhere between the piston and the shock (first-type solution). Below we show that thisassumption results in a consistent solution. For ω ≥ n , the mass near the piston diverges and there are no consistent self similarsolutions of the entire flow. We discuss second-type self-similar solutions for this regime in Section A.6. In the first-type case,the shock must propagate at a constant velocity, as can be derived from a few arguments. Most simply, by using dimensionalanalysis, R sh can only be constructed from v t . Another argument concerns the mass element adjacent to the piston. Such anelement is part of the self-similar flow and therefore its normalized position, ξ p = v t/R sh ( t ) , must be constant, implying that δ = 0 . Moreover, since its position coincides with a C characteristic of the self-similar solution, given by Equation (18), wemust have U ( ξ p ) = 1 . By using Equation (7), we can derive the shock velocity, ˙ R sh = v /ξ p , which also equals the constant A .A more physical argument for the constant velocity of the shock can be made by considering the total energy of the flow. Theenergy contained in the self-similar solution is E s ( R sh ) = f ( n ) Z R sh ξ p R sh drr n − (cid:18) ρu + 1 γ − p (cid:19) = f ( n ) A R δ + n − ω sh K [ I k ( ξ p ) + I i ( ξ p )] , (A9)with I k ( ξ ) = Z ξ dξ ′ ξ ′ n +1 G U , I i ( ξ ) = Z ξ dξ ′ ξ ′ n +1 G γ ( γ − C , (A10)and f ( n ) = for n = 1 , π for n = 2 , π for n = 3 . (A11)The I k and I i terms describe the kinetic and internal energy of the gas, respectively. Since at any given shock position R sh theenergy of the gas must be finite (nonzero), I k and I i cannot diverge and at least one of them is nonzero. Therefore, E s ( R sh ) diverges as R δ + n − ω sh when R sh → ∞ . The energy of the gas is supplied from the work done on it by the piston, W ( R sh ) ∝ t ( R sh ) Z dtv ( v t ) n − R δ − ω sh P ( ξ p ) ∝ P ( ξ p ) R δ + n/α − ω sh . (A12)In order for the work done by the piston to diverge in accordance with the energy of the gas, we must have α = 1 ( δ = 0 ). For thework done by the piston to be finite (nonzero), P ( ξ p ) must be finite (nonzero), and therefore the pressure near the piston behavesas p ( t ) ∝ t − ω . Note that if δ = 0 for n ≥ ω , then the energy of the gas is not increasing as R sh → ∞ , which is not physicalgiven that the piston is performing work on the gas. This immediately shows that for n ≥ ω the self-similar solution cannotbe valid near the piston, and only second-type self-similar solutions are possible, with the possibility that the shock accelerates, δ > (see Section A.6).It is straightforward to generalize Equation (4) to the general case ω < n , by noting that ρ ( t = 0 , r ) = Kr − ω ∝ m ω/ ( ω − n ) .The density and the speed of sound profiles near the contact scale as ρ ( t, m ) ∝ t − ω/γ m ω ( γ − γ ( ω − n ) ,c s ( t, m ) ∝ t − ω ( γ − γ m − ω ( γ − γ ( ω − n ) . (A13)These results can be verified directly from the asymptotic behavior near the piston, derived below. The density near the pistonvanishes (diverges) and the speed of sound there diverges (vanishes) for ω < ( ω > ). For ω = 0 both density and speed ofsound are finite (non-zero) near the piston (for n = 1 this is a simple Riemann problem, while for n = 2 , the solutions werederived by Sedov 1945; Taylor 1946; Sedov 1959).In Sections A.3, A.4, and A.5 we show for ω < , ω = 0 , and < ω < n , respectively, that I k and I i are finite, such that E s ( R sh ) diverges as R n − ω sh when R sh → ∞ , in accordance with the work done by the piston on the gas. This showes that thefirst-type self-similar solutions are consistent. The derived shock velocity is shown for some values of γ and ω in Figures 4 and 7for the n = 1 and n = 2 , cases, respectively.0 U C ω < ω = 0 ( n > < ω < n Sonic lineStrong shock point(
U, C ) = (1 , C p ) F IG . 6.— Different types of the C ( U ) curves for the solutions of the piston problem, obtained for ω values in the regimes ω < (solid), ω = 0 (dashed), and < ω < n (dot-dashed). Also shown as a dotted line is the sonic line, U + C = 1 . The circle marks the strong shock point, Equations (13), and the squaremarks the point ( U, C ) = (1 , C p ) . A.3. Increasing density profile ( ω < ) Numerical integrations of Equations (A4) starting at the strong shock point, Equations (13), indicate that C diverges as U approaches 1 (see Figure 6 for a representative case). Analysis of the the behavior of the solution near ( U, C ) = (1 , ∞ ) showsthat f ∝ C ν , where ν = 2 γ ( n − ω/γ ) ω ( γ − < . (A14)The quadrature, Equation (A5), gives to leading order in ff = (cid:18) n − ωγ (cid:19) ln (cid:18) ξξ p (cid:19) . (A15)Using these results and Equation (A6) we find G ∝ f ω (1 − γ ) / ( nγ − ω ) . (A16)Therefore, I k and I i are finite. A.4. Constant density profile ( ω = 0 ) For planner symmetry ( n = 1 ) this is a simple Riemann problem, and in what follows we consider the cylindrical and spherical( n = 2 , , respectively) cases (Sedov 1945; Taylor 1946; Sedov 1959). Numerical integration of Equations (A4) shows thatfor solutions starting at the strong shock point, C reaches a finite (nonzero) value C p as U approaches 1 (see Figure 6 for arepresentative case). Note that ( U, C ) = (1 , C p ) is not a singular point of Equations (A4), but nevertheless the integrationends there, as the boundary condition for the piston is U ( ξ p ) = 1 . We next analyze the the behavior of the solution near ( U, C ) = (1 , C p ) . To leading order in f , Equation (A4) is given by dfdC = − nC p ⇒ f = n (cid:18) − CC p (cid:19) . (A17)1The quadrature, Equation (A5), gives to leading order in ff = n ln (cid:18) ξξ p (cid:19) . (A18)Using these results and Equation (A6) we find that G has some finite (nonzero) value near the piston. Therefore, I k and I i arefinite. A.5. Moderately decreasing density profile such that the mass near the piston is finite ( < ω < n ) Numerical integration of Equations (A4) shows that for solutions starting at the strong shock point, C reaches zero as U approaches 1. We next analyze the the behavior of the solution near ( U, C ) = (1 , . We consider the planar case seperately( n = 1 ) since its analysis is somewhat different than that of the cylindrical and spherical cases ( n = 2 , ). A.5.1 Planar symmetry ( n = 1 ) Equation (A4) is given, to leading order in f , by d ln fd ln C = − f − C (cid:16) − ωγ (cid:17) f + ω ( γ − γ C . (A19)Equation (A19) implies that lim f → d ln fd ln C = ν, (A20)where either ν = 1 or ν = 2( γ − ω ) /ω ( γ − > . The quadrature, Equation (A5), gives to leading order in ff = θ ln (cid:18) ξξ p (cid:19) , θ = ( (cid:16) − ωγ (cid:17) , ν > , ( γ − γ +1 , ν = 1 . (A21)Using these results and Equation (A6) we find G ∝ f − µ , µ = (cid:26) ( γ − ω/ ( γ − ω ) , ν > , (( γ + 1) ω − / ( γ − , ν = 1 . (A22)Numerical integration of Equations (A4) and (A5) shows that solutions starting at the strong shock point approach ( U, C ) = (1 , along a ν > curve. However, for both ν = 1 and ν > , I k and I i are finite.As can be seen in Figure 4, the ratio between the shock velocity and v diverges in the limit where ω → . We note that P ( ξ p ) approaches zero in this limit. A.5.2. Cylindrical and spherical symmetries ( n = 2 , ) In this case, equation (A4) is given, to leading order in f , by d ln fd ln C = − f − C (cid:16) n − ωγ (cid:17) − ( γ − n − f + ω ( γ − γ C . (A23)The solution of Equation (A23) must also be of the form given by Equation (A20). Assuming f tends to 0 slower than C , i.e., ν < , leads to a contradiction since Equation (A23) gives ν = 2 / ( γ − n − , which is larger than for γ < and n = 2(3) . Therefore, ν must satisfy ν ≥ . For ν > , Equation (A23) gives ν = 2( nγ − ω ) ω ( γ − , (A24)which satisfies ν > for ω < n . For ν = 1 , Equation (A23) gives f = 2 nγ − ( γ + 1) ω γ − γ ( γ − n − C . (A25)The solution of the quadrature, Equation (A5), gives f = θ ln (cid:18) ξξ p (cid:19) , θ = ( (cid:16) n − ωγ (cid:17) , ν > , n ( γ − γ +1 , ν = 1 . (A26)Using these results and Equation (A6) we find G ∝ f − µ , µ = (cid:26) ( γ − ω/ ( nγ − ω ) , ν > , (( γ + 1) ω − n ) /n ( γ − , ν = 1 . (A27)2Numerical integration of Equations (A4) and (A5) shows that solutions starting at the strong shock point approach ( U, C ) = (1 , along a ν > curve. However, µ < for both ν = 1 and ν > , and thus I k and I i are finite.As can be seen in panel (d) of Figure 7, ˙ R sh /v is some finite value (which depends on γ ) at ω = 2(3) for n = 2(3) (unlikethe planner case, for which ˙ R sh /v diverges as ω → ). A.6. Steeply decreasing density profile such that the mass near the piston diverges ( ω ≥ n ) We already commented, based on the energy contained in the first-type self-similar solution, that for ω ≥ n the self-similarsolution cannot be valid near the piston, and only second-type self-similar solutions are possible. In what follows, we derive thisdirectly from the properties of the self-similar solutions, assuming δ = 0 .Numerical integration of Equations (A4) for n = 1 shows that for solutions starting at the strong shock point, Equations (13), C ( U ) crosses the sonic line ( U + C = 1 ) at a non-singular point (either ∆ = and/or ∆ = 0 ). Therefore, δ = 0 results ina non-physical solution. For n = 2 , , the integration shows that in the range n ≤ ω ≤ ω g ( γ, n ) ( ω g is increasing with γ , ω g (1 , n ) = n and ω g (5 / , n ) ≃ . . for n = 2(3) , see Waxman & Shvarts 1993, for the n = 3 case), known as the“gap” region, C reaches zero as U approaches 1. For ω values above the “gap”, C ( U ) crosses the sonic line at a non-singularpoint. Thus, we only need to show that in the “gap”, δ = 0 results in a non-physical solution. This is achieved by examiningEquation (A27), which shows that µ > for both ν = 1 and ν > , and therefore I k diverges in the limit ξ → ξ p (this was shownfor n = 3 in Kushnir & Waxman 2010).The second-type self-similar solutions were already found for the strong explosion problem, in which a large amount of energyis deposited within a small region at the center of an initially cold gas with an initial density ρ = Kr − ω , for ≤ ω ≤ ω g ( γ, by Kushnir & Waxman (2010) and for ω g ( γ, < ω by Waxman & Shvarts (1993, 2010). The equations for the strong explosionproblem and for the piston problem are identical for second-type self-similar solutions (since δ is uniquely determined by theexistence of such a solution, Kushnir & Waxman 2010). Given that the self-similar part of the flow does not depend on thedetails of the non-self-similar part of the flow (one case involving an explosion, following which each mass element moves atan asymptotically constant velocity and the second case involving a piston moving at a constant velocity), these second-typeself-similar solutions are also the solutions for the piston problem.Generalizing the results of Waxman & Shvarts (1993, 2010) to n = 1 , and of Kushnir & Waxman (2010) to n = 2 is straightforward, and the details are not given here. The results are that in the “gap” the shock propagates at a constant velocity ( δ = 0 )and ξ c ( R sh ) is a C characteristic. Note that in this case the solution complies with the global conservation laws, since nearthe piston the flow deviates from the self-similar flow and the arguments given above for the divergence of I k do not hold. For ω > ω g ( γ, n ) the shock accelrates ( δ > ), ξ c ( R sh ) is a C + characteristic, and the solution passes thorough a singular point. Theself-similar exponent δ for γ = 4 / and γ = 5 / as function of ω and n is plotted in Figure 7 for both the piston problem andthe strong explosion problem.The behavior of the solution in the “gap” region shows some interesting properties. Since δ = 0 and ξ c ( R sh ) approachesthe piston (see Equation (18)), we can infer the shock velocity ˙ R sh = v /ξ p (there is an analog for v for the strong explosionproblem, since for ω > n most of the mass is near the explosion center and acquires a typical velocity). As can be seen in panel(d) of Figure 7, ˙ R sh /v is some finite value (which depends on γ ) at ω = ω g ( γ, n ) . In other words, the transition to acceleratingshocks happens from a finite ratio of ˙ R sh /v . Analyzing the solutions near the transition, we find that the transition happens asthe causal connection between the piston and the shock is lost. This is natural, since for δ > there cannot be such a causalconnection (as explained in the beginning of the appendix), and indeed these solutions include a sonic point. The loss of thecausal connection happens differently for different geometries. For n = 1 , the pressure at the piston P ( ξ p ) approaches zero as ω approaches ω g ( γ, n ) . For n = 3 the pressure at the piston approaches some finite (nonzero) value at these limit, and the lossof connection is due to some ξ a ( γ, n ) > ξ p , for which U ( ξ a ) + C ( ξ a ) = 1 (and therefore the time for a sound wave to cross thispoint is infinite). B. NUMERICAL 1D PLANAR LAGRANGIAN SCHEME
The cells are separated into two groups based on their proximity to the contact surface. Cells i = 1 , , ...i pr closer than achosen Lagrangian point m pr ,which is chosen as one of the nodes (node i pr between cell i pr and i pr + 1 ), are assumed to have auniform pressure p pr and are treated separately from the rest of the cells. Other cells, i > i pr , are advanced by a standard scheme.During each time step, properties of the cells i ≤ i pr are calculated as follows: for a given (unknown) next step value of thepressure p pr ( t + ∆ t ) , the new densities ρ i ( t + ∆ t ) can be calculated for each cell using the equation of state. Using the constantcell masses, m i , and new densities ρ i ( t + ∆ t ) , the total length of the region can be calculated. By comparing this length to thatobtained from the velocity and force equations on the border node i pr , an algebraic equation is obtained for p pr ( t + ∆ t ) , whichis solved iteratively. We next write down this equation explicitly. Fore simplicity, we ignore nuclear burning which is straightforward to incorporate.We assume that the thermodynamic variables ρ i ( t ) and c s,i ( t ) are defined at the cells’ centers. The force equation for node i pr can be solved similarly as to a regular node ( i > i pr ), and therefore r i pr ( t + ∆ t ) is known. The task is to solve for all othervariables of the uniform pressure region at the time t + ∆ t . The length of the region (or equivalently the position of the node i pr )at t + ∆ t is given by z pr ( t + ∆ t ) = i pr X i =1 ∆ z i ( t + ∆ t ) = i pr X i =1 m i ρ i ( t + ∆ t ) , (B1)3 ω δ n = 1 Firsttype Secondtype δ ST δ piston γ = 5 / γ = 4 / (a) ω δ n = 2 Firsttype Secondtype δ ST δ piston ω g (4 / , ≃ . ω g (5 / , ≃ . γ = 5 / γ = 4 / (b) ω δ n = 3 Firsttype Secondtype δ ST δ piston ω g (4 / , ≃ . ω g (5 / , ≃ . γ = 5 / γ = 4 / (c) −10 −8 −6 −4 −2 0 1 2 3 410 −2 −1 ω ˙ R s h / v − = ξ − p − n = 2 n = 3 γ = 5 / γ = 4 / (d)F IG . 7.— Panels (a)-(c): the self-similar acceleration exponent δ (where ˙ R sh ∝ R δ sh ) as a function of the density profile index ω (where ρ ( t = 0 , r ) ∝ r − ω )for γ = 5 / (solid line) and γ = 4 / (dashed line). Panel (a): planar symmetry ( n = 1 ). Panel (b): Cylindrical symmetry ( n = 2 ). Panel (c): Sphericalsymmetry ( n = 3 ). For the strong explosion problem and ω < n , the self-similar exponent satisfies δ ST = ( ω − n ) / (Sedov 1946; von Neumann 1947; Taylor1950). For n ≤ ω ≤ ω g ( γ, n ) , and n = 2 , , the self-similar exponent satisfies δ = 0 (Kushnir & Waxman 2010), and for ω g ( γ, n ) < ω the self-similarexponent was found by Waxman & Shvarts (1993). Panel (d): the shock velocity (normalized by v ) as function of ω for γ = 5 / (solid line) and γ = 4 / (dashed line) for n = 2 , . The inset zooms in around the “gap” region to make it clear that the shock velocity is finite there (and crosses smoothly the point ω = n , shown as dotted line). and similarly i pr X i =1 m i ρ i ( t ) = r i pr ( t ) , (B2)leading to r i pr ( t + ∆ t ) − r i pr ( t ) = i pr X i =1 m i (cid:18) ρ i ( t + ∆ t ) − ρ i ( t ) (cid:19) ≃ − i pr X i =1 m i ( ρ ) i ( t + ∆ t/
2) ( ρ i ( t + ∆ t ) − ρ i ( t )) , (B3)In the uniform pressure region the flow is isentropic, and therefore ρ i ( t + ∆ t ) − ρ i ( t ) = p ( t + ∆ t ) − p ( t )( c s ) i ( t + ∆ t/ (B4)4holds for every cell there. Using this with Equation (B3) we get p ( t + ∆ t ) = p ( t ) − r i pr ( t + ∆ t ) − r i pr ( t ) P i pr i =1 m i ( ρ c s ) i ( t +∆ t/ . (B5)Equations (B5) and (B4) are solved with iterations p k ( t + ∆ t ) = p ( t ) − r i pr ( t + ∆ t ) − r i pr ( t ) P i pr i =1 m i ( ρ c s ) i,k − ( t +∆ t/ ,ρ i,k ( t + ∆ t ) = ρ i ( t ) + p k ( t + ∆ t ) − p ( t )( c s ) i,k − ( t + ∆ t/ , (B6)supplemented by the equation of state c s,i,k ( t + ∆ t ) = c s ( p k ( t + ∆ t ) , ρ i,k ( t + ∆ t )) (B7)and by ρ i,k ( t + ∆ t/
2) = 12 ( ρ i,k ( t + ∆ t ) + ρ i ( t )) , (cid:0) c s (cid:1) i,k ( t + ∆ t/
2) = 12 h(cid:0) c s (cid:1) i,k ( t + ∆ t ) + (cid:0) c s (cid:1) i ( t ) i , (cid:0) ρ c s (cid:1) i,k ( t + ∆ t/
2) = 12 h(cid:0) ρ c s (cid:1) i,k ( t + ∆ t ) + (cid:0) ρ c s (cid:1) i ( t ) i , (B8)where k is the number of the iteration and for the initial guess, k = 0 , we use ρ i, ( t + ∆ t/
2) = ρ i ( t ) , (cid:0) c s (cid:1) i, ( t + ∆ t/
2) = (cid:0) c s (cid:1) i ( t ) , (cid:0) ρ c s (cid:1) i, ( t + ∆ t/
2) = (cid:0) ρ c s (cid:1) i ( t ) . (B9)Typically, a few iterations are sufficient for convergence.We implemented this scheme in the 1D, Lagrangian version of the VULCAN code (for details, see Livne 1993).(B9)Typically, a few iterations are sufficient for convergence.We implemented this scheme in the 1D, Lagrangian version of the VULCAN code (for details, see Livne 1993).