Numerical Simulations of Mass Transfer in Binaries with Bipolytropic Components
Abstract
We present the first selfconsistent, three dimensional study of hydrodynamic simulations of mass transfer in binary systems with bipolytropic (composite polytropic) components. In certain systems, such as contact binaries or during the common envelope phase, the coreenvelope structure of the stars plays an important role in binary interactions. In this paper, we compare mass transfer simulations of bipolytropic binary systems in order to test the suitability of our numerical tools for investigating the dynamical behaviour of such systems. The initial, equilibrium binary models possess a coreenvelope structure and are obtained using the bipolytropic selfconsistent field technique. We conduct mass transfer simulations using two independent, fully threedimensional, Eulerian codes  FlowER and Octotiger. These hydrodynamic codes are compared across binary systems undergoing unstable as well as stable mass transfer, and the former at two resolutions. The initial conditions for each simulation and for each code are chosen to match closely so that the simulations can be used as benchmarks. Although there are some key differences, the detailed comparison of the simulations suggests that there is remarkable agreement between the results obtained using the two codes. This study puts our numerical tools on a secure footing, and enables us to reliably simulate specific mass transfer scenarios of binary systems involving components with a coreenvelope structure.
keywords:
binaries: close – hydrodynamics – methods: numerical1 Introduction
About half of the stars in the sky are in reality binaries or multiple star systems (Trimble, 1983), and with an expected surge in detection of transient events due to largescale surveys, it becomes more relevant to understand the nature of transient phenomena involving binary interactions. In recent years there has been considerable interest in modelling dynamical interactions of binary systems involving main sequence as well as post main sequence stars. For example, Nandez et al. (2014) conducted a smooth particle hydrodynamics (SPH) study to simulate the evolution of the progenitor binary of V1309 Scorpii, which was a contact binary merger event, through the common envelope phase with a red giant primary and a degenerate as well as a main sequence secondary. This project is motivated by our interest in the understanding of the mass transfer events in binary systems, and in particular, the evolutionary status, internal structure as well as stability criteria of contact binary systems. The contact binaries in equilibrium exhibit many apparently contradictory properties such as Kuiper’s paradox (conformity to Roche geometry, Kuiper (1941)), Lucy’s paradox (equal surface temperature of unequal mass components, Lucy (1968)) and secular stability over billions of years (Li et al., 2007). Stȩpień (2009) postulates a model for energy transported from the primary to the lower mass component by means of a steadily circulating stream in the equatorial region, in order to resolve Lucy’s paradox, while also maintaining thermal stability of the system. The dynamic stability and exact nature of the largescale circulations has not yet been verified from first principles with models which are fully selfconsistent. Modelling the components with a coreenvelope structure may help elucidate the dynamical stability of these flows in contact binary systems.
With the above issues in the theoretical understanding of contact binary systems in mind, in this paper we describe a set of gridbased numerical tools suitable for studying the hydrodynamic evolution of mass transferring binary systems with each component having a core and an envelope. One way to model the internal structure of such a star is using bipolytropes (sometimes also referred to as composite polytropes), where the core and the envelope of the star are allowed to have different polytropic indices (Henrich & Chandrasekhar, 1941). We are investigating the dynamical evolution of binaries involving stars having a coreenvelope structure, and, in particular, contact binaries with hydrodynamic simulations in a series of three papers. In Kadam et al. (2016) (henceforth paper I) we established the viability of the selfconsistent field method to produce single rotating bipolytropic structures in equilibrium. In the current paper (paper II) we present a set of Eulerian numerical tools that we have developed for investigating the dynamical phases of mass transfer in bipolytropic binary systems. A suite of hydrodynamic simulations of bipolytropic binary systems is conducted in order to test the suitability of our numerical methods for investigating such scenarios. We plan to focus on contact binary systems, the evolution of the progenitor of V1309 Scorpii, in paper III.
There are two key steps in order to conduct fully three dimensional hydrodynamic simulations of binary systems. In the first step, we build the initial models using a modified version of Hachisu’s selfconsistent field technique (Hachisu, 1986b), which we call the bipolytropic selfconsistent field (BSCF) technique. A bipolytrope is a structure with two polytropes stacked on top of each other, so that the core and the envelope have two different equations of state. The composition difference between the core and the envelope can be mimicked by specifying a jump in the molecular weight in the BSCF method. In the next step, we use the bipolytropic binary systems as quiet initial conditions and study their dynamic evolution using suitable hydrodynamic codes. We used two independent, explicit, Eulerian codes – FlowER (Motl et al., 2002; D’Souza et al., 2006) with a fixed cylindrical grid and Octotiger (Marcello et al., 2017) with a Cartesian AMR (adaptive mesh refinement) grid. The binary systems were chosen to match closely across the two codes, so that the results could be compared directly and could be used as a benchmark for each code. In the first set of simulations the initial mass ratio of the binary systems was set to 0.7, so that with a convective envelope of the donor, the mass transfer would be unstable. These simulations were conducted until merger at two grid resolutions with both the codes. In order to compare behaviour of the codes for binaries undergoing stable mass transfer, we also conducted simulations with an initial mass ratio of 0.6 for a reasonably large number of orbits. We examined the time evolution of key parameters, including separation, mass transfer rate, mass ratio, orbital and spin angular momenta, as well as density cross sections. The dynamical behaviour of the systems agreed remarkably well and the differences can be accounted for by considering the differences between the two codes, especially the inclusion of shock heating in Octotiger and the limited resolution of the stellar cores in FlowER. This study proved the reliability of our numerical techniques and highlighted the suitability of Octotiger for investigating dynamical binary interactions with bipolytropic components.
The structure of this paper is as follows. Section 2 describes the BSCF method that is used for generating the initial models. In section 3 we describe the two hydrodynamics codes that are used for the dynamical evolution of the binary systems. Section 4 explains the setup and properties of the suite of binary systems used as initial conditions for the simulations. In section 5 we describe how the binary parameters are calculated, and then compare and analyse the simulation results for both unstable as well as stable mass transfer scenarios. In section 6 we summarize our findings and discuss the applicability of our numerical tools for future investigations.
2 Bipolytropic Binaries
First introduced by Ostriker & Mark (1968), the selfconsistent field (SCF) technique is an iterative method of finding the equilibrium solutions of a rotating, gravitationally bound fluid. This method was further improved in Hachisu (1986a, b) to produce selfgravitating, rapidly rotating stars as well as systems without axial symmetry, such as binary or multiple star systems. This improved method of calculating equilibrium structures is called Hachisu’s selfconsistent field (HSCF) technique. In the past, the HSCF technique has been used to construct the initial models of detached, close binary systems with a wide range of mass ratios and Roche lobe filling factors ( Motl et al. (2002); D’Souza et al. (2006); Even & Tohline (2009)). The binary systems that were generated had either a single polytropic or a zerotemperature white dwarf equation of state for both the stars, and they were subsequently used as initial conditions for the study of their dynamical behaviour. In this section we describe a computational method to obtain the equilibrium configuration of close binary systems having a bipolytropic structure for each component. This is achieved by using the bipolytropic selfconsistent field (BSCF) method. This technique is similar to the method described in paper I, which is used for constructing rapidly rotating single bipolytropic spheroidal as well as toroidal structures in equilibrium.
2.1 Implementation
We make the following assumptions while using the BSCF technique for a binary system. The stars are assumed to be rotating synchronously in a circular orbit with a uniform angular velocity () such that the system is in hydrostatic equilibrium in the frame of reference corotating with the system. The basic equations describing the internal structure of each bipolytropic star are similar to those described in paper I.
The polytropic equation of state assumes a power law relationship between the pressure and the density of the gas
(1) 
where is the polytropic index, is the polytropic constant and is the polytropic exponent. A bipolytropic star has one such polytropic index for the core () and one for the envelope (). Since the pressure and temperature are continuous functions across the coreenvelope interface, this implies
(2) 
where , are the densities at the interface and is the ratio of average molecular weights between the core and the envelope  and respectively. In the corotating frame of reference, the steadystate equation of motion for each stellar component, , in a binary can be written as
(3) 
where is the specific enthalpy of the gas, is the gravitational potential and is the coordinatedependent part of the centrifugal potential. Here are the integration constants which are analogous to the Bernoulli constant in classical fluid flows. These constants are different for each polytrope, in general, hence in a bipolytropic binary configuration there are four such constants. The density and hence the enthalpy is zero at the surface of each star.
The BSCF method for obtaining a binary system requires seven boundary conditions in order to solve the hydrostatic equations for the two stars along with Poisson’s equation. These are depicted in Figure 1. The boundary conditions are chosen to be the maximum densities of the two stars , , the location of the outer boundary of star 1 (A), the location of the inner boundary of the same star (B) and the location of the inner boundary of star 2 (C). Here subscript denotes the star on the right in Figure 1, so that two of its boundaries are specified. The three points are collinear and lie along the line joining the centre of mass of the two stars. In order to obtain a bipolytropic binary, the locations of the boundary points of the core for both the stars (points D and E) need to be specified. In addition to the above, the polytropic indices of the core and the envelope, as well as the ratio of molecular weights between them for both the stars need to be specified, the choices of which depend on the type of stars one wants to represent.
An iterative scheme to solve Equation 3, along with Poisson’s equation in the BSCF method works as follows. Note that the integration constants can be different for the core and the envelope of each star. An initial density distribution satisfying all the outer boundary conditions must be guessed for both the stars. The gravitational potential can be calculated by solving Poisson’s equation. The numerical technique to solve this equation is identical to that used during the hydrodynamical evolution, and is described in section 3 for each code. Since the density, and hence the enthalpy, at the surface is forced to vanish, the surface boundary conditions at point A, B and C are
(4) 
The envelope integration constant is the same at points A and B, which allows us to calculate the angular velocity of the binary,
(5) 
When we put this value of angular velocity back in Equation 3 for points A and C, the envelope integration constants for both the stars, and , can be found.
Now we are at the same point as Equation 18 in paper I, where the rotational velocity of the system and the envelope integration constant for each star are known. The new density distribution for each star can be found by following similar steps. In the next iteration cycle, the axis of rotation is allowed to move to the location of the centre of mass of the system along the line joining the centre of mass of each star. The centrifugal potential is calculated with respect to this new axis of rotation. The computation is performed for two nonoverlapping regions for star 1 and 2, such that each star lies completely inside only one of the regions. The iteration procedure is repeated until the relative changes in , , , and are smaller than a prescribed convergence criterion, .
The quality of a converged solution is quantified using the virial error
(6) 
where T is the total kinetic energy, W is the total gravitational energy and is the volume integral of the pressure,
(7) 
VE is a global measure of the accuracy of the equilibrium of a nonlinear dynamical system and should ideally be zero. The calculations are performed in dimensionless form of the variables, or “code units", where the maximum density on the grid (), the gravitational constant () and a suitable length, typically the size of the computational grid, are set to unity.
Similar to the HSCF method, the BSCF technique suffers from one major drawback. Although close binaries with a wide range of mass ratios and core mass fractions can be generated using the BSCF technique, physical quantities of interest such as the mass ratio, Roche lobe filling factors and the core mass fractions cannot be specified a priori or controlled precisely. The control parameters , and locations of points A through E need to be adjusted through trial and error, to produce systems which have desired properties. Note that there are no analytical solutions for binary bipolytropes (or polytropes) to compare the converged solutions with. However, one way to test the viability of the solutions obtained is to evolve these models with a hydrodynamic code and compare the results with our expectations from well known analytical approximations.
3 Numerical Algorithms
Once an equilibrium binary system is generated using the BSCF method, it can be used as a quiet initial condition and can be evolved numerically using a suitable hydrodynamic code. We conduct the binary evolution using two explicit, gridbased (Eulerian), fully threedimensional, hydrodynamic codes – FlowER and Octotiger. This section briefly describes these two computational fluid dynamics algorithms.
3.1 The FlowER Code
FlowER is an explicit, conservative, finitevolume, threedimensional Eulerian code, which is second order accurate in both time and space. The code is similar to the ZEUS code (Stone et al., 1992) and is implemented on a uniform cylindrical grid. The code tracks the following fluid quantities: mass density (), angular momentum density (), vertical momentum density () and radial momentum density (). In addition, the entropy tracer
(8) 
where is the internal energy per unit mass and is the ratio of specific heats, is also tracked with an assumption of adiabatic evolution. The set of differential equations solved correspond to the conservation laws of above four quantities, in addition to the internal energy equation for in a reference frame rotating with a constant angular velocity.
The gravitational potential at each timestep is obtained by numerically solving Poisson’s equation with the same technique as described in paper I. An ideal gas equation of state is used to close the system of equations
(9) 
Here, is the polytropic exponent used while evolving the system in time in the hydrodynamic code. All simulations in this paper are performed using a constant value of , so that the fluid is convectively stable. The outer boundary conditions are outflow and are implemented in such a way that the material is free to move out of the grid, but it is not allowed to move back in from the outermost cells. At the axis of the cylindrical grid, the innermost annulus of cells is averaged in the azimuthal direction at a given layer in the vertical direction. For a detailed description of the treatment of the advection terms, source terms, time centreing and numerical tests, see Motl et al. (2002) and D’Souza et al. (2006).
With the BSCF method, the centre of mass is allowed to drift along the line joining the two stars. We transport the density field from the BSCF grid to the FlowER grid such that the centre of mass coincides with the coordinate axis of the cylindrical grid as closely as possible. The translation is implemented on the pressure field and the density is recalculated at the end, in order to preserve the discontinuity at the interface of the core and the envelope. This translation is performed using a bicubic spline interpolation scheme described in Press et al. (1986).
Consider a star placed in a cylindrical computational grid. If the star is placed such that its centre lies at the axis of the grid, the resolution in the azimuthal direction equals that of the grid. However, if the star is placed away from the axis, as a component of a binary system, its azimuthal resolution decreases sharply and keeps getting worse as it moves farther from the axis. The core in a bipolytropic star is especially affected by this degradation of the resolution because of its relatively small size. Unlike a single bipolytrope placed at the centre of the cylindrical grid, or a polytropic binary, we noticed that hydrodynamic evolution of bipolytropic binary systems in FlowER showed a numerical artefact– the central density (defined as the maximum density of a star) of each bipolytropic component showed a fall of the order of a few per cent on a dynamical timescale. To investigate this issue further, we conducted a set of simulations with noninteracting, bipolytropic binaries having a mass ratio of unity. In conclusion, we found that after a threshold, the fall in the central density was roughly proportional to the absolute slope of the density near its maximum. This implied that if the cores of bipolytropic binaries were not sufficiently resolved, the finite volume method as implemented in FlowER with a cylindrical geometry could not preserve steep density gradients near the centre of a bipolytropic star. This resulted in numerical diffusion of bipolytropic cores, and was reflected by the fall of the central densities of both the binary components. The numerical diffusion can rearrange gas inside a star and affect the dynamical behaviour of the binary. Increasing the resolution of the FlowER simulations, however, is computationally very expensive. Although this resolutiondependent numerical effect appears to be a major drawback, the study of mass transferring bipolytropic binaries, as discussed in sections 5.2 and 5.3, shows that the overall behaviour of a binary system remains unaffected.
3.2 The Octotiger Code
Octotiger is a state of the art, explicit, conservative, hydrodynamic code with a Cartesian AMR grid, developed for modelling mass transfer in binary systems (Marcello et al., 2017). The AMR uses an octree data structure. Refined nodes have eight regularly spaced children, providing twice the resolution in their respective domains. The selfsimilarity in octree structures provides a relatively simple to implement AMR framework, as compared to other alternatives such as patchbased AMR. The equations that are numerically solved in Octotiger, conservation of mass and momentum as well as the internal energy equation (evolution of ) are similar to those in FlowER. The code also implements conservation of total energy explicitly, in a manner similar to Marcello & Tohline (2012). The ideal gas equation of state used in Octotiger is identical to that implemented in the FlowER code. However, the energy expression used for the treatment of shocks at high velocities, where the kinetic energy dominates the internal energy of the fluid, is different. In general, the specific internal energy is calculated as , where is the specific total energy inside a cell and is the velocity. In shocks, the very high velocities involved could make calculations of inaccurate, as . Accurate computation of the temperature is impossible in this case, since the pressure is proportional to the internal energy. Octotiger uses the dual energy formalism of Bryan et al. (1995) to address this issue. The treatment of the boundary in Octotiger similar to that in FlowER and the material is allowed to flow out but not back into the grid from the outermost boundary points. Thus FlowER evolves fluids adiabatically with an ideal gas equation of state, while Octotiger takes shock heating into account and uses the dual energy formalism with an ideal gas equation of state.
Gridbased hydrodynamic codes using finite volume methods typically conserve linear momentum or angular momentum, but not both. Octotiger solves Poisson’s equation using a modified Fast Multipole Method (FMM) with complexity where the potential is computed as a multipole expansion to third order. Octotiger implements a modified version of the FMM of Dehnen (2000) which conserves both linear and angular momenta to machine precision (Marcello, 2017). The hydrodynamic module also conserves linear and angular momenta to machine precision through the novel technique of Després & Labourasse (2014). The conservation of mass and total energy as well as both linear and angular momenta to machine precision is important for investigating long term, dynamical stability of an evolving binary system.
Octotiger is implemented in and parallelized using the High Performance ParalleX (HPX) runtime system developed in collaboration with the Stear^{1}^{1}1http://stellar.cct.lsu.edu/about/ (Systems Technology, Emergent Parallelism, and Algorithm Research) group at the Center for Computation & Technology (CCT) at the Louisiana State University (Kaiser et al., 2015). HPX allows efficient use of computational resources through better scaling and performance as compared to more conventional programming models such as MPI (Message Passing Interface).
4 The Initial Models
In this section we describe the initial models for the suite of bipolytropic binary simulations conducted using the two codes, FlowER and Octotiger. We conducted a total of six benchmark simulations (see Table 1). A set of four simulations started with an initial mass ratio , are called Q0.7; while two additional simulations started with an initial mass ratio are called Q0.6. The properties of the initial models were matched closely for a given mass ratio.
The choice of the internal structure of the initial models was affected by the limitations of the code FlowER, as well as our intention of testing the numerical tools with bipolytropic binaries, as opposed to simulating a real evolutionary scenario. Both Q0.7 and Q0.6 models were constructed such that the lower mass, secondary component will overflow its Roche lobe. As explained in section 3.1, in the case of FlowER simulations the steep density gradients near a star’s centre result in a proportional numerical diffusion of the core over a dynamical time. Thus for the purpose of benchmarking, we chose moderate values for the core polytropic indices of the bipolytropic components. The core and envelope polytropic indices for all stars in the simulations were . These values are suitable for a radiative core and a convective envelope. Since the gradient of the central density also increases with the increasing , we selected a value of for all stars. This value is close to the ratio of molecular weights between fully ionized helium core and a solar composition envelope (). We would like to emphasize that this internal structure of the bipolytropes is not fully consistent with what is expected from a real evolutionary scenario. However, the detailed study of simulations with bipolytropic binary components having deep convective envelopes is useful for code testing and comparison. Such benchmarking is essential so that the novel behaviour of simulations, with stars having realistic internal structure, can be trusted in the future.
Table 1 lists the resolutions used to carry out the benchmark simulations. The prefix “cyl" or “oct" denotes the code used (FlowER or Octotiger respectively) and postfix “LR" or “HR" denotes the use of low or high resolution respectively. The cylLR model had a grid resolution of . In the cylHR model, the resolution of each star was doubled in the direction, and approximately doubled in the other two directions. The total resolution of cylHR model was . In the FlowER grid, the mesh interval in the direction was set to that in the direction, . The domain of cylHR simulation was smaller in the radial direction as compared to that of cylLR, and contained only about 58% of the volume. The AMR grid in Octotiger changes dynamically according to the matter distribution, hence we specify only the initial effective resolution of each simulation. The resolution criteria in Octotiger are based on the total number of refinement levels specified as well as the density distribution. The core of both the stars were resolved to the finest refinement level, and both the envelopes are resolved to the next level. The rest of the grid was refined such that each cell has approximately an equal amount of mass. The initial effective resolution for an Octotiger simulation was obtained by multiplying the total number of nodes by the number of subgrids (), which equals the total number of cells on the grid. The model in octHR had two more levels of refinement as compared to octLR and its computational domain was also twice as large (eight times more volume), which effectively doubled the resolution of the cores of both the donor and the accretor. For comparison, the core resolutions are listed in Table 1. In the case of the FlowER code, the core resolution is specified as the number of cells in the radial and the azimuthal direction, and in the case of Octotiger, the core resolution is specified as the number of cells in the X and Y direction. Distances in all four simulations were normalized with respect to the radius of the BSCF grid in the cylLR model. Note that the size of the computational grid in a simulation can be different from that of the corresponding BSCF model. The sizes of the computational domains for the Q0.7 simulations are also listed in Table 1 in code units.
Benchmark  Model ID  Resolution ➀  Domain Size ➁  Core Resolution ➂  Notes 

