Tidal interaction in compact binaries: a postNewtonian affine framework
Abstract
We develop a semianalytical approach, based on the postNewtonian expansion and on the affine approximation, to model the tidal deformation of neutron stars in the coalescence of black holeneutron star or neutron starneutron star binaries. Our equations describe, in a unified framework, both the system orbital evolution, and the neutron star deformations. These are driven by the tidal tensor, which we expand at postNewtonian order, including spin terms. We test the theoretical framework by simulating black holeneutron star coalescence up to the onset of mass shedding, which we determine by comparing the shape of the star with the Roche lobe. We validate our approach by comparing our results with those of fully relativistic, numerical simulations.
pacs:
04.25.Nx, 04.30.Dg, 04.25.dkI Introduction
Coalescing binaries composed of neutron stars (NS) and/or black holes (BH), are among the most promising sources of gravitational waves (GWs) to be detected by gravitational wave interferometers like Virgo and LIGO virgoligo (). These systems are also interesting since they are thought to be related to short gammaray burst NPP09 ().
The process of coalescence has been studied mainly using postNewtonian (PN) and fully relativistic, numerical simulations. PN expansion B06 () has the advantage of providing a semianalytic description of the evolution of the system, but it is poorly convergent in the strong field limit; therefore it is appropriate to study the inspiral phase only. These limitations have been overcome by PN resummed formulations like EOB EOB (), but these formulations are not able to describe the dynamical features of the stellar deformation. Moreover, in the standard PN expansion the compact objects are treated as pointlike up to the 4.5 (included) postNewtonian order. Finite size effects are, formally, of order 5PN; however their contribution is larger than what a naive counting of PN orders may suggest MW04 (). Tidal deformations have recently been included in the PN framework through the “Love number” approach FH08 (); H08 (), which assumes that the tidal field is proportional to the quadrupole momentum (see below). The effects of the tidal deformation on the orbital motion have been studied in VFH11 (); EOBtidal ().
Fully relativistic codes are the most powerful tool to investigate the latest phases of the inspiral and merger (see F09 () for a review on the subject). They are, however, not exempt from drawbacks: their computational cost is high, therefore the parameter space cannot be explored at large; furthermore, initial data solvers are still unable to provide accurate initial data for binaries with nonaligned spins, and may introduce spurious numerical effects which, if not appropriately cured, affect the subsequent evolution of the system. These problems are of particular relevance in BHNS binaries, where the lack of symmetry makes more difficult to follow the entire process of coalescence by fully relativistic simulations. For these reasons, the process has been studied in the literature using some simplifying assumptions, or for a restricted set of parameters. For instance, in TBFS07 (); S09q (); FGP09 (); FGP10 () the inspiral is modeled as a sequence of quasiequilibrium circular orbits with decreasing radius; in TBFS07 (); S09q () the process is studied by fully relativistic simulations, whereas FGP09 (); FGP10 () use the affine approach (see below). In F06a (); F06b () Einstein’s equations are evolved assuming that the black hole is non rotating, and for large values of the mass ratio , whereas in S07 (); S09 (); S10 (); E09 () takes values ; in E09 (); S11 () the black hole is assumed to rotate with spin parallel to the orbital angular momentum, and different values have been considered. For a recent review on fully relativistic simulations of BHNS binaries see ST11 () (the literature on NSNS coalescing binaries is much more extended, and we do not report it here).
In this paper we develop a semianalytic approach to study BHNS and NSNS coalescence, by merging two different frameworks: the PN approach, which accurately describes the system orbital motion, and the affine model CL85 (); LM85 (); WL00 (); CFS06 (); FGP09 (); FGP10 (), which describes the stellar deformations induced by the tidal field. To this aim, we compute the tidal tensor associated to the PN metric of a twobody system, defined in terms of the PN Riemann tensor and of the local tetrad of the deformed body
(1) 
up to . This tensor was derived with a different approach in DSX92 (); RF05 (); VF09 () up to ; our expression coincides with that of VF09 () and also includes terms, associated to the spins of the compact objects.
In the affine model the NS is described as a deformable ellipsoid, subject to its selfgravity, to internal pressure forces and to the tidal field of the companion. In the original formulation of this approach, the NS structure was considered at a Newtonian order, assuming a polytropic equation of state (EOS) CL85 (); LM85 (); WL00 (); CFS06 (). A first improvement was introduced in FGP09 (); FGP10 (), where general relativity was taken into account in the description of the stellar structure, and nonpolytropic EOSs were considered. This approach was used to study quasiequilibrium configuration sequences of BHNS systems, in order to estimate the critical distance at which the NS is disrupted by the tidal interaction FGP09 (), to determine the corresponding cutoff frequency in the emitted gravitational wave signal FGP10 (); PROR11 (), and to estimate the mass of the torus which forms after the NS is disrupted PTR11 ().
In this paper we further improve the affine model introducing a more accurate description of the orbital motion and of the tidal interaction. Our approach differs from existing work on NS tidal deformation in compact binaries in the following aspects.

