Orbital Parameters of Binary Radio Pulsars : Revealing Their Structure, Formation, Evolution and Dynamic History
Orbital parameters of binary radio pulsars reveal the history of the pulsars’ formation and evolution including dynamic interactions with other objects. Advanced technology has enabled us to determine these orbital parameters accurately in most of the cases. Determination of post-Keplerian parameters of double neutron star binaries (especially of the double pulsar system PSR J0737-3039) provide clean tests of the general theory of relativity and in the future may lead us to constrain the dense matter Equations of State. For binary pulsars with main-sequence or white dwarf companions, knowledge about the values of the orbital parameters as well as of the spin periods and the masses of the pulsars and the companions might be useful to understand the evolutionary history of the systems. As accreting neutron star binaries lead to orbit circularization due to the tidal coupling during the accretion phase, their descendants binary millisecond pulsars are expected to be in circular orbits. On the other hand, dense stellar environments inside globular clusters cause different types of interactions of single stars with pulsar binaries like “fly-by”, “exchange”, and “merger”. All these interactions impart high eccentricities to the pulsar binaries. So it is quite common to get eccentric millisecond pulsar binaries in globular clusters and we find that “fly-by” causes intermediate values of eccentricities while “exchange” or “merger” causes high values of eccentricities. We also show that “ionization” is not much effective in the present stage of globular clusters. Even in the absence of such kinds of stellar interactions, a millisecond pulsar can have an eccentric orbit as a result of Kozai resonance if the pulsar binary is a member of a hierarchical triple system. PSR J1903+0327 is the only one eccentric millisecond pulsar binary in the galactic disk where stellar interactions are negligible. The possibility of this system to be a member of a hierarchical triple system or past association of a globular cluster have been studied and found to be less likely.
Inter University Centre for Astronomy and Astrophysics, Pune 411 007, India
Radio pulsars are rotating, magnetized neutron stars (NSs) emitting radio beams from magnetic poles. Their masses are in the range of , radii km, densities and surface values of the dipolar magnetic field are around gauss. Presently there are 1864 radio pulsars with their spin periods covering a wide range 1.39 ms 11.78 s
Pulsars are believed to be born with spin period s which increases very slowly with time (at the rate of ) as their rotational energy is emitted in the form of electromagnetic energy. Millisecond spin periods are explained as a result of spin-up process, which occurs when a neutron star has a normal star as its binary companion and accretes matter from the companion (at its red giant phase) causing transfer of angular momentum. Depending upon the evolution of the companion, the resultant millisecond pulsar can be either isolated or a member of a binary system. In this article we concentrate on binary pulsars and discuss how the knowledge of their orbital parameters help us to understand the structure, formation, evolution and dynamic history of pulsars. As an example, binary millisecond pulsars are expected to be in circular orbits as tidal coupling during the accretion phase lead to orbit circularization. But if there is no history of accretion phase, the orbit of the binary pulsar can be eccentric because of the asymmetric kick imparted to the pulsar during its birth in a supernova. A millisecond pulsar binary can also have an eccentric orbit if it experiences interactions with other near-by stars (which is possible in dense stellar environments like globular clusters) or if it is a member of a hierarchical triple.
It is noteworthy to mention here that in addition to radio wavelengths, pulsars emit in x-rays and/or gamma rays too. But radio pulsars are the most useful astronomical objects firstly because the underlying neutron stars are more stable than in the case of X-ray pulsars which have episodically varying accretion torques; secondly high sensitivities of radio telescopes and data analysis software make it possible to determine their spin periods very accurately. Moreover, for radio pulsars in binary systems ( of presently known radio pulsars are in binary systems), it is possible to determine the mass of the pulsar and its companion, size and shape of the orbit by modelling Keplerian and post-Keplerian parameters. Such accurate determination of those parameters are difficult for X-ray pulsars at present. That is why we concentrate only on binary radio pulsars when we try to use the knowledge about their orbital parameters to understand the properties of pulsars. The first binary pulsar PSR B1913+16 was discovered in 1974 at Arecibo observatory by R. Hulse and J. H. Taylor  which is a 59 ms pulsar having orbital period 0.323 days and orbital eccentricities 0.671. Interestingly, its mass companion star is most probably another neutron star the first discovered binary pulsar is also the first discovered double neutron star (candidate) system. Afterwords nine other double neutron star (DNS) systems have been discovered. Among these ten DNSs, seven systems are confirmed to be DNSs and the rest are most likely to be DNSs . These DNSs are useful tools to test the validity of the general theory of relativity  and even in principle to constrain the dense matter Equations of State . In Section 2 we discuss how the study of DNSs can lead to a better understanding of the structure of neutron stars ( the Equation of State), the physics of the pulsar magnetosphere and the formation scenario for the special case of the double pulsar system.
Most of the binary pulsars have either main sequence (MS) or white dwarf (WD) companions. Knowledge about the values of the orbital parameters as well as of the spin periods and the masses of these pulsars (and their companions) might be useful to understand the evolutionary history of these systems. The possible causes of eccentricities of millisecond pulsar binaries have been studied in this article. When these systems are in dense stellar environments as in globular clusters, the eccentricities can arise as a result of different types of interactions of single stars with pulsar binaries like “fly-by”, “exchange” and “merger” (Section 3). The eccentricity of the only eccentric millisecond pulsar PSR J in the galactic disk is difficult to explain as the hierarchical triple scenario faces difficulties to satisfy observed parameters of this system and the globular cluster origin also appears to be unlikely (section 4). All the findings of the present work have been summarized in section 5.
2 Double Neutron Star Binaries
After performing pulsar observations, one determines “time of arrivals” (TOAs) of pulses at the telescope. Then various clock corrections are applied to the TOAs, followed by some other corrections “Römer delay”, “Shapiro delay” and “Einstein delay” to calculate the TOAs with respect to the solar system barycenter. Then the correction for the frequency dependent delay (dispersion) by the interstellar medium is added. Finally, the pulse phase is fitted with time as a Taylor series expansion like where and are reference phase and time, and are the spin frequency and its first time derivative. In case of a binary pulsar, one needs to model the orbital motion of the pulsar, and before doing that same type of delays for the binary system ( another set of “Römer delay”, “Shapiro delay” and “Einstein delay”) should be taken into account. For a non-relativistic binary pulsar, five Keplerian parameters are required to describe the orbit, these are the orbital period (), projected semi-major axis of the pulsar orbit (, where is the inclination of the orbit with respect to the sky plane), orbital eccentricity (), longitude of periastron () and the epoch of the periastron passage (). The semi major axes of the pulsar orbit and the companion orbit can be written as and respectively where is the semi major axis of the relative orbit coming in the Kepler’s law . Here is the mass of the pulsar and is the mass of the companion. Pulsar timing analysis means first getting an initial fit for the orbital parameters by inspecting the variation of the pulse period (which is equal to the spin period of the pulsar) and then improving the values of the parameters through a phase connected timing solution. But for relativistic binaries ( when the companion of the pulsar is another neutron star or a white dwarf), five post-Keplerian (PK) parameters are needed, which are the relativistic advance of periastron (), redshift parameter (), Shapiro parameters (range “” and shape “”) and the rate of change of orbital period () due to gravitational wave emission. Details of all the above mentioned parameters and how they are used in the timing models are described very nicely in “Handbook of Pulsar Astronomy” by D. R. Lorimer and M. Kramer .
Post-Keplerian descriptions of the orbits are essential for “double neutron star” (DNS) binaries which include PSR J0737-3039(A, B), PSR B1534+12, PSR J1756-2251, PSR J1811-173, PSR J1906+0746, PSR B1913+16, PSR B2127+11C (confirmed DNSs), PSR J1518+4904, PSR B1820-11, PSR J1829+2456 (candidate DNSs - orbital parameters hint DNS nature but companion mass informations are insufficient to confirm). Only one of the DNSs, namely PSR J0737-3039 is a double pulsar system both the neutron stars emit radio pulses. This feature has made this system a unique astrophysical laboratory and we shall discuss about it in more details in subsections 2.1, 2.2 and 2.3. Moreover, there are three relativistic NS-WD systems - PSR J0751+1807, PSR J1757-5322, PSR J1141-6545 for which also post-Keplerian descriptions are unavoidable. Let us now try to understand how DNSs can help us in better understanding of the structure of neutron stars, physics of the pulsar magnetosphere and formation scenario (for the special case of the double pulsar system) of neutron stars.
2.1 Constraining Dense Matter Equations of State
After the first proposal of neutron star as the ultimate fate of stars (of mass ) by Baade and Zwicky in 1934  and the immediate prediction of mass and radius by Oppenheimer and Volkoff in 1939 , many people worked to understand the nature of matter at such high densities, resulting a large number of Equations of States (EsoS) including the EsoS for strange quark matter. Different EsoS result different mass-radius relations - see Fig 2 of  to realize how different the mass-radius curves can be. So the question arises, how to test the validity of those EsoS. Presently there are a number of efforts to constrain the EsoS through astronomical observations in X-rays. The usual approach is to determine the mass and the radius of the stars with the help of various observational features like gravitational redshifts (z) from spectral lines, cooling characteristics, kHz quasi-periodic oscillations (QPO) etc [9, 10, 11, 12, 13, 14]. But these methods are not foolproof, the value of z used in Özel’s  analysis of EXO 0748676 could not be reproduced afterwords . Moreover, to constrain EsoS from QPO observations, one needs to believe in a particular model of QPO which is again a subject of debate. An alternative method to constrain the dense matter EsoS might be the measurement of the moment of inertia of a neutron star in a relativistic binary, especially for a DNS system and to match this observationally determined value with the theoretically predicted one. The Equation of state (EoS) which gives the best agreement between observed and the theoretical values of the moment of inertia of the NS can be considered as the best description of NS structure.
The moment of inertia of a neutron star of known mass can be theoretically calculated for any particular EoS either using Hartle’s approximations  or by performing exact calculations [17, 18] for the rotating axisymmetric neutron star. It has been checked that the values of moment of inertia for a particular NS obtained by this two approaches do not differ more than 1% for a star like PSR J0737-3039A [17, 19, 18]. We found the moment of inertia of PSR J0737-3039A to be in the range for a set of EsoS having different stiffness .
Let us now discuss about observational determination of the moment of a neutron star by accurate measurement of the periastron advance () of the orbit of a relativistic binary pulsar (preferably a DNS) through timing analysis. The theoretical expression of including the effect of spin-orbit coupling is known to be as follow :
where the first term inside the square bracket is the 1PN term, the second term is the 2PN term, and the third term represents the spin-orbit coupling (SO). The parameters used in the above expression are defined as :
where is the gravitational constant, is the moment of inertia of the body (), is its angular velocity of rotation (here means modulo 2, 2+1=1) and . is the spin period of the star and is the orbital period of the binary. and . is the orbital angular momentum vector, is the spin vector, is the line of sight vector perpendicular to the sky and is the inclination of the orbit with respect to the sky implying that the angle between and is also .
In the expression of , the moment of inertia of both members of the DNS system are present (through and in Eqn. 1) which means that we have only one observable and two unknowns making Eqn. (1) unsolvable. But fortunately, is negligible for ms. So if ms, equation for periastron advance becomes
So DNS binaries having one fast neutron star ( ms) and another slow neutron star ( ms) are needed to extract the moment of inertia of the faster one by equating the observed with the prediction of Eqn. (6) using measured values of , , , , and . It is clear that the angle between and plays a vital to in the expression of . If , the angle between and becomes and we get . But if this is not the case the angle between and has a non-zero value, will be maximum if is parallel to the vector giving , for other orientation of , will be between to . So the knowledge of the angle between and is needed for determination of the moment of inertia of the faster member of a DNS system.
At present, it seems that the use of Eqn. (6) to determine the moment of inertia of a neutron star is possible only for the double pulsar system PSR J0737-3039 [21, 19], as for this system, timing analysis of the pulses from both the stars has already enabled us to determine , , , , and [22, 23]. The two pulsars which are named as PSR J0737-3039A and PSR J0737-3039B have masses and spin periods as , , ms and s with orbital parameters as , days, , . So it is clear that we can easily neglect the spin-orbit term due to pulsar B. Although the angles between the vectors and are not known which are needed to be determined first to measure the moment of inertia of the pulsar A accurately. Just for the sake of completeness, we wish to mention that for PSR J0737-3039, the value of is 1.2176 whereas the value of is 0.60876. The accuracies of other relevant parameters are also needed to be improved.
In general, it is clear that the term in Eqn. (6) will have higher value for (i) smaller value of , (ii) smaller value of and (iii) higher value of . Unfortunately, the eccentricity of PSR J0737-3039 is small in comparison to other DNSs. Future discovery of any double pulsar system satisfying these three conditions better than PSR J0737-3039 will be useful for the purpose of moment of inertia measurement (in addition to the necessary condition ms).
2.2 Interplay between Two Members of the Double Pulsar
Among many interesting features observed from the double pulsar system PSR J0737-3039, one is the variation of the flux density of the slower pulsar B at its different orbital phases due to the distortion of B’s magnetosphere by A’s wind [23, 24]. On the other hand, the flux density of pulsar A is constant for most of its orbit except s near its superior conjuction when the radio beam from the pulsar A is eclipsed by the magnetosphere of the pulsar B . The eclipse is asymmetric, with the flux density decreasing more slowly during ingress than it increases during egress. The duration of the eclipse is very mildly dependent on the observing frequencies, with the eclipse lasting slightly longer at lower frequencies . Theoretical models explain these eclipse properties in the context of synchrotron absorption of A’s radio emission by the plasma in the closed field line region of B’s magnetosphere [26, 27, 28]. By modeling the eclipse, the precession of B’s spin axis about the orbital momentum vector has been found to match with the prediction of the general theory of relativity within the measurement uncertainty . It is noteworthy to mention here that till date, the eclipse has been studied in the observing frequency range 427-1400 MHz. The study of PSR J0737-3039 (both A and B) at a much lower frequency of 326 MHz is being performed using Ooty Radio telescope (RAC-TIFR, Ooty, India). The details of eclipse characteristic at this low frequency will be useful to test the validity of the eclipse models and better understanding of the pulsar magnetosphere.
2.3 Formation of the Double Pulsar
Using observed parameters of the double pulsar system, different formation models were proposed [30, 31, 32]. Initially, it was proposed that this system has originated from a close helium star plus a neutron star binary in which at the onset of the evolution the helium star had a mass in the range 4.0 6.5 and an orbital period in the range 0.1 0.2 day . In this model, a kick velocity in the range 70 230 must have been imparted to the second neutron star at its birth. Afterwords an electron capture supernova scenario was proposed . At the same time, Piran & Shaviv  estimated that it is most likely that the progenitor of the pulsar B had a mass of with a low kick velocity 30 . On the other hand, a study based on the population synthesis method  favored the standard formation scenario with reasonably high kick velocities (70-180 ). Later, it has been suggested  that the pulsar B is a strange star rather than a neutron star as the baryon loss is smaller for the first case supporting the models with small kick.
3 Binary Pulsars in Globular Clusters
Globular clusters (GCs) are dense spherical collection of stars orbiting around the galactic center containing a significant number of binary stars. The first radio pulsar discovered in a globular cluster is PSR B1821-24A in M28  in 1987, 20 years after the discovery of first pulsar by Jocelyn Bell. This pulsar is at present an isolated pulsar but its very low spin period (3 millisecond) hints to a binary origin and spin-up due to accretion . How old and evolved stellar systems like GCs can contain apparently young, active objects like pulsars have led to suggestions, based on lifetime constraints from orbital decay by gravitational radiation , that some of these systems have been formed recently, or are forming even today. Interaction between binary stars and single stars in GCs is believed to be an important dynamical process as the resulting binary systems may provide a substantial source of energy for the host GC, as the binding energies of a few, very close binaries (like neutron star binaries) can approach that of a moderately massive GC [36, 37]. Observable parameters of the pulsars and their binary companions in GCs, such as spin, orbital period and eccentricity, projected radial position in the cluster, companion mass and their distributions provide a tracer of the past history of dynamical interactions of the binary NSs in individual GCs. These parameters can even provide a valuable test-bed to examine the theoretical scenarios of formation and evolution of recycled pulsars.
3.1 Archival Data of Binary Pulsars in Globular Clusters
Till today, a total of 140 pulsars in 26 GCs
In Fig 1, we plot Vs for these 73 binary pulsars in GCs with known orbital solutions. A logarithmic scale in both eccentricity and orbital periods are chosen, as the enormous range of both variables and the regions occupied by observed pulsars are less obvious in linear scales. Downward arrows have been drawn for the systems for which only upper limits of eccentricities (see table 2) are known. Non-millisecond pulsars are enclosed by s and the DNS system PSR B212711C in M15 is enclosed by a green . In the present work, we define millisecond pulsars as the pulsars having ms, we don’t go for other definitions of millisecond pulsars based on both and (see  for one such definition) as for GC binaries, the values of are effected significantly by the cluster potential.
Observed binaries in GCs can be categorized in three groups : (I) large eccentricity pulsars (), 21 pulsars (22 if we include M30B), (II) moderate eccentricity pulsars (), 20 pulsars and (III) small eccentricity pulsars (), 32 pulsars. In the database, orbital eccentricities of group (III) pulsars have been listed as zero, and in this logarithmic plot we assign them an arbitrarily small value of .
In addition we plot on the right of Fig 1, dashed lines corresponding to the freeze-out eccentricity - orbital period relation predicted by Phinney  on the basis of the fluctuation dissipation theorem for convective eddies in the erstwhile red giant envelopes surrounding the white-dwarf cores which end up as companions of the neutron stars. This relation is obtained from the following expression:
where is the orbital angular rotation frequency, is the semi-major axis of the binary, is the orbital eccentricity, L is the luminosity of the companion, is the radius of the companion, is the envelop mass of the companion, is the pulsar mass, is the companion mass and is the reduced mass of the system. Phinney used when the mass transfer ceases and freezes. Then the above equation reduces to where is a function of the effective temperature (), and . is in days. Phinney calculated (see Eqn. 7.35 of ) using average values of and for population I stars ()  which gives the upper dashed green line in the right side of Fig 1. However for globular clusters, population II stars are more appropriate. So we recalculated for population II red giant stars  and obtained . This line is plotted as the lower dashed green line in Fig. 1.
Earlier, it had been found that while the low mass X-ray binaries (LMXBs) are predominantly found in very dense GCs, the radio pulsars tend to be found more evenly distributed among low and high density clusters . We also don’t find any obvious correlation of the number of binary radio pulsars with any GC parameters like , , or (where is the number density () of single stars in units of and is the velocity dispersion () in units of 10 in GCs and is the distance of the GCs from sun in kpc see Table 1). But if we take only GCs with more than 3 observed binary pulsars, we find the number of binary pulsars to increase with , but there are only five such GCs (Table 1).
A comparison of binary pulsars in globular clusters and binary pulsars in the galactic disk
There are total 83 binary pulsars in the galactic disk (and 1652 isolated pulsars) among which orbital period and/or eccentricity values are unknown for four binaries. So we exclude them from our study and we consider the rest 78 binary pulsars. In the database, eccentricities of three binaries have assigned the value “0” for which we put an arbitrary small value (as in the case of GC binary pulsars).
In Fig 2, we plot Vs for all binary pulsars with known orbital solutions. Blue s are for GC binaries and red s are for disk binaries. Three disk binaries with days have not been plotted here as we have selected the range of up to 1000 days in this plot. Here we take upper values of eccentricities as the actual values. DNSs are enclosed by green s. Non-millisecond pulsars are enclosed by s. It is believed that millisecond pulsars in binary systems formed from mass and angular momentum transfer due to Roche lobe overflow and the resultant tidal effects should appear in . Since many highly eccentric binary millisecond pulsars are found in GCs, this indicates that stellar interactions are important for inducing higher eccentricities which will be discussed later in this section. On the other hand, PSR J is the only one eccentric millisecond pulsar in the disk and there are lots of works in the recent past to explain the eccentricity of this system [43, 44, 45]. We shall discuss about this system in somewhat details in section 4.
Observational selection effects may be influencing the distribution seen in this figure for the two sets , the most important selection effect being operative towards the left of the diagram, namely, it is more difficult to detect millisecond pulsars in very short orbital period and highly eccentric binaries; distance to the pulsars also is an important selection effect, since at large distance only the brightest pulsars can be observed. Nevertheless, it is clear that there is a large abundance of pulsars with long orbital periods and intermediate eccentricity near the eddy “fluctuation dissipation” lines among the galactic disk pulsars which are absent in the GC pulsar population, where the important selection effects are unlikely to play a major role. The observed preponderance of DNSs in the disk with the very short orbital periods and large eccentricities may be related to selection effects, but there is at least one DNS in GC - PSR J2127+1105C in M15 in the same phase-space. The disk DNSs are by and large at small distances while M15 and other host GCs are usually far away.
In the top panel of Fig. 3 we compare the distribution of s for the GC and the disk binaries and found that there are more low mass companions for GC binaries. In the bottom panel of Fig. 3 we compare the distribution of s for the GC and the disk binaries and found that in GCs, most are millisecond pulsars while for disk pulsars, are more uniformly distributed (not completely, more pulsars with lower ) over a wide range.
In the left panel of Fig. 4, we plot , and for all binaries while in the right panel we plot , and . In both cases, blue points represent GC binaries and red points represent disk binaries. Three disk binaries with (see table 3) have been omitted from the plot in the left panel, even then we can see that companion masses of disk binaries are usually higher than those of GC binaries. Nearly circular binary pulsars have low mass companions.
It was proposed that orbital parameters of binary pulsars with white-dwarf companions satisfy the following relation :
In the upper panel of Fig. 5, we plot against for binary pulsars in GCs and also in the galactic disk. DNS binaries and the binaries with (Chandrasekhar limit) have been excluded. GC binaries are shown with blue s and disk binaries are shown with red s. The points correspond to median values of taking while the upper and the lower limits of s correspond to and respectively. Theoretical prediction following Eqn. (8) is the middle solid green line. The upper and lower lines are obtained by multiplying the right hand side of Eqn. (5) by 2.5 (see ). Note that the model used to derive Eqn. (8) is not valid when , though we have plotted the binaries with (located left of the vertical black dotted line) just to see how is related with for these systems. For recycled (millisecond) pulsars, where significant amount of mass had been transferred to the pulsars in past, one should expect lower values of and . for these cases the relation becomes
The deviation of GC binary pulsars from the theoretical prediction of Eqn. (8) is well understood as their eccentricities have been effected by stellar interactions. So for the millisecond pulsar binaries in GCs (which naturally exclude the DNS in M15), we plot Eqn. (9) in the lower panel of Fig. (5). Even then we found significant discrepancy between the theoretical prediction and the observed values. Moreover, for binary pulsars in the disk (upper panel of Fig. 5) where eccentricities and orbital periods should bear the clear signature of stellar evolutions, we also found disagreement between the model and the observation.
We summarize the relevant properties of the host GCs and in Table (1). The properties of binary pulsars in globular clusters and in the galactic disk are given in tables (2) and (3) respectively. Interested readers are requested to see original sources for better precision of the given parameters and for the values of other parameters.
Through a combination of two-body tidal captures and exchange interactions where an isolated neutron star replaces a normal star in a binary, numerous neutron star binaries are formed in GCs . This explains why the ratio of binary pulsars to isolated pulsars is so high in GCs in comparison to that in the galactic disk. Tidal capture leads to tight binaries whereas exchange usually leads to wider binaries (see Eqns. 14). If the mass transfer in the tidally captured neutron star binary happen in a stable manner, it can spin up the neutron star to a millisecond spin period. Thus the observed preponderance of millisecond pulsars in GC binaries can be attributed to the tidal capture process (details will be discussed in the next subsection). We explain the eccentricities of millisecond pulsar binaries in GCs as a result of different types of interactions of single stars with pulsar binaries using present properties of the GCs. Some of the results can be found in our earlier papers [48, 49]. For the sake of simplicity, we don’t go for a full population synthesis study of the evolution of the globular clusters including all kinds of stellar interactions binary-binary, binary-single star and single star - single star interactions. We consider only binary-single star interactions faced by binary pulsars preceded by a discussion on the formation of binary pulsars through tidal capture.
3.2 Orbital Parameter Changes Due to Stellar Interactions and Gravitational Radiation
Globular clusters contain a large number of LMXBs compared to the galactic disk; this had led to the suggestion  that a binary is formed by the tidal capture of a non-compact star by a neutron star in the dense stellar environment of the GC cores. If mass and angular momentum transfer from the companion to the neutron star took place in a stable manner, this could lead to recycled pulsars in binaries or as single millisecond pulsars [51, 52, 53]. However large energies may be deposited in tides after the tidal capture of a neutron star by a low mass MS star in a GC. The resultant structural readjustments of the star in response to the dissipation of the modes could be very significant in stars with either convective or radiative damping zones [54, 55] and the companion star can undergo a size “inflation” due to its high tidal luminosity which may exceed that induced by nuclear reactions in the core, by a large factor. Efficiency of viscous dissipation and orbit evolution is crucial to the subsequent evolution of the system as viscosity regulates the growth of oscillations and the degree to which the extended star is bloated and shed. The evolution of the system could lead to either a merger (leading to a Thorne Zytkow object) or a neutron star surrounded by a massive disk comprising of the stellar debris of its erstwhile companion or, in other cases, an LMXB or even a detached binary following envelope ejection. Tidal oscillation energy can return to the orbit (i.e. the tide orbit coupling is a dynamical effect) thus affecting the orbit evolution or the extent of dissipative heating in the less dense companion. The evolutionary transition from initial tidal dynamics to any likely final, quiescent binary system is thus regulated by viscosity. Tide orbit coupling can even lead to chaotic evolution of the orbit in some cases  but in the presence of dissipation with non-liner damping, the chaotic phase may last only for about 200 yrs after which the binary may undergo a periodic circularization and after about 5 Myrs finish circularizing . Of the stars that do not directly lead to mergers, a roughly equal fraction of the encounters lead to binaries that either become unbound as a result of de-excitation or heating from other stars in the vicinity or they are scattered into orbits with large pericenters (compared to the size of the non-compact star) due to angular momentum transfer from other stars . Thus “intermediate pericenter” encounters can lead to widened orbit (even 100 days orbital period in an eccentric orbit) in tidal encounters involving main sequence stars and neutron stars, which in the standard scenario are attributed to encounters between giants and neutron stars. The complexities of the dynamics of tidal capture and subsequent evolution are manifold and attempts to numerically simulate collisions of neutron stars with red giants have been made with 3D Smoothed Particle Hydrodynamics (SPH) codes [59, 60]. We note in this context that Podsiadlowski  have stated “it is not only premature to rule out tidal capture as a formation scenario for LMXBs, but that LMXBs in globular clusters with well determined orbital periods actually provide observational evidence in its favor”. The different formation channels of pulsars in binaries in GCs and the birthrate problems of millisecond pulsar binaries vs LMXBs were considered in [61, 62].
Until the early eighties there was no evidence of a substantial population of primordial binaries in any globular clusters, and it was even thought that GCs are significantly deficient in binaries compared to a younger galactic population . Theoretical modeling of GCs was often started off as if all stars were singles. During the 1980s several observational techniques began to yield a rich population of binary stars in GCs [64, 65, 66]. When the local binary fraction is substantial, the single star - binary interaction can exceed the encounter rate between single stars by a large factor . The existence of a significant population of primordial binary population in GCs indicated that three body processes have to be accounted for in any dynamical study of binaries involving compact stars. An encounter between an isolated star and a binary may lead to a change of state of the latter, e.g.: i) the original binary may undergo a change of eccentricity and orbital period but otherwise remain intact – a “preservation” or a “fly-by” process; ii) a member of the binary may be exchanged with the incoming field star, forming a new binary – an “exchange” process; iii) two of the stars may collide and merge into a single object, and may or may not remain bound to the third star - a “merger” process; iv) three of the stars may collide and merge into a single object - a “triple merger” process; v) all three stars may become unbound - an “ionization” process. The ionization is always a “prompt” process, whereas the others can be either prompt or “resonant” process. In resonant processes, the stars form a temporarily bound triple system, with a negative total energy but which decays into an escaping star and a binary, typically after 10-100 crossing times.
The value of the binding energy of the binary and the velocity of the incoming star determine the type of interaction the binary will encounter. The critical velocity of the incoming star, for which the total energy (kinetic plus potential) of the three body system is zero, can be defined as
where and are the masses of the binary members, is the mass of the incoming star, is the semi-major axis of the binary and G is the gravitational constant.
In case of a binary-single star interaction, semi major axis of the final binary () is related to the semi major axis of the initial binary () as follows  :
where and are masses of the members of the initial binary, and are masses of the members of the final binary. is the fractional change of binary binding energy . The binding energies of the initial and final binaries are
For fly-by, , giving
For exchange, , (star 2 is being replaced by star 3) giving
For merger, , (star 2 is being merged with star 3) giving
The actual value of is not known. Putting simplifies equations (11), (13), (14) and (15) giving for fly-by interactions. Similar expressions can be derived for other cases star 1 is being replaced by star 3 or star 1 being merged with star 3. But in the present study, we always assume star 1 to be the neutron star, so these processes are irrelevant. Merger interactions which do not result in binary systems are also irrelevant for this work.
Interactions between binaries and single stars can enhance the eccentricity of the binary orbit and this will be discussed in detail in the subsequent sections. On the other hand, binary pulsars emit gravitational waves which reduce both the sizes and eccentricities of the orbit leading to mergers of the binary members on a timescale of which can be calculated as follows  :
which is almost same as the expression of calculated from .
In this work, we use eq. (16) to calculate .
Figs. 1, 2 and Table 2 show that most of the binary pulsars in GCs are millisecond pulsars (69 pulsars out of 73 have ms, one has ms and the other three having values as 80.34 ms, 111.61 ms and 1 s). Theoretically one expects that spun-up, millisecond pulsars in binary systems formed from mass and angular momentum transfer due to Roche lobe overflow and the resultant tidal effects should appear in low eccentricity orbits . Since many highly eccentric binary millisecond pulsars are found in globular clusters (see the right panel of Fig 4), this indicates that stellar interactions are important in GCs for inducing higher eccentricities. Rasio and Heggie [69, 70] studied the change of orbital eccentricity ( where is the initial eccentricity and is the final eccentricity which is observable at present) of an initially circular binary () following a distant encounter with a third star in a parabolic orbit. They used secular perturbation theory, i.e. averaging over the orbital motion of the binary for sufficiently large values of the pericenter distance where the encounter is quasi-adiabatic and used non secular perturbation theory for smaller values of where the encounter is non-adiabatic. In the first case varies as a power law with and in the second case varies exponentially with . The power law dominates for wheras the exponential law dominates for . The timescales for these processes are  (where we write in place of for the sake of simplicity) :
where is in days and is in years. The values of are different for different GCs (table 1) and we grouped them into six different groups according to their values of and calculated for mean values of for each group. The unit of is which we will not mention explicitly anymore.
Cluster binary pulsars in the plane with contours
of yrs (solid lines) and yrs (dashed lines)
for each group are shown in Fig 6.
Contours of yrs (solid lines) and yrs (dashed lines) for binaries with and or are also shown. Pulsars with projected positions inside the cluster cores are marked with , those outside the cluster cores are marked with and the pulsars with unknown position offsets are marked with . Solid red line corresponding to years for is outside the range of axes used in our plots (away from the upper left corner). Pulsars in a particular group are marked with the same color as of the contours for that group. The upper limits of eccentricities are taken as the actual values of eccentricities. If a pulsar is located on the upper left half of the corresponding years line in the region where years, then its eccentricity cannot be due to fly-by interactions. There are two such binaries,
one is PSR B1718-19
We have seen that there are many zero eccentricity pulsars which lie in the region where years. So it is of interest to understand why they still appear in circular orbits despite the possibility of fly-by interactions. We will discuss about this point in subsection 3.4.
Exchange and merger interactions
Since exchange and merger interactions have not been amenable to
analytic treatment one needs to use numerical techniques. We performed numerical scattering experiments using the STARLAB
where is the cross section and is the critical velocity as defined earlier (eq. 10). So from the reported values of , one can calculate and then the interaction time scale as where is the number density of single stars. We took different sets of stellar parameters () to understand the effect of each of them and for each set of stellar parameters, we varied over a wide range in such a way that the initial orbital periods (, calculated using Kepler’s law) are in the range of days.
For each set of parameters (), we take the maximum trial density () as 5000 which is uniformly distributed in the impact parameter () over the range where simply corresponds to a periastron separation of . The impact parameter range is then systematically expanded to cover successive annuli of outer radii with trials each until no interesting interaction takes place in the outermost zone . Thus the total number of trials become significantly large. As an example, for the parameters we had a sample total of 15570 scatterings. Out of these 10960 were fly-bys, 3722 exchanges, 982 two mergers and 6 three-mergers.
We performed simulations over widely different values of 11.76 , 13.15 , 7.79 and 3.29 chosen in such a way that almost the entire range of ( respectively) has been covered. The variation in is also significant. But we mainly concentrate on the runs for as in Ter 5, the host of the largest number of known binary radio pulsars.
In Fig. (7) we plot the variation of for exchange (), merger () and ionization (, whenever significant) with and check the dependencies on different parameters - (i) In the first figure, we set , , and take two values of 13.16 (red) and 3.27 (green) which are the highest and the lowest values of among our choice (see Table 1). (ii) In the second figure, we set , and take different values of (red), (green) and (blue). Ionization starts for at days. (iii) In the third figure, we set , and take different values of (red), (green). (iv) In the fourth figure, we set , and take the star to be either an MS (red) or a WD (green). The results can be summarized in a tabular form :
|change of variable||result|
|decrease in the value of||and|
|increase slightly at higher|
|decrease in the value of||and|
|increase significantly at lower .|
|Ionization is significant at very low|
|values of and higher values of|
|increase in the value of||and|
|nature of the incoming star||is slightly higher at lower when the|
|(of mass )||incoming star is a WD rather than a MS.|
|is significantly lower when the|
|incoming star is a WD rather than a MS.|
Eqn. (23) reveals that for any fixed set of stellar parameters () which becomes after using the expression of as in Eqn. (10). So the interaction timescale can be written as . As it is clear from the top-left panel of Fig. (7), that for a fixed set of , the variation of is negligible for any given . So we can write . Remember, the same relation has been used for in Eqns. (22a, 22b).
STARLAB also gives the properties of final states, eccentricities and semi major axes of the final binaries for each set of inputs (). So for each value of , we get a number of values of and . For each value of , we calculated percentile, median and percentile values of as well as the values from analytic expressions given in Eqns. (14) and (15) putting . As all these values are very close to each other, we use the relation corresponding to throughout (eq. 11). As an example, in Fig. 8, we have plotted Vs as reported by STARLAB (scatter plots) for both merger and exchange with , , . The line for 25 percentile, median, 75 percentile along with the analytic relation of with (eq. 11) are shown in the same plot.
In each panel of Figs. (9) and (10), we plot of the initial binary (comprising of stars and ) along the top x-axis and along the bottom x-axis. The left y-axis gives the final eccentricities while the right y-axis gives the time scales of the interactions. Purple points in the left panels are for exchange interactions while the green points in the right panels are for the merger interactions. Interaction time scales are plotted with black s. The vertical orange lines form the boundaries of the allowed orbital period regions where interaction time scales are less than years. It is clear from the scatter plots (Figs. 9, 10) that the final binaries will most probably have if they undergo either exchange or merger events. The observed binary pulsars with in Ter 5 can be found in Table 2 where the companion masses () are for inclination angle . These binaries are also shown in all panels of Figs. 9 and 10 (red in color, named in lower panels of Fig. 10). In Fig 9, , and we vary from to and . We have also performed the similar kind of numerical simulations by taking and found that the outcomes of the simulations i.e. the appearance of the scatter plots (Fig 10) does not change much in comparison to the change with the variation in the value of (for fixed ).
In Section 3.2.1, we have seen that PSR B1638+36B (in M13) has such a combination of that it can not be explained by fly-by interactions. Moreover, the scatter plots (Figs. 9, 10) also reveal that its eccentricity can not be caused by exchange or mergers as these interactions imparts eccentricities larger than whereas the upper limit of the eccentricity of PSR B1638+36B is 0.001. One possibility is that the true eccentricity of this system is significantly lower than the upper limit so that it falls in the yrs region. Another possibility is that it could have undergone eccentricity pumping of the inner orbit due to the presence of a third star in a wide outer orbit, i.e. in a hierarchical triplet like the PSR B1620-26 in M4 system [72, 73, 74, 75].
If an observed eccentric binary pulsar lies in the region where time scale for a particular interaction is greater than years, then that interaction can not be responsible for its eccentricity. Moreover, an exchange interaction to be the origin of the eccentricity of a particular pulsar, we should have and for merger