Set  Donor Accretor  
Q0.7  cylLR  M  
cylHR  M  
octLR  M  8 levels of refinement  
octHR  M  10 levels of refinement  
Q0.6  cylEE  M  
octEE  M  8 levels of refinement 

➀ The resolution in the cylindrical grid is specified as . The initial resolution in the AMR grid is specified as total number of grid cells obtained by multiplying the total number of nodes by subgrids.
➁ The sizes of the computational domains are specified in code units. For the FlowER simulations the sizes are in and directions and for the Octotiger the sizes are in , and direction.
➂ The core resolution for cylindrical grid is in terms of the number of cells in and direction, and that for AMR grid is the number of cells in and direction.
The setup for the Q0.6 simulations was similar to the Q0.7 simulations at the low resolution, the “LR" models. The postfix “EE" in Table 1 denotes our intention of an extended evolution, in order to demonstrate that stable mass transfer can be handled adequately by our codes. Thus cylEE model had a grid resolution of in FlowER, while the corresponding model octEE had 8 levels of refinement. The internal structure of both the components was also similar to the Q0.7 models, and .
The initial values of the binary system parameters for the benchmark simulations are listed in Table 2 in code units. Here a subscript denotes the donor, denotes the accretor, denotes the core and denotes the envelope. The quantities listed in this table are– mass (), effective radius (), maximum density (), core mass fraction , rotational (spin) angular momentum , polytropic constant of the core () and the envelope () for each of the components, along with the mass ratio (), the Roche lobe filling factor by volume (Rlff), separation (), period (), orbital angular momentum , virial error (VE, see Equation 6) and the minimum density on the grid at the beginning of simulation (). The Rlff is defined in terms of the ratio of the volume of the component above a density threshold of to the volume of its Roche lobe. Although for a given mass ratio the constructed initial models closely match each other, one may notice that they are not identical. The accuracy of the converged solution depends on the resolution and the departure from equilibrium should tend to converge to zero with the increase in resolution. Thus even if the input parameters are precisely specified via the locations of the five boundary points specified in the BSCF method, obtaining identical models for this study is not possible. Based on our previous experience with mass transferring binary systems with polytropic components (Motl et al., 2017), this slight difference in the initial models should not affect the results. Assuming a primary with solar mass and radius, the properties of the initial binary can be calculated in solar units, which are listed in Table 3. Such binary systems involving lowmass stars (approximately ) are thought to remove angular momentum through magnetic braking in order to form a semidetached or contact configuration (Kawaler, 1988).
Parameter ➀ / Model ID  cylLR  cylHR  octLR  octHR  cylEE  octEE 