In WL00 (); CFS06 (); FGP09 (); FGP10 (); OK96 (); PTR11 () the affine model was used assuming that the NS follows a timelike geodesic of Kerr’s spacetime; this approximation fails when the mass ratio is low. In addition, timedilation factors were neglected.
Furthermore, most works employing the affine model WL00 (); CFS06 (); FGP09 (); FGP10 () do not evolve the dynamical equations of stellar deformation. Rather, they find, at each value of the orbital radius, the corresponding stationary configuration describing the deformed star. In OK96 (); PTR11 () the dynamical equations were solved; however, while the orbital evolution was described in PN coordinates, the BH tidal field was expressed in the BoyerLindquist coordinates of the Kerr metric describing a single BH; neglecting the difference between these coordinate systems yields a loss of accuracy of the model.
These problems are solved in the fully consistent approach presented in this paper, where the BHNS or NSNS systems are described by a twobody PN metric, which holds for any value of the massratio . The tidal tensor itself is expressed in the PN coordinates, and the proper time of each compact body is expressed in terms of the PN time coordinate, through the appropriate Lorentz factor. Our approach is valid up to the onset of mass shedding, which occurs when the deformed star crosses the Roche lobe; after that, it can no longer be applied, since the assumption that the star is a deformed ellipsoid is significantly violated. We describe the orbital motion of the compact objects using the PN equations for pointlike objects, with nexttoleadingorder tidal corrections; the NS tidal deformation is driven by the tidal tensor of the PN metric. The dynamical equations are a system of (nonlinear) ordinary differential equations in time.

