Radio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 28 September 2018 (MN L A TEX style file v2.2)
Radio Flares of Compact Binary Mergers: the Effect ofNon-Trivial Outflow Geometry
Ben Margalit , & Tsvi Piran Racah Institute of Physics, Hebrew University of Jerusalem, Israel Physics Department, Columbia University, 538 West 120th Street New York, NY 10027
ABSTRACT
The next generation gravitational waves (GW) detectors are most sensitive toGW emitted by compact (neutron star/black hole) binary mergers. If one of those isa neutron star the merger will also emit electromagnetic radiation via three possiblechannels: Gamma-ray bursts and their (possibly orphan) afterglows (Eichler et al.1989), Li-Paczynski Macronovae (Li & Paczy´nski 1998) and radio flares (Nakar & Pi-ran 2011). This accompanying electromagnetic radiation is vitally important in con-firming the GW detections (Kochanek & Piran 1993). It could also reveal a wealthof information regarding the merger and will open a window towards multi-messengerastronomy. Identifying and characterizing these counterparts is therefore of utmostimportance. In this work we explore late time radio flares emitted by the dynami-cally ejected outflows. We build upon previous work and consider the effect of theoutflow’s non-trivial geometry. Using an approximate method we estimate the radiolight-curves for several ejected matter distributions obtained in numerical simulations.Our method provides an upper limit to the effect of non-sphericity. Together with thespherical estimates the resulting light curves bound the actual signal. We find thatwhile non-spherical geometries can in principle lead to an enhanced emission, in mostcases they result in an increase in the timescale compared with a corresponding spher-ical configuration. This would weaken somewhat these signals and might decrease thedetection prospects.
Compact binary mergers (hereafter called simply mergers)have recently been the focus of extensive research in the ef-forts to detect gravitational waves (GW) emitted in thesemerger events (see Faber & Rasio 2012, for a recent re-view). The decreasing orbital period of the binary pulsars(e.g. Hulse & Taylor 1975; Taylor & Weisberg 1989; Krameret al. 2006; Weisberg et al. 2010) provides so far the bestevidence for the existence of GW. Still, as yet, GW have notbeen directly detected. While gravitational radiation shouldbe produced in various different scenarios, current GW de-tectors are most sensitive to radiation emitted in mergers,and hence the quest of detecting GW is closely linked withan understanding of these events.A detection of electromagnetic counterparts to thesemerger events is of utmost importance (Kochanek & Piran1993). Such electromagnetic signals, if observed, could helppave the way to confirm the first GW detection and wouldreveal a wealth of information about the merger process (e.g.Nakar & Piran 2011). Since mergers are considered as likelysources of short Gamma-Ray Bursts (GRBs) (Eichler et al.1989), these have often been considered as possible electro-magnetic counterparts to mergers. However GRBs are highlybeamed and are only observed when their relativistic jet points towards us. This decreases significantly the chancesof a coincident GW-GRB detection. Although off-axis GRBscannot be observed in gamma-rays, they are followed by alate time isotropic radio signal, namely the orphan afterglow,which may possibly be observed. These orphan afterglowsare more promising electromagnetic counterparts in this re-spect, yet their signals are expected to be weaker than theradio flares from dynamical ejecta from the mergers thatwe discuss here (Nakar & Piran 2011; Hotokezaka & Piran2015).In addition to GRBs, two main electromagnetic coun-terparts have been suggested. These are the “Macronova”(or alternatively “Kilonova”, as preferred by a few authors)which results from the radioactive decay of freshly synthe-sized r-process material and thus resembles a supernova (Li& Paczy´nski 1998; Kulkarni 2005; Metzger et al. 2010; Piranet al. 2013; Kasen et al. 2013; Barnes & Kasen 2013; Tanaka& Hotokezaka 2013; Grossman et al. 2014), and radio flareswhich arise from the interactions between mergers’ outflowsand the surrounding ISM and are analogous to Radio Su-pernovae and GRB afterglows (Nakar & Piran 2011; Piranet al. 2013). The latter are produced via synchrotron radia-tion of non-thermal electrons accelerated at the outflow-ISMshock front. Owing to the typical velocities and masses in-volved in this ejecta ( ∼ . c and ∼ − M (cid:12) respectively), c (cid:13) a r X i v : . [ a s t r o - ph . H E ] M a r B. Margalit & T. Piran this afterglow component should peak in the radio band ona timescale of a year after the merger. So far, radio flareshave been considered only for spherically symmetric models(Nakar & Piran 2011; Piran et al. 2013). Our goal here isto characterize the effect of departure from spherical sym-metry on the resulting radio signals. This is an importantrefinement of previous results, since merger simulations sug-gests that the outflows produced by these mergers are sig-nificantly ashperical. For a related discussion on the effect ofnon-sphericity on the Macronova signal, see Grossman et al.(2014).In this work we treat only the sub/mildly relativisticoutflow components ejected dynamically by mergers, andfocus on the late time radio signal produced by these com-ponents. We neglect ultra-relativistic jet components thatmay produce GRBs or other short-lived signals, as well asnon-dynamically ejected outflows such as those that mayarise from strong neutrino driven winds (see Hotokezaka &Piran 2015, for a review of these components and their ex-pected radio flares). Our method follows Piran et al. (2013),hereafter PNR13, which assumed a spherically symmetricoutflow, and address the impact of nontrivial geometry onthe resulting radio signals and it’s repercussions for the de-tectability of such signals using current and future radiotelescopes.The use of an approximate method is justified in lightof the large astrophysical uncertainties regarding these sys-tems. In particular, the circum-binary density is a key un-known that will depend on the binary’s location within it’shost galaxy. This density is crucial in reproducing the exactdynamics of the system, and variations in this density willaffect the timescale and luminosity of the radio light-curvesdramatically. Additionally, the exact outflows expected frommergers are still rather ill constrained. The dynamicallyejected outflow depends on the exact merger scenario as wellas on unknown microphysics and in particular on the neu-tron star’s equation of state. We use the dynamical outflowsfound by Rosswog et al. (2014), yet other groups using dif-ferent numerical methods and equations of state find slightlydifferent outflow characteristics (e.g Hotokezaka et al. 2013).Furthermore, other mass outflows arise in mergers besidesthe dynamical ejecta considered here (see e.g. Hotokezaka& Piran 2015). In view of these many uncertainties in therelevant physical parameters, we choose to use an approxi-mate analysis instead of a full fledged numerical simulation.As we discuss later, our approximate method maximizes theimpact of the non-sphericity and as such it provides an up-per limit to the full effect. Our motivation is to estimatethe detectability of these late time radio signals in com-parison with other suggested electromagnetic counterpartssuch as GRB orphan afterglows. We therefore focus here onevaluating only the peak timescale and peak flux of thesesignals (these will suffice for estimating the detectability),and we do not presume to reproduce the exact radio flarelight-curves. We expect that the real peak time scale andpeak flux are somewhere between those estimtated here andthose estimated assuming spherical symmetry.This paper is organized as follows - we begin in § § §
4. In § § §
5. Finally, we conclude in § Even though we are interested in non-spherical outflows inthis work, our method is based on results for spherical sys-tems (see § E and an initialvelocity v decelerates on a time scale, t dec : t dec = R dec v = (cid:18) E πnm p (cid:19) / v − / , (1)where n is the circum-binary ISM number density, and m p the proton mass. The homologously expanding outflow ischaracterized by the energy distribution E ( (cid:62) v ). At earliertimes the ejecta expands freely, whereas at later times theoutflow follows the Sedov-Taylor self-similar solution (Se-dov 1946; Taylor 1950). At times between t dec of the fastestmoving shell in the ejecta, and t dec of the slowest movingshell we describe the approximate dynamics of the contactdiscontinuity (between the shocked ejecta and shocked ISM)by solving:4 π nm p R v = M ( R ) v = E ( (cid:62) v ) , (2)to obtain the radius R ( t ) and velocity v ( t ) of the contactdiscontinuity. The leftmost equality of Eq. 2 assumes spher-ical geometry and a constant ISM density. It is correct upto a factor of order unity, which will depend on the exactstructure profile of the ejecta. The dynamics implied by thisequation is governed solely by the outflow’s angular energydensity, dE/d Ω, which is just E/ π for the spherical case. Inthe following sections we generalize this for a non-sphericaloutflow confined to a solid angle Ω, and we use this term( dE/d Ω) here for convenience.The collisionless shock wave accelerates electrons to apower-law distribution dN/dγ ∼ γ − p . This holds for γ > γ m ,where γ m is the minimal Lorentz factor of the accelerated c (cid:13) , 000–000 adio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry electrons, and 2 (cid:54) p (cid:54) p ≈ . − . p ≈ . − (cid:15) e and (cid:15) B of the total energy density at the shock front are given to theelectron and the magnetic energy respectively. Synchrotronradiation is emitted by these accelerated electrons. FollowingPNR13, the typical synchrotron frequency ν m and typicalflux F m can be written as: ν m ( t ) ≈ . (cid:16) (cid:15) B . (cid:17) / (cid:16) (cid:15) e . (cid:17) n / β ( t ) , (3) F m ( t ) ≈ . (cid:18) R ( t )10 cm (cid:19) n / (cid:16) (cid:15) B . (cid:17) / (cid:18) d cm (cid:19) − β ( t ) , (4)where d is the distance between the source and the observer(for simplicity we neglecting any cosmological effects).The synchrotron emission is self absorbed below ν a , thesynchrotron self-absorption frequency. For simplicity we fo-cus on a single standard scenario in which ν a , ν m < ν obs ,where ν obs is the observed frequency. This generally holdsfor ν obs = 1 . ν m and ν a will normally besignificantly smaller than 1 GHz (see Eq. 3 and equation 5of PNR13 plugging in typical velocities of the merger’s dy-namic ejecta, β ≈ . − . t contributeto the signal, the light-curve scales as F ν obs ( t ) = F m ( t ) (cid:18) ν obs ν m ( t ) (cid:19) − p − ∝ E ( (cid:62) v ( t )) v ( t ) − p − , (5)where the last proportionality is given by relating R ( t ) tothe ejecta’s energy distribution via Eq. 2. It describes thefact that only the energy of shells initially moving fasterthan v ( t ) contribute. The complete light-curve is thereforegiven by solving Eq. 2 for the shock radius and velocity andsubstituting these into Eq. 5.In the following we discuss some properties of light-curves arising from this model. In particular we examinethe peak flux and peak timescale for these signals whichare important for the rest of this work, but also quantifythe shape of the light-curve by defining the peak width. Webegin by examining an idealized toy model for the outflowenergy distribution E ( (cid:62) v ), namely the single-velocity shell,and continue by showing that any arbitrary energy distribu-tion can be approximated by this toy model for the purposesof evaluating the peak light-curve properties.A spherical “single-velocity shell” is a spherically sym-metric outflow with energy distribution dE/dv = E δ ( v ).While this distribution does not well describe the merger’soutflows, we will show later on that results for single-velocityshells can be very useful in approximating the peak light-curve properties (see § t peak = t dec ∝ (cid:18) dMd Ω (cid:19) / (cid:18) dEd Ω (cid:19) − / . (6) The corresponding peak flux can similarly be expressed as: F peak ∝ E (cid:18) dEdM (cid:19) − p − . (7)The light-curve of a single-velocity shell scales as F ν obs ( t ) ∝ (cid:40) t ; t (cid:54) t dec t − p − ; t dec < t . (8)Using this result we find the times at which the flux reacheshalf it’s peak value ( t ( − )1 / and t (+)1 / ): t ( − )1 / = 2 − / t peak ≈ . t peak , (9) t (+)1 / = 2 / (15 p − t peak ≈ . t peak , (10)where in the last equality in Eq. 10 we have substituted p = 2 .
5. With these expressions, we define the peak width,∆ t , as the time period during which the flux exceeds halfit’s peak value. For a spherical single-velocity shell, we find∆ t = t (+)1 / − t ( − )1 / ≈ . t peak . We use this measure for thepeak width to quantify the extent with which two distinctlight-curves overlap.An important note regarding ∆ t given above, is thatthe light-curve scaling laws on which it is based are self-similar solutions towards which (or from which) the systemsteadily evolves, and are only applicable for a single veloc-ity shell. For more ‘realistic’ shells which have an arbitrary E ( (cid:62) v ), the light-curve around the peak flux differs signifi-cantly from this approximation, changing smoothly betweenthese scaling laws which only hold at very early times or latetimes. This effect significantly increases the peak width fornon single-velocity shells. Still, ∆ t is useful in that it maybe used as a (rather loose) lower bound on the peak widthof an outflow.We turn now to examine an arbitrary energy distribu-tion E ( (cid:62) v ) outflow. The velocity v peak at which the light-curve peaks for such an outflow is easily found by differen-tiating Eq. 5 with respect to v . This peak velocity satisfies: d log E ( (cid:62) v ) d log v (cid:12)(cid:12)(cid:12)(cid:12) v peak = − p − . (11)This result is straightforward once Eq. 5 is established -at earlier times than t peak the shock velocity is larger than v peak , yet the amount of energy available E ( (cid:62) v ) decreasesfaster than v − (5 p − / , hence the overall flux is smaller. Atlater times, the shock velocity is smaller than v peak and theenergy E ( (cid:62) v ) increases, yet this increase is slower than v − (5 p − / , hence the decrease in velocity dominates andthe total flux is once again smaller. Figure 2 depicts a rep-resentative energy distribution for one of the merger’s dy-namic outflow configurations discussed below (ns14ns14),and shows that the velocity v peak at which the light-curvepeaks corresponds to the value expected from Eq. 11.We find that the peak flux of a radially structured out-flow can therefore be characterized using only two parame-ters v peak , and E ( (cid:62) v peak ). It does not depend on the entiredistribution E ( (cid:62) v ). Unfortunately we cannot characterizethe time at which the light-curve peaks using only these twovariables, for that the detailed energy distribution must be c (cid:13)000
5. With these expressions, we define the peak width,∆ t , as the time period during which the flux exceeds halfit’s peak value. For a spherical single-velocity shell, we find∆ t = t (+)1 / − t ( − )1 / ≈ . t peak . We use this measure for thepeak width to quantify the extent with which two distinctlight-curves overlap.An important note regarding ∆ t given above, is thatthe light-curve scaling laws on which it is based are self-similar solutions towards which (or from which) the systemsteadily evolves, and are only applicable for a single veloc-ity shell. For more ‘realistic’ shells which have an arbitrary E ( (cid:62) v ), the light-curve around the peak flux differs signifi-cantly from this approximation, changing smoothly betweenthese scaling laws which only hold at very early times or latetimes. This effect significantly increases the peak width fornon single-velocity shells. Still, ∆ t is useful in that it maybe used as a (rather loose) lower bound on the peak widthof an outflow.We turn now to examine an arbitrary energy distribu-tion E ( (cid:62) v ) outflow. The velocity v peak at which the light-curve peaks for such an outflow is easily found by differen-tiating Eq. 5 with respect to v . This peak velocity satisfies: d log E ( (cid:62) v ) d log v (cid:12)(cid:12)(cid:12)(cid:12) v peak = − p − . (11)This result is straightforward once Eq. 5 is established -at earlier times than t peak the shock velocity is larger than v peak , yet the amount of energy available E ( (cid:62) v ) decreasesfaster than v − (5 p − / , hence the overall flux is smaller. Atlater times, the shock velocity is smaller than v peak and theenergy E ( (cid:62) v ) increases, yet this increase is slower than v − (5 p − / , hence the decrease in velocity dominates andthe total flux is once again smaller. Figure 2 depicts a rep-resentative energy distribution for one of the merger’s dy-namic outflow configurations discussed below (ns14ns14),and shows that the velocity v peak at which the light-curvepeaks corresponds to the value expected from Eq. 11.We find that the peak flux of a radially structured out-flow can therefore be characterized using only two parame-ters v peak , and E ( (cid:62) v peak ). It does not depend on the entiredistribution E ( (cid:62) v ). Unfortunately we cannot characterizethe time at which the light-curve peaks using only these twovariables, for that the detailed energy distribution must be c (cid:13)000 , 000–000 B. Margalit & T. Piran −1 v [ c ] E ( ≥ v ) [ e r g ] v peak , E ( ≥ v peak ) ∼ v − p − Figure 1.
Cumulative energy distribution E ( (cid:62) v ). The red crosssignifies the velocity at which the calculated flux peaks. Thedashed red line shows that the logarithmic slope at this peakvelocity is − (5 p − /
2, as expected by Eq. 11. The cumulative dis-tribution decreases faster than this slope above v peak , but not assteeply as a step function (which would describe a single-velocityshell). accounted for. Still, we may find t peak up to a factor of orderunity using only these two parameters.We can change the time at which the v peak velocityshell decelerates while keeping E ( (cid:62) v peak ) fixed by allocat-ing more/less energy into higher velocity shells. If a largefraction of E ( (cid:62) v peak ) is put into very fast shells (whichdecelerate quickly), only a small fraction is left in the v peak shell, and by Eq. 1, this velocity shell decelerates earlier. Anupper bound on this deceleration time is therefore given bythe limiting case in which all the matter above this velocity, M ( (cid:62) v peak ), moves with v peak (as a single-velocity shell),and the distribution E ( (cid:62) v ) decreases infinitely fast above v peak , as a step function. In reality the distribution will de-crease gradually, some of the energy will go into faster shellswith short deceleration times, and the shell moving with v peak will be affected earlier. Thus this approximation willbe altered by a numerical factor < t peak (cid:46) t dec = (cid:18) E ( (cid:62) v peak )4 πnm p (cid:19) / v − / peak . (12)A lower bound on this numerical factor (and thus on thetimescale) may be obtained by considering the case where E ( (cid:62) v ) falls off as slowly as possible above v peak . Eq. 11 inconjunction with the fact that both E ( (cid:62) v ) and d log E/d log v decrease monotonically, implies that E ( (cid:62) v ) must falloffsteeper than v − (5 p − / . Hence we may find the peak timein the limiting case where E ( (cid:62) v ) decays as a power lawabove v peak with a slope − (5 p − /
2. By solving Eq. 2 incase of a general power law profile E ( (cid:62) v ) ∝ v − ( k − , wecan write the expression for the shock velocity (as given byequation 16 of PNR13). We are interested in the case of k = (5 p + 3) /
2. The peak time is thus bound from belowby the time at which v ( t ) for this power law profile equals v peak . Solving this condition by substituting the appropriate k and replacing v min by v peak in equation 16 of PNR13, wefind that t peak > k − k (cid:18) E ( (cid:62) v peak )4 πnm p (cid:19) / v − / peak = 5 p − p + 3 t dec . (13)Since both the upper and lower bounds scale equally(as t dec ), we have pinpointed the peak time to within a fac-tor of (5 p − / (5 p + 3) ≈ .
6. We can therefore reasonablycharacterize a complicated outflow distribution by a char-acteristic velocity v peak (which is defined by Eq. 11), and acharacteristic energy E ( (cid:62) v peak ), and reduce realistic radi-ally structured ejecta into roughly equivalent single-velocityshell building blocks. These ‘equivalent’ single-velocity shells(characterized by velocity v peak and energy E ( (cid:62) v peak )) willnot reproduce identical light-curves, but will reproduce thecorrect peak flux, and approximately the correct peak time.Since these are the focus of this present work, such a reduc-tion will suffice when necessary. So far we have discussed spherical outflows. Turning nowto non-spherical outflows we use a Piecewise Spherical Ap-proximation in which we do not calculate the full three-dimensional (3D) hydrodynamics of the system in detail.Instead, we decompose the 3D system into many effectiveone-dimensional (1D) radial problems which are solved sep-arately. We then combine the one dimensional results to es-timate the overall outflow evolution and the correspondingsynchrotron emission.The method follows schematically the following stages:(i) We divide the system into solid angle elements and ex-tract the cumulative energy distribution E ( (cid:62) v ) / Ω for eachelement. This distribution determines the dynamics of eachelement. (ii) We calculate the spherical equivalent dynam-ics of each solid angle using E ( (cid:62) v ) / Ω according to Eq. 2.(iii) We obtain the resulting synchrotron flux for each ele-ment. (iv) We estimate the system’s overall light-curve bysumming the contribution of all elements.We call this approximation “piecewise-spherical” in thesense that we decompose the aspherical outflow into severalsolid angle ‘pieces’ (denoted by index i ) and treat each solidangle component i as if it is a part of a spherically sym-metric outflow with an energy distribution dE ( (cid:62) v ) /d Ω = E i ( (cid:62) v ) / Ω i . Using the results of § c (cid:13) , 000–000 adio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry Ω Ω , Ω , … , Ω 𝑁 𝐹 𝜈1 (𝑡) 𝑡 𝐸 (≥ 𝑣)/Ω 𝑣 𝑟 𝑣(𝑡)𝑅(𝑡) Ω 𝑟 𝐸 (≥ 𝑣)/Ω 𝑣 Ω 𝑣(𝑡) 𝑅(𝑡) 𝐹 𝜈2 (𝑡) 𝑡 𝐹 𝜈 (𝑡) 𝑡 (𝑖𝑣)(𝑖) ... 𝑁 .. . 𝑁 (𝑖𝑖) (𝑖𝑖𝑖) Figure 2.
A cartoon sketch illustrating our method. In step (i) we divide the spherical system into N solid angle bins Ω , ..., Ω N , andextract the cumulative energy distribution E ( (cid:62) v ) from each. In step (ii) we use each bin’s energy density distribution to calculate thetemporal evolution of the spherical equivalent system via Eq. 2. Using these dynamics, we calculate in (iii) the flux emitted by each solidangle. Finally we arrive at the total light-curve by summing up the N contributions in (iv). We turn now to discuss the important issue of the limitationsof this approximate method. We begin by stating clearly theunderlying assumptions of the piecewise-spherical approxi-mation:(1) Each solid angle element evolves according to it’sspherical equivalent system, and therefore the material ineach solid angle expands only in the radial direction.(2) Each solid angle element evolves separately, and is de-coupled from any interactions with other elements. In partic-ular, neighboring solid angles elements do not transfer massor energy from one to another and they don’t affect eachother’s dynamics.Both assumptions will hold well if the angular variationsin energy density and velocity are small, so that the systemis nearly spherically symmetric. To illustrate this considerthe limiting case of a jet confined within a solid angle Ω. Inthis case the energy density and velocity change abruptlyand dramatically on the boundary and assumption (1) willbreak down because the jet will expand azimuthally. Thisjet expansion occurs as the shocked region is heated andthe hot shocked material expands into the surrounding coldISM at roughly the speed of sound. If, on the other hand,we consider an angular region Ω which is a part of a nearlyspherical system, such a sideways expansion is hindered.Assumptions (1) and (2) are closely related to eachother, as any azimuthal expansion of one bin essentially in-volves mass transfer to the adjacent bins. This mass trans-fer would in turn affect both bins’ evolution and assump-tion (2) will break down. However, we expect that even incase the system is far from spherical, both these assump-tions will hold reasonably well until the time of peak flux .This is expected because the peak flux time in this regime( ν m , ν a < ν obs ) is the deceleration time, t dec , of the out-flow. Because t dec is by definition the time at which the surrounding ISM begins affecting the ejecta’s dynamics, itis reasonable to expect that the ISM-ejecta interactions willnot result in any significant sideways expansion of the ejecta before this point.Faster moving ejecta components will expand az-imuthally before their slower surroundings, and energywill transfer from high to low velocity regions, effectivelysmoothing out the anisotropy in the initial distribution.Thus, our approximation overestimates the effect of as-phericity, and we expect the piecewise-spherical method toprovide an upper limit, such that the true light-curve liessomewhere in between our present results and those calcu-lated using the spherical approximation. We begin by using the piecewise spherical approach de-scribed in § M tot andenergy E tot uniformly distributed over a solid angle of 4 π can be morphed into an aspherical outflow by redistributingthe mass/energy unevenly between different solid angle ele-ments. Under the piecewise spherical approach, each of thesesolid angle elements can then be evolved independently, al-lowing an estimate of the overall resulting light-curve. Toobtain an analytic intuition to the behavior of non-sphericalsystems we consider a simple model in which we approxi-mate each solid angle element as a single-velocity shell, as c (cid:13)000
A cartoon sketch illustrating our method. In step (i) we divide the spherical system into N solid angle bins Ω , ..., Ω N , andextract the cumulative energy distribution E ( (cid:62) v ) from each. In step (ii) we use each bin’s energy density distribution to calculate thetemporal evolution of the spherical equivalent system via Eq. 2. Using these dynamics, we calculate in (iii) the flux emitted by each solidangle. Finally we arrive at the total light-curve by summing up the N contributions in (iv). We turn now to discuss the important issue of the limitationsof this approximate method. We begin by stating clearly theunderlying assumptions of the piecewise-spherical approxi-mation:(1) Each solid angle element evolves according to it’sspherical equivalent system, and therefore the material ineach solid angle expands only in the radial direction.(2) Each solid angle element evolves separately, and is de-coupled from any interactions with other elements. In partic-ular, neighboring solid angles elements do not transfer massor energy from one to another and they don’t affect eachother’s dynamics.Both assumptions will hold well if the angular variationsin energy density and velocity are small, so that the systemis nearly spherically symmetric. To illustrate this considerthe limiting case of a jet confined within a solid angle Ω. Inthis case the energy density and velocity change abruptlyand dramatically on the boundary and assumption (1) willbreak down because the jet will expand azimuthally. Thisjet expansion occurs as the shocked region is heated andthe hot shocked material expands into the surrounding coldISM at roughly the speed of sound. If, on the other hand,we consider an angular region Ω which is a part of a nearlyspherical system, such a sideways expansion is hindered.Assumptions (1) and (2) are closely related to eachother, as any azimuthal expansion of one bin essentially in-volves mass transfer to the adjacent bins. This mass trans-fer would in turn affect both bins’ evolution and assump-tion (2) will break down. However, we expect that even incase the system is far from spherical, both these assump-tions will hold reasonably well until the time of peak flux .This is expected because the peak flux time in this regime( ν m , ν a < ν obs ) is the deceleration time, t dec , of the out-flow. Because t dec is by definition the time at which the surrounding ISM begins affecting the ejecta’s dynamics, itis reasonable to expect that the ISM-ejecta interactions willnot result in any significant sideways expansion of the ejecta before this point.Faster moving ejecta components will expand az-imuthally before their slower surroundings, and energywill transfer from high to low velocity regions, effectivelysmoothing out the anisotropy in the initial distribution.Thus, our approximation overestimates the effect of as-phericity, and we expect the piecewise-spherical method toprovide an upper limit, such that the true light-curve liessomewhere in between our present results and those calcu-lated using the spherical approximation. We begin by using the piecewise spherical approach de-scribed in § M tot andenergy E tot uniformly distributed over a solid angle of 4 π can be morphed into an aspherical outflow by redistributingthe mass/energy unevenly between different solid angle ele-ments. Under the piecewise spherical approach, each of thesesolid angle elements can then be evolved independently, al-lowing an estimate of the overall resulting light-curve. Toobtain an analytic intuition to the behavior of non-sphericalsystems we consider a simple model in which we approxi-mate each solid angle element as a single-velocity shell, as c (cid:13)000 , 000–000 B. Margalit & T. Piran opposed to treating the full energy distributions E ( (cid:62) v ) foreach solid angle element.Consider two simple examples: (a) The entire mass andthe entire energy of the system are contained within half thesky (2 π rad ). (b) The mass is distributed spherically but theentire energy is within 2 π rad . In both cases, half the spherecontains no energy and therefore it does not contribute tothe signal whatsoever. In case (a) the specific energy (andhence the velocity) of the relevant half sphere is the sameas in the spherical case, yet the energy density has doubled.Therefore, the timescale increases by 2 / , and the flux re-mains the same (the specific energy as well as the mass inhalf the sphere are identical to the spherical specific energyand total mass). On the other hand, in case (b) the relevanthalf spehere (the part which will contribute to the overallsignal) contains the same mass density as in the sphericallysymmetric case, but twice the energy density. In this casethe timescale decreases by 2 − / , and the flux increases by2 (5 p − / . This illustrates the fact that the timescale andthe flux of an asymmetric system may either increase or de-crease with respect to it’s spherical equivalent depending onthe details of the mass and energy redistribution.We now turn to apply the same considerations to a set of N solid angle elements characterized by dM i , dE i , d Ω i wherethe index i enumerates the various solid angle elements.These parameters satisfy: (cid:80) i dM i = M tot , (cid:80) i dE i = E tot and (cid:80) i d Ω i = 4 π . Since we are primarily interested in thetotal light curve’s peak timescale, we can estimate this underthe piecewise spherical approximation by averaging over thepeak flux times of all constituent solid angle elements. Themean peak time is a sensible estimate for the overall peaktimescale as long as the temporal density is large enough.In other words, if the constituent light-curves of each solidangle overlap significantly around their peak fluxes, this esti-mate will be valid. This can be quantified by stating that thepeak times of all solid angle components should be within ∼ ∆ t of each other (see Eqs. 9,10).We define the flux weighted mean peak time normalizedby the original peak time of the spherical equivalent system t ( sphr ) peak as a dimensionless function τ : τ ( dΩ , dE , dM ) = (cid:104) t peak i (cid:105) t ( sphr ) peak = (cid:80) Ni =1 d Ω − / i dE − / i dM / i × (cid:32) dE p − i dM − p − i (cid:80) Nj =1 dE p − j dM − p − j (cid:33) (4 π ) − / E − / tot M / tot = (cid:80) ˜ d Ω i − / ˜ dE i p − ˜ dM i − p − (cid:80) ˜ dE j p − ˜ dM j − p − , (14)where in the last equality we have introduced the normalizedvariables ˜ x i ≡ x i /x tot (so that 0 < ˜ x i < (cid:80) ˜ x i = 1).This constraint eliminates three variables so that we are leftwith a total of 3 × ( N −
1) free parameters.Particularly interesting choices of parameters are “equaldistribution” points which retain the system’s spherical sym-metry (yielding τ = 1). Equal distribution points (e.d.) ex-ist besides the trivial case of ˜ dx k = 1 /N , as long as themass and energy densities remain identical to the spherical distribution, i.e. ˜ dE k / ˜ d Ω k = ˜ dM k / ˜ d Ω k = 1. These equaldistribution points are critical points of τ , since ∂τ ˜ ∂ Ω k (cid:12)(cid:12)(cid:12)(cid:12) e.d = ∂τ ˜ ∂E k (cid:12)(cid:12)(cid:12)(cid:12) e.d = ∂τ ˜ ∂M k (cid:12)(cid:12)(cid:12)(cid:12) e.d = 0 . (15)This extremum is a saddle point of τ and hence with appro-priate choice of variables one can shorten or lengthen thetimescale as desired ( τ < < τ respectively). We il-lustrate this fact graphically in Fig. 3 by plotting contoursof τ for the case of N = 2.We focus now on a particular scenario in which massand energy are redistributed under the additional constraintthat dE/dM is kept constant. Under this constraint, all masselements have equal velocity. In this case the equal distribu-tion point discussed previously becomes a global minimumof τ . The timescale of a constant specific energy redistribu-tion can therefore only increase with respect to the sphericaldistribution: τ | dE/dM =1 (cid:62) N = 2 in Fig. 3.The dashed black line in these figures depicts the curve˜ dE/ ˜ dM = 1. Under the constant specific energy constraint, τ should be along this curve. It can be seen that the curvealways passes through the equal distribution point at which τ = 1. As one moves away from this point along the con-straining curve (in either direction), the timescale only in-creases ( τ > F peak F ( sphr ) peak (cid:46) N (cid:88) i =1 (cid:18) ˜ dE ˜ dM (cid:19) − p − i ˜ dE i = 1 . (16)Of course this is a rough order of magnitude estimate thatcan be used as an upper bound for F peak , since the under-lying peak fluxes should not perfectly overlap.We have therefore shown that in case the specific en-ergy is kept constant in a redistribution of mass and energy,the optically thin synchrotron peak time (as estimated bythe weighted mean τ ) will necessarily increase, and the peakflux will not change with respect to the spherical equivalentlight-curve. Of course one can always redistribute the massand energy of the system creating a fast moving jet compo-nent. Such a jet will peak at earlier times than the sphericalequivalent system and may even lower the average peak time(so that τ < τ , breaks down as an es-timate of the overall timescale, as the jet component’s peakflux may exhibit itself as a completely separate componentin the light-curve (in other words the two components will c (cid:13) , 000–000 adio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry ˜ dM ˜ d E τ ( ˜ d Ω = 0 . , ˜ dM, ˜ dE ) (a) ˜ dM ˜ d E τ ( ˜ d Ω = 0 . , ˜ dM, ˜ dE ) (b) ˜ dM ˜ d E τ ( ˜ d Ω = 0 . , ˜ dM, ˜ dE ) (c) ˜ dM ˜ d E τ ( ˜ d Ω = 0 . , ˜ dM, ˜ dE ) (d) Figure 3.
The weighted mean peak time (normalized with respect to the spherical peak time) - τ , as a function of the redistributionvariables ˜ d Ω , ˜ dE, ˜ dM , for N = 2 angular divisions. Plotted are contours of τ in the ( ˜ dE, ˜ dM ) plane, where successive figures (from top leftto bottom right) depict decreasing values of ˜ d Ω. The equal distribution point can be spotted in Fig. 3(a) at ( ˜ d Ω , ˜ dE, ˜ dM ) = (0 . , . , . d Ω , ˜ dE, ˜ dM ) = (0 . , . , . τ at these is always 1. The dashed black line corresponds to the constraint˜ dE/ ˜ dM = 1, which physically means that the specific energy (and hence the velocity) is kept fixed. This line always passes through theequal distribution point, which serves as a minimum along this line. Therefore, with the additional constraint keeping ˜ dE/ ˜ dM fixed, thetimescale necessarily increases ( τ (cid:62) not significantly overlap). As long as the dynamic ejecta re-tains a roughly spherical distribution, one should not expectsuch a jet component or any other strong variations in spe-cific energy. In this case our results stand and we expecta longer overall light-curve timescale. In particular we willshow in § We turn now to apply our method to the results of neutronstar merger outflow simulations by Rosswog et al. (2014),from which we obtain the 3D hydrodynamic configuration of the dynamical ejecta. These SPH simulations are a contin-uation of previous works (Rosswog 2013; Piran et al. 2013)which simulated neutron star mergers till several millisec-onds after the merger event, and find that roughly 10 − M (cid:12) is dynamically ejected from the mergers, with a strong de-pendence on any asymmetry in the binary masses (equalmass mergers do not tidally disrupt as much mass as non-equal binaries). Rosswog et al. (2014) followed up by simu-lating the long term evolution of the dynamical ejecta, tak-ing into account the effects of nuclear heating on the outflowdynamics. The major hydrodynamic effect of this nuclearenergy deposition is in accelerating slightly the outflow anddiminishing significant irregularities in the outflow, essen-tially smoothing out the ejecta’s angular distribution (seeRosswog et al. 2014, for more details). Since our calcula-tions attempt to characterize the affects of asymmetry onthe ejecta’s afterglow signal, it is important we take the c (cid:13)000
The weighted mean peak time (normalized with respect to the spherical peak time) - τ , as a function of the redistributionvariables ˜ d Ω , ˜ dE, ˜ dM , for N = 2 angular divisions. Plotted are contours of τ in the ( ˜ dE, ˜ dM ) plane, where successive figures (from top leftto bottom right) depict decreasing values of ˜ d Ω. The equal distribution point can be spotted in Fig. 3(a) at ( ˜ d Ω , ˜ dE, ˜ dM ) = (0 . , . , . d Ω , ˜ dE, ˜ dM ) = (0 . , . , . τ at these is always 1. The dashed black line corresponds to the constraint˜ dE/ ˜ dM = 1, which physically means that the specific energy (and hence the velocity) is kept fixed. This line always passes through theequal distribution point, which serves as a minimum along this line. Therefore, with the additional constraint keeping ˜ dE/ ˜ dM fixed, thetimescale necessarily increases ( τ (cid:62) not significantly overlap). As long as the dynamic ejecta re-tains a roughly spherical distribution, one should not expectsuch a jet component or any other strong variations in spe-cific energy. In this case our results stand and we expecta longer overall light-curve timescale. In particular we willshow in § We turn now to apply our method to the results of neutronstar merger outflow simulations by Rosswog et al. (2014),from which we obtain the 3D hydrodynamic configuration of the dynamical ejecta. These SPH simulations are a contin-uation of previous works (Rosswog 2013; Piran et al. 2013)which simulated neutron star mergers till several millisec-onds after the merger event, and find that roughly 10 − M (cid:12) is dynamically ejected from the mergers, with a strong de-pendence on any asymmetry in the binary masses (equalmass mergers do not tidally disrupt as much mass as non-equal binaries). Rosswog et al. (2014) followed up by simu-lating the long term evolution of the dynamical ejecta, tak-ing into account the effects of nuclear heating on the outflowdynamics. The major hydrodynamic effect of this nuclearenergy deposition is in accelerating slightly the outflow anddiminishing significant irregularities in the outflow, essen-tially smoothing out the ejecta’s angular distribution (seeRosswog et al. 2014, for more details). Since our calcula-tions attempt to characterize the affects of asymmetry onthe ejecta’s afterglow signal, it is important we take the c (cid:13)000 , 000–000 B. Margalit & T. Piran most realistic initial angular distribution as conditions forour calculations.We use snapshots of the dynamical ejecta, taken fromGrossman et al. (2014) ∼
10 s after the merger event, as ini-tial conditions for our present calculations. We choose snap-shots taken several seconds after the merger, as at thesetimes the ejecta’s interaction with the surrounding ISMis negligible and hence the calculations of Rosswog et al.(2014), which did not account for these interactions, arevalid. On the other hand, these times are late enough sothat the nuclear decays have deposited practically all theirenergy into the outflow and affected it’s dynamics. Rosswoget al. (2014) show that from this point onward the ejectacontinues expanding homologeously (neglecting ISM inter-actions), and therefore the nuclear decays have little effecton the outflow at later times.We examine outflows from simulations of three partic-ular systems, distinguished by the neutron star masses: acanonical equal mass 1 . M (cid:12) -1 . M (cid:12) binary, a slightly asym-metric binary with masses 1 . M (cid:12) -1 . M (cid:12) , and an extremelyasymmetric system with masses 1 . M (cid:12) -1 . M (cid:12) . We denotethese systems ‘ns14ns14’, ‘ns14ns13’, and ‘ns16ns12’ respec-tively. We additionally restrict our focus in this work to thecase where ν obs = 1 . ν m , ν a < ν obs , regime. Through-out this work we take the microphysical system parametersas: (cid:15) e = 0 . (cid:15) B = 0 .
1, and p = 2 .
5, and use an ISM num-ber density of n = 1 cm − , but these parameters are easilychanged with qualitatively well known effects.This section is divided into two subsections loosely cor-responding to stages (i) and (ii-iv) of our method respec-tively (see Fig. 2). In § § § Using spherical coordinates centered on the merger location,we divide each system’s outflow into N = 20 ×
20 solid an-gles, keeping the number of SPH particles in each solid anglefixed. In Appendix A we show that this arbitrary binningallocation does not affect our results.Figure 4 shows mollweide projections of the character-istic energy density dE ( (cid:62) v peak ) /d Ω and characteristic ve-locity v peak of each angular bin for the three systems. Theseare a reduction of each bin’s cumulative distribution E ( (cid:62) v )into a single-velocity shell equivalent which reproduces thesame peak flux and the approximate peak time (see § θ = 0)so that the outflow is somewhat ‘flattened’ (see also Fig. 6of Rosswog et al. 2014). The qualitative difference betweenequal and non-equal mass mergers is mostly apparent inthat non-equal mergers show a stronger asymmetry insidethe equatorial (merger) plane. This breakdown of the axialsymmetry is a trace of the tidally disrupted tails producedduring the merger events. As expected unequal mass binariestidally disrupt each other more significantly than equal masssystems.The azimuthal symmetry breaking is especially promi-nent in the case of the ns16ns12 merger shown in Fig. 4(b).Here a large jump in the characteristic velocity is apparentat ϕ ≈ π/
3, with no coinciding jump in the energy density.In fact, while the ϕ ≈ π/ θ ≈ ϕ further from this point, we find a large gap region contain-ing no ejecta. This gap region is contained within the first(0 < ϕ (cid:46) π/
3) bin, hence the sudden drop in v peak . There-fore this unsmooth transition has it’s roots in the asymmetryof the underlying equatorial distribution.To illustrate this point further we plot the ns16ns12SPH particle distribution projected onto the equatorialplane in Fig. 5(a). The particles’ velocities at this stage arevery nearly proportional to their radii (as they are freelyexpanding) so that the distance from the center is a tracerof the velocities. The angular binning in ϕ is also plotted inthis figure, allowing comparison with the mollweide projec-tion of Fig. 4(b). The red lines indicate the two leftmost binsin Fig. 4(b) over which the velocity jump occurs. This figureclearly shows that one bin is dominated by the end of thetidally disrupted tail, while the other bin is dominated byslower moving ejecta and encompasses the large gap regionin the distribution (hence the significant angular size of thisbin).We show that the gap region in the equatorial distri-bution does not affect the resulting light-curve significantly,and therefore does not introduce any binning bias. This isshown by excluding the gap region from the analysis andcomparing the results with and without this region. Figure5(b) illustrates the SPH particle distribution in the equato-rial plane for the renewed analysis. The shaded region cor-responds to the excluded angle encompassing the gap in thedistribution, and the radially extending grey curves showthe renewed binning of the system. The system is rotatedclockwise by 1 .
95 rad with respect to the original distribu-tion (depicted in Fig. 5(a)) so that the excluded region endsat ϕ = 2 π and the first bin begins at ϕ = 0. Figure 6(a)continues by plotting the characteristic velocity map of thissystem, where the black region corresponds to the gap whichhas been excluded from the analysis. The velocity in thiscase varies smoothly and monotonically from right to leftas the bins are dominated by faster moving ejecta. Finally,Fig. 6(b) depicts the light-curves for the original distribu-tion, and the distribution from which the gap region has c (cid:13) , 000–000 adio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry dE ( ≥ v peak ) d Ω ergrad $ (a) ns16ns12 - characteristic energy density v peak [ c ] (b) ns16ns12 - characteristic velocity dE ( ≥ v peak ) d Ω ergrad $ (c) ns14ns13 - characteristic energy density v peak [ c ] (d) ns14ns13 - characteristic velocity dE ( ≥ v peak ) d Ω ergrad $ (e) ns14ns14 - characteristic energy density v peak [ c ] (f) ns14ns14 - characteristic velocity Figure 4.
Mollweide projections of the velocity and energy density angular distributions for ns16n12 (a,b), ns14ns13 (c,d) and ns14ns14(e,f). The energy of the ejecta is concentrated in all three cases into a region close to the binary orbital plane ( θ = 0). Additionally, theejecta energy is concentrated in smaller ϕ angles for larger binary mass ratios. This is a tracer of the tidally disrupted tail, which is moreprominent with increasing binary mass ratio. While the characteristic velocity varies by a modest factor of ∼
2, the energy density variesby nearly 3 orders of magnitude and it is therefore the decisive factor in determining the peak timescale.c (cid:13)000
2, the energy density variesby nearly 3 orders of magnitude and it is therefore the decisive factor in determining the peak timescale.c (cid:13)000 , 000–000 B. Margalit & T. Piran [h] −10 −8 −6 −4 −2 0 2 4 6 8x 10 −8−6−4−20246810 x 10 x [cm] y [ c m ] (a)(b) Figure 5.
The SPH particle distribution for ns16ns12 projectedon the equatorial (x-y) plane. Ejecta velocities are proportional totheir radius. The radially extending grey lines correspond to theangular binning of the system in ϕ . (a) The original distribution.The red lines emphasize the bins between which a significant ve-locity jump occurs (see Fig. 4(b)). This plot illustrates the factthat the velocity jump is due to the transition between the end ofthe tidally disrupted tail and a region containing slower movingejecta. (b) The rotated distribution for the renewed analysis. Theshaded black area indicates the gap region excluded from thisanalysis. v peak [ c ] (a) characteristic velocity map −1 t (yr) F ν ( m J y ) Gap region excluded from distributionOriginal distribution (b) light-curve comparison
Figure 6. (a) A velocity map of the ns16ns12 system, excludingthe gap region. The black area corresponds the excluded region.This distribution is rotated with respect to Fig. 4(b) (see Fig.5(b)). The figure illustrates that the velocity varies smoothly andmonotonically from right to left apart from the jump across theexcluded black region. (b) A comparison of the ns16ns12 light-curves produced by the original distribution and by the distri-bution excluding the gap region. The two light-curves are nearlyidentical, indicating that the gap region does not affect the re-sults. been excluded. It is evident from this figure that the gap re-gion is inconsequential to the results, and that including it inthe original analysis does not affect or bias the results. Thisis clear as the two light-curves presented in Fig. 6(b) overlapsuch that they are nearly indistinguishable. Additionally, weshow more generally in Appendix A that the arbitrarinessof any binning choice does not affect the results for any ofthe three merger scenarios.
After examining the angular binning and the distribution ofcharacteristic variables we calculate the dynamic evolutionof each angular bin according to it’s isotropic equivalent evo-lution, and the resulting synchrotron light-curves. The light c (cid:13) , 000–000 adio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry [h] −2 −1 t (yr) F ν ( m Jy ) ns16ns12 Light−Curve spherical equivalent light−curvecalculated light−curve (a) ns16ns12 −2 −1 t (yr) F ν ( m Jy ) ns14ns13 Light−Curve spherical equivalent light−curvecalculated light−curve (b) ns14ns13 −2 −1 t (yr) F ν ( m Jy ) ns14ns14 Light−Curve spherical equivalent light−curvecalculated light−curve (c) ns14ns14 Figure 7.
The light-curves for the three binary merger systems. The flux is calculated at ν obs = 1 . ν m , ν a < ν obs for all bin components) and taking a typical distance of d = 10 cm. The dashed grey lines correspond to the sphericalequivalent light-curves for these systems, which we compare to calculations of the present work (red lines). In all three cases the light-curves peak at later times than their spherical equivalents, and attain roughly the same peak fluxes. These two effects increase in theirsignificance for larger binary mass ratios, which may be explained since these produce less symmetric outflows. Both these results arewell explained theoretically (see §
4) if the characteristic velocities do not vary significantly between solid angle bins (see Fig. 4 and 8).Quantitative results for the peak time and flux are presented in table 1, along with these signals detectability prospects (see § curves are compared to their spherical equivalents in Fig. 7.It is not surprising that the less symmetric the underlyingsystem is, the greater is the differences between the resultsand their spherical equivalents. In particular the peak timeincreases dramatically in the ns16ns12 system when com-pared with the timescale increase of the ns14ns14 merger.The basic reason for this result is the fact that nearly allof the ns16ns12 outflow is confined to the binary plane,and moreover it is significantly unevenly distributed withinthis plane, as most of the energy is concentrated in a smalltidally disrupted tail. Thus, most of the system’s energy isdistributed over a small fraction of the sky, leading to asignificantly larger energy density dE/d Ω when comparedwith the spherical equivalent situation (in which the energyis distributed evenly over 4 π rad ). As the peak time is pro-portional to ( dE/d Ω) / , this immediately translates into anincrease in the light-curve’s timescale. This reasoning worksonly as long as the characteristic velocity does not vary sig- nificantly, so that the timescale is determined almost entirelyby the energy density.Indeed, when looking back at figures 4(a) and 4(b) wenotice that while the energy density fluctuates over threeorders of magnitude, the velocity varies by a modest factor oftwo. Thus, to a zeroth order, the main effect contributing tothe timescale delay is the compression of the ejected matterinto a narrow solid angle which gives rise to a dramaticincrease in the energy density dE/d Ω (or equivalently tothe mass density dM/d
Ω).In order to verify that the timescale is mainly influ-enced by the concentration of the outflow into smaller solidangles, we plot in Fig. 8 the peak time of each solid an-gle’s signal versus the two independent variables contribut-ing to the peak time, dM ( (cid:62) v peak ) /d Ω and v peak . Here, wechoose to parameterize the peak time using the mass in-stead of energy density, since the energy density correlates c (cid:13)000
Ω).In order to verify that the timescale is mainly influ-enced by the concentration of the outflow into smaller solidangles, we plot in Fig. 8 the peak time of each solid an-gle’s signal versus the two independent variables contribut-ing to the peak time, dM ( (cid:62) v peak ) /d Ω and v peak . Here, wechoose to parameterize the peak time using the mass in-stead of energy density, since the energy density correlates c (cid:13)000 , 000–000 B. Margalit & T. Piran by definition with the characteristic velocity, whereas themass should be independent of it. From Eq. 12 and 13 wefind that the peak time scales with the deceleration time,i.e. t peak ∝ ( dM ( (cid:62) v peak ) /d Ω) / v − peak . Figure 8(a) showsclearly that this relationship holds quite strongly, where thesmall deviations arise from the fact that an additional degreeof freedom (which depends on the detailed energy distribu-tion E ( (cid:62) v peak )) determines the exact value of t peak . In anycase it is evident from this figure that t peak lies within thelower and upper bounds given by Eq. 13 and 12 which areillustrated as dashed grey lines in Fig. 8(a).Figures 8(b) and 8(c) show the dependence of the peaktime on each one of the parameters independently. Alsoshown is the Pearson Correlation between the two variables.Clearly, t peak depends strongly on the mass density andonly very weakly on the velocity distribution. This supportsour earlier conclusion that the timescale is governed by theejecta’s concentration into smaller portions of the sky, andnot by dramatic changes in the characteristic velocities.Finally, we compare the numerical results with the re-sults of §
4, in which we have shown that the flux weightedaverage of each solid angle’s peak time must increase for anyredistribution in which the specific energy is kept fixed. Fig-ure 8 shows that the peak time correlates only very weaklywith the characteristic velocity, so that the specific energyis inconsequential and may be treated as fixed. This jus-tifies a comparison with the fixed specific energy scenariodiscussed in §
4. Still, in order to complete the comparisonwe must ascertain that the flux weighted mean peak time, τ , is a valid estimate of the total timescale. In § τ is a reasonable estimate as long as the differentcontributing light-curves overlapped significantly. This wasquantified by estimating the necessary overlap width, ∆ t ,for single-velocity shells. Non-single-velocity shells give riseto larger light-curve widths, yet we can still use the single-velocity shell results as lower bounds on ∆ t in this case.Figure 9 depicts the flux weighted histograms of eachbin’s peak time for the three systems , as well as the ad-ditional relevant timescales of the problem. Also illustratedare arrows showing the width, ∆ t , of of the underlying light-curves around the weighted mean. Short/dark arrows showthe lower bounds on this width (as given by Eq. 9 and10 for single-velocity shells), and the longer/lighter arrowsshow the actual measured width of a typical component’slight-curve. The distribution of widths is centered around t ( − )1 / ∼ . × t peak , and t (+)1 / ∼ . × t peak . It is clear thatmost of the components significantly overlap and hence weshould expect the weighted mean to serve as a reasonableestimate of the total timescale. This is also apparent directlyby comparing the calculated weighted mean (depicted as adashed red line in these figures) with the actual total light-curve peak time (depicted as dashed green lines). In all threecases these curves are in proximity of each other when com-pared with the spherical equivalent peak time (dashed blueline). Another point worth mentioning is the fact that theactual peak time is always greater than the weighted meanestimate. This is most likely caused by the asymmetry ofthe light-curve around the peak (the light-curve varies atlate/early times as t − . / t respectively). As this asymme-try is unaccounted for by the weighted mean peak time, it [] t p e a k ( y r ) t dec ∝ dM ( ≥ v peak ) d Ω / v − peak (yr) Correlation = 0.999 t dec p − p +3 t dec (a) t p e a k ( y r ) dM ( ≥ v peak ) d Ω / (arbitrary units) Correlation = 0.789 (b) t p e a k ( y r ) v − peak (arbitrary units) Correlation = 0.136 (c)
Figure 8.
The correlation of each solid angle bin’s peak time withits mass and velocity distributions for the ns16ns12 system. Fig.8(a) clearly shows that the peak time scales as t dec . Figures 8(b)and 8(c) illustrate that the main factor determining the timescalefor this system is the mass density rather than velocity. This isquantified via the Pearson Correlation between the variables ineach case. The fact the peak time is determined mainly by themass density allows a comparison with the theoretical expecta-tions in case the specific energy is fixed (see § t peak , asgiven by Eq. 12 and 13. The points are all well within these limitswith very little scatter, indicating that the functional form of thecumlative energy distribution is nearly identical for all bins.c (cid:13) , 000–000 adio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry (a) ns16ns12 (b) ns14ns13(c) ns14ns14 histogramWeighted h t peak i Total light-curve peakSpherical equivalent peaktypical ∆ t lower bound on ∆ t Figure 9.
The flux weighted distributions of peak times for each solid angle component for the three systems. The red dashed linesindicate the weighted mean of these components. The green dashed lines depict the actual peak time of the total light-curve. The bluedashed lines show the peak time in the spherical equivalent case. Arrows depict the light-curves’ width around the weighted mean, definedas the times at which the flux diminishes to half it’s peak value. Dark/short arrows show the lower bounds on this width (Eq. 9 and10), whereas light/long arrows show the typical widths for our systems. The figures show that the weighted mean peak time is a goodestimate of the total light-curve’s peak time. This is consistent with the fact that the underlying components overlap significantly (theirwidth’s overlap over some region). can be expected that the actual timescale be delayed slightlyfurther than this estimate would suggest.Analyzing ‘realistic’ merger dynamical outflows usingour piecewise spherical approximation, we have found anincrease in the timescale of the merger radio light-curves incomparison with their spherical equivalent estimates studiesin previous works (Nakar & Piran 2011; Piran et al. 2013).This shifts the peak timescales to the order of ∼
10 years. Incontrast, the peak flux of the light-curves is comparable toit’s spherical equivalent and it does not change significantlywhen taking into account the nontrivial outflow geometry.Both these results are consistent with the analytic estimatespresented in § We turn now to discuss the detectability of compact binarymergers’ radio flares based on these results. The number ofradio flares observable by an instrument with detection level F lim in an all sky snapshot is N all − sky ≈ R V ∆ t, (17)where R is the merger event rate, V is the detectable vol-ume and ∆ t is the time that the flux is above the detectionthreshold F lim . We can improve on Eq. 17 by taking intoaccount the dependence of ∆ t on the distance to the mergerevent as follows: N all − sky = (cid:90) ∞ Θ ( F ( r ) − F lim ) R ∆ t ( r ) 4 πr dr = (18)= R π (cid:18) L πF lim (cid:19) / (cid:90) ∞
32 ∆ t ( χ ) χ − / dχ ≡ R V (cid:104) ∆ t (cid:105) . c (cid:13)000
32 ∆ t ( χ ) χ − / dχ ≡ R V (cid:104) ∆ t (cid:105) . c (cid:13)000 , 000–000 B. Margalit & T. Piran
Here L is the merger’s peak luminosity, χ ≡ F ( r ) /F lim ,and we assumed isotropy and neglected any cosmologicalcorrections to the merger rate R ( r ) and peak flux F ( r ) = L/ (4 πr ). ∆ t ( χ ) is the time the flux is above the detectionlimit, and it depends on the exact shape of the radio flarelight-curve as well as on the ratio between the peak flux F and the detection threshold F lim . In the simple case of asingle-velocity shell we can write down ∆ t ( χ ) explicitly us-ing the light-curve scaling relations given by Eq. 8 (assumingagain that ν m , ν a < ν obs ),∆ t ( χ ) = (cid:16) χ p − − χ − / (cid:17) t dec . (19)Using this result we continue in finding the relevant time-scale (cid:104) ∆ t (cid:105) : (cid:104) ∆ t (cid:105) = 18(5 p + 3)11(45 p − t dec ≈ . t dec (20)Plugging this expression into Eq. 18 and simplifyingthe result gives an estimate for the number of merger rem-nants observable over the whole sky. This result is similar tothat given by Nakar & Piran (2011), which estimated that∆ t ≈ t dec , but is larger due to a round-off error in theircalculations. Expressing t dec and L in terms of the single-velocity outflow energy E and velocity βc using Eq. 4 and 1,we find that the the number of spherically symmetric single-velocity shell radio remnants observable at 1 . N . all − sky ≈ (cid:18) E erg (cid:19) / (cid:16) (cid:15) B . (cid:17) p +1)8 (cid:16) (cid:15) e . (cid:17) p − × n p +124 (cid:18) β . (cid:19) p − (cid:18) F lim . (cid:19) − / × (cid:18) R − yr − (cid:19) (cid:18) (cid:104) ∆ t (cid:105) /t peak . (cid:19) (21)We use our own canonical values in evaluating this ex-pression, which are somewhat different than those chosen byPNR13. We notice that the expression depends sensitivelyon the outflow velocity: N all − sky ∝ β . which we take to beslightly lower than 0 . c . N all − sky depends also on the ejectaenergy which should be in the range of E ∼ − erg forthe outflows considered. A larger value than our canonical10 erg will cause an increase in N all − sky compared with ourcanonical value. Additionally, these results are sensitive tothe ISM density ( N all − sky ∝ n . ) so that binaries mergingin less tenuous surroundings will decrease the detectabilityas well. Based on current observation of Galactic neutronstar binaries we expect that at least a significant fraction ofthem be located in the disk of their host galaxy, where thetypical ISM densities are ∼ − (Draine 2011). A caveatto this point is that it stems from the observed Galactic bi-nary neutron star population, which is both small, and con-fined within Milky-Way alone, so that other values of n arequite possible and this parameter is left ill-constrained. Forexample, Fong et al. (2014) find n ≈ × − −
30 by evaluat-ing the 130603B short GRB afterglow. This case illustratesthe large uncertainty in ascertaining the circum-binary den-sity using sGRB afterglow observations, and although kickscould relocate neutron star binaries outside their host galaxywhere n would be quite small, there is no clear evidence torule out the galactic merger scenario. Using these light-curves we can estimate the detectabil-ity for each merger case by evaluating Eq. 18 numerically.The results are given in table 1, which also states thepeak flux and peak time found for the different mergerlight-curves. Additionally the table includes the values of (cid:104) ∆ t (cid:105) /t peak which are the corrections to the approximation∆ t ≈ t peak taking into account the specific light-curve shape.These calculated values (which are always of order unity)should be taken with a grain of salt because these light-curves may not necessarily reproduce the exact light-curvesdue to approximations used (which may break down afterthe ejecta begins it’s deceleration at t peak ). In any case, aconservative lower bound on N all − sky is given by replacing (cid:104) ∆ t (cid:105) /t peak with the single velocity shell value of 0 .
86 (seeEq. 20). This substitution is straightforward since N all − sky scales linearly with (cid:104) ∆ t (cid:105) /t peak .We also provide in table 1 the outflow velocity at t peak and the energy of the outflow above this velocity for thespherical equivalent cases. As explained in §
2, these are thecharacteristic energy and velocity which allow us to estimatethe peak flux and time using single-velocity shell results (byexchanging
E, v with E ( (cid:62) v peak ) , v peak in the relevant equa-tions). Substituting these characteristic values into Eq. 21for instance, yields similar results to N all − sky calculated nu-merically. The only reason the results do not coincide exactlyis the uncertainty in evaluating the peak time using the char-acteristic velocity and energy. This time can be found usingthese two variables only up to a factor of order unity (andconfined between 0 . −
1, see Eq. 12, 13).Our results for N all − sky are significantly larger thanthose estimated by previous works (Nakar & Piran 2011;Piran et al. 2013), even for the spherical equivalent cases.There are several reasons for this increase. The first is acalculation round-off error in these papers, which leads toa factor of ∼ ∼
20 as estimated inNakar & Piran (2011) as opposed to ∼
110 given in thiswork, taking (cid:104) ∆ t (cid:105) = t peak and the canonical variables ofNakar & Piran (2011) for the comparison). The second hasto due with the fact that the typical energies and velocitiesof the merger outflows are different than those evaluatedNakar & Piran (2011). Specifically, for typical energies of 2 × erg and velocities of 0 . c we get a further increase of (2 × / ) / (0 . / . ≈ . ∼ N all − sky determines thenumber of signals above the detection threshold of a givenradio facility, it is an entirely different question whetherthese signals can be identified as merger radio flares. Themain difficulty our present results pose to the identificationprocess is the increase in time-scale compared with previ-ous works (Piran et al. 2013). The time-scales we estimatein this work are on the order of ∼
10 yr. This long periodwill severely influence the ability to identify these signalsas transients and therefore associate them with merger out-flows.In order to identify the source signal as a transient, a c (cid:13) , 000–000 adio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry F peakν [mJy] t peak [yr] N all − sky (cid:104) ∆ t (cid:105) /t peak v peak [ c ] E (cid:0) (cid:62) v peak (cid:1) [erg]ns16ns12 0 .
53 14 9 . × . .
53 6 . . × .
59 0 .
187 7 . × ns14ns13 0 .
090 10 . . .
108 4 . .
32 0 .
174 1 . × ns14ns14 0 .
084 8 . . .
093 5 . .
27 0 .
153 2 . × Table 1.
Numerical results illustrating the detectability of the systems considered. Each row designates one of the three systems and isdivided into a top portion stating the results of the current analysis, and bottom portion stating the results for the spherical equivalentsystems. The results are calculated at an observed frequency of ν obs = 1 . d = 10 cm, ISM densityof n = 1cm − , and microphysical parameters (cid:15) e = (cid:15) B = 0 . p = 2 .
5. The detectability is estimated by the number of observablesources in an all sky snapshot, N all − sky , assuming a merger rate of R = 300Gpc − yr − and a detection threshold of F lim = 0 . significant change in it’s flux must be measured. This is usu-ally achieved by returning and pointing the radio telescopetowards the source at several different times, but on the scaleof a decade (on which the signal changes significantly), mostinstruments will have been upgraded or replaced by moreadvanced facilities, making any recurring measurements dif-ficult to conduct and contaminated by instrumental changes.An alternative method of classifying a transient event is de-tecting smaller variations in the signal which would occurover smaller time periods. The problem with this method isthat a typical detection should be very close to the detec-tion threshold at a low S/N ratio, and would therefore notbe sensitive enough to small variations in the flux. Thus,the only way such an approach could work were if an un-likely strong signal were detected, the chances of which aredrastically diminished. We conclude this section by comparing the detectability ofmerger radio flares with the detectability of GRB orphan af-terglows, which are often considered as more promising elec-tromagnetic counterparts to mergers. The well constrainedobservables of GRBs are the isotropic equivalent energy ofthe burst emitted in gamma-rays E γ, iso , and the rate of ob-served GRBs R obsGRB . Orphan afterglows on the other handare produced after the beamed ultra-relativistic jet has ex-panded into a roughly spherical shell, emitting isotropic ra-diation. Afterglow luminosities and rates should thereforedepend on the total energy of the GRB and overall GRB rate(not the observed rate which only include bursts pointing to-wards Earth). Assuming that the GRB emission mechanismis very efficient, the isotropic energy of the burst should becomparable to the isotropic energy emitted in gamma-rays, E iso ∼ E γ, iso (Nakar 2007). The total burst energy and over-all rate can then be calculated based on the well constrained E iso and R obsGRB , if we assume that the initial GRB jet coversa fraction f b of the sky: E tot = E iso f b ; R GRB = R obsGRB f − b . (22)We use these results in estimating the number of GRBorphan afterglows detectable at 1 . N orphanall − sky , similarlyto the calculation for merger radio flares. The observed rateof GRBs is enery dependent, and can be approximated as apower-law over the energy range of 10 − erg: R obsGRB ≈
10 Gpc − yr − (cid:0) E iso / erg (cid:1) − α where α ≈ . − N orphanall − sky ≈ . f / b (cid:18) E iso erg (cid:19) / − α (cid:16) (cid:15) B . (cid:17) p +1)8 (23) × (cid:16) (cid:15) e . (cid:17) p − n p +124 (cid:18) F lim . (cid:19) − / , where we take (cid:104) ∆ t (cid:105) /t peak = 1, and β = 1 since the GRBoutflow is ultra-relativistic.Using an expression similar to Eq. 23, Levinson et al.(2002) showed that although the GRB afterglow rate in-creases with smaller beaming factors, the overall detectabil-ity will decrease due to the decrease in the afterglow energy(see Eq. 22). Most of the variables in Eq. 23 are well con-strained, with the exception of the beaming factor f b , whichcould be on the order of f b = 1 /
20 (by definition f b < § ∼
10 yr, GRB afterglows expand relativisticallyand will most likely be dominated by low energy bursts of ∼ erg (because of the observed luminosity function) sothat they will peak on much shorter time-scales of ∼ We review here our main results and discuss their impli-cations in a broader context. In particular we address thequestion of whether the radio flares found here show promise c (cid:13)000
10 yr, GRB afterglows expand relativisticallyand will most likely be dominated by low energy bursts of ∼ erg (because of the observed luminosity function) sothat they will peak on much shorter time-scales of ∼ We review here our main results and discuss their impli-cations in a broader context. In particular we address thequestion of whether the radio flares found here show promise c (cid:13)000 , 000–000 B. Margalit & T. Piran as detectable electromagnetic signals from compact binarymergers. We finish by suggesting further research directionsinto this and related problems.Throughout this work we have studied radio signals re-sulting from the deceleration of merger dynamic outflows asthey interact with the surrounding ISM. We have expandedon previous works (Nakar & Piran 2011; Piran et al. 2013) bytaking into account the asymmetry in these outflows (whichhave previously been approximated as spherically symmet-ric). We focus on observed frequency ν obs = 1 . ∼ ∼ n = 1cm − . The ISM densityis an ill constrained parameter and could be orders of mag-nitude lower if the merger occurs outside it’s host galaxydisk. In such cases the timescale will be even longer than es-timated above, as t peak ∝ n − / . We additionally find thatthe light-curve’s peak flux remains roughly the same as inthe approximated spherical treatment.The implications of our findings are important in assess-ing the detectability prospects of these signals. The longer timescale is a double edged sword in addressing this issue.On one hand, longer timescales make the identification moredifficult. This is important for both all sky searches and forattempts to identify radio flares following identification ofa GW candidate. On the other hand, the longer timescalesmean a larger number of observable sources in an all skysnapshot. The combined effect will depend on the telescopeused and on the observation strategy as well as the unknownbackground field of radio transients that might compete withthese events.Support for this work was provided by an ERC ad-vanced grant “GRBs”, by the I-CORE Program of thePlanning and Budgeting Committee and The Israel ScienceFoundation grant No 1829/12 and by ISA grant 3-10417 . REFERENCES
Barnes J., Kasen D., 2013, ApJ, 775, 18Draine B. T., 2011, Physics of the Interstellar and Inter-galactic MediumEichler D., Livio M., Piran T., Schramm D. N., 1989, Na-ture, 340, 126Faber J. A., Rasio F. A., 2012, Living Reviews in Relativity,15, 8Fong W., Berger E., Metzger B. D., Margutti R., ChornockR., Migliori G., Foley R. J., Zauderer B. A., Lunnan R.,Laskar T., Desch S. J., Meech K. J., Sonnett S., DickeyC., Hedlund A., Harding P., 2014, ApJ, 780, 118Grossman D., Korobkin O., Rosswog S., Piran T., 2014,MNRAS, 439, 757Guetta D., Piran T., 2006, A&A, 453, 823Hotokezaka K., Kiuchi K., Kyutoku K., Okawa H.,Sekiguchi Y.-i., Shibata M., Taniguchi K., 2013,Phys. Rev. D, 87, 024001Hotokezaka K., Piran T., 2015, ArXiv e-printsHulse R. A., Taylor J. H., 1975, ApJ, 195, L51Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25Kochanek C. S., Piran T., 1993, ApJ, 417, L17Kramer M., Stairs I. H., Manchester R. N., McLaughlinM. A., Lyne A. G., Ferdman R. D., Burgay M., LorimerD. R., Possenti A., D’Amico N., Sarkissian J. M., HobbsG. B., Reynolds J. E., Freire P. C. C., Camilo F., 2006,Science, 314, 97Kulkarni S. R., 2005, ArXiv Astrophysics e-printsLevinson A., Ofek E. O., Waxman E., Gal-Yam A., 2002,ApJ, 576, 923Li L.-X., Paczy´nski B., 1998, ApJ, 507, L59Metzger B. D., Mart´ınez-Pinedo G., Darbha S., QuataertE., Arcones A., Kasen D., Thomas R., Nugent P., PanovI. V., Zinner N. T., 2010, MNRAS, 406, 2650Nakar E., 2007, Phys. Rep., 442, 166Nakar E., Gal-Yam A., Fox D. B., 2006, ApJ, 650, 281Nakar E., Piran T., 2011, Nature, 478, 82Piran T., Nakar E., Rosswog S., 2013, MNRAS, 430, 2121Rosswog S., 2013, Royal Society of London PhilosophicalTransactions Series A, 371, 20272Rosswog S., Korobkin O., Arcones A., Thielemann F.-K.,Piran T., 2014, MNRAS, 439, 744Sedov L. I., 1946, Prikl. Mat. Mekh., 10, 241Tanaka M., Hotokezaka K., 2013, ApJ, 775, 113 c (cid:13) , 000–000 adio Flares of Compact Binary Mergers: the Effect of Non-Trivial Outflow Geometry Taylor G., 1950, Royal Society of London Proceedings Se-ries A, 201, 159Taylor J. H., Weisberg J. M., 1989, ApJ, 345, 434Weisberg J. M., Nice D. J., Taylor J. H., 2010, ApJ, 722,1030
APPENDIX A: BINNING ERRORS
In the following, we estimate the errors associated with thearbitrary choice of binning, and show that they are insignif-icant.To estimate this, we calculate the light-curve for differ-ent divisions of solid angles keeping the number of bins fixedat N = 20 ×
20. We enlarge and decrease randomly the sizeof each solid angle element and with it the number of SPHparticles it contains. As we need to keep a minimal numberof SPH particles in each bin for the analysis to work, we re-stricted the random noise so that the smallest possible solidangle bin would still satisfy this requirement. We addition-ally rotated the system randomly about the ϕ axis so as tobegin the division at different points (essentially so that thefirst bin will not be constrained to the same location everytime).Since the main variable we have studied in this work isthe peak time of such light-curves, we extract this parameterand asses it’s statistical variability as a measure of the errorsassociated with the arbitrary binning process. Fig. A1 showshistograms for the peak timescale extracted from analyzingthe same data using different angular discretizations. Thestandard deviations are of order ∼
13 13.5 14 14.5 15051015202530354045 t peak (yr) C o un t s ns16ns12 t peak Distribution histogramgaussian fi t: µ = 14 . σ = 0 . (a) t peak (yr) C o un t s ns14ns13 t peak Distribution histogramgaussian fi t: µ = 10 . σ = 0 . (b) t peak (yr) C o un t s ns14ns14 t peak Distribution histogramgaussian fi t: µ = 8 . σ = 0 . (c) Figure A1.
Distributions of the total light-curve’s peak time, t peak , calculated using an ensemble of randomized N = 20 × ∼ § (cid:13)000