➀ The subscripts D, A, c and e stand for donor, accretor, core and envelope respectively.
Model ID  cylLR  cylHR  octLR  octHR  cylEE  octEE 

12.99 
Before analysing the results, let us consider the outcome that we may expect from analytical approximations for the two sets of simulations, Q0.7 and Q0.6. There are two critical mass ratios which determine the ultimate outcome of a conservative mass transfer event. The first critical mass ratio is for the separation; below this value the separation increases on mass transfer and above this value it decreases. The next critical mass ratio, , arises from the adiabatic response of the donor as well as the response of its Roche lobe to mass loss. This interplay is quantified by comparing the massradius exponent, , of the Roche lobe () to that of the adiabatic value for the Roche lobe filling star () (Webbink, 1985). Dynamically unstable mass transfer occurs when , and it continues in a quasisteady state, on a thermal or nuclear timescale otherwise. For a conservative mass transfer, (Tout et al., 1997). The donor in the benchmark simulations has a deep convective envelope with polytropic index, . Thus for our analysis, we can assume an approximate value of , corresponding to pure polytrope. This leads to the critical mass ratio, . For the benchmark simulations Q0.7 and Q0.6, the equals about 0.18 and 0.39 respectively. Thus mass transfer in Q0.7 simulation should be dynamically unstable, while that in the Q0.6 should be stable. Note that simplifying assumptions are made in the above analytical treatment, such as point mass potentials and constant mass ratio during mass transfer. The behaviour of a real system depends on fully three dimensional calculations. The initial mass ratio of 0.7 offers an interesting mass transfer case, as we can test behaviour of the simulations through the merger. The system with an initial mass ratio of 0.6 allows us to test if our tools are able to cope with simulations which undergo stable mass transfer. We shall see in the next section that both the results of the benchmark simulations are consistent with the above expectations. Conversely, consistent hydrodynamical evolution confirms that the initial bipolytropic binary models obtained through the BSCF method are of sufficient quality to be useful as quiet initial conditions for such investigations.
5 Mass Transfer Simulations
In this section we describe the Q0.7 and Q0.6 set of simulations, as well as the framework for analysing the results. Since the initial conditions were chosen to match closely, for the given mass ratio, we can expect agreement in the dynamical behaviour of these systems as well as the ultimate outcome. As explained in section 3, the code FlowER evolves adiabatically with an ideal gas equation of state, while Octotiger takes into account shockheating and uses the dual energy formalism with an ideal gas equation of state. The flow dynamics of the binary systems during mass transfer are not expected to be affected by these differences, as we shall also confirm in sections 5.2 and 5.3.
The BSCF method produces binary models which are detached in general. A driving mechanism was incorporated in each hydrodynamic code such that a certain fraction of specific angular momentum was artificially removed from the system per orbit. This in turn reduced the separation thus eventually initiating mass transfer. The lower mass star initially filled a larger fraction of its Roche lobe, as compared to its companion, in all of our models. Hence the mass transfer proceeded from the lower to the higher mass component. The Q0.7 set of simulations was driven for 4 initial orbits by removing 2% of the angular momentum per orbit, while the Q0.6 set was driven for 3 orbits at the rate of 1.5% of the angular momentum per orbit. Note that as a consequence of the driving, the binary components were rotating nonsynchronously at the onset of mass transfer. According to our experience of studying binaries with polytropic components, ( Motl et al. (2002), Marcello & Tohline (2012), Motl et al. (2017)) the driving phase does not affect the simulation results. The subsequent hydrodynamical evolution of the binary systems was treated without any special assumptions or symmetries.
5.1 Diagnostics
In order to study the evolution of the binaries in the hydrodynamic simulations, we analysed the temporal evolution of the following global parameters, called diagnostic quantities:

Binary separation,

Mass transfer rate normalized to the initial mass and period, , for the donor

Mass ratio, and normalized mass ratio,

Normalized orbital angular momentum,

Normalized spin angular momentum of the donor and accretor, ,

Maximum density (also called central density) of the donor and accretor, and
Here a subscript “0" denotes the initial value, “D" denotes the donor and “A" denotes the accretor. The diagnostic parameters were calculated in a similar way for the simulations conducted using both hydrodynamic codes. The donor and accretor material were separated by a dividing plane, normal to the line joining the maximum density of the two stars at the point. During the evolution, the point was considered to be the point of the maximum effective gravitational potential along the line between the two centres of mass. The separation was calculated as the distance between the centre of mass of the donor and the accretor. The centre of mass, as well as the mass for each component, was calculated for the total mass on the either side of the dividing plane. The spin angular momentum of each component ( and ) was calculated in the inertial frame with respect to its centre of mass. The orbital angular momentum was calculated by subtracting the spin angular momentum of each component from the total angular momentum of the system, all calculated in the inertial frame. In addition to the above diagnostic quantities, we studied the density distributions in the equatorial as well as meridional plane. A data file containing the density field, called a “frame", was written 120 times per initial orbit in FlowER and 100 times per initial orbit in Octotiger.
The time in our simulations is always normalized with respect to the initial orbital period for the respective model. For Q0.7 simulations, we synchronized the simulations with respect to a single event during the evolution, the merger of the binary occurring at time . We defined for the octHR simulation such that the the accretor spin angular momentum achieves the maximum value during the merger phase. We then compare the frame corresponding to for the octHR to the other simulations, by eye. The time corresponding to this closest matching frame is identified as for that particular simulation (see Figure 9). Once has been identified for each simulation, we could align all the timedependent plots for the Q0.7 simulations. The particular choice of the accretor spin angular momentum is discussed later in section 5.2. The choice of the temporal “zero point" is arbitrary and we chose the beginning of cylHR simulation as this point, since it completed the largest number of orbits. We shifted the diagnostic plots for the rest of the simulations such that the for each of them coincides. This shifted time was given by
(10) 
where
(11) 
for each simulation. We terminated the diagnostic plots soon after the .
In the case of Q0.6 simulations the binary systems exhibit a stable mass transfer throughout the evolution. Thus we do not have a clear event for synchronization and the diagnostic quantities are simply plotted from the beginning of the two simulations.
5.2 Hydrodynamic Evolution of Q0.7 Simulations
In this section we describe the evolution of Q0.7 set of binaries and discuss the results. As explained in section 4, the initial mass ratio lies above the critical value for the system, , hence we expect an unstable mass transfer over a dynamical timescale. The benchmark for these simulations was to evolve the initial bipolytropic binary to a fixed point in time, to merger. Table 4 describes the summary of the Q0.7 simulations along with some quantities of interest, and Table 5 lists the movies of the temporal evolution of the equatorial cross sections of the density fields. The timedependent behaviour of key diagnostic parameters are plotted in Figures 27, and Figures 16 and 17. We found that all four binary models were unstable to mass transfer, resulting ultimately in a merger. The end states as well as the intermediate behaviour agree very well with the simple analytical approximations and also among the simulations. Although there were a few key differences among the simulations, all of these can be accounted for by our current understanding of the numerical tools.
Model ID  cylLR  cylHR  octLR  octHR 

Equation of state ➀  Ideal gas  Ideal gas  Ideal gas with  Ideal gas with 
Dual Energy  Dual Energy  
Formalism  Formalism  
Total orbits ()  
➁  
➂  
Timesteps  M  M  k  k 
➃  
Number of cores  
Cost (CPUhr/orbit)  k  k  k  k 
Cost (calculations/orbit)  
Machine  QueenBee  QueenBee  SuperMIC  SuperMIC 
(LONI)  (LONI)  (LSU HPC)  (LSU HPC)  
Processors  2.33GHz 4Core  2.33GHz 4Core  2.8GHz 10Core  2.8GHz 10Core 
Xeon 64bit  Xeon 64bit  Ivy BridgeEP  Ivy BridgeEP  
Xeon 64bit  Xeon 64bit 

➀ Equation of state used in the simulation for determining the pressure. ➁ The time of merger (as determined in section 5.1). ➂ The difference between for cylHR simulation and the current one (Equation 6.2). ➃ The time needed to complete the simulation.
Model ID  Movies ➀  

cylLR  cylLR  
cylHR  cylHR  
octLR  octLRzoomed ➁  octLRwide ➂ 
octHR  octHRzoomed  octHRwide 