NS tidal deformations have also been studied in a series of paper FH08 (); H08 (); HLLR10 (); DN09 (); BP09 (); DN10 (), where the deformation properties have been encoded in a set of numbers, the Love numbers, which relate the quadrupole tensor (or, more generally, the multipole moments of the star) to the tidal tensor. This approach is grounded on the adiabatic approximation, i.e., on the assumption that the orbital evolution timescale is much larger than the timescale needed for the star to set into a stationary configuration. In this approximation, the quadrupole tensor is proportional to the tidal tensor:
(2) with constant. The Love number can be computed by studying the response of a single star to an external tidal tensor FH08 (); H08 (); DN09 (); BP09 (). This model has been employed to determine the effect of tidal deformation on the orbital motion of a NS in a binary system HLLR10 (); VF09 (); DN10 (); VFH11 ().
We also compute the Love number (see Section III), without assuming the adiabatic approximation: the stellar deformation is found by solving dynamical equations.
To test the accuracy of our approach, we compare the results with the existing literature on BHNS binaries. As a preliminary check, we verify that our PN description of the orbital motion accurately reproduces the fully relativistic results S09 (); E09 (); privcomm (). Then, we verify that the onset of mass shedding we determine, is consistent with the results of fully relativistic simulations S09 (); E09 (); privcomm (). Finally, we check that the stellar deformations predicted by our model are consistent with existing computations of the Love number FH08 (); H08 (); DN09 (); BP09 ().
This paper focuses mainly on the theoretical framework, and on its validation by comparison with the existing literature, where available. The tool we develop will be used in future works to study the dynamics of compact binaries.
Ii The model
We use notations and conventions introduced in FGP09 (), where the affine model (partially improved with respect to the original formulation CL85 (); LM85 ()) is widely discussed. In the following are the masses of the two compact objects; we shall consider the tidal deformation of the NS with mass and radius ; the companion, with mass , can either be a BH or a NS. Furthermore, we define and , and the mass ratio .
ii.1 Improved affine model
The basic assumption of the affine model is that the NS, deformed by the tidal field, maintain an ellipsoidal shape; it is an type Riemann ellipsoid, i.e., its spin and vorticity are parallel, and their ratio is constant EFE (). The deformation equations are written in the star principal frame, i.e., the frame comoving with the star, whose axes coincide with the ellipsoid principal axes. In what follows, is the axis which points toward the companion; and are the axes orthogonal to , with lying in the orbital plane; the indices 1,2,3 label the corresponding directions. Surfaces of constant density inside the star form selfsimilar ellipsoids and the velocity of a fluid element is a linear function of the coordinates in the principal frame. Under these assumptions, the infinite degrees of freedom of the stellar fluid motion can be reduced to five CL85 (); LM85 (); WL00 (), and are associated to dynamical variables governed by a set of nonlinear differential equations, which describe the evolution of the stellar deformation. These variables are the three principal axes of the ellipsoid and two angles, , , defined as
(3) 
where is the NS proper time, and is the ellipsoid angular velocity, measured in the parallel transported frame associated with the star center of mass. is the angle between the principal frame and the parallel transported frame. is defined as follows:
(4) 
where is the vorticity along the axis in the principal frame. The NS internal dynamics is described in terms of the Lagrangian
(5) 
where is the fluid kinetic energy, is the internal energy and is the star selfgravity. In the original approach introduced by Carter and Luminet, these are defined in a Newtonian framework
(6)  
(7)  
(8) 
where is the baryonic mass, the mass density, the Newtonian energy density, the gravitational potential; all these quantities are solutions of the Newtonian equations of stellar structure. Furthermore, and satisfies the virial theorem, which states that, in the spherical configuration, , where
(9) 
The variation of the Lagrangian (5) gives the equations of motion for the five dynamical variables .
In CL85 (); LM85 () it was shown that, under the affine hypothesis, the integrals in Eqns. (68) and their variations, can be expressed in terms of integrals on the fluid variables computed for the spherical configuration (), and of functions of the dynamical variables. With this simplification, the equations of motion can easily be found. In the following, a superscript hat will denote quantities computed for the spherical star. and expressed in terms of the “hatted” quantities and of the dynamical variables are
(10)  
(11) 
where is the scalar quadrupole moment
(12) 
(the subscript means that the integration is performed on the spherical star) and is the selfgravity potential of the spherical star.
The procedure to make explicit the dependence of the internal energy on the dynamical variables is more subtle. The internal energy variation can be written as
(13) 
The pressure integral given by Eq. (9) can not be factorized in a spherical integral and a function of the axes; however, it can be expressed as
(14) 
where is the fluid equation of state, and is the rescaled mass density
(15) 
with mass density in the spherical configuration. When , the pressure integral reduces to the spherical pressure integral .
A first improvement to this approach was introduced in FGP09 (), where the Newtonian description of the NS equilibrium configuration was replaced by the relativistic equations of stellar structure (TOV)
(16) 
Here is the relativistic massenergy density in the spherical configuration, is the radial coordinate in a Schwarzschild frame associated to the non rotating NS, and is the gravitational mass enclosed in a sphere of radius . We remark that this is a major change, since the relativistic radius is smaller than the Newtonian radius by . We also remark that the Schwarzschild coordinate is different from the radial coordinate in the Newtonian frame . The selfgravity integral was changed accordingly, as
(17) 
where is an effective relativistic gravitational potential of the spherical star, defined in terms of the TOV equations as follows
(18) 
With this definition the virial theorem
(19) 
is still satisfied, with solution of the TOV equations (16) and baryon mass density. As shown in Section II.5, some terms in the dynamical equations cancel in the spherical limit, only if the virial theorem is satisfied. A nonexact cancellation of these terms would lead to strong instabilities.
A further improvement, which we introduce in this paper, consists in a careful treatment of the coordinate frames. To describe the integrals in the spherical configuration, the relevant coordinate systems are: (i) the Schwarzschild frame, with radial coordinate , in which the TOV equations (16) are expressed; (ii) the Newtonian frame for a spherical star , which we now replace with the corresponding PN postNewtonian frame B06 (), with isotropic radial coordinate and metric (for a single star)
(20) 
where . Following Chan65 () the transformation between the postNewtonian isotropic radial coordinate and the Schwarzschild coordinate inside the star is given by
(21) 
The scalar quantity can be expressed, in the Schwarzschild frame, in terms of the corresponding spatial threemetric :
(22)  
The integrand in the quadrupole moment (12) is expressed in the postNewtonian coordinates, i.e., it is . The integrals then take the form
ii.2 The postNewtonian metric
To derive the equations describing the orbital motion of the binary and the tidal tensor, we shall use a PN metric written in harmonic coordinates ():
(24)  
(25)  
(26) 
where the potentials ,,, are defined in terms of retarded integrals over the source densities BFP98 (); FBA06 (). We stress that the potential appearing in the metric of the twobody system (24)(26), is different from the potential in Eq. (20), which is the metric of a single star. Since these potentials are written as expansions of powers of , in the following we shall identify the order of expansion with a superscript index. Thus, defines the scalar potential of order in , is the term and so on.
ii.3 The orbital motion
Following BAF11 (), we assume that the orbit evolves as a slow adiabatic inspiral of a quasicircular orbit, i.e., the energy lost through gravitational waves is balanced by a change of the total binding energy of the system
(27) 
where and the GW flux are expressed in terms of the PN variable
(28) 
being the orbital frequency. Eq.(27) yields
(29) 
We neglect the orbital eccentricity because, due to gravitational wave emission, the orbit circularizes well before the latest stages of the inspiral which we are studying P64 (). We use the approach named “Taylor T4 approximant” DIST01 (); SOA10 (), in which the righthand side of eq.(29) is expanded to PN order including spin terms. We also include the effects of the NS tidal deformation on the orbital motion, up to nexttoleadingorder VFH11 (). The orbital phase and the orbital frequency are computed by numerically integrating the following ODEs
(30)  
(31) 
where the pointparticle contribution reads
(32) 
with the coefficient given in Appendix A, and the tidal term is given by
(33)  
where is the Love number of the body , and means the same terms but whit the label and exchanged. As we discuss in Section III.2, the values of the Love number for different stellar models can be computed with our approach, and agree with the values obtained in the literature H08 ().
The orbital separation is evaluated through the PN expression for , which is known up to order PN, including spin terms Fav11 (), and is found solving the equation
(34)  
where and the spin variables are defined as follows
(35)  
(36) 
are the spin angular momenta of bodies , with dimensionless spin parameters and unit direction vectors .
It is important to remark that the adiabatic inspiral of the orbital motion and the “adiabatic approximation” for the Love number, are two different approximations: the first assumes that the orbital timescale is much larger than that associated to the gravitational wave energy loss (orbital adiabatic approximation); the second assumes, as mentioned in the Introduction and in Section III.2, that the orbital timescale is much larger than the timescale associated to the NS internal dynamics (Love number adiabatic approximation). In this paper, we use the orbital adiabatic approximation, but we drop the Love number adiabatic approximation.
ii.4 PostNewtonian tidal deformations
Tidal interactions in binary systems have been studied by many authors in the framework of general relativity (see for instance M75 (); M83 ()). They are described by the equation of geodesic deviation:
(37) 
where is the Riemann tensor, , and, in the present case, is the velocity of the star center and is the separation vector between and a generic fluid element. By introducing an orthonormal tetrad field associated with the frame centered in , parallel transported along its motion, and such that , Eq. (37) can be cast in the form P56 ()
(38) 
where the , and are the components of the relativistic tidal tensor, defined in Eq. (1). In the affine approach, Eq. (38) applies with replaced by .
In the following subsections, starting from the PN metric given in Eqns. (24)(26), we write the explicit expression of the parallel transported tetrad, and compute the tidal tensor assuming equatorial motion.
The parallel transported tetrad
The orthonormal tetrad associated to the PN metric (24)(26), satisfies the FermiWalker transport equations (expressed in terms of the coordinate time , rather than of the proper time of ) F88 ():
(39) 
where
(40) 
is the acceleration of and its coordinate velocity, with . The tetrad vectors are F88 ():
where
(41)  
and
(42) 
is the angle describing geodesic precession and frame dragging, given in terms of the antisymmetric matrix Q defined as
(43) 
is the post Newtonian potential, associated with the components of the PN metric (25), and is an arbitrary integration constant.
The tidal tensor
Having defined the tetrad field, we have explicitly computed (with the help of the symbolic manipulation software maple and the package GRTensor) the Riemann and the tidal tensors, up to order . The general structure of the tidal tensor components in terms of the derivatives of the PN potentials, is given in Appendix B; here we show, as an example, the component :
(44)  
(45)  
To hereafter we omit the parentheses to indicate the tetrad components of the tidal tensor. In Eqns. (45) means the same term but with the labels and exchanged; and , where x is the field point and the trajectory of (and similarly for ) ; is the coordinate velocity, the relative displacement between the two masses, and the relative velocity. is the dimensional antisymmetric LeviCivita symbol with , and is the spinvector.
The steps needed to evaluate the tidal tensor at the center of the NS (i.e., at location 1 in our conventions) are the following:

