A Comparison of Grid-Based and SPH Binary Mass-Transfer and Merger Simulations
Patrick M Motl, Juhan Frank, Jan Staff, Geoffrey C Clayton, Christopher L Fryer, Wesley Even, Steven Diehl, Joel E Tohline
AA Comparison of Grid-Based and SPH Binary Mass-Transfer and MergerSimulations
Patrick M. Motl , Juhan Frank , Jan Staff , Geoffrey C. Clayton , Christopher L. Fryer , WesleyEven , Steven Diehl , , and Joel E. Tohline ABSTRACT
Currently there is great interest in the outcomes and astrophysical implications ofmergers of double degenerate binaries. In a commonly adopted approximation, thecomponents of such binaries are represented by polytropes with an index n = 3 /
2. Wepresent detailed comparisons of stellar mass-transfer and merger simulations of poly-tropic binaries that have been carried out using two very different numerical algorithms— a finite-volume “grid” code and a smoothed-particle hydrodynamics (SPH) code. Wefind that there is agreement in both the ultimate outcomes of the evolutions and theintermediate stages if the initial conditions for each code are chosen to match as closelyas possible. We find that even with closely matching initial setups, the time it takesto reach a concordant evolution differs between the two codes because the initial depthof contact cannot be matched exactly. There is a general tendency for SPH to yieldhigher mass transfer rates and faster evolution to the final outcome. We also presentcomparisons of simulations calculated from two different energy equations: in one serieswe assume a polytropic equation of state and in the other series an ideal gas equationof state. In the latter series of simulations an atmosphere forms around the accretorwhich can exchange angular momentum and cause a more rapid loss of orbital angularmomentum. In the simulations presented here, the effect of the ideal equation of stateis to de-stabilize the binary in both SPH and grid simulations, but the effect is morepronounced in the grid code.
Subject headings: binaries: close, hydrodynamics, methods: numerical, white dwarfs,AM CVn, Type Ia supernovae, R CrB Indiana University Kokomo, School of Sciences, P.O. Box 9003, Kokomo, IN, 46903-9004, USA; [email protected] Louisiana State University, Department of Physics & Astronomy, 202 Nicholson Hall, Baton Rouge, LA 70803-4001, USA Department of Astronomy, University of Florida, Gainesville, FL, 32611, USA Computational Methods Group (CCS-2), Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM,87545, USA Theoretical Astrophysics Group (T-6)/Computational Methods Group (CCS-2), Los Alamos National Labora-tory, P.O. Box 1663, Los Alamos, NM 87545, USA a r X i v : . [ a s t r o - ph . S R ] F e b
1. Introduction
Dynamically unstable mass transfer and stellar mergers involve complex hydrodynamical flowsof material through a varying gravitational potential. Investigating these rapid phases of stellarevolution requires the solution of Poisson’s equation (for the self-gravity of the binary) coupled tothe equations of hydrodynamics to predict the outcome of the mass transfer event. The complexityof these coupled equations necessitate a numerical solution, of which several possible computationaltechniques are available – each with their own strengths, weaknesses and numerical artifacts.In this paper we investigate mass transfer, tidal disruption and merger events in binariesusing two very different numerical codes and we compare the results obtained. We follow thecommonly adopted approximation in which the two stars are represented by n = 3 / (cid:46) § § § § §
2. Numerical Algorithms
From a computational science point of view, our separate groups have shared the goal ofdeveloping computational fluid dynamic (CFD) algorithms that can be used to accurately simulatebinary mass-transfer events in a variety of astrophysical settings. We treat both stars as continuumfluid structures whose Newtonian gravitational fields and tidal interactions are modeled in a fullyself-consistent fashion. We also employ numerical schemes that are designed to spatially resolve theinterior structures of both stars as well as faithfully represent the complex, nonlinear dynamicalflows that result from mass exchange between the two stars. The SPH technique serves as thefoundation of the numerical code designed by the “LANL” group; the numerical code designed bythe “LSU” group — hereafter, referred to as the “grid” code — employs an Eulerian grid and usesa finite-volume discretization technique. The distinguishing attributes of both codes are describedin § P is derived from the mass density ρ via the relation, P = κ i ρ γ , (1)where the adopted adiabatic exponent is γ = and the polytropic constant κ i , which specifiesthe specific entropy of a fluid element, can have one of two values depending on whether the fluidelement is originally part of the donor star ( i = D) or part of the accretor ( i = A). (Values of κ A and κ D , are reported in Table 1 for the initial models considered here.) For “ideal gas” simulations,we use the relation, P = ( γ − (cid:15)ρ , (2)where γ = and where, at each integration time step, the specific internal energy of the gas, (cid:15) , isdetermined from the PDE that governs energy conservation. As noted in Table 1, this equation ofstate allows us to simulate low mass (compared to the Chandrasekhar mass) double white binarieswhen both binary components are on the same adiabat. In other cases, the binary components willhave the mass – radius relation appropriate to low-mass, fully convective stars.Both codes employ mechanisms for handling the formation and movement of shock discontinu-ities. Via the energy equation, for example, some kinetic energy is converted to heat when the fluidpasses through a shock. Hence, the specific entropy of a fluid element should also increase whenthe fluid passes through a shock. In practice, this additional heat will be realized by an increase in (cid:15) that is above what would be expected through adiabatic compression alone. When the ideal gas 5 –EOS is used, the gas pressure also will reflect this heat gain because, via Eq. (2), P is proportionalto (cid:15) . However, when the polytropic EOS is used, this additional heat is effectively lost from thesystem because the specific entropy of each fluid element is forced to remain constant, independentof whether or not the fluid has passed through a shock. By evolving otherwise identical systems,first, with a polytropic (P) EOS then, second, with an ideal gas (I) EOS, the set of simulationspresented here not only permits us to compare results from different algorithmic techniques butalso to isolate the effects that shock heating has on binary mass-transfer evolutions.Results from CFD simulations can depend sensitively on initial conditions. With this in mind,we have made a concerted effort to start each simulation from a well-defined, equilibrium binaryconfiguration, and to start each pair of comparison simulations from initial configurations that arenearly identical. The initial binary configuration for each of the 23 simulations identified in Table2 was drawn from one of the equilibrium configurations listed in Table 1. By design, the modelsspan a range of initial binary mass ratios, 1 . ≥ q = M D /M A ≥ . q > q crit ) to smaller initial mass ratios where masstransfer may be stable depending on the fate of the angular momentum of the accretion stream. Forour chosen polytropic exponent of γ = , the standard critical mass ratio for dynamical stability is q crit = (Rappaport et al. 1982; Soberman et al. 1997), but could be significantly lower for directimpact accretion (e.g. q crit = 0 .
22 Marsh et al. 2004; Gokhale et al. 2007). Table 1 includes theradii of each binary component as well as the accretor’s radii in units of the circularization radius, R circ (where R circ is calculated from equation 4.20 of Frank et al. 2002), showing that all of ourinitial models will experience direct impact accretion.The initial models identified in Table 1 as Q1 .
3, Q0 .
5, and Q0 . a were selected for this in-vestigation because their polytropic simulations had already been studied using the grid code (seeDMTF06 and MFTD07). In each of these cases the detailed structure of the initial equilibriumconfiguration was derived using an SCF technique as described in § § q = 0 . q crit = 2 /
3, and filled agap in parameter values that would otherwise exist between models Q1 . .
5. The standardSCF technique was used to construct the model identified as Q0 . a in Table 1 as the additional ini-tial model for the grid code. But, in this case, a separate relaxation technique was used to generatethe model identified in Table 1 as Q0 . b for insertion into the SPH code (See § . a and Q0 . b are very similar but not identical to one another; for example κ A (cid:54) = κ D inQ0 . a . An additional SCF model, Q0 . b , at higher numerical resolution, was also constructed. Asdetailed in Table 1, the polytropic models Q0 . a and Q0 . b differ slightly in their parameters dueto the difference in numerical resolution used for the SCF code.As constructed, the donor star did not quite fill its Roche lobe in each initial model. Afterintroducing an initial model into one of the two CFD codes, we “drove” the donor into contact 6 –with its Roche lobe by explicitly extracting angular momentum from the binary system at rate ofone percent per orbit for a specified “driving duration.” (The precise manner by which angularmomentum was extracted from the binary system in the two separate CFD codes is described in § We use the SCF technique developed by Hachisu (Hachisu 1986; Hachisu et al. 1986) to generateour initial models. This approach allows us to form semi-detached equilibrium structures for thebinary components that are synchronously rotating assuming a polytropic equation of state. Theseequilibrium structures naturally incorporate the effects of tidal distortion in their self-consistentRoche geometry as well as rotational flattening of the stars.In a frame of reference co-rotating with the binary, the fluid equations reduce to a singlescalar equation for each binary component that equates the sum of the gravitational potential,centrifugal potential and enthalpy to an integration constant. Three boundary points are set alongthe line of centers where the mass density is forced to vanish and these boundary points allow usto algebraically solve for the integration constant for each star and the angular frequency of therotating frame.We iteratively solve the scalar equations starting from spherical Gaussian density distribu-tions and the Keplerian estimate for the orbital frequency. Upon solving Poisson’s equation (on acylindrical mesh identical to the one used in the grid-based CFD code) we have an estimate of theenthalpy field. For our chosen polytropic equation of state, the enthalpy determines the next guessat the matter density where we rescale the density to maintain prescribed maximum density valuesfor each star and an improved estimate of the centrifugal potential is computed about the currentcenter of mass of the binary. Poisson’s equation is solved again and new values for the integrationconstants and orbital frequency are computed from the boundary conditions.The iteration proceeds until the structure changes by an amount smaller than a given tolerancein each subsequent iteration. The converged value of the orbital frequency, Ω , is used to initializethe velocity field in our binary simulations and to define the angular frequency of the grid code’srotating frame of reference. The equilibrium solution is not exact however and the orbit will not beprecisely circular but rather exhibit epicyclic oscillations about the mean binary separation evenif artificial driving is not imposed. For the models constructed in this paper, the amplitude of ˙ aa ,where time is measured in the initial orbital period of the binary, ranges from ∼ .
02 to almost0 . P ≡ π/ Ω ,as quantified in Table 1, is used to normalize evolutionary times. In the grid-based simulations,the converged SCF model is shifted slightly within the grid, via an interpolation scheme, before it 7 –is introduced into the CFD code in order to better ensure that the center of mass of the binarysystem initially coincides with the origin of the cylindrical coordinate system. For SPH simulations, the particle distribution has to be chosen a priori . Since we want topay particular attention to resolving the Roche lobe overflow accretion stream, we need sufficientresolution (i.e. particles) in the outer layers of the donor. Since the SCF method yields asymmetricstellar configurations, we developed a new setup technique for SPH initial conditions that canhandle arbitrary spatial configurations. This technique is based on weighted Voronoi tessellation(WVT) and is discussed in greater detail in Diehl et al. (2012). For the sake of completeness, wehere discuss the essentials as well as details specific to the SCF setup problem.The WVT initial condition setup is based on an approach developed by Diehl & Statler (2006)to adaptively bin two-dimensional astronomical X-ray images to a constant signal-to-noise ratio.We extend this technique to work in arbitrary dimensions, and optimize the algorithm to workwith an underlying, arbitrary probability distribution instead of images. The algorithm starts outwith an initial guess for the particle distribution, which we choose to be random positions drawnfrom the probability distribution. The WVT setup then iteratively applies artificial short-range(typically within 2-3 smoothing lengths) forces to the particles until they settle into a smoothparticle distribution. The strengths of these forces are determined by the desired input resolutionat the particle positions. This process generally converges rather quickly ( <
40 iterations) andresults in a setup with very low particle-noise, and the desired particle distribution. For uniformdistributions, the configuration resembles that of gravitational glasses, which are often employedin cosmological simulations (e.g. Diehl et al. 2012; Wang & White 2007).In the context of this investigation, we apply a uniform particle distribution for the accretor. Ina high-resolution SPH simulation, each accretor is modeled with a total of 150,000 SPH particles.In order to properly resolve the accretion stream and beat down the noise in the mass transferrate due to a small number of SPH particles, we decided to resolve the outer layers of the donorstar with much higher precision. The resolution at the center of the donor star is chosen tomatch approximately that of the accretor. However, the smoothing lengths in the outer layers areapproximately 4 times smaller, putting effectively ∼
64 more resolution elements in these layers,adding up to a total of 850,000 SPH particles in the donor star and, hence, approximately 1,000,000particles for each full simulation. Figure 1 shows a slice through the meridional plane of the WVTSPH setup for initial model Q0 . a , as defined in Table 1. For clarity, only particles within onesmoothing length of the meridional plane are plotted.The SPH density interpolation of the particle masses in the WVT method generally reproducesthe SCF grid setup to better than 1%, as shown in Figure 2 where we show mass profiles forthe stellar components in the initial states. Despite the fact that the stars are not spherically 8 –symmetric, we construct the mass profile as the mass within spherical shells centered on the centerof mass of each stellar component. Integral quantities for the initial states including the totalangular momentum, orbital angular momentum and the spin angular momentum of each binarycomponent agree to a level of 1% or better for models Q0 . a , Q0 . .
3. For model Q0 . b ,a separate, simpler relaxation method was used starting with two perfect spheres placed at theexpected orbital separation. Centrifugal and tidal distortions were then introduced adiabaticallyallowing the components to relax to equilibrium. Our finite-volume “grid” code is an explicit code that is second order accurate in both time andspace and has been used extensively to simulate self-gravitating astrophysical flows (Motl et al. 2002;D’Souza et al. 2006; Motl et al. 2007; Staff et al. 2012). The code utilizes cylindrical coordinatesand evolves equations for the conservation of mass, radial, vertical, and angular momentum and anenergy equation, as described in more detail below. With our choice of coordinates and formulationof the fluid equations we are guaranteed that the total angular momentum of the system will notchange due to advection or due to gravitational or pressure forces.The hydrodynamics treatment is similar to that in the Zeus code (Stone & Norman 1992);we use advection operators based on Van Leer slopes and a similar numerical viscosity to stabilizethe fluid equations in the presence of shocks. The fluid equations are formulated for a rotatingreference frame (including source terms for the associated Coriolis and centrifugal forces) and weset the code’s reference frame to rotate at the initial orbital frequency of the binary to minimizethe effects of numerical diffusion through the grid. Further details of the hydrodynamics algorithmare available in Motl et al. (2002) and D’Souza et al. (2006).For polytropic simulations our energy equation is cast in terms of an entropy tracer τ = ( (cid:15)ρ ) /γ (3)where (cid:15) is the specific internal energy of the fluid and ρ is the mass density. The entropy tracerevolves with an advection equation similar to the continuity equation for polytropic simulations.For simulations with an ideal gas equation of state we evolve the total energy density, E , with ∂E∂t + ∇ · v E = −∇ · p v − ρ v · ∇ Φ eff (4)where Φ eff includes both the gravitational and centrifugal potentials. We compute the internalenergy of the fluid by subtracting the kinetic energy density of the fluid from the total energydensity.The fluid equations evolve densities of conserved quantities like mass density or momentumdensity but updating the fluid state requires additional quantities that don’t obey conservation 9 –laws like the velocity of the fluid. We are thus continually dividing by the mass density of the fluidand therefore no cell in the grid can have a vanishing mass density. To prevent division by zero,we floor the mass density at every timestep and for the runs presented here the floor density is ρ min = 1 × − compared to a maximum stellar density of ∼ − . The entropy tracer, τ , is floored similarly to ρ in polytropicevolutions to remain above a value τ min = (cid:18) κ min γ − (cid:19) γ ρ min (5)where κ min is the smaller of the polytropic constants in the binary system. For simulations withan ideal gas equation of state, the total energy density E is floored to remain above the smallestinternal energy density, τ γ min .Poisson’s equation is solved by performing a Fourier analysis along the periodic azimuthaldirection and an alternating-direction implicit technique is used for the radial and vertical directions.The boundary value of the gravitational potential is computed via a numerical integration of themass density with cylindrical Greens functions based on half-integer degree Legendre polynomialsof the second kind (Cohl & Tohline 1999; Cohl et al. 2001). We employ outflow boundary conditionsat the outer edge of the computational domain and over the course of a merger simulation, up to ∼ .
3% of a system’s total mass can be lost from the active grid.As was alluded to above, each binary simulation was begun with a phase of artificial “driving”that causes the orbit to slowly decay and forces the donor star into contact with its Roche lobe.This establishes a mass transfer stream that can be well resolved by the code, given the chosennumerical resolution. In the grid code the phase of driving was accomplished by removing a fractionof the specific angular momentum of each fluid element each time step of the simulation at a ratesuch that one percent of the initial angular momentum is removed per initial orbital period, P .A typical numerical simulation with the grid code (at a numerical resolution of 162 radialzones by 98 vertical zones by 256 azimuthal zones for a total of approximately 4 million cells aslisted in Table 2) required about 24 days of run time on 64 cores of the Louisiana Optical NetworkInitiative (LONI) machine QueenBee, assuming a polytropic equation of state. With the ideal gasequation of state, the sound speed is higher which further reduces the Courant-Friedrichs-Lewy(CFL) limited time step and increases the typical run time to approximately 150 days. The highersound speed for ideal gas equation of state simulations results from a hot, low density envelopeof shock heated material that forms around the accreting star. Note that, since we are using acylindrical coordinate grid, the CFL limiting constraint on the time step is often reached near thecoordinate axis and is thus significantly smaller than would be the case with a Cartesian coordinategrid. 10 – The SPH code used in this work is called SNSPH, and is described in much greater detail inFryer et al. (2006). Here we will give a global overview of the code, focusing on the aspects mostimportant to this project.SNSPH is a tree-based, massively parallel SPH code that has undergone rigorous testing andhas been applied to a large variety of astrophysical problems (e.g. Fryer & Warren 2002; Young et al.2006; Diehl et al. 2008). The gravitational forces are evaluated in an order N log N gravitationaltree, where long-range forces are replaced by monopole terms. Our multipole acceptance criterionis described in Warren & Salmon (1995) and is designed to provide a direct bound on the errorbudget associated with the monopole approximation.There are several flavors of SNSPH, and each integrates the hydrodynamics equations in aslightly different manner. For this investigation we have chosen the version that integrates thespecific internal energy in order to facilitate a direct determination of temperatures from the SPHsimulations. For our SNSPH calculations, we use the simple viscosity prescription utilizing bothbulk ( α ) and von Neumann-Richtmyer ( β ) viscosity coefficients (Benz 1990; Monaghan 1992, 2005;Fryer et al. 2006). We do not evolve these constants, and most of our calculations used standardvalues: α = 0 . β = 2 .
0. However, we did run a test calculation using α = 0 .
125 and β = 0 . Coyote located at Los Alamos NationalLaboratory.
Our analysis of each binary simulation has been carried out primarily using two diagnostictechniques: (1) Throughout each simulation, a variety of global binary parameters were evaluatedby performing integrals over the computational volume as detailed in Appendix A, and then curveswere plotted to show the variation with time of these parameters. (2) For each simulation, a moviewas generated to show the time-variation of the matter distribution. We determined the degree towhich a simulation performed with the SPH code resembled the simulation performed with the gridcode by comparing the behavior of various parameter curves and/or by comparing correspondingframes from movies generated with the two separate codes.Most of the figures shown in subsequent sections of this paper display results from either oneor both of these diagnostics. For example, the curves drawn in Figure 3a show the time-dependentbehavior of six binary system parameters derived from three separate Q1.3P simulations — orbitalseparation, a , binary mass ratio, q , log of the donor’s mass-transfer rate, | ˙ M D | , orbital angularmomentum, J orb , and spin angular momentum of the accretor, J A , and donor, J D . The colorcontour plots shown in Figure 3b are movie frames drawn from two separate points in time duringtwo of these Q1.3P simulations. Since the two numerical codes use different frames of reference,an inertial reference frame in the case of the SPH simulations and a frame of reference rotating atthe initial frequency of the binary in the grid simulations, there will be a trivial phase angle offsetbetween the two codes. For the purpose of aiding comparisons between the two codes, the gridcode data are rotated to match the phase of the SPH simulation for the color contour plots of theprojected mass density. Throughout the paper, such figures are used to illustrate key findings fromour comparison simulations. We also note that, in general, the frames taken from our grid codesimulations show an outline of the outer edge of the cylindrical grid at low density levels. This is anumerical artifact resulting from imperfect outflow boundary conditions applied at the outer edgeof the computational domain.Initially we agreed that the group that developed the SPH code and the group that developedthe grid code would perform these diagnostic analyses independently. Eventually we realized thatsome differences that appeared when comparing diagnostics from different simulations arose from 12 –subtle differences that appeared in the volume integrals that were being used by the two differentgroups to define various global binary parameters. There are, for example, different ways to definethe orbital angular momentum or the spins of the two stars when some of the matter resides in anenvelope that is shared by the two stars. We ultimately decided, for diagnostic purposes only, tomap results from the SPH simulation onto a 3D grid and to use the diagnostic software developedfor the grid code to also analyze results from each SPH simulation. The diagnostic figures presentedthroughout the remaining sections of this paper were derived in this manner. Appendix A explainshow each volume integral was defined in order to evaluate the various global binary parameters.The initial data for the Q0.7 binary were developed independently by the two groups resultingin the models Q0.7a and Q0.7b described in Table 1. We made use of the freedom of the chosenpolytropic equation of state to re-scale the data from the grid evolutions to match the total massand initial separation of the SPH model Q0.7b for the purposes of analyzing the simulation results.Both thee line plots and images of the projected mass density from the grid code have been rescaledaccordingly.
3. Evolutions
As is cataloged in Table 3, we are reporting on results from 23 separate binary mass-transfersimulations. Each of the simulations was followed through at least eight, usually more than adozen, and up to fifty orbits with each individual orbit requiring tens of thousands of Courant-limited integration time steps (for the grid code simulations). Hence, our critical examination ofthe behavior of these systems and, most importantly, the comparison of behavior between differentsimulations, has required the analysis of a tremendous amount of simulation data. In what follows,our formal discussion of the results is broken down into a series of subsections that, generallyspeaking, march through the models in the order in which they are listed in Table 2.Before diving into the details, the reader may find it advantageous to skim through the figuresand/or some of the accompanying movies (21 in total, as labeled in column 5 of Table 3) in orderto sample the variety of physical phenomena that arise in these types of systems. Examples are: • Binary merger – see video01 or video02, and Figure 3b; • Tidal disruption of the donor – see video04, video16, or video21, and Figure 4b or Figure 14; • The hot, tenuous common envelope that forms in “ideal gas,” as opposed to “polytropic,”simulations – compare video07 to video17, and see Figure 9; • Nonaxisymmetric modes having four-sided (box-like) and three-sided (triangular) shapes ex-cited in an accretion disk – see video07 or video12, and Figure 8; • Binary detachment following a long phase of relatively stable mass transfer – see video07,video09, or video11. 13 –
Here we discuss and compare results from simulations that were carried out using a polytropicequation of state. The relevant simulations are labeled as Q1.3P, Q0.7P, Q0.5P, and Q0.4P alongwith a two-character model identification suffix (for example, simulation Q1.3P G
1) in the first 13rows of Table 2. In each case the structure of the binary system at time t = 0 was prescribedby one of the five, initially uniformly rotating, n = polytropic binary models identified in Table1. At the start of a simulation, the donor star was driven into contact with its Roche lobe byartificially draining angular momentum from the system at a rate of one percent per orbit for the“driving duration” specified in column 5 of Table 2. Thereafter the system was evolved in a fullyself-consistent fashion without further external influences. As is documented in column 4 of Table2, some SPH simulations were carried out using 100 thousand (100 k ) particles and others used 1million (1 M ) particles; most simulations performed with the grid code used a uniform cylindricalgrid containing approximately four million grid cells but we include one grid simulation of theQ0.4P system at a higher resolution of approximately 47 million (47 M ) cells.To aid in our analysis, we have established the time of an event common to each mergersimulation, the merger time, t merge , as a means of synchronizing these simulations. While thechoice of the exact merger time is in some sense arbitrary, we have found empirically that the firstdata set where the merged object has only one density concentration could be identified relativelyeasily by eye. We found this particular choice of merger time could be identified in a repeatablemanner for both grid and SPH data sets from the animations corresponding to each simulation.Merger times are recorded in the third column of Table 3 for those simulations that merge. Witha common event identified in all merger simulations, we can measure time relative to t merge toeliminate trivial phase differences and variations in instability growth rates in different simulationsto compare the highly nonlinear state of the evolving fluid in a clear way. From t merge we haveestablished a zero point time relative to a particular grid-based simulation as,( t zpt ) ID = ( t merge ) G − ( t merge ) ID . (6)Time sequences of data from a particular simulation can be compared then by plotting the behaviorof various system parameters using the shifted time, t shift ≡ t + t zpt , (7)to label the temporal axis. For reference, column 4 of Table 3 contains t zpt for each mergersimulation. Q1.3P
Model Simulations
We carried out three separate adiabatic simulations of the binary polytrope that had an initialmass ratio q = 1 . . P . As is indicated by the note referenced in the last column of Table2, some results from the grid-code simulation Q1.3P G1 already have been presented as simulationQ1.3-D in § G
1) or the SPH code (models S S q crit . As Figure 3a illustrates, throughout each model simulation, thebinary separation ( a ) and system mass-ratio ( q ) steadily decreased while the donor mass-transferrate ( | ˙ M D | ) steadily increased. As is shown in the accompanying movies – video01 from simulationQ1.3P G
1, video02 from simulation Q1.3P S
1, and video03 from simulation Q1.3P S G G t merge = 12 . P . The lower right-hand image inFigure 3b is the frame picked from the movie of SPH-code simulation Q1.3P S t merge = 7 . P for simulation Q1.3P S
1. In a similar manner we have deduced that t merge = 8 . P for the lower resolution SPH-code simulation Q1.3P S G S G S
1. These two images alsodisplay a high degree of quantitative similarity. They illustrate what a careful viewing of video01,video02, and video03 reveals: the two different numerical schemes produce very similar dynamicalresults, especially as the simulations each evolve through the time of merger. Note that the SPHsimulations are conducted in an inertial frame whereas the grid-code simulations are computed in 15 –a frame initially corotating with the binary. Consequently, in all the figures showing comparisonsof equatorial density distributions obtained with both codes, we have rotated the density contoursobtained with the grid code to match the orientation of those obtained with the SPH code.We cannot ignore the fact that the precise time of merger was different in the three simula-tions. The decidedly different merger times seem to suggest that different code algorithms or modelresolutions produce different results. But a more careful examination reveals that if an appropri-ate zero-point time-shift, t zpt , is added to the recorded evolutionary time of each simulation, thesimulations show a high degree of agreement. For this particular set of models, the merger timeshave been used to define t zpt as shown in column 4 of Table 3. In Figure 3a the time-evolutionarybehavior of various system parameters has been displayed using the shifted time that synchronizesmerger between the simulations. Viewed from this perspective, we conclude again that the threeseparate Q1.3P simulations show a remarkable level of agreement.Consider first the middle panel on the right-hand-side of Figure 3a, which displays the time-dependent behavior of the system’s orbital angular momentum. During the first two orbits, J orb steadily decreases at a rate of 1% per orbit due to the imposed driving, then it remains nearlyconstant until the merger occurs. Over the last ∼ J orb . After the zero-point shift in time has been takeninto account, all three curves that describe the time-evolutionary behavior of the binary separation(see the top panel on the left-hand side of Figure 3a) also lie on top of one another over the last ∼ P .All three simulations show that, once mass-transfer commences, the rate increases by morethan two orders of magnitude in less than 10 orbits; see the middle panel on the left-hand sideof Figure 3a. It is gratifying to see that, after t zpt has been taken into account, all three | ˙ M D | diagnostic curves converge immediately prior to merger. However, at earlier times the slope of themass-transfer curve produced by the grid-code simulation does not match the slope of the curvesproduced by the SPH simulations. At the end of the initial phase of driving, the value of | ˙ M D | in the grid-code simulation is also somewhat lower than the value of | ˙ M D | in the SPH simulation.This is the first indication that, at relatively low mass-transfer rates, the simulations conductedwith our two different numerical techniques generate results that are quantifiably different fromone another. Figure 3a also shows that, while the blue and red curves illustrating the behavior of q ( t ) and J A ( t ) overlap at late times, they do not lie precisely on top of each other at earlier times.For both of these binary parameters, this behavior probably reflects the discrepancy already notedin mass-transfer rates at early times.The second indication that simulations conducted with the two different numerical codes gen-erate results that are quantifiably different from one another comes from the bottom panel on the 16 –right-hand side of Figure 3a, which illustrates the behavior of the spin of the donor star. In the gridcode, J D slowly but steadily decreases up until the moment of merger, whereas in the SPH code (atboth resolutions) J D initially increases, then decreases for a while, before increasing dramaticallyat the moment of merger. Overall, the Q1.3P simulation does not appear to be seriously impactedby this discrepancy in the behavior of J D ( t ). This is likely because the donor stores less than 10%of the system’s total angular momentum and, because q >
1, the system is doomed to merge,independent of the precise treatment of spin-orbit coupling.
Q0.7P
Model Simulations
As is recorded in Table 2, we carried out three separate polytropic simulations of binary systemshaving an initial mass ratio of q = 0 .
7. Simulations conducted with the grid code (models G G
2) were started from the initial model identified in Table 1 as Q0.7a while the simulationconducted with the SPH code (model S
1) was started from the initial model identified as Q0.7b.These initial models were very similar, but not identical to one another; for example, the two starsin model Q0.7b were constructed as n = 3 / n = of q crit = . Simulations G G G . P while G . P . Model S . P and was conducted with100 k SPH particles.As Figure 4a illustrates, the time-variation of key binary parameters exhibited by the variousQ0.7P simulations was very similar and resembled the time-variations seen in simulation Q1.3P.For example, throughout each simulation a and q steadily decreased while | ˙ M D | steadily increased.We conclude that, as with the Q1.3P simulation, the Q0.7P simulation was dynamically unstableto mass transfer. It is not entirely appropriate to describe the Q0.7P simulation as a merger,however. As is shown in the accompanying movies – video04 from simulation Q0.7P G
1, video05from simulation Q0.7P G
2, and video06 from simulation Q0.7P S G t merge = 9 . P for simulation Q0.7P G
1. Thebottom right-hand image in Figure 4b is the frame picked from the movie of SPH-code simulationQ0.7P S t merge = 11 . P for simulation Q0.7P S
1. In a similarmanner we have deduced that t merge = 21 . P for grid-code simulation Q0.7P G
2. These merger 17 –times are recorded in the third column of Table 3. The “instant of merger” frames selected fromsimulations Q0.7P G
1, Q0.7P G
2, and Q0.7P S k particles.The left-hand images in Figure 4b present frames from the movies of simulations Q0.7P G G
2, and Q0.7P S t = − . P ) prior to t merge . These images also displaya high degree of quantitative similarity. The pairs of images in Figure 4b illustrate what a carefulviewing of the relevant movies reveals: the two different numerical schemes produce very similardynamical results, especially as the simulations each evolve through the phase of tidal disruptionof the donor. Furthermore, a comparison of the top pair of images with the middle pair of imagesshows that the outcome of this particular simulation is relatively insensitive to the precise amountof driving that is applied in order to initially bring the surface of the donor into contact with itsRoche lobe.As was seen in the context of the Q1.3P simulations, the time of merger is different in thethree Q0.7P simulations. It is not surprising that tidal disruption of the donor occurs earlier ingrid-code model G G G
1, the “driving duration” was longerand, hence, the donor was initially driven into deeper contact with its Roche lobe. As Figure4a illustrates, all three Q0.7P simulations show very good agreement when the time-evolutionarybehavior of various system parameters is displayed using the shifted time defined by equations (6)and (7).Again, it is worth pointing out ways in which results from the grid-code simulations differ fromresults obtained from the SPH-code simulation. As Figure 4a shows, throughout the simulation themass-transfer rate increases somewhat more rapidly in simulation Q0.7P S G G q ( t ) curves do not lieprecisely on top of one another at early times. The curves that trace the time-dependent behaviorof J orb and J A in the SPH-code simulation (blue curves) lie slightly below the (red) curves drawnfrom the grid-code simulations. This is almost certainly because, in this case, the initial modelswere not identical. While the blue curve describing the time-evolutionary behavior of J D does notlie precisely on top of the red curves, the general shape of the curves derived from the two differentnumerical techniques appears to be more similar than was seen in Figure 3a for simulation Q1.3P.We note that the short period ( ∼ P ) oscillations that are visible in all a ( t ) curves — seethe upper left-hand panels of Figures 3a and 4a — arise because, in practice, none of the binaryorbits is perfectly circular. The amplitude of the resulting “epicyclic oscillations” is determined, inpart, by the duration of driving that was initially imposed on each model simulation. The absolutephases of the epicyclic oscillations generally do not — and should not be expected to — matchwhen the a ( t ) curves are plotted as a function of the shifted time. 18 – Q0.5P and
Q0.4PWe have simulated polytropic mass-transfer in three binaries with initial mass ratios q < q crit :Simulations Q0.5P with q = 0 .
500 and polytropic simulations Q0.4a with q = 0 . q = 0 . G S . P . The evolutionary behavior of model Q0.5P G S G § G
1, Q0.4P G
2, and Q0.4P S
1, which were initiallydriven into contact at a rate of 1% per orbit for 1 . P . Simulations Q0.4P G G G M grid cells while simulationQ0.4P G M grid cells, partly in an effort to establish convergence of resultswith the grid code and to measure the lowest initial mass transfer rate that can be establishedin a numerical simulation of a semi-detached binary. Movies showing the evolutionary behaviorof these three models are labeled video09, video10, and video12, respectively. The left panel ofFigure 5 displays the time-evolutionary behavior of a , q , | ˙ M D | , J orb , J A , and J D from simulationsQ0.5P G S
1. Analogous information is displayed in the right panel of Figure 5 fromsimulations Q0.4P G G S
1. In all of the Figure5 plots, red curves present results from grid-code simulations and blue curves present results fromthe 1 M -particle SPH-code simulations.The Q0.5P and Q0.4P model simulations differ from the simulations discussed earlier (Q1.3Pand Q0.7P) in one especially significant way: The binary separation, a , does not decrease mono-tonically with time. Instead, a decreases for only a brief period of time associated with the initiallyimposed episode of driving, then the system reaches some minimum separation, a min , and a increasesmonotonically thereafter. After a begins to increase, the qualitative behavior of each system’s sub-sequent evolution appears to depend on the depth of contact that was established between thedonor and its surrounding Roche lobe during the phase of driving, as well as the rate of masstransfer, | ˙ M D | , exhibited at the end of the phase of driving, which correlates with the depth ofcontact. Simulations in which sufficiently deep contact is established display a mass-transfer ratethat continues to increase with time even as the system separates; hereafter, these will be referredto as “Case A” evolutions. If, on the other hand, the donor established mass-transfer rate growsfor only a short period of time, reaches a maximum rate | ˙ M max | , then levels off or declines as thesystem separates; these will be referred to as “Case B” evolutions. 19 –The dashed curves sketched in Figure 6 illustrate in a qualitative fashion how other key param-eters of the Q0.5P and Q0.4P binary systems respond to a steadily increasing mass-transfer rate in“Case A” evolutions; the solid curves illustrate how these key parameters respond to the “Case B”behavior. In Case A evolutions ( e.g. , simulation Q0.4P S e.g. , simulation Q0.4P G q ( t ), J orb ( t ), and J D ( t ) each developan inflection point, indicating a slowing of the rate at which q , J orb , and J D are declining. In someCase B evolutions, angular momentum is eventually returned to the donor.Comparing the two grid-based evolutions Q0 . G . G
2, we find that the higher-resolution simulation Q0 . G . G . G q model simulations, we again used both numerical codes to conduct polytropic simulationsstarting from initial model Q0.4 (see Table 1) but subjected the system to 1% per orbit drivingfor only 1 . P . These are identified as simulations Q0.4P G S G § . P . A movie showing the first 1 . P of this “baseline simulation” was published by MTFD07;here video11 carries the same simulation through 47 P ). Key results from these two simulationsare displayed as dashed curves in Figure 7 — the red and blue curves are, respectively, from thegrid-code and SPH-code simulations. For comparison purposes, the solid red and blue curvesdisplayed in Figure 5b have been redrawn in Figure 7. Both of the SPH-code simulations displaycharacteristics of a “Case A” evolution whereas both of the grid-code simulations (the dashed andsolid red curves) display characteristics of a “Case B” evolution.Based on the trends seen in Figure 7, it seems likely that, by subjecting the initial model to aphase of driving that is even shorter than 1 . P , the SPH code would produce an evolution with“Case B” characteristics. Alternatively, it seems likely that, by subjecting the initial model to aphase of driving that is even longer than 1 . P , the grid code would produce an evolution with“Case A” characteristics. Via a set of grid-code simulations, MTFD07 have previously provided 20 –evidence (see especially their Figure 4) that, as the donor star is initially driven into deeper anddeeper contact with its Roche lobe, a transition can be made from simulations that exhibit “CaseB” characteristics to simulations that display “Case A” characteristics.Our analysis of Q0.4P and Q0.5P model simulations also suggests the following. If a simulationof binary mass transfer traverses a “Case B” evolutionary trajectory, the donor may detach from itsRoche lobe and, in so doing, avoid merger or tidal disruption. Three of the grid-code simulationsactually show the mass-transfer stream detaching from the accretor late in the simulation: Video07shows the binary detaching in simulation Q0.5P G t ≈ P ; as video11 shows, starting fromthe same initial state but with a shorter episode of driving, the stream detaches later, at 46 . P ; andvideo09 shows the binary detaching in simulation Q0.4P G t ≈ P . The corresponding SPHsimulations (Q0.5P S S As has already been mentioned, MFTD07 used the same grid code as we have used here toanalyze the behavior of Q0.4P simulations that were subjected to a variety of driving rates anddriving durations. In § ≥ m ≥ t = 15 . P in MFTD07’s Q0 . D simulation, the disk surroundingthe accretor has a triangular shape. As the image in the upper-left-hand corner of our Figure 8shows, at time t = 13 . P in the same MFTD07 simulation, the disk surrounding the accretor hasa box shape; the time-evolutionary transition between the box and triangular shapes is shown inthe mpeg animation of simulation Q0.4D that accompanies the MFTD07 publication.The bottom two images in Figure 8 have been drawn from our video12. They illustrate thatanalogous box- and triangular-shaped distortions developed around the accretor in SPH simulationQ0.4P S
1. The box-shaped distortion appears at time t = 16 . P and the triangular-shapeddistortion appears at time t = 17 . P . We were gratified to see these azimuthal distortions developspontaneously in simulation Q0.4P S
1. It demonstrates that our two entirely different numericaltechniques for simulating binary mass-transfer can encounter and follow with comparable fidelitythe development of unexpected nonlinear-amplitude structures in the midst of a complex binaryinteraction.Box- and triangular-shaped distortions are not apparent in our grid-code simulation (Q0.4P G S
1. Conversely, a careful viewing of video07 shows development of box- and triangular- 21 –shaped distortions in the grid-code simulation of model Q0.5P, whereas these distortions are notapparent in the corresponding SPH simulation (model Q0.5P S
1; see video08). We suspect – asdid MFTD07 – that “standing wave” azimuthal distortions of this type are routinely excited inpolytropic simulations as a result of dynamical interactions between the mass-transfer stream andthe accretion disk. But the distortions do not grow to nonlinear amplitude, and therefore are notvisible to the eye, unless | ˙ M D | is sufficiently large. In this light, it is not surprising that our eyesare unable to detect box- or triangular-shaped distortions in video09 even though they are presentin video12 because, for model Q0.4P, the mass-transfer rate is always higher in the “Case A” SPHsimulation than in the “Case B” grid-code simulation (see the Figure 5b plot of the time-dependentbehavior of | ˙ M D | ). And it is not surprising that our eyes are unable to detect nonlinear amplitude,standing-wave distortions in video08 even though they are present in video07 because, for modelQ0.5P, the mass-transfer rate remains relatively low in the SPH simulation whereas it grows to mod-erately high levels in the long, grid-code simulation (see the Figure 5 plot of the time-dependentbehavior of | ˙ M D | ). Here we discuss and compare results from simulations that were carried out using an ideal gas(“I”) equation of state. The relevant simulations are labeled as Q1.3I, Q0.7I, Q0.5I, and Q0.4I alongwith a two-digit model identification suffix in the last 10 rows of Table 2. As with the polytropicsimulations, in each case the structure of the binary system at time t = 0 was prescribed by one ofthe initially uniformly rotating, n = 3 / M ) particles; one model (Q0.7I S
1) used 100 thousand (100 k ) particles. Every simulationperformed with the grid code used a uniform cylindrical grid containing approximately 4 milliongrid cells.In all of the ideal gas simulations, one structural feature develops that distinguishes eachsimulation from its polytropic counterpart: During the phase of mass transfer, a relatively hot,tenuous “envelope” develops around the accretor. This happens because the entropy of the materialincreases as the donor material passes through the shock front that is established where the accretionstream impacts the surface of the accretor. (The shock front also develops in all of the polytropicsimulations but, by design, the entropy of the gas was not allowed to change.) This low-densityenvelope is easy to spot in the movie that accompanies each ideal gas simulation, especially if it isviewed alongside the movie from its polytropic counterpart; the lowest density, blue contour regionis significantly extended in each ideal gas simulation. The accretion stream is also noticeably thicker 22 –(hotter) in each ideal gas simulation, relative to its polytropic counterpart. And in some cases thehot tenuous envelope migrates over to the donor, establishing a classic “common envelope” in whichthe two stars are embedded. The image displayed in the top, right-hand corner of Figure 9 hasbeen extracted at time t = 17 . P from video17, which shows the ideal-gas simulation of modelQ0.5 conducted with the grid code – specifically, simulation Q0.5I G
1. The common envelopestructure is very apparent when compared with a frame (top, left corner of Figure 9) extractedfrom the middle ( t = 14 . P ) of video07, which shows the polytropic simulation conducted withthe grid code that started from the same initial model (Q0.5) – simulation Q0.5P G
1. Similarly,the bottom two panels of Figure 9 display, for comparison, images from the polytropic simulation(video08; model Q0.5P S
1) and the ideal gas simulation (video18; model Q0.5I S
1) that startedfrom initial model Q0.5 and were conducted using the SPH code.In ideal gas simulations performed with the grid code, there is a tendency for a moderatelyextended “blue” atmosphere to develop around the binary components within the first orbit. Thishappens because, in the grid code, the region surrounding the initial binary system is not actuallya vacuum but, rather, contains a very low density gas at a defined “floor” level. Once a simulationbegins, this low-density gas free-falls down onto the donor, crashes into the surface supersonically,and heats up as a result. (The same dynamics occur in each corresponding polytropic simulationbut a hot donor atmosphere does not result because, again by design, the entropy of the gas wasnot allowed to increase as it crashed into the stellar surface.) By contrast, an extended atmospheredoes not typically develop around the donor during the first orbit of most SPH simulations becausethe region surrounding the binary system initially contains no material.
Q1.3I
Model Simulations
We carried out two separate ideal-gas simulations of binaries that had an initial mass ratio q = 1 . . P in order to initiatemass transfer. Whether the simulation was conducted using the grid code (model G
1) or the SPHcode (model S a ) and system mass-ratio( q ) steadily decreased while the donor mass-transfer rate ( | ˙ M D | ) steadily increased. As is shownin the accompanying movies – video13 from simulation Q1.3I G
1, and video14 from simulationQ1.3I S G S S S G G P . This is reflected in, for example, a longer period of time for the masstransfer rate to grow in the grid evolution shown in Figure 10a compared to the polytropic evolutionin Figure 3a. This delay is likely related to the moderately extended, tenuous “blue” atmospherethat developed in the star’s Roche lobes in the grid-code, but not the SPH, simulation. A tenuousatmosphere appears early in video13 but is absent in the polytropic simulation shown in video01 andSPH clearly does not suffer from this artifact in either polytropic (video02) or ideal gas (video14)evolutions. With the relatively long driving time required to bring the donor in to contact with itsRoche lobe, the simulation Q1 . I G . P G Q0.7I
Model Simulations
As part of our comparison, we conducted two simulations of the binary with an initial massratio of q = 0 . k particles and was driven intocontact for only 1 . P while the grid simulation was conducted with approximately 4 million gridcells and was driven into a deep contact with a driving period of 2 . P .The effects of a lower resolution in the SPH simulation Q0.7I S S G S t zpt = 2 . P both numerical techniques yield a very similarsolution for the tidal instability that destroys the donor star. Images drawn from these simulationsduring tidal instability are shown in Figure 11b.As in the Q1 . I discussion, the SPH code gives nearly identical results regardless of the assumedequation of state. For the grid code, the results are more similar to the polytropic evolutions andnoticeably, there is no delay in the time to merger and the mass transfer rate does not stall in theideal gas simulation. One difference that can be discerned in the grid-based evolution Q0 . I G . P G . P G
2, the donor’s spin angular momentumdeclines gently before rising steeply as the donor star is disrupted.
Q0.5I and
Q0.4ITo complete our testbed of simulations, we have repeated simulations of the two low massratio models with an ideal gas equation of state yielding the six simulations Q0.5I G
1, Q0.5I S G
1, Q0.4I G
2, Q0.4I S S
2. The simulations of the Q0.5 model offer a goodrepresentation of a DD system as both components in the initial data share the same value of thepolytropic constant (see Table 1). As the donor star was relatively far from contact initially, thesystem was driven for 2 . P in both the SPH and the grid-code simulations.The SCF model Q0.4a is somewhat more complicated as the polytropic constants for the donorand accretor are not equal and the donor star is too large (compared to what would be the casewith a zero-temperature white dwarf equation of state). The relatively large volume of the donor’sRoche lobe that is occupied by the star has important consequences for the outcome of the idealgas simulations of this model. We have run this model four times, twice with each code where wehave varied the initial depth of contact and hence the initial mass transfer rate. The simulationsQ0.4I G S . P to establish contact while the simulations Q0.4I G S . P to establish a more shallow initial contact.Including the ideal gas equation of state introduces a new phenomenon into simulations ofbinaries with initial mass ratios below the cutoff for dynamically unstable mass transfer, q crit .Namely, the accreting star will expand to the point that it also fills its Roche lobe forming acommon envelope about both stars with the possibility that mass is lost from the binary throughthe outer Lagrange points (see the movies video17 and video18 as well as the still images in thebottom row of Figure 14). The incorporation of shock heating in the ideal gas equation of stateintroduces an additional, spurious effect in grid-based simulations. In an Eulerian treatment, nogrid cell is allowed to be completely empty and a floor level of the mass density is maintainedartificially even in regions that should clearly be vacant. At the start of a simulation, this materialwill begin to free-fall onto the binary system and shock heat as it crashes into the stars, heatingtheir outer layer. As can be seen in the plot of mass transfer rate for the Q0.5I simulations in theleft hand panel of Figure 12, the grid simulation Q0.5I G S
1. The hotter envelope for the donor star in the grid-based simulation means that thecommon envelope phase where the separation varies about a constant mean value is longer than inthe SPH simulation Q0.5I S G
1. In the grid-based simulation a common envelopephase is established where the mass transfer rate grows with relatively little change in the orbitalseparation or the distribution of angular momentum between the orbital motion and the spin of eachcomponent. This phase is followed by a rapid tidal instability where the orbital angular momentumplummets while the spin angular momentum of each star grows rapidly. In the SPH evolution ofthis model there is still a common envelope phase but the evolution seems to be approaching thetidal disruption of the donor star, not a tidal instability.A similar pattern can be seen in the ideal gas simulations of the Q0.4a initial data. The righthand panel of Figure 12 shows data for the two simulations with a greater initial depth of contact(corresponding to a longer driving period of 1 . P ). The grid simulation Q0.4I G . P in simulation Q0.4I S G S G
1. While the initial mass transfer rates and subsequent time scales aredifferent between the two runs, there is detailed agreement between the structure of the commonenvelope, the tidal disruption of the donor star and the resulting disk of debris wrapped aroundthe accretor.
4. Summary and Conclusions
For reference, we have summarized our conclusion for the fate of each simulation (when itcould be determined given existing data) in the seventh column of Table 3. Possible outcomesrange from the donor star detaching from its Roche lobe, the destruction of the binary throughtidal instability (the spin angular momentum of each component grow rapidly while the orbital 26 –angular momentum crashes), and the tidal disruption of the donor star (where the spin angularmomentum of the donor remains relatively constant while the donor is tidally shredded). Theseoutcomes are in some instances modified by the formation of a common envelope phase which canbe identified as relatively steady orbital separation while the mass transfer rate grows.The standard linear analysis of the dynamical stability of mass transfer in binaries with spher-ical components, in the Roche approximation, and driven by angular momentum losses (AML)(Rappaport et al. 1982; Soberman et al. 1997; Marsh et al. 2004; Gokhale et al. 2007), reveals twocritical values of the mass ratio q = M D /M A : 1) if q > q a , the effect of the mass transfer by itselfwill tend to reduce the binary separation, whereas it will tend to expand the binary otherwise; 2) if q < q crit , the mass transfer is stable and will tend toward a secular value proportional to the AMLrate, whereas the mass transfer will be unstable otherwise and grow dynamically for q > q crit .In the simplest case of conservative mass transfer (no loss of mass from the system) andneglecting the effects of tidal distortions and spin-orbit coupling, q a = 1, and as already statedearlier, for donors obeying a polytropic EOS with n = 3 / q crit = 2 / < q a . Therefore, basedon these critical mass ratios, the expected evolutionary behavior of interacting, polytropic binarysystems can be divided into three broad categories, depending principally upon whether the system’sinitial mass ratio q falls above q a , in between q crit and q a , or below q crit . If q > q a , both AMLand mass transfer contribute to shrinking the binary separation and driving the donor into deepercontact. Thus we expect the mass-transfer rate to increase monotonically while the separationdecreases monotonically, probably leading to a merger. If q crit < q < q a , AML will shrink thebinary initially, but once contact is established, the dynamically unstable mass transfer is expectedto grow, halt the contraction and then rapidly increase the separation. This behavior resemblesclosely the “Case A” evolutions of § q < q crit , afterAML shrinks the binary into contact, the mass transfer will settle toward the expected secularequilibrium value. If, as in our case, q crit < q a , the expectation is that the orbit will settle into asecular expansion, with the donor, and thus the binary, surviving. This behavior resembles in partthe “Case B” evolutions of § q > q a ), following the onset of mass transfer, are dynamically unstable from the 27 –outset. As the binary separation decreases, the mass-transfer rate monotonically increases and,in a finite time, the stellar components either merge or the donor tidally disrupts. On the otherhand, once mass-transfer has been initiated in systems that have a relatively small q < q a , thebinary separation will generally increase with time. Our two distinctly different CFD algorithmsagree on the evolutionary outcome of these dynamically unstable systems. For example, in everyone of our Q1.3 simulations, the cores of the two stars merge; whereas, in every Q0.7 simulation,the donor tidally disrupts and material from the donor disperses to form a rotationally flatteneddisk around the accretor. As is illustrated by Figures 3b and 4b (and their associated animationsequences), this agreement is more than qualitative: detailed structures that arise during mergeror tidal disruption events are present in simulations that are conducted with both codes.Highly non-linear, transient structures such as the box- and triangle-shaped resonances in theQ0 . a polytropic evolutions have emerged in simulations with our two disparate codes as highlightedin Figure 8. These transitory equatorial distortions indicate that both techniques are yieldingsolutions remarkably close to one another in the parameter space of possible numerical solutionsdespite their different error terms, numerical artifacts, and fundamental representations of a fluidstate. These resonance features represent a non-trivial verification of our codes.For a given model, our two different CFD codes generally do not agree on the precise time ofmerger (or of tidal disruption), t merge . But as is summarized in Table 3, even the same code willproduce different values of t merge when a simulation is repeated at a different code resolution or witha different driving duration. These apparent discrepancies arise because the precise value of t merge is sensitive to the depth of contact that the Roche lobe makes with the donor star following theinitial phase of driving. All other things being equal, a relatively brief phase of driving will result inrelatively shallow contact, in which case the binary simulation will begin with a low mass-transferrate and it will take the system a relatively long time to merge. If the phase of driving lasts longer,the simulation will begin with a higher mass-transfer rate and, as a result, the time to merger willbe shorter. For a given driving duration, we have found that the time to merger is shorter in anSPH-code simulation than in a grid-code simulation. In effect, this implies that for a specifieddriving duration the Roche lobe establishes deeper contact with the donor star in the SPH codethan in the grid code. We must emphasize, however, that the ultimate outcome of these simulationsis not sensitive to code resolution or to driving duration. Also, as stated earlier, the outcome doesnot depend on the selected CFD algorithm: If the chosen binary system is dynamically unstabletoward merger (or to tidal disruption of the donor), all simulations ultimately converge to the samesolution.Our collection of simulations demonstrates that the evolutionary outcome of systems in thethird category (relatively small q < q crit ) is sensitive to the depth of contact that the Rochelobe makes with the donor star during an initial phase of driving. (See especially the discussionassociated with Figure 6.) If, for a given q , the depth of contact is sufficiently deep, then themass-transfer rate will monotonically increase even as the binary separation is increasing and,ultimately, the donor star will tidally disrupt. This is what has been referred to in § q and driving duration —at least within the range of parameters we have explored in this low q category, using a polytropicequation of state — simulations conducted with the SPH code tend to result in “Case A” simulationswhereas simulations conducted with the grid code tend to result in “Case B” simulations. As aconsequence, we find that the donor star tends to tidally disrupt in simulations conducted withthe SPH code whereas the donor tends to survive in simulations conducted with the grid code, atleast for the length of time we’ve been able to follow the simulations. A given system can thereforeappear to be either unstable or stable, depending on the CFD algorithm that has been selected toperform the simulation. This is not a particularly desirable result. But it is understandable and itis consistent with the behavior deduced from our analysis of relatively large q simulations. For aspecified driving duration the Roche lobe typically establishes deeper contact with the donor starin the SPH code than in the grid code. In turn, deeper contact will tend to result in a “Case A”rather than a “Case B” evolution.When we compared polytropic and ideal gas simulations of the same initial model we noted theformation of a low-density hot envelope and the greater thickness of the stream both as a result ofthe higher temperature of the gas in the ideal gas case. However, while the higher temperature ofthe envelope resulting from the accretion shock appears both in the SPH and grid simulations andis likely a true physical effect, the dissipation of energy due to “vacuum” material settling on thebinary during the first orbit is only present in grid simulations and should be considered a spuriouseffect. The spurious heating generates early on a hotter, thicker atmosphere, which allows a widerstream and enhances the mass transfer. This tends to speed up the simulation, sets it on a moreunstable pathway and makes it difficult for the system to return to stability when the donor is nearcontact initially. This has little or no effect on the ultimate fate of the Q1.3I and Q0.7I models,but it does affect the simulation of the models Q0.5I and Q0.4I. In particular, the higher initialmass transfer rate in the grid simulations for the Q0.4I model results in a faster tidal disruptionfor the grid simulations than in the SPH simulations. This is in contrast to most of the other casespresented here, where the mass transfer rate in SPH simulations is higher yielding also a faster finaloutcome.An obvious difficulty in the interpretation of our low q simulations arises in connection with theuncertainty in the effective value of q crit because of the noted differences between the assumptionsmade in analytic treatments and the “reality” of the simulations. Since direct impact occurs in oursimulations, but also disk-like flows arise around the accretor, we may guess that the effective valueis somewhere in the range 0 . < q crit < /
3. However, tidal effects are obvious in the images andmovies presented here, and our simulations suggest that the ultimate fate of the low q systems maywell depend on the depth of contact achieved in the initial driving into contact. The levels of masstransfer in the simulations exceed the likely levels in real systems by many orders of magnitude.This may not be important for systems that are expected to be dynamically unstable, but caution 29 –should be exercised when interpreting results for marginally unstable systems.In summary, we conclude that aside from understandable shifts in the times taken by the samemodels evolved with SPH and grid hydrodynamics to reach the same evolutionary stages, and thesplit into Case A/B evolutions for low q systems, the agreement between the results produced byour codes is remarkably close in the final outcomes and during intermediate evolutionary stages.We acknowledge valuable interactions that we have had with Dominic Marcello, KundanKadam, and Zachary Byerly. In addition, we would like to thank the anonymous referee fortheir very careful and thoughtful review of our manuscript. This work has been supported, inpart, by grants AST-0708551, NSF CREATIV grant AST-1240655, and DGE-0504507 from theU.S. National Science Foundation and, in part, by grant NNX07AG84G from NASA’s ATP pro-gram. This research also has been made possible by grants of high-performance computing timeon the TeraGrid (MCA98N043), at LSU, and across LONI (Louisiana Optical Network Initiative).This material is based upon work supported by the National Science Foundation under Grant No.ACI-0338618l, OCI-0451237, OCI-0535258, and OCI-0504075. This research was supported in partby the Indiana METACyt Initiative. The Indiana METACyt Initiative of Indiana University issupported in part by Lilly Endowment, Inc. This work was supported in part by Shared UniversityResearch grants from IBM, Inc., to Indiana University. This material is based upon work supportedby the National Science Foundation under Grant No. CNS-0521433. 30 – A. Diagnostics
In an effort to quantitatively compare results from our various simulations of binary mass-transfer events, we have tracked the time-evolutionary behavior of a set of six key binary systemparameters: Orbital separation, a , binary mass ratio, q , mass-transfer rate of the donor, ˙ M D , orbitalangular momentum, J orb , and spin angular momentum of the accretor, J A , and donor, J D . Foreach set of model simulations — for example, the simulations labeled Q1.3P that are discussed in § q is naturally dimensionless in both codes; a , J orb , J D , and J A are in code units (see Table 1); time is normalized to the binary system’s original orbital period, P ; and ˙ M D has been normalized by a mass transfer rate corresponding to the initial mass of thedonor in the time span of the initial orbital period, that is,˙ M norm ≡ M D ( t = 0) P . (A1)For each model, the value of this normalization quantity also can be determined from the informa-tion provided in Table 1.The simulation data from both numerical techniques were analyzed with the same tools in aneffort to ensure that results from both numerical schemes were viewed through the same quantities.We chose to use diagnostic tools that were initially developed to analyze data sets from the cylin-drical coordinate-based grid-code simulations. In order to enable analysis of the SPH simulationswith this same set of diagnostic tools, data sets from the SPH simulations for mass, velocity, etc.were mapped onto a Cartesian grid (using a nearest grid point algorithm). Time series data forthe six key variables were then computed for each simulation as follows.The maximum density location for each star is identified and the line of centers of the binaryis assumed to run through the two density maxima. A plane is then constructed parametricallythat is perpendicular to the line of centers and divides the computational domain into two regions(donor and accretor) near the inner Lagrange point to enable identification of fluid elements thatbelong to the donor or accretor. The mass of each component and the binary separation are thencomputed from simple volume integrals of the density and density-weighted position of each star.The volume integrals extend over the entire grid region ascribed as donor or accretor based onthe dividing plane. The momentum and mass within a given grid cell are combined to give thethree components of velocity for the material in each cell. The orbital angular momentum is thencalculated from the center of mass position and velocity of the two stellar components while thespin angular momentum of each component is computed by referencing velocities to the relevantcenter of mass. The mass transfer rate is computed as the centered finite difference of the mass ofthe donor star with respect to time and is quite noisy. The mass transfer rates shown in the figuresare smoothed using a sliding box car average with a width equal to three times the initial orbitalperiod, P . Hence the first data point for the mass transfer rate shown is at 1 . P after the startof the simulation. 31 – REFERENCES
Amaro-Seoane, P., Aoudia, S., Babak, S., Bin´etruy, P., Berti, E., Boh´e, A., Caprini, C., Colpi, M.,Cornish, N. J., Danzmann, K., Dufaux, J.-F., Gair, J., Hinder, I., Jennrich, O., Jetzer, P.,Klein, A., Lang, R. N., Lobo, A., Littenberg, T., McWilliams, S. T., Nelemans, G., Petiteau,A., Porter, E. K., Schutz, B. F., Sesana, A., Stebbins, R., Sumner, T., Vallisneri, M., Vitale,S., Volonteri, M., Ward, H., & Wardell, B. 2013, GW Notes, Vol. 6, p. 4-110, 6, 4Badenes, C., & Maoz, D. 2012, ApJ, 749, L11Benz, W. 1990, in Numerical Modelling of Nonlinear Stellar Pulsations Problems and Prospects,ed. J. R. Buchler, 269–288Brown, W. R., Kilic, M., Hermes, J. J., Allende Prieto, C., Kenyon, S. J., & Winget, D. E. 2011,ApJ, 737, L23Cohl, H. S., Rau, A. R. P., Tohline, J. E., Browne, D. A., Cazes, J. E., & Barnes, E. I. 2001,Phys. Rev. A, 64, 052509Cohl, H. S., & Tohline, J. E. 1999, ApJ, 527, 86Dan, M., Rosswog, S., Br¨uggen, M., & Podsiadlowski, P. 2014, MNRAS, 438, 14Dan, M., Rosswog, S., Guillochon, J., & Ramirez-Ruiz, E. 2012, MNRAS, 422, 2417Diehl, S., Fryer, C., & Herwig, F. 2008, in Astronomical Society of the Pacific Conference Series,Vol. 391, Hydrogen-Deficient Stars, ed. A. Werner & T. Rauch, 221Diehl, S., Rockefeller, G., Fryer, C. L., Riethmiller, D., & Statler, T. S. 2012, ArXiv e-printsDiehl, S., & Statler, T. S. 2006, MNRAS, 368, 497D’Souza, M. C. R., Motl, P. M., Tohline, J. E., & Frank, J. 2006, ApJ, 643, 381Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition, 398Fryer, C. L., Rockefeller, G., & Warren, M. S. 2006, ApJ, 643, 292Fryer, C. L., Ruiter, A. J., Belczynski, K., Brown, P. J., Bufano, F., Diehl, S., Fontes, C. J., Frey,L. H., Holland, S. T., Hungerford, A. L., Immler, S., Mazzali, P., Meakin, C., Milne, P. A.,Raskin, C., & Timmes, F. X. 2010, ApJ, 725, 296Fryer, C. L., & Warren, M. S. 2002, ApJ, 574, L65Gokhale, V., Peng, X. M., & Frank, J. 2007, ApJ, 655, 1010Hachisu, I. 1986, ApJS, 62, 461Hachisu, I., Eriguchi, Y., & Nomoto, K. 1986, ApJ, 311, 214 32 –Han, Z. 1998, MNRAS, 296, 1019Herant, M. 1994, Memorie della Societa Astronomica Italiana, 65, 1013Iben, Jr., I., & Tutukov, A. V. 1984, ApJS, 54, 335Kilic, M., Brown, W. R., Allende Prieto, C., Kenyon, S. J., & Panei, J. A. 2010, ApJ, 716, 122Kilic, M., Brown, W. R., Kenyon, S. J., Allende Prieto, C., Andrews, J., Kleinman, S. J., Winget,K. I., Winget, D. E., & Hermes, J. J. 2011, MNRAS, 413, L101Marsh, T. R., Nelemans, G., & Steeghs, D. 2004, MNRAS, 350, 113Monaghan, J. J. 1992, ARA&A, 30, 543—. 2005, Reports on Progress in Physics, 68, 1703Motl, P. M., Frank, J., Tohline, J. E., & D’Souza, M. C. R. 2007, ApJ, 670, 1314Motl, P. M., Tohline, J. E., & Frank, J. 2002, ApJS, 138, 121Nelemans et al. 2005, A&A, 440, 1087Nelson, A. F. 2006, MNRAS, 373, 1039Pakmor, R., Edelmann, P., R¨opke, F. K., & Hillebrandt, W. 2012, MNRAS, 424, 2222Pakmor, R., Kromer, M., Taubenberger, S., & Springel, V. 2013, ApJ, 770, L8Pakmor, R., Springel, V., Bauer, A., Mocz, P., Munoz, D. J., Ohlmann, S. T., Schaal, K., & Zhu,C. 2016, MNRAS, 455, 1134Parsons, S. G., Marsh, T. R., G¨ansicke, B. T., Drake, A. J., & Koester, D. 2011, ApJ, 735, L30Postnov, K. A., & Yungelson, L. R. 2006, Living Reviews in Relativity, 9, 6Rappaport, S., Joss, P. C., & Webbink, R. F. 1982, ApJ, 254, 616Raskin, C., Scannapieco, E., Fryer, C., Rockefeller, G., & Timmes, F. X. 2012, ApJ, 746, 62Raskin, C., Scannapieco, E., Rockefeller, G., Fryer, C., Diehl, S., & Timmes, F. X. 2010, ApJ, 724,111Raskin, C., Timmes, F. X., Scannapieco, E., Diehl, S., & Fryer, C. 2009, MNRAS, 399, L156Soberman, G. E., Phinney, E. S., & van den Heuvel, E. P. J. 1997, A&A, 327, 620Staff, J. E., Menon, A., Herwig, F., Even, W., Fryer, C. L., Motl, P. M., Geballe, T., Pignatari,M., Clayton, G. C., & Tohline, J. E. 2012, ApJ, 757, 76 33 –Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753Thomas, P. A., & Couchman, H. M. P. 1992, MNRAS, 257, 11Wang, J., & White, S. D. M. 2007, MNRAS, 380, 93Warren, M. S., & Salmon, J. K. 1995, Computer Physics Communications, 87, 266Webbink, R. F. 1984, ApJ, 277, 355Young, P. A., Fryer, C. L., Hungerford, A., Arnett, D., Rockefeller, G., Timmes, F. X., Voit, B.,Meakin, C., & Eriksen, K. A. 2006, ApJ, 640, 891Yungelson, L. R., Livio, M., Tutukov, A. V., & Saffer, R. A. 1994, ApJ, 420, 336Zhu, C., Chang, P., van Kerkwijk, M. H., & Wadsley, J. 2013, ApJ, 767, 164
This preprint was prepared with the AAS L A TEX macros v5.2.
34 –Table 1. Initial (Polytropic) Binary System Parameters
Model q κ D ρ , D R D κ A ρ , A R A R A R circ M tot a sep P = π Ω Analog a Notes b Q1.3 1 .
323 0 . .
600 0 . . .
000 0 . .
602 3 . × − . .
74 MS (1)Q0.7a 0 .
700 0 . .
608 0 . . .
000 0 . .
302 2 . × − . .
20 WD · · ·
Q0.7b 0 . .
02 0 . . .
02 1 . . .
584 1 . × − .
75 31 .
70 WD · · ·
Q0.5 0 .
500 0 .
016 0 .
235 0 . .
016 1 .
000 0 . .
496 9 . × − . .
52 WD (2)Q0.4a 0 . . .
710 0 . . .
000 0 . .
768 2 . × − . .
75 WD (3)Q0.4b 0 . . .
740 0 . . .
000 0 . .
811 2 . × − . .
05 WD · · · a Stellar analog of polytropic binary where MS means fully-convective, low mass, main sequence stars and WD means a white dwarf where bothcomponents are low in mass compared to the Chandrasekhar limit. b Information drawn from: (1) Table 4 of DMTF06; (2) Table 5 of DMTF06; (3) Table 1 of MFTD07.
35 –Table 2. Simulation Prescriptions.Simulation a Model Code Resolution b Driving Notes c ID Duration(1) (2) (3) (4) (5) (6)Q1.3P G . P (1) S . P S . P Q0.7P G . P G . P S . P Q0.5P G . P (2) S . P Q0.4P G . P G . P G . P (3) S . P S . P Q1.3I G . P S . P Q0.7I G . P S . P Q0.5I G . P S . P Q0.4I G . P G . P S . P S . P Simulations labelled with “P” were evolved with a polytropicequation of state while simulations labelled with “I” utilized an idealgas equation of state. b For SPH simulations “resolution” means number of particles. Forgrid simulations, “resolution” means number of uniform cylindricalgrid cells; ( N R , N θ , N Z ) = (162 , , .
36 – c Identical to model simulation (1) Q1.3-D, first discussed in § § § . P .
37 –Table 3. Simulation Diagnostics.Simulation Model t merge t zpt Movie Outcome a ElaborationID P P ID(1) (2) (3) (4) (5) (6) (7)Q1.3P G .
38 0 .
00 video01 TI Binary merges at 12 . P S .
57 +4 .
81 video02 TI Binary merges at 7 . P S .
17 +4 .
21 video03 TI Binary merges at 8 . P Q0.7P G .
70 0 .
00 video04 TI Donor disrupts at 9 . P G . − .
41 video05 TI Donor disrupts at 21 . P S . − .
11 video06 TI Donor disrupts 11 . P Q0.5P G · · · .
00 video07 D Binary detaches at 34 P S · · · .
00 video08 ? Undetermined after 14 . P Q0.4P G · · · .
00 video09 D Binary detaches at 40 P G · · · .
00 video10 ? Undetermined after 50 . P G · · · .
00 video11 D Binary detaches at 46 . P S · · · .
00 video12 ? Undetermined after 18 . P S · · · . · · · ? (no movie)Q1.3I G .
54 0 .
00 video13 TI Binary merges at 15 . P S .
69 +7 .
85 video14 TI Binary merges at 7 . P Q0.7I G .
37 0 .
00 video15 TI Donor disrupts at 12 . P S .
04 +2 .
33 video16 TI Donor disrupts at 10 . P Q0.5I G · · · .
00 video17 CE + TI Donor disrupts at 19 . P S · · · .
00 video18 CE + TD? Undetermined after 15 . P Q0.4I G · · · .
00 video19 CE + TI Donor disrupts at 8 . P G · · · .
00 video20 CE + TI Donor disrupts at 12 . P S · · · .
00 video21 CE + TI Donor disrupts at 21 . P S · · · . · · · ? (no movie) a TI = tidal instability; TD = tidal disruption; D = detaches; CE = common envelope 38 –Fig. 1.— Slice through the meridional plane of the Q0 . a SPH setup, only particles within onesmoothing length of the plane are shown. Note how the outer layers of the donor (right side) aremuch better resolved with a smaller particle spacing, whereas the accretor (left side) is set up witha uniform particle density. 39 –Fig. 2.— The mass profile of each stellar component taken from spherical shells centered on thecenter of mass of the accreting star (left column) and donor star (right column) for the four initialbinary configurations Q1.3, Q0.7a and Q0.7b, Q0.5, and Q0.4a. The SPH initial data is plottedin blue while the Eulerian initial data is shown in red. Both codes are thus proceeding from verysimilar initial states and in the Eulerian representation, the donor stars (which nearly reach theirRoche lobe) are resolved to a somewhat lower density level. 40 –Fig. 3a.— Simulation Q1.3P. Curves show the time-dependent behavior of six system parameters— orbital separation, a , binary mass ratio, q , log of the donor mass-transfer rate, | ˙ M D | , orbitalangular momentum, J orb , and spin angular momentum of the accretor, J A , and donor, J D — asderived from models Q1.3P G S S q = 0 . t zpt , as detailed in column 4 of Table 3; units usedalong the vertical axis for each plot are defined in Appendix A. 41 –Fig. 3b.— Simulation Q1.3P. Images show contours of equatorial-plane column densities frommodel Q1.3P G S left-most images : t G = 12 . P and t S = 7 . P ; right-mostimages : t G = 12 . P and t S = 7 . P . The top two images have been extracted from video01and the bottom two have been extracted from video02. 42 –Fig. 4a.— Simulation Q0.7P. Curves show the time-dependent behavior of six system parameters— a , q , log | ˙ M D | , J orb , J A and J D — as derived from models Q0.7P G G S q = 0 . G S t zpt , as detailed in column 4 of Table 3; units used along the vertical axis for each plot aredefined in Appendix A. 43 –Fig. 4b.— Simulation Q0.7P. Images show contours of equatorial-plane column densities fromsimulation Q0.7P G G S left-most images : t G1 = 12 . P , t G2 = 20 . P , and t S1 = 11 . P ; right-most images : t G1 = 13 . P , t G2 = 21 . P , and t S1 =11 . P . The images displayed for simulations Q0 . G . G . S Q . Q . F i g . . — S i m u l a t i o n s Q . P a nd Q . P . C u r v e ss h o w t h e t i m e - d e p e nd e n t b e h a v i o r o f s i x b i n a r y s y s t e m p a r a m e t e r s : a , q , l og | ˙ M D | , J o r b , J A a nd J D . T h e un i t s u s e d t o n o r m a li zee a c h q u a n t i t y a r e d e fin e d i n A pp e nd i x A . L e f t P a n e l: R e s u l t s d e r i v e d f r o mm o d e l s Q . P G ( r e d c u r v e s ) a nd Q . P S ( b l u ec u r v e s ) . R i g h t P a n e l: R e s u l t s d e r i v e d f r o mm o d e l s Q . P G ( r e d s o li d c u r v e s ) , Q . P G ( d a s h e d r e d c u r v e s ) a nd Q . P S ( b l u ec u r v e s ) .
45 –Fig. 6.— Accompanying the § . P ; for comparison,dashed curves present results from simulations Q0.4P G S . P . 47 –Fig. 8.— Equatorial-plane surface-density contours from (bottom) model Q0.4P S left-most images : t Q04D = 13 . P and t S1 = 16 . P ; right-most images : t Q04D = 15 . P and t S1 = 18 . P . The bottom two images havebeen extracted from our video08; a movie showing the development of the structures displayed inthe top two images has been published in MFTD07. 48 –Fig. 9.— The images presented exemplify the qualitative difference between polytropic and idealgas simulations. A hot tenuous envelope develops around one or both stars due to shock heatingin the ideal gas simulations. Four selected frames from four simulations starting with initial modelQ0.5 are shown: top row grid results, bottom row SPH results, left column polytropic simulations,right column ideal gas simulations. Since t zpt is undefined in this case, the four frames were selectedby eye to represent a similar stage in their time evolution. 49 –Fig. 10a.— Simulation Q1.3I. Curves show the time-dependent behavior of six system parametersas derived from models Q1.3I G S t zpt = 7 .
85. 50 –Fig. 10b.— Simulation Q1.3I. Images show contours of equatorial-plane column densities frommodel Q1.3I G S left-most images : t G = 15 . P and t S = 7 . P ; right-mostimages : t G = 15 . P and t S = 7 . P . 51 –Fig. 11a.— Simulation Q0.7I. Curves show the time-dependent behavior of six system parametersas derived from models Q0.7I G S t zpt = − .
7. 52 –Fig. 11b.— Simulation Q0.7I. Images show contours of equatorial-plane column densities frommodel Q0.7I G S left-most images : t G = 11 . P and t S = 9 . P ; right-most images : t G = 12 . P and t S = 10 . P . 53 – Q . Q . F i g . . — S i m u l a t i o n s Q . I a nd Q . I . S a m e a s F i g u r e e x ce p t f o r i d e a l - ga s e q u a t i o n o f s t a t e . L e f t P a n e l: R e s u l t s d e r i v e d f r o mm o d e l s Q . I G ( r e d c u r v e s ) a nd Q . I S ( b l u ec u r v e s ) . R i g h t P a n e l: R e s u l t s d e r i v e d f r o mm o d e l s Q . I G ( r e d c u r v e s ) a nd Q . I S ( b l u ec u r v e s ) .
54 –Fig. 13.— As in Figure 7, curves show the time-dependent behavior of six binary system parametersfrom model simulations Q0.4I. Blue curves are derived from SPH-code simulations and red curvesare derived from grid-code simulations. Solid curves exactly reproduce information provided in theright panel of Figure 12 from simulations that were driven into contact for 1 . P ; for comparison,dashed curves present results from simulations Q0.4I G S . P . 55 – F i g . . — S e v e n i m ag e s e x t r a c t e d f r o m v i d e o19a nd v i d e o21 t h a t ill u s t r a t e (t o p r o w ) s i m u l a t i o n Q . I S v e r t h e t i m e p e r i o d . ≤ t / P ≤ . nd ( b o tt o m r o w ) s i m u l a t i o n Q . I G v e r t h e t i m e p e r i o d . ≤ t / P ≤ . , du r i n g t h ee p i s o d e o f t i d a l d i s r up t i o n o f t h e d o n o r . I nb o t h c a s e s , a d j a ce n t i m ag e s a r e s e p a r a t e d i n t i m e b y a pp r o x i m a t e l y . P a nd t h e i m ag e s e q u e n cec o v e r s a t o t a l o f . P0