➀ The movies are available for viewing in the online supplemental material. ➁ Zoomed movies show equatorial cross sections at the length scale of the cylLR grid. ➂ Wide movies show the full octHR grid.
Consider the separation between the two components, plotted in Figure 2. Oscillations with an amplitude of approximately can be noticed throughout the evolution up to merger for all four simulations. In general, an initial binary model obtained through the BSCF method is not in perfect equilibrium and the orbits are not perfect circles. The Q0.7 binaries were further perturbed by the artificial removal of the angular momentum during the initial driving phase, exacerbating the situation. The slight deviations from equilibrium resulted in binary orbits with small eccentricity, in turn causing the oscillations of the components about the mean separation. The final amplitude of the oscillations also depends on their exact phase when the driving is stopped, hence the absolute phase of these oscillations should not be expected to match for any two simulations. The next prominent feature is an initial, roughly linear fall in the separation for the first 4 orbits, which was a direct result of the artificial removal of angular momentum, as each binary was driven together. The logarithmic derivative of the separation, , is proportional to twice that of the angular momentum, for a conservative mass transfer (Frank et al., 2002). A total of decrease in the total angular momentum thus resulted in about decrease in the separation, as observed in Figure 2. We made sure that a steady mass transfer stream was established before the driving was stopped. The behaviour of the Q0.7 systems agreed remarkably well after this point, after . The separation decreased gradually over about six orbits, until the binary merges. The merger was rapid and occurred within only about two orbits.
Figure 3 shows the normalized and smoothed mass transfer rate of the donor (). Note that the ordinate is in a logarithmic scale. We needed to use the smoothed rates, since the plots for the instantaneous mass transfer rate turned out to be very noisy because of the subtraction of two nearly equal quantities both in the numerator as well as the denominator. The smoothing was performed using a moving boxcar average with a width of three initial orbital periods. The relatively low amplitude oscillations in the before the sharp rise near the end are caused by the change in the mass transfer rate due to the epicyclic oscillations in the orbital separation. Once the driving ended and a steady mass transfer began at , increased steadily as the binaries evolved. The average rate of the normalized mass transfer after was of the order of 0.1 per orbit, which explains the lifetime of the binary with the merger occurring within about 8 orbits. The cylHR simulation lasted the longest because at the end of its driving, , while octLR simulation lasted the shortest time as the mass transfer rate at the end of its driving phase was the highest amongst the Q0.7 simulations, . The shifted plots for matched very well for all simulations after the end of the initial driving phase at , until about two orbits before the merger.
As shown in Table 4, the time taken from the beginning of the simulation to the merger, , varied by a large factor across the simulations. For the FlowER simulations, was longer as compared to the Octotiger simulations, and it was also longer at the higher resolutions for both the codes. We can explain the observed trend in as follows. The time of merger in a simulation depends on the average mass transfer rate, which in turn depends on the depth of initial contact as well as the resolution. First consider the trend across the two codes, since they have a key difference – Octotiger evolves the total energy equation explicitly which allows for shock heating of the gas. Thus in Octotiger simulations, when the accretion stream impacts the accretor, the supersonic velocities shockheat the outer atmosphere of the stars, forming a puffed up common envelopelike structure. This condition is exacerbated by the fact that radiative transport is not implemented in the code, which prevents the material from cooling. With a hot atmosphere, the mass transfer rate should be higher in Octotiger simulations, as it should scale with the sound speed of the fluid near the point. This common envelope can also possibly cause numerical drag on the binary, leading to a relatively rapid merger even at the same level of initial contact depth. As seen in Figure 3, all simulations showed a sharp increase in near the end, which indicates a rapid merger. However, during the final phase of the merger, there was a noticeable difference between the FlowER and the Octotiger simulations. We conjecture that the higher mass transfer rate during the merger process in the Octotiger simulations was also caused by the hotter atmosphere of the donor. For a given code, the higher resolution simulations lasted longer because a finer grid can resolve a lower mass transfer rate. We will show that the trends in the diagnostic plots support this hypothesis.
Figure 4 shows the normalized mass ratio of the Q0.7 simulations. The four curves in this plot essentially represent the history of mass transfer in each simulation. As the binaries evolved, decreased monotonically, and it continued to decrease even after driving was stopped, indicating a continuing mass transfer from the lower mass component to its companion. The plots for all Q0.7 simulations, except for the cylHR, show a high degree of agreement. The discrepancy between the evolution of the mass ratio of the cylHR simulation and the other simulations combined can be explained as follows. As a simulation progresses, if we assume the same mass transfer rate, the mass ratio depends on the initial masses of both the components in a binary as well as the amount of mass transferred from the donor to the accretor. Consider Figure 3 again, which shows that the mass transfer in the case of the cylHR simulation lasted almost twice as long as the other three simulations. Before we can expect an agreement between the codes, before , the donor in the cylHR simulation was already transferring mass for about 12.5 orbits. Over this period, the normalized mass transfer rate can be considered as approximately per orbit. Thus the donor had transferred about 2.5 % of its mass before , making the mass ratio smaller by the same factor, as observed in Figure 4. Notice that over the last two orbits, the Octotiger simulations show a steeper slope of as compared to the FlowER simulations which indicates a larger rate of mass transfer, as we have already seen in Figure 3.
The normalized orbital angular momenta are plotted in Figure 5, for the Q0.7 simulations. Since all binaries were driven for four orbits at the rate of 2% per orbit, the orbital angular momentum for each simulation initially showed a linear decline of a total of 8%. Figures 6 and 7 show the plots for the normalized spin angular momenta of the donor and accretor respectively. First consider the qualitative behaviour of these plots during the mass transfer. For all Q0.7 simulations, decreased steadily even after the driving was stopped and just before the merger, it fell precipitously over the last two orbits. We see that as decreased right before the merger, both and peaked sharply, indicating that the orbital angular momentum was being used to spin up each of the binary components. This behaviour of the angular momenta was a clear indication of tidal instability of the binary system. With this instability, each component tried to maintain synchronicity, which extracted angular momentum from the orbit, resulting in a merger. In certain cases the binary can develop a masstransfer instability, wherein the donor keeps transferring mass to the accretor in a dynamically unstable manner, such that the separation decreases and the donor essentially ends up “falling" on the accretor. Thus the spin angular momentum of the donor remains relatively constant until the merger, unlike the results of Q0.7 simulations. We determined on the basis of (section 5.1) because in all cases, the accretor angular momentum grows rapidly during the final stages of binary merger.
All three angular momenta (orbital and spin for both the components) generally agree across the Q0.7 simulations, however, there are two differences. The accretor angular momentum of the cylHR simulation starts to deviate early on from the corresponding curves of the other simulations, and the two codes behave differently during the merger over the last two orbits before . Consider the evolution of in the cylHR simulation. As seen earlier, the donor had transferred about 2.5% of its mass before . This mass also carried angular momentum, which resulted in spinning up of the accretor. The component of the velocity of the accreting material at the point, perpendicular to the line joining the two stars can be given as , where is the Roche lobe radius of the accretor. Using the values in Table 2, we can estimate that the % mass lost by the donor transferred about % of its angular momentum to the accretor before . The accretor also had sufficient time to interact with, and gain angular momentum from, . The combination of these two factors can explain the observed discrepancy of in Figure 7. During the final two orbits before the merger, the behaviour of the spin angular momentum of the donor as well as the accretor was different between the two codes. The sharper rise in the slope of both and in the case of the Octotiger simulations suggests that the orbital angular momentum was extracted more efficiently during the merger. This could be facilitated by a higher mass transfer rate near the end. Thus the difference in the spin angular momenta between the FlowER and Octotiger simulations is consistent with our analysis of the mass transfer rate earlier in this section. The agreement between the diagnostic plots indicates that the initial Q0.7 binary models generated using the BSCF methods were sufficiently closely matched.
The spin up of the stars at the expense of the orbital angular momentum was rapid and occurred over a dynamical timescale. From the initial values of the angular momenta in Table 2, it is clear that the spin angular momenta remained much smaller than the orbital angular momentum until the merger. Thus the Darwin (1880) instability can not play a role in the dynamical evolution. The tidal synchronization timescale for a noninteracting binary with parameters matching the initial systems, as estimated by Zahn (1977), is of the order of years. However, this is in the case of a detached binary system. The binaries in Q0.7 simulations undergo unusually high rate of mass transfer. At the end of initial driving phase, the rate of mass transfer was about orbit as compared to about yr for observed Algollike systems. This resulted in a corresponding shortening in the response of the individual spins and a relatively rapid spin evolution.
Figures 810 show equatorial cross sections of the density fields during the evolution of the Q0.7 simulations. We used an open source visualization tool VisIt (Childs et al., 2012) to render the cross sectional frames for the simulations. The distance scale used was the same for each simulation and each frame was 4 code lengthunits wide. The density field ranged from to in code units and is plotted on a logarithmic scale, with the intention of visualizing the low density accretion stream, while also depicting the coreenvelope structure. Each figure corresponds to the same during the binary simulation, thus representing the same stage of the dynamical evolution. Each figure is also rotated through a trivial angular phase about the axis of rotation to match the octHR simulation, so that the key features can be directly compared. Figure 8 corresponds to when a steady, direct impact mass transfer stream was established. The time is close to the end of the driving phase of the octLR simulation (see Figure 5). The behaviour of the Q0.7 systems agreed remarkably well after this point, with the exceptions discussed above. The resolutiondependent “numerical wind" can be noticed in the FlowER simulations. This figure also shows the formation of a puffed up, common envelope in the Octotiger simulations, due to the shockheated atmosphere. Note that a real outflow of material may occur in binaries from these outer Lagrangian points during mass transfer.
Figure 9 was used to determine for all simulations by matching as closely as possible with the octHR density. All four simulations showed the inspiral of the Q0.7 binary systems during the merger and the material lost through the outer Lagrange points, and , forming a double pinwheel pattern. The contrast between the core and the envelope in the FlowER simulations was noticeably less sharp as compared to the Octotiger ones, especially in case of the donor star. This is because in the FlowER simulations, the contact discontinuity at the coreenvelope boundary dissolved with time due to insufficient resolution.
Figure 10 shows the density in the equatorial plane right after the merger when disruption of the donor was complete according to the diagnostic plots. However, the remnants of the former core of the donor can be seen being wrapped around the accretor as a long red streak of inspiralling material above a density threshold of about 0.4. This former core material of the donor was poorly resolved in the FlowER simulations. The two cores did not collide in a headon manner. The donor was tidally broken apart and, due to its high specific angular momentum, its material dispersed as a thick disk around the final merged object. Figure 11 shows the meridional cross section of the density after the merger. The central, nonspherical merged object has the same basic structure in all four simulations. The remnant of the secondary core can be seen in cross section of the Octotiger simulations, whereas the core was much more diffuse and poorly resolved in the FlowER simulations. The major difference between the FlowER and the Octotiger simulations was the thickness of the disk in the latter case, due to high temperatures caused by the shock heating during the merger. The shockheated disk was almost twice as thick with Octotiger, as compared to the FlowER simulations. Notice that all the frames in Figures 811 qualitatively agree very well across all four Q0.7 simulations.
Figures 12 and 13 are identical to Figures 9 and 10 for Octotiger simulations, with the difference that the former pair of frames shows the density slices at the length scale of the entire computational domain of the octHR model. The extended, spiral outflow structure, about an order of magnitude larger than the radius of the accretor, was formed in both the Octotiger simulations. A simulation with the FlowER code will be computationally very expensive if we achieve similar size of the grid which resolves the circumstellar disk in its entirety, while also trying to maintain enough resolution of the core. Figure 14 shows a meridional cross section of the density field right after the merger, again at the length scale of the octHR simulation. The full thick disk is visible in this view in the case of the octHR simulation and most of the material remained in the equatorial plane. These frames demonstrate the much better capability of Octotiger to resolve the circumbinary disk and the extended outflow structure using AMR. Due to computational constraints, we terminated the simulations within about an orbit after the merger. Figure 15 shows the density distribution of octLR simulation, which was the only simulation continued, at three orbits after . Only about 1.6 % mass and 8.4% of the total angular momentum of the system was lost at this point.
We now discuss some salient differences among the simulations. Figures 16 and 17 show plots of the central density of the donor and accretor respectively as the Q0.7 binaries evolve. Consider the evolution of the central density of the FlowER simulations first. The FlowER simulations showed a decrease in the central density of both the stars over a dynamical time scale. As explained in section 3.1, this implies that the cores were not wellresolved. Table 1 lists the core resolution of the donor and the accretor for the initial Q0.7 models for reference. We expect this decrease to be a function of resolution, since increased resolution can better resolve the steep density gradients near the core. We can confirm this prediction by comparing the plots for cylLR and cylHR simulations. For both the donor and the accretor, as the resolution increased, the central density showed less deviation from the initial value of unity. Note that the mass transfer began just before the driving phase finished after the first four orbits in each simulation. In Figures 16 and 17 we can observe that as soon as the mass transfer started, the behaviour of the central density changed for both the donor and the accretor in the FlowER simulations. A significant change in the central density of the donor or the accretor did not occur in the case of Octotiger simulations before the end of the initial driving phase. As the mass transfer continued, the central density of the donor in the Octotiger simulations decreased, while that of the accretor started to increase. The difference between the two Octotiger simulations may again be a resolution effect. The decrease in the central density in the FlowER simulations appears to be a major numerical artefact. However, from our comparative analysis of the diagnostic plots as well as the density cross sections we can conclude that this did not significantly affect the dynamical evolution. The agreement among the results also implies that the simulations with both the hydrodynamic codes were concordant.
5.3 Hydrodynamic Evolution of Q0.6 Simulations
In this section we describe the evolution of Q0.6 set of simulations cylEE and octEE. As explained in section 4, the initial mass ratio lies below the critical value (), hence we expect a stable mass transfer. The benchmark for these simulations was to evolve the initial bipolytropic binary for a sufficiently large number of orbits, so that we can demonstrate that the binary systems remain separated, without showing signs of an imminent merger over a dynamical timescale. These simulations were not repeated at a higher resolution because we have already studied resolution dependance with Q0.7 set of simulations, and it is computationally expensive to run a large number of orbits at a high resolution. Table 6 describes the summary of Q0.6 simulations along with some quantities of interest, and Table 7 lists the movies of the temporal evolution of the equatorial cross sections of the density fields. The timedependent behaviour of key diagnostic parameters is plotted in Figures 1923. We found that Q0.6 binary models were stable on mass transfer with both the codes up to about 60 orbits, and at the end of the simulations the mass transfer rate stabilized to a constant value. Note that since cylEE simulation was conducted at an overall low resolution, only a general agreement during the evolution can be expected, while the ultimate outcome of achieving a stable mass transfer with both codes is paramount. After accounting for the resolution limits for the cylEE simulation as well as the difference in the stable mass transfer rates, the diagnostic quantities show a reasonable agreement among the two simulations.
Consider the the normalized and smoothed mass transfer rate of the donor for Q0.6 set of simulations, plotted in Figure 18. The smoothing was done with a moving boxcar average method, identical to that used for Q0.7 simulations. Each system was driven for 3 orbits at the rate of per orbit, and we ensured that a small but steady mass transfer stream was established before the driving was stopped. The mass transfer rate initially increased and then eventually plateaued for both the simulations. Note that this behaviour is qualitatively different from Q0.7 set of simulations. Thus as we anticipated, the change in the radius of the donor can keep up with that of its Roche lobe on mass transfer. However, the stable normalized mass transfer rate of cylEE simulation, , was about 22 times larger than that of octEE simulation, . We discuss the numerical reasons behind this difference later in this section. The stable mass transfer rate in real binaries depends on how the angular momentum is removed from the system, as well as the thermal or nuclear evolution of the components. The cylEE simulation can be considered to have progressed at a faster rate as compared to octEE, and the diagnostic parameters of the former reflected this rapid evolution.
Figure 19 shows the normalized separation between the two binary components of Q0.6 set of simulations. Oscillations with an amplitude of approximately can be noticed throughout the evolution due to small eccentricity of the orbits. The roughly linear fall in the separation for the first 3 orbits was caused by the artificial removal of angular momentum. A total of of the angular momentum was removed, resulting in the initial approximately decrease in the separation. The binary systems remained separated in both the simulations, which is consistent with our analytical expectations. In the case of cylEE simulation after the driving phase ended, the separation decreased over the next 20 orbits, and then the binary started separating, while the separation for octEE simulation shows little change after the initial driving phase. Near the end of the simulations, the cylEE binary also ended up at a larger separation as compared to octEE. These differences can be understood on the basis of the different mass transfer rates. Due to over an order of magnitude larger stable mass transfer rate, by the end of the simulations, the donor in cylEE transferred about 13% of its mass to its companion, while the donor in octEE transferred only 0.7% of its mass. This can also be inferred from Figure 20, which shows the plot of normalized mass ratio for the Q0.6 simulations. Since the mass ratio of the Q0.6 binaries was less than unity, the analytical equations predict that the separation should increase on mass transfer. The timescale over which the separation increases would be proportional to the rate of mass transfer. Thus the change in separation in case of octEE was much smaller than that in cylEE over the length of the simulations. The magnitude of the increase in separation would be qualitatively proportional to the amount of mass transferred; thus the binary in cylEE ended up at a greater separation then that in octEE. At the end of the simulations, the the donor in cylEE was less massive, as compared to the donor in octEE. Since the convective envelope with expands on mass loss, the donor in Q0.6 set of simulations would expand to a larger size. This expansion of the donor explains the fact that at , the binary in cylEE ended up at a separation wider than at the beginning of the mass transfer phase, while the donor was still capable of transferring mass.
Model ID  cylEE  octEE 