estimate the various derivatives of the PN potentials ;

compute the tidal tensor at the source location , and apply a regularization procedure;

express all quantities in the twobody center of mass frame;

switch to the star principal frame defined in Section II.1.
The second point requires further clarifications. Since the PN metric refers to pointlike sources, the PN potentials (45) diverge when computed at the source locations and . To remove this divergence, we apply the Hadamard regularization procedure BFP98 (), which we briefly describe. Let be a function depending on the field point x and on the two source locations ,, and admitting, when x approaches , an expansion of the form
(46) 
The regularized value of at the point is the Hadamard part finie, which is the average, with respect to the direction , of the term in the sum (46):
(47) 
We use this procedure to evaluate and their derivatives at the source point. We remark that the regularization procedure should be applied separately on the Riemann tensor and on the orthonormal tetrad. However, at the PN order we are considering, it is perfectly equivalent to apply the regularization procedure directly to the tidal tensor .
We now express the point particle positions in the system center of mass frame, by the following coordinate transformation BI03 ()
(48)  
where
(49) 
Finally, we express in the principal frame using the rotation matrix
(50) 
where is defined by Eq. (3)
(51) 
The complete form of the tidal tensor is given by:
(55)  
where the lag angle describes the misalignment between the axis and the line between the two objects. In the tidal tensor components the dot indicates differentiation with respect to the coordinate time . In the principal frame, the geodesic deviation equation for the tidal deformation can be written as
(56) 
It should be stressed that, as noted in B06 (), if the system is in quasicircular inspiral, the radial motion is due only to gravitational backreaction; consequently