Equation of state  Ideal gas  Ideal gas with 
Dual Energy  
Formalism  
Total orbits ()  
Timesteps  M  M 
Number of cores  
Cost (CPUhr/orbit)  k  k 
Cost (calculations/orbit)  
Machine  QueenBee  QueenBee 
(LONI)  (LONI)  
Processors  2.33GHz 4Core  2.33GHz 4Core 
Xeon 64bit  Xeon 64bit 
Model ID  Movies 

cylEE  cylEE 
octEE  octEEzoomed 
A comparison of the behaviour of the angular momenta for the Q0.6 simulations requires us to account for small numerical diffusion effects in the cylEE simulation which are not present in the corresponding octEE simulation. Figure 21 shows normalized orbital angular momenta for the Q0.6 simulations. The initial driving can be noticed in the plot with a linear decline of about 4.5% over the first three orbits. Figures 22 and 23 show the plots for the normalized spin angular momenta of the donor and accretor respectively. Note that all three angular momenta plots for Q0.6 showed changes which are much smaller in magnitude as compared to Q0.7 set of simulations. After the end of the driving phase, decreased and the accretor was spun up in both the cases. Figure 22 shows that the donor angular momenta in both the cylEE and octEE simulation levelled out to a similar value. Thus the spin up of the accretor was at the cost of the orbital angular momentum through the mass transfer. The change in as well as in these simulations was qualitatively proportional to the average mass transfer rate of the respective simulation. However, a sign of discrepancy was seen in the early behaviour of in cylEE simulation. The binary systems in Q0.6 set of simulations were driven for the first three orbits, which removed AM from the system as a whole. Thus all three angular momenta – , and – should show an initial decrease before the mass transfer could be established. As soon as cylEE simulation started and before any mass transfer could begin, however, the spin angular momentum of the donor started increasing. This behaviour of in cylEE simulation was anomalous, since it was neither expected nor seen in the Octotiger counterpart. Note, however, that this anomalous increase in was very small, approximately two orders of magnitude below the angular momentum removed during the driving, and thus unlikely to affect the evolution significantly.
In previous numerical experiments with single polytropic components the FlowER code showed excellent behaviour (Motl et al., 2002; D’Souza et al., 2006; Even & Tohline, 2009). Thus, most likely the early behaviour of in cylEE before the onset of mass transfer is due to the numerical diffusion of the bipolytropic cores in the FlowER code. We have tested this hypothesis by running three supplementary simulations without driving: 1) FlowER with bipolytropic components; 2) FlowER with polytropic components; and 3) the same binary with bipolytropic components using Octotiger. Only the first simulation displayed the small anomalous increases in and accompanied by diffusion in the cores while the total angular momentum remained constant within code limits. The orbital angular momentum and the separation decreased slightly to compensate for the anomalous increase of the spin angular momenta. In the interest of conciseness and the primary focus of this paper, we do not present the results of these three simulations in detail. In spite of the drawbacks of the FlowER code, we stress that there was an excellent qualitative agreement between cylEE and octEE simulations. We conclude that the simulations with both the hydrodynamic codes were concordant in the case of a stable mass transfer as well over a dynamical timescale.
Figure 24 shows the equatorial cross sections of the density fields at , when the binary separation was approximately equal for both the Q0.6 simulations. The cross section for cylEE is rotated through a trivial angle for ease of comparison. The mass transfer rate in both the simulations was stabilized at this point and a more substantial mass transfer stream can be clearly observed in the case of cylEE. As explained in section 5.2, since Octotiger code allowed for shock heating in combination with the lack of radiative cooling, the components in octEE simulation showed a puffed up atmosphere. Note that the minimum value of the density threshold is chosen at a higher value than before in order to make the mass transfer stream more clearly noticeable in octEE simulation. At this threshold, the numerical wind in cylEE simulation and most of the puffed up atmosphere in octEE simulation (similar to Figure 8 (a) and (c)) is not visible. The differences between the two codes in terms of the treatment of energy equation did not affect the ultimate outcome of Q0.6 set of simulations.
We calculated the amount of mass and angular momentum lost by the binary system as a whole at the end of Q0.6 simulations. The mass is considered lost if it lied outside the point, and the angular momentum was also measured within this region after accounting for the initial driving phase. The binary in cylEE simulation lost a total of about of its mass over the duration of the simulation, however it gained approximately of the total angular momentum. These quantities were dominated by two competing numerical effects in FlowER simulations. First, there was numerical wind which carried away mass and angular momentum from the binary. Secondly, the implemented density floor reset added a finite amount of mass after each timestep. This mass added angular momentum to the system as well, since the calculations were performed in a rotational frame of reference. Thus cylEE simulation could not accurately track the lost mass and angular momentum over such a large number of orbits. The Octotiger simulation can give better estimates of the lost mass and angular momentum, since there was no density floor reset or numerical wind. The binary in octEE simulation lost only about of its total mass, which carried away of the total angular momentum. The primary reason for this mass loss through the point was the swelling of the atmosphere of the donor due to the treatment of shock heating. The lost material in Octotiger simulations can be better tracked with a larger size of the computational domain.
The numerical diffusion and associated artefacts in FlowER simulations could be ameliorated by increasing the resolution, however this is computationally prohibitively expensive. When the grid resolution is doubled for a unigrid simulation, the increase in the computational cost increases by a factor of because of the increase in the size of the data. However, since the general benchmark is evolving to a fixed point in time, the Courant et al. (1928) condition needs to be obeyed. Thus the timestep in Cartesian geometry decreases by a factor of two, which makes the computations 16 times more expensive. In cylindrical geometry the penalty of Courantlimited timestep is higher, making the computation up to 32 times more expensive in FlowER simulations. Consider Tables 4 and 6 where the cost of performing each simulation in terms of computational time is described. The cylHR simulation was about 21 times more expensive per orbit, as compared to the cylLR simulation, while the octHR simulation was only about 8 times more expensive than octLR per orbit. Thus in terms of scaling with increasing problem size, FlowER scaled worse than a factor of 16, while Octotiger fared much better. Octotiger can also utilize a much larger number of cores because of its use of the HPX runtime system. We would like to emphasize that Octotiger code as well as the HPX runtime system both are under active development and are being continuously improved. We used a more optimized version of Octotiger for octEE simulation, which incorporated adaptive angular velocity during evolution. The binary axis was kept fixed along the Xaxis of the computational domain during evolution, thus reducing the need of constant recalculation of the AMR grid (see movie octEEzoomed listed in Table 7). Compared to octLR simulation, which was conducted with a similar effective resolution, octEE simulation was over 8 times more efficient. Thus the code Octotiger is much better suited for conducting numerical simulations of close binary interactions as compared to FlowER.
6 Summary and Discussions
In this paper we established the reliability of the set of gridbased numerical tools – the BSCF method and Octotiger – for simulating the evolution of bipolytropic binary systems on a dynamical scale, orbits. We achieved this through a comparative analysis with an independent and well tested hydrodynamic code FlowER. With the BSCF method, binary systems with a general coreenvelope structure for both the stars can be obtained, which were used as quiet initial conditions for the simulations. The initial models in the suite of simulations closely matched each other, so that the results can be used as benchmarks. The codes were compared across binary systems with two mass ratios, and , the former at two resolutions. The code FlowER evolved a binary adiabatically on a cylindrical mesh with an ideal gas equation of state. Octotiger used a dual energy formalism (which allowed for shock heating via solving total energy equation explicitly) with ideal gas equation of state and evolved the binary on a Cartesian AMR grid. The binary components were initially driven together until a steady accretion stream was formed, and the subsequent dynamical behaviour was analysed without any special assumptions.
We summarize the results obtained from the suite of simulations as follows–

Q0.7 set of simulations resulted in an unstable mass transfer and an eventual merger of the systems due to tidal disruption of the donor. This was consistent with our theoretical expectations, since the initial mass ratio was greater than the approximate critical value for such system.

Q0.6 set of simulations resulted in a stable mass transfer, which was observed for about 60 orbits. This behaviour was consistent with our expectations as well, since the initial mass ratio was less than the critical value.

In the case of Q0.7 set of simulations, a discrepancy between the two codes was observed with respect to the behaviour of mass transfer rate and angular momenta just before the merger. This can be accounted for by the higher mass transfer rate in Octotiger during the merger due to shock heating, as compared to the FlowER simulations.

The FlowER simulations showed an anomalous decrease in the central densities for both the donor and the accretor, which was caused by numerical diffusion due to relatively poor resolution of the cores. Neither the overall behaviour nor the ultimate fate of the binary systems were affected by this numerical artefact. In the case of Q0.6 set of simulations, indirect effects on the mass transfer rate and angular momenta can be inferred. The effects can be attributed to the change in the moment of inertia of individual components due to the internal redistribution of the mass.

The intermediate stages, in terms of diagnostic plots and density cross sections, closely matched for the Q0.7 simulations at two resolutions. The stable behaviour of Q0.6 simulations also reasonably matched for the two simulations. This established the concordance of simulations for both the codes, while also proving that the codes can reliably handle stable as well as unstable mass transfer scenario.
Further investigations using Octotiger and BSCF method could shed light on the complex evolutionary history of contact binary systems and, especially on dynamical aspects that are not wellunderstood, such as largescale circulations, internal structure and stability criteria. In particular, in paper III these tools will be used to investigate the merger that led to V1309 Scorpii outburst. The progenitor contact binary had a rather extreme mass ratio of 0.11, and it is proposed that the Darwin (1880) instability of the progenitor contact binary resulted in the observed merger (Stȩpień, 2011). These tools can also be useful in studying the interesting case of KIC 9832227, which may be the first prediction of a luminous red nova event Molnar et al. (2017). With the ability of Octotiger to conserve key physical quantities up to machine precision, as well as the increase in efficiency due to the HPX runtime system, we can conduct detailed, high resolution simulations that were previously not feasible. In conclusion, this study is a step forward in the direction of investigating dynamical interactions of mass transferring binary systems. In the future, the continuation of this project will allow us to numerically study and better understand the evolution of close and contact binary systems, as well as related astrophysical phenomena.
7 Acknowledgements
We thank the referee, Prof. Christopher Tout, for his insightful comments and his suggestion to include Q0.6 set of simulations. We wish to acknowledge the support from the National Science Foundation through CREATIV grant AST1240655. The numerical work was carried out using the computational resources of the Louisiana Optical Network Initiative (LONI) and Louisiana State University’s High Performance Computing (LSU HPC).
References
 Bryan et al. (1995) Bryan G. L., Norman M. L., Stone J. M., Cen R., Ostriker J. P., 1995, Computer Physics Communications, 89, 149
 Childs et al. (2012) Childs H., et al., 2012, in , High Performance Visualization–Enabling ExtremeScale Scientific Insight. pp 357–372
 Courant et al. (1928) Courant R., Friedrichs K., Lewy H., 1928, Mathematische Annalen, 100, 32
 D’Souza et al. (2006) D’Souza M. C. R., Motl P. M., Tohline J. E., Frank J., 2006, ApJ, 643, 381
 Darwin (1880) Darwin G. H., 1880, Astronomische Nachrichten, 96, 217
 Dehnen (2000) Dehnen W., 2000, ApJ, 536, L39
 Després & Labourasse (2014) Després B., Labourasse E., 2014, Technical report, Angular Momentum preserving cellcentered Lagrangian and Eulerian schemes on arbitrary grids, Sep 2014, hal01065105
 Even & Tohline (2009) Even W., Tohline J. E., 2009, ApJS, 184, 248
 Frank et al. (2002) Frank J., King A., Raine D., 2002, AccretionPower in Astrophysics, 3rd edn. Cambridge
 Hachisu (1986a) Hachisu I., 1986a, ApJS, 61, 479
 Hachisu (1986b) Hachisu I., 1986b, ApJS, 62, 461
 Henrich & Chandrasekhar (1941) Henrich L. R., Chandrasekhar S., 1941, ApJ, 94, 525
 Kadam et al. (2016) Kadam K., Motl P. M., Frank J., Clayton G. C., Marcello D. C., 2016, MNRAS, 462, 2237
 Kaiser et al. (2015) Kaiser H., Heller T., Bourgeois D., Fey D., 2015, Proceedings of the First International Workshop on Extreme Scale Programming Models and Middleware, ESPM ‘15, (New York, NY, USA), pp 29–37
 Kawaler (1988) Kawaler S. D., 1988, ApJ, 333, 236
 Kuiper (1941) Kuiper G. P., 1941, ApJ, 93, 133
 Li et al. (2007) Li L., Zhang F., Han Z., Jiang D., 2007, ApJ, 662, 596
 Lucy (1968) Lucy L. B., 1968, ApJ, 151, 1123
 Marcello (2017) Marcello D. C., 2017, ApJS, p. under review
 Marcello & Tohline (2012) Marcello D. C., Tohline J. E., 2012, ApJS, 199, 35
 Marcello et al. (2017) Marcello D. C., Kadam K., Clayton G. C., Frank J., Kaiser H., Motl P. M., 2017, Proceedings of Science, Under review
 Molnar et al. (2017) Molnar L. A., et al., 2017, in American Astronomical Society Meeting Abstracts. p. 417.04
 Motl et al. (2002) Motl P. M., Tohline J. E., Frank J., 2002, ApJS, 138, 121
 Motl et al. (2017) Motl P. M., Frank J., Staff J., Clayton G. C., Fryer C. L., Even W., Diehl S., Tohline J. E., 2017, preprint, (arXiv:1702.03562)
 Nandez et al. (2014) Nandez J. L. A., Ivanova N., Lombardi Jr. J. C., 2014, ApJ, 786, 39
 Ostriker & Mark (1968) Ostriker J. P., Mark J. W.K., 1968, ApJ, 151, 1075
 Press et al. (1986) Press W., Teukolsky S., Vetterling W., Flannery B., 1986, Numerical Recipes in Fortran, 2nd edn. Cambridge University Press
 Stȩpień (2009) Stȩpień K., 2009, MNRAS, 397, 857
 Stȩpień (2011) Stȩpień K., 2011, A&A, 531, A18
 Stone et al. (1992) Stone J. M., Mihalas D., Norman M. L., 1992, ApJS, 80, 819
 Tout et al. (1997) Tout C. A., Aarseth S. J., Pols O. R., Eggleton P. P., 1997, MNRAS, 291, 732
 Trimble (1983) Trimble V., 1983, Nature, 303, 137
 Webbink (1985) Webbink R., 1985, Interacting Binary Stars. Cambridge University Press, p. 39
 Zahn (1977) Zahn J.P., 1977, A&A, 57, 383