Luc Blanchet

Alessandra Buonanno

Guillaume Faye

Tail-induced spin-orbit effect in the gravitational radiation of compact binaries

Keywords: numbers: 04, 30, -w, 04, 25, -g

Gravitational waves contain tail effects which are due to the back-scattering of linear waves in the curved space-time geometry around the source. In this paper we improve the knowledge and accuracy of the two-body inspiraling post-Newtonian (PN) dynamics and gravitational-wave signal by computing the spin-orbit terms induced by tail effects. Notably, we derive those terms at 3PN order in the gravitational-wave energy flux, and 2.5PN and 3PN orders in the wave polarizations. This is then used to derive the spin-orbit tail effects in the phasing through 3PN order. Our results can be employed to carry out more accurate comparisons with numerical-relativity simulations and to improve the accuracy of analytical templates aimed at describing the whole process of inspiral, merger and ringdown.

I. INTRODUCTION A. Motivation

During the last ten years a network of groundbased laser-interferometer gravitational-wave detectors has been built and has taken data at design sensitivity. It is a worldwide network composed of the Laser Interferometer Gravitational wave Observatory (LIGO), Virgo, GEO-600, and TAMA and it has operated in the frequency range 10-10 3 Hz. Coalescing binary systems composed of black holes and/or neutron stars are among the most promising sources for those detectors. By 2016 the gravitational-wave detectors will be upgraded to a sensitivity such that event rates for coalescing binary systems will increase by approximately a factor one thousand, thus making likely the first detection of gravitational waves from those systems. In the future, space-based detectors like LISA should detect supermassive black-hole binary systems in the low frequency band 10 -4 -10 -2 Hz.

The search for gravitational waves from coalescing binary systems and the extraction of source parameters are based on the matched-filtering technique, which requires a rather accurate knowledge of the waveform of the incoming signal. In particular, the detection and subsequent data analysis are made by using a bank of templates modeling the gravitational wave emitted by the source. The need of a faithful template bank has driven the development of accurate templates over the last thirty years.

The post-Newtonian (PN) expansion is the most powerful approximation scheme in analytical relativity capable of describing the two-body dynamics and gravitational-wave emission of inspiraling compact binary systems [1]. The PN expansion is an expansion in the ratio of the characteristic orbital velocity of the binary v to the speed of light c. However, as the black holes approach each other toward merger, we expect the PN expansion to lose accuracy because the velocity of the holes approaches the speed of light. At that point, numerical relativity (NR) [2][3][4] plays a crucial role providing us with the dynamics and gravitational-wave emission of the last cycles of inspiral, followed by the merger and ringdown phases. Furthermore, by properly combining PN predictions and NR results, it is possible to describe analytically and/or numerically with high accuracy, the full gravitational-wave signal [5][6][7][8].

Black holes in binary systems can carry spin, and when spins are not aligned with the orbital angular momentum, spins induce precession of the orbital plane (see e.g. Ref. [9]). This adds substantial complexity to the gravitational waveforms, making it indispensable to include spin effects in templates used for the search. Moreover, as found long time ago [10][11][12][13][14][15][16][17][18], gravitational waves contain tail effects which are due to the back-scattering of linear waves in the curved space-time geometry around the source (and primarily generated by its mass). This causes the gravitational-wave signal to depend on the entire history of the binary system.

In this paper we improve the knowledge and accuracy of the two-body inspiraling dynamics and gravitationalwave signal by computing the spin-orbit (SO) terms induced by tail effects. This is the continuation of our previous work on spins [19,20] where we obtained the next-to-leading 2.5PN SO contributions in the equations of motion and gravitational-wave energy flux. Here, we derive those SO terms at 3PN order in the gravitationalwave energy flux, where they are entirely due to tails. Furthermore we obtain the SO terms at 2.5PN and 3PN orders in the wave polarizations that are specifically due to tails, leaving aside other SO terms at these orders that come from instantaneous (non-tail) linear contributions and which will not be computed here.

We obtain the energy flux in two independent ways, first directly using the radiative multipole moments, and second by differentiating and squaring the gravitationalwave polarizations. To compute the SO tail effects in the wave polarizations we solve the two-body dynamics taking into account spin precession. Assuming quasicircular adiabatic inspiral, we also compute the 3PN SO terms induced by tails in the gravitational-wave phasing. These results can be used to improve the accuracy of inspiraling templates, to carry out comparison with numerical-relativity predictions, and to improve the accuracy of effective-one-body and phenomenological templates [5][6][7][8].

As an important check of our results we obtain the 3PN SO tail terms in the energy flux in the test-particle limit, and find perfect agreement with earlier PN computations based on black-hole perturbation induced by the motion of a test particle around a massive black hole [21].

This paper is organized as follows. In Sec. II we review the post-Newtonian multipole moment formalism and discuss relevant properties of tails. In Sec. III, we describe how spin effects are included in the PN formalism and derive the binary's evolution equations when black holes carry spins. In Sec. IV we obtain the time evolution of the moving triad and solve the precessing dynamics at the relevant PN order. In Sec. V we compute the 2.5PN and 3PN SO tail effects in the gravitational waveform and polarizations. Restricting ourselves to quasi-circular adiabatic inspiral, we derive in Sec. VI the 3PN SO tail effects in the energy flux and in the gravitational phasing. Section VII contains our main conclusions. In Appendix A we check that the equations of motion do not contain physical SO effects at 3PN order. In Appendix B we investigate the secular evolution of the spins, and finally in Appendix C we give the explicit gravitationalwave polarizations.

B. Notation

In this paper we adopt the following notations. G is the Newton constant and c is the speed of light. As usual we refer to nPN as the post-Newtonian terms with formal order O(c -2n ) relative to the Newtonian acceleration in the equations of motion, or to the lowestorder quadrupole-moment formalism for the radiation field. Greek indices are space-time indices, and Latin are space indices. The quantity

L = i 1 • • • i ℓ denotes a multi- index composed of ℓ multipolar spatial indices i 1 , • • • , i ℓ ranging from 1 to 3. Similarly L -1 = i 1 • • • i ℓ-1 and kL -2 = ki 1 • • • i ℓ-2 ; N L = N i1 • • • N i ℓ is the product of ℓ spatial vectors N i (similarly for x L = x i1 • • • x i ℓ ). We indicate with ∂ L = ∂ i1 • • • ∂ i ℓ and ∂ kL-2 = ∂ k ∂ i1 • • • ∂ i ℓ-2
the product of partial derivatives ∂ i = ∂/∂x i . In the case of summed-up (dummy) multi-indices L, we do not write the ℓ summations from 1 to 3 over their indices. Furthermore, the symmetric-trace-free (STF) projection is indicated using brackets,

T L = STF[T L ]; thus U L = U L and V L = V L for STF moments. For in- stance we write x i v j = 1 2 (x i v j + x j v i ) -1 3 δ ij x • v.
The transverse-traceless (TT) projection operator is denoted P TT ijkl = P ik P jl -1 2 P ij P kl where P ij = δ ij -N i N j is the projector orthogonal to the unit direction N = X/R of a radiative coordinate system X µ = (c T, X). The quantity ε ijk is the Levi-Civita antisymmetric symbol such that ε 123 = 1. Finally, we denote time derivatives with a superscript (n), and we indicate the symmetrization operation with round parentheses.

II. WAVE GENERATION FORMALISM

The gravitational waveform h TT ij , generated by an isolated source described by a stress-energy tensor T µν with compact support, and propagating in the asymptotic regions of the source, is the TT projection of the metric deviation at the leading-order 1/R in the distance to the source. It is parametrized by STF mass-type moments U L and current-type ones V L , which constitute the observables of the waveform at infinity from the source and are called the radiative moments [22]. The general expression of the TT waveform, in a suitable radiative coordinate system X µ = (c T, X), reads, when neglecting terms of the order of 1/R 2 or higher,

h TT ij = 4G c 2 R P TT ijkl +∞ ℓ=2 N L-2 c ℓ ℓ! U klL-2 - 2ℓ c(ℓ + 1) N m ε mn(k V l)nL-2 . (2.1)
Here the radiative moments U L and V L are functions of the retarded time T R ≡ T -R/c in the radiative coordinate system (we denote R = |X|). The integer ℓ refers to the multipolar order, and N = X/R = (N i ) is the unit vector pointing from the source to the far away detector.

The TT projection operator P TT ijkl and other notations are defined in Sec. I B. With P = (P i ) and Q = (Q i ) denoting two unit polarization vectors, orthogonal and transverse to the direction of propagation N , the two "plus" and "cross" polarization states of the waveform read as

h + = P i P j -Q i Q j 2 h TT ij , (2.2a)
h × = P i Q j + P j Q i 2 h TT ij . (2.2b)
Our convention for the choice of the polarization vectors P and Q in the case of binary systems will be specified in Fig. 2. Plugging Eq. (2.1) into the standard expression for the gravitational-wave energy flux we get [22]

F = +∞ ℓ=2 G c 2ℓ+1 (ℓ + 1)(ℓ + 2) (ℓ -1)ℓ ℓ!(2ℓ + 1)!! U (1) L U (1) L + 4ℓ(ℓ + 2) c 2 (ℓ -1)(ℓ + 1)!(2ℓ + 1)!! V (1) L V (1) L . (2.3)
A. Expression of the radiative moments

In the multipolar-post-Minkowskian formalism [13][14][15], the radiative moments are expressed in terms of two other sets of moments, referred to as the "canonical" moments M L , S L , and which are relevant to the description of the source's near zone. The relation between the radiative moments U L , V L and the canonical ones M L , S L encodes all the non-linearities in the wave propagation between the source and the detector [15]. Those relations may be re-expanded in a PN way and are then seen to contain, at the leading 1.5PN order, the contribution of the so-called gravitational-wave tails, due to backscattering of linear waves onto the space-time curvature associated with the total mass of the source itself. The explicit expressions at 1.5PN order are [15,23]

U L (T R ) = M (ℓ) L + 2GM c 3 TR -∞ dt M (ℓ+2) L (t) ln T R -t 2τ 0 + κ ℓ + O 1 c 3 non-tail , (2.4a) V L (T R ) = S (ℓ) L + 2GM c 3 TR -∞ dt S (ℓ+2) L (t) ln T R -t 2τ 0 + π ℓ + O 1 c 3 non-tail , (2.4b)
where M is the Arnowitt-Deser-Misner (ADM) mass associated with the source. It also coincides with the mass monopole moment. The remainders O(c -3 ) non-tail in Eqs. (2.4) denote some correction terms which are at least of order 1.5PN and are instantaneous or contain the non-linear memory effect which will not be considered in the present computation. Here κ ℓ and π ℓ denote some numerical constants given by [23]

κ ℓ = 2ℓ 2 + 5ℓ + 4 ℓ(ℓ + 1)(ℓ + 2) + ℓ-2 k=1 1 k , (2.5a)
π ℓ = ℓ -1 ℓ(ℓ + 1) + ℓ-1 k=1 1 k . (2.5b)
The constant τ 0 in Eqs. (2.4) is a freely specifiable time scale entering the relation between the radiative time T R and the corresponding retarded time in harmonic coordinates.

The canonical moments M L and S L are themselves linked to six sets of multipole moments characterizing the source, collectively called the "source" moments and denoted I L , J L , W L , X L , Y L , Z L . The point is that those source moments are known as explicit integrals extending over the pseudo-energy-momentum tensor of the matter fields and the gravitational field of the source [23,24]. In the following we shall essentially need I L and J L which represent the main mass and current moments of the source. The other moments W L , X L , Y L and Z L play a little role because they are associated with a possible gauge transformation performed at linear order. It turns out that the difference between the canonical moments M L , S L and the source moments I L , J L arises only at the 2.5PN order:

M L = I L + O 1 c 5 , (2.6a)
S L = J L + O 1 c 5 . (2.6b)
For the present application it will be sufficient to consider the source mass moments I L at 1PN order and the current ones J L at Newtonian order (see the discussion in Sec. III B). These are given by [25]

I L = d 3 x xL σ + 1 2c 2 (2ℓ + 3) xL |x| 2 σ (2) - 4(2ℓ + 1) c 2 (ℓ + 1)(2ℓ + 3) xiL σ (1) i + O 1 c 4 , (2.7a)
J L = ε ij i ℓ d 3 x xL-1 i σ j + O 1 c 2 . (2.7b)
The other moments we shall need are the mass monopole moment M and the monopole of the moment W L , which are given by

M = d 3 x σ + O 1 c 2 , (2.8a)
W = 1 3 d 3 x x i σ i + O 1 c 2 .
(2.8b)

The mass, current and tensor densities σ, σ i , σ ij in Eqs. (2.7)-(2.8), are defined as (where

T ii ≡ δ ij T ij ) σ = T 00 + T ii c 2 , σ i = T 0i c , σ ij = T ij . (2.9)
We recall that, e.g., σ (n) in Eqs. (2.7) means taking ntime derivatives. The spin parts of the source moments in Eqs. (2.7) will come from the model we adopt for the stress-energy tensor T µν appropriate to spinning compact binaries (see details in Sec. III). Importantly, we notice that for the accuracy required by our calculation of the spin effects due to tails all integrands in Eqs. (2.7) have compact support. This is in contrast with the spin effects at 2.5PN order which necessitate non-compact supported higherorder terms in the source moments [20]. (We also find that the second term in I L can be ignored in the present application to spins.)

B. Computing the tail integrals

The tail integrals in Eqs. (2.4) extend over the entire past of the evolving source and it is a priori a non trivial task to compute them. Here we recall, based on Refs. [15,16], that the tails are actually very weakly sensitive (in a post-Newtonian sense) to the past history of the source, and can essentially be computed using the current dynamics, i.e. at current time T R , of the source.

We have to compute, e.g., the integral appearing in the radiative mass multipole moment (2.4a), in which we can replace, following (2.6a), the canonical moment M L by the source moment I L . Thus,

U L (T R ) = TR -∞ dt I (ℓ+2) L (t) ln T R -t 2τ ′ 0 , (2.10)
where we pose τ ′ 0 = τ 0 e -κ ℓ . Let us introduce a constant time interval T to split the integral (2.10) into some contribution coming from the "recent past", and extending from the current time T R to T R -T , and the remaining contribution called the "remote past", from T R -T to -∞ in the past. The recent past can be thought of as corresponding to the most recent orbital period of a compact binary system, while the remote past will include the details (eventually unknown) of the formation and early past evolution of the compact binary. However we shall prove that our result is independent of the chosen time scale T .

To control the convergence of the tail integral (2.10) in the past we make a physical assumption regarding the behavior of the multipole moment I L (t) when t → -∞. We assume that at very early times the source was formed from a bunch of freely falling particles initially moving on some hyperbolic-like orbits, and forming at a later time a gravitationally bound system by emission of gravitational radiation. The gravitational motion of initially free particles is given by x i (t) = V i t+W i ln(-t)+X i +o(1), where V i and X i denote constant vectors, and W i = GmV i /V 3 (see Ref. [26] for a proof; to simplify we consider the relative motion of two particles with total mass m). Here the Landau remainder o-symbol satisfies ∂ n o(1)/∂t n = o(1/t n ). From that physical assumption we find that the multipole moment behaves when t → -∞ like

I L (t) = A L t ℓ +B L t ℓ-1 ln(-t)+C L t ℓ-1 +o(t ℓ-1 ) , (2.11)
where A L , B L and C L are constant tensors. The time derivatives of the moment appearing in Eq. (2.10) are therefore dominantly like

I (ℓ+2) L (t) = D L t -3 + o(t -3 ) , (2.12)
which ensures that the integral (2.10) is convergent.

Next we integrate the "remote-past" integral (from T R -T to -∞) by parts and make use of our assumption (2.11)-(2.12) to arrive at (posing

t = T R -T x) U L (T R ) = I (ℓ+1) L (T R ) ln T 2τ ′ 0 +T 1 0 dx ln x I (ℓ+2) L (T R -T x) + +∞ 1 dx x I (ℓ+1) L (T R -T x) . (2.13)
At this stage it is convenient to perform a Fourier decomposition of the multipole moment, i.e.

I L (t) = +∞ -∞ dΩ 2π I L (Ω) e -iΩt . (2.14)
(The Fourier coefficients satisfy I * L (Ω) = I L (-Ω) since the moment is real.) Inserting (2.14) into (2.13) we obtain a closed-form result in the Fourier domain thanks to the mathematical formula [START_REF] Gradshteyn | Table of Integrals[END_REF]

λ 1 0 dx ln x e iλx + i +∞ 1 dx x e iλx = - π 2 s(λ) -i ln |λ| + γ E , (2.15)
where λ = ΩT , with s(λ) and |λ| denoting the sign and the absolute value, and γ E being the Euler constant. Finally the result reads

U L (T R ) = i +∞ -∞ dΩ 2π (-iΩ) ℓ+1 I L (Ω) e -i Ω TR × π 2 s(Ω) + i ln(2|Ω|τ ′ 0 ) + γ E . (2.16)
We observe that the arbitrary time scale T has cancelled from this result.

Later we shall apply this result to the computation of the waveform and energy flux of a spinning compact binary system. A priori, since the tail integral (2.10) depends on all the past history of the binary (with the current binary's dynamics being the result of its long evolution by gravitational radiation emission), we expect that the binary's continuous spectrum of frequencies Ω should contain all orbital frequencies at any epoch in the past, say ω(t) with t ≤ T R , besides the current orbital frequency ω(T R ). However, it has been shown in the Appendix of Ref. [16] that one can actually compute the tail integral by considering only the current frequency ω(T R ). Indeed the error made by this procedure is small in a post-Newtonian sense, being of the order of O(ξ ln ξ), where ξ = ω/ω 2 denotes the adiabatic parameter associated with the gravitational radiation emission, and evaluated at the current time T R . In a PN expansion we have ξ(T R ) = O(1/c 5 ) so the error made by replacing the past dynamics by the current one is of the order of O(ln c/c 5 ) and can be neglected. The proof given in Ref. [16] is based on a simple model of binary evolution in the past, where an always circular orbit is decaying by radiation following the lowest order quadrupole formula, and spins are neglected. In this paper we shall assume that this result remains valid for spinning binaries.

III. APPLICATIONS TO SPINNING BINARIES A. Spin vectors for point-like objects

Following our previous work [19,20] we base our calculations on the model of point-particles with spins [START_REF] Papapetrou | [END_REF][29][30][31][32][33][34][35][START_REF] Dixon | Isolated systems in general relativity[END_REF][START_REF] Bailey | [END_REF][38][39][40]. The stress-energy tensor T µν of a system of spinning particles is the sum of a monopolar piece, made of Dirac delta-functions, plus the dipolar or spin piece, made of gradients of delta-functions:

T µν =c 2 A +∞ -∞ dτ A p (µ A u ν) A δ (4) (x -y A ) √ -g A - 1 c ∇ ρ Sρ(µ A u ν) A δ (4) (x -y A ) √ -g A , (3.1)
where δ (4) is the four-dimensional Dirac function, x µ is the field point, y µ A is the world-line of particle A, u µ A = dy µ A /(cdτ A ) is the four-velocity, such that g A µν u µ A u ν A = -1 where g A µν ≡ g µν (y A ) denotes the metric at the particle's location, p µ A is the linear momentum of the particle, and Sµν

A denotes its antisymmetric spin angular momentum. Our notation and conventions are the same as in Refs. [19,20] which provide more details, except that here we shall denote using an overbar (i.e. Sµν A ) the original spin variable used in [19,20]. Note that with our convention the spin variable has the dimension of an angular momentum times c.

In order to fix unphysical degrees of freedom associated with an arbitrariness in the definition of Sµν in the case of point particles (and associated with the freedom in the choice for the location of the center-of-mass worldline of extended bodies), we adopt the covariant supplementary spin condition also called Tulczyjew condition [33,34]:

Sµν A p A ν = 0 , (3.2)
which allows the natural definition of the spin four-vector SA µ in such a way that

Sµν A = - 1 √ -g A ε µνρσ p A ρ m A c SA σ , (3.3)
where ε µνρσ is the four-dimensional antisymmetric Levi-Civita symbol (such that ε 0123 = 1). For the spin vector SA µ itself, we choose a four-vector that is purely spatial in the particle's instantaneous rest frame, which means that in any frame

SA µ u µ A = 0 . (3.4)
This choice is also adopted in Refs. [41][42][43][44]. As a consequence of the condition (3.2), we can check that the mass defined by m 2 A c 2 = -p µ A p A µ is constant along the trajectories, i.e. dm A /dτ A = 0.

Important simplifications occur in the case of SO interactions, which are linear in the spins. Neglecting quadratic (spin-spin) interactions, the linear momentum is simply linked to the four velocity as p µ A = m A cu µ A , so the supplementary spin condition (3.2) reduces to Sµν A u A ν = 0, and the equation of evolution of the spins is given by

D SA µ dτ A = 0 , (3.5)
which means that the spin is parallely transported along the particle's trajectory. Following [19,20] we adopt in a first stage as the vector spin variable the contravariant components of the vector Si A , which are obtained by raising the index on SA i by means of the spatial metric γ ij A , denoting the inverse of the covariant spatial metric γ A ij ≡ g A ij evaluated at point

A (i.e. such that γ ik A γ A kj = δ i j ). Hence our initial spin variable is Si A = γ ij A SA j . (3.6)
This definition of the spin vector SA = ( Si A ) agrees with the choice already made in Refs. [43,44].

At the leading SO approximation, the contravariant spin variables SA defined by Eq. (3.6) coincide with the spin variables with constant magnitude broadly used in the literature (see, e.g., Ref. [42]). At the next-to-leading order, the variables SA differ from constant magnitude spins and their relationship has been worked out at 2PN order in Eq. (7.4) of Ref. [20]. In the present paper we shall denote the constant magnitude spins by S A (although they were denoted S c A in Refs. [19,20]). We know that it is actually better when presenting final results to switch to the constant magnitude spins S A since they have a simpler precession equation (and turn out to be secularly conserved, i.e., over a radiation-reaction time scale; see Ref. [45] and Appendix B below).

For two bodies (A = 1, 2) the relationship between the constant magnitude spins and the original spin variables up to 1PN order is:

S 1 = S1 + 1 c 2 - 1 2 (v 1 S1 ) v 1 + Gm 2 r 12 S1 + O 1 c 4 , (3.7)
together with the relation for the other particle obtained by exchanging all particle labels. We denote by v A = dy A /dt the coordinate velocity of the particle A (with mass m A ) and by r 12 = |y 1y 2 | the relative distance. See Ref. [20] for more accurate formulas extending (3.7) to 2PN order. In the present paper we shall consistently work only with the constant magnitude spins S A .

In the case of binary systems it is convenient to pose

S = S 1 + S 2 , (3.8a) Σ = S 2 X 2 - S 1 X 1 , (3.8b)
where

X 1 = m 1 /m and X 2 = m 2 /m (with m = m 1 + m 2
). In addition we find it useful to occasionally use the dimensionless spin variables

s = S G m 2 , σ = Σ G m 2 .
(3.9)

B. Multipole moments with spin-orbit effects

The matter-source densities (2.9) depend on the components of the stress-energy tensor. At the leading PN order, the spin contribution therein (indicated by the subscript S) reduce to

σ S = - 2 c 3 ε ijk v i 1 S j 1 ∂ k δ 1 + 1 ↔ 2 + O 1 c 5 , (3.10a) σ S i = - 1 2c ε ijk S j 1 ∂ k δ 1 + 1 ↔ 2 + O 1 c 3 , (3.10b) σ S ij = - 1 c ε kl(i v j) 1 S k 1 ∂ l δ 1 + 1 ↔ 2 + O 1 c 3 , (3.10c)
where δ 1 (x, t) = δ[x-y 1 (t)] means the three-dimensional Dirac delta-function evaluated on the particle 1, and 1 ↔ 2 means the same quantity but corresponding to the particle 2. In Ref. [20] the SO terms have been computed in the source mass quadrupole moment I ij up to next-to-leading 2.5PN order and the source current quadrupole moment J ij up to next-to-leading 1.5PN order. All the other source moments were computed at the leading SO order. Those results are sufficient for our purpose. Actually, to compute the specific contributions of tails we need only the moments at leading SO order, given for general ℓ by

I S L = 2ℓ c 3 (ℓ + 1) ℓv i 1 S j 1 ε ij i ℓ y L-1 1 (3.11a) -(ℓ -1)y i 1 S j 1 ε ij i ℓ v i ℓ-1 1 y L-2 1 + 1 ↔ 2 + O 1 c 5 , J S L = ℓ + 1 2c y L-1 1 S i ℓ 1 + 1 ↔ 2 + O 1 c 3 . (3.11b)
Because the leading SO terms scale as O(1/c 3 ) in the mass source moments, and as O(1/c) in the current source moments, the number of non-linear terms needed in the radiative moments [Eqs. (5.1) below] is small. We refer to Sec. V of [20] for higher-order expressions of SO contributions of the source quadrupole moments.

C. Equations of motion with spin-orbit effects

Here we investigate the case where the binary's orbit is nearly circular, i.e., whose radius is constant apart from small perturbations induced by the spins (as usual we neglect the gravitational radiation damping at 2.5PN order). We denote by x = y 1y 2 the relative position of the particles (and v = dx/dt). Following Ref. [42] we introduce an orthonormal moving triad {n, λ, ℓ} defined by n = x/r as before, ℓ = L N /|L N | where L N ≡ mν x × v with ν = X 1 X 2 denotes the Newtonian orbital angular momentum and ν the symmetric mass ratio, and λ = ℓ×n. Those vectors are represented on Fig. 1, which shows the geometry of the system. The orbital frequency ω is defined for general, not necessarily circular orbits, by v = ṙn + rωλ where ṙ represents the derivative of r with FIG. 1: We show (i) the source frame defined by the orthonormal basis (x, y, z), (ii) the instantaneous orbital plane which is described by the orthonormal basis (x ℓ , y ℓ , ℓ), (iii) the moving triad (n, λ, ℓ), and (iv) the direction of the total angular momentum J (agreeing by definition with the z-direction). Dashed lines show projections into the x-y plane.

respect to the coordinate time t. It is also equal to the scalar product of n and v which we denote as (nv) = ṙ. The components of the acceleration a = dv/dt along the basis {n, λ, ℓ} are then given by

n • a = r -rω 2 , (3.12a) λ • a = r ω + 2 ṙω , (3.12b)
ℓ • a = -rω λ • dℓ dt . (3.12c)
We project out the spins on this orthonormal basis, defining S = S n n + S λ λ + S ℓ ℓ and similarly for Σ. Next we impose the restriction to quasi-circular precessing orbits which is defined by the conditions r = 0 = ṙ so that v 2 = r 2 ω 2 (neglecting radiation reaction damping terms). In this way we find [19] that the equations of the relative motion in the frame of the center-of-mass are

dv dt = -ω r ω n + ω prec ℓ + O 1 c 6 . (3.13)
There is no component of the acceleration along λ. Comparing with Eqs. (3.12) in the case of circular orbits, we see that ω is indeed the orbital frequency, while what we call the "precessional frequency" ω prec = λ • dℓ/dt is proportional to the variation of ℓ in the direction of the velocity v = rωλ. We know that ω 2 is given by

ω 2 = G m r 3 1 + γ (-3 + ν) + γ 3/2 (-5s ℓ -3δσ ℓ ) + O 1 c 4 , (3.14)
where we denote δ ≡ X 1 -X 2 and s ℓ ≡ (sℓ) = s • ℓ, where the spin variables are defined by Eq. (3.9). The PN parameter is γ ≡ Gm/(rc 2 ) and we have included only the 1PN non-spin term and the leading SO correction at 1.5PN order. On the other hand, we get [19]

ω prec = -ω γ 3/2 7s n + 3δσ n + O 1 c 4 , (3.15)
where s n ≡ (sn) = s • n. At the leading 1.5PN SO order the orbital frequency (3.14), as well as ω prec , remain unchanged if we were to substitute some other variables to the spins S, Σ. However, when working at a higher PN approximation, it is more convenient to use the spin variables S, Σ, built from the constant magnitude spins.

The main advantage of the spins S A is that they satisfy the usual-looking precession equations

dS A dt = Ω A × S A , (3.16)
showing that the spins precess around the direction of Ω A , and at the rate Ω A = |Ω A |. The equation (3.16) could in principle be extended to any PN order (at the linear SO level). The precession's angular-frequency vectors Ω A have been computed up to the 2PN order for circular orbits in Ref. [20]. Here, we shall only need the 1PN leading order:

Ω 1 = ω γ 3 4 + ν 2 - 3 4 δ ℓ + O 1 c 4 . (3.17)
To obtain Ω 2 we simply have to change δ into -δ. Both precession frequencies are constant in magnitude and independent of the spins in the 1.5PN dynamics.

The equations of motion (3.13) and the precession equations (3.16) together leave invariant the total angular momentum,

J = L + 1 c S , dJ dt = 0 , (3.18)
where L denotes the orbital angular momentum. For future reference we give the components of L along the triad basis at 1PN order for non-spin effects and at the leading 1.5PN order for spin ones [19,20]:

L ℓ = Gm 2 ν c x -1/2 1 + 3 2 + ν 6 x , (3.19a)
L n = ν x c 1 2 S n + 1 2 δΣ n , (3.19b)
L λ = ν x c -3S λ -δΣ λ . (3.19c)
Note that the components L n and L λ are due to spin effects arising at order O(c -3 ). See Eq. (7.10) of Ref. [20]. Equivalently, introducing the orbital precession vector ω = ω ℓ and spin precession vector ω prec = -ω prec n, these equations can be written as

dn dt = ω × n , (4.2a)
dλ dt = ω + ω prec × λ , (4.2b
)

dℓ dt = ω prec × ℓ . (4.2c)
We recall that the spin precession frequency is given by Eq. (3.15) or equivalently

ω prec = -ω x 3/2 7s n + 3δσ n + O 1 c 4 , (4.3)
where we pose x ≡ (Gmω/c 3 ) 2/3 which defines a gaugeinvariant PN parameter agreeing with γ at leading order.

We shall now solve the evolution equations for the moving triad {n, λ, ℓ} at the 1.5PN order in the conservative dynamics (i.e., neglecting radiation reaction). We recall that the spin variables we use in this calculation are those with constant magnitude, denoted here S A .

It is convenient to introduce a fixed (inertial) orthonormal basis {x, y, z}; see Figs. 1 and2. For a given value of the total angular momentum J (a constant vector), and of the direction N = X/R of the detector as seen from the source, a canonical choice of the basis vectors is as follows: (i) z is defined to be the normalized value of J, namely J/J; (ii) y is orthogonal to the plane spanned by N and z and points to the direction that corresponds to the positive orientation of the acute angle z, N , i.e. y = z × N /|z × N |; (iii) x completes the triad. We see that x, z and N are coplanar by construction. Then, we introduce the standard spherical coordinates with the inclination angle measured from the zenith direction z and the azimuthal angle measured from x. The spherical coordinates of N and ℓ are denoted as (θ, ϕ) and (α, ι) respectively, and since N lies in the same plane as x and z, we have ϕ = 0 (see Fig. 2). Since ι is the angle between the total and orbital angular momenta, we have

sin ι = |J × ℓ| J . (4.4)
The angles (α, ι) are referred to as the precession angles. We now derive the time evolution of our triad vectors from that of the precession angles (α, ι), and of an appropriate phase Φ that specifies the position of n with respect to some reference direction. Following Ref. [46], we introduce the unit vectors

x ℓ = J × ℓ |J × ℓ| , y ℓ = ℓ × x ℓ , (4.5)
such that {x ℓ , y ℓ , ℓ} is an orthonormal basis. The phase angle Φ is defined by (see Fig. 1):

Φ = x ℓ , n = y ℓ , λ . (4.6)
The rotation takes place in the instantaneous orbital plane spanned by n and λ, and we have

n = cos Φ x ℓ + sin Φ y ℓ , (4.7a)
λ = -sin Φ x ℓ + cos Φ y ℓ , (4.7b)
from which we deduce

e -i Φ = x ℓ • n + iλ = J λ -i J n J 2 n + J 2 λ . (4.8)
Combining (4.8) with (4.4) we also get

sin ι e -i Φ = J λ -i J n J . (4.9)
By identifying the right-hand sides of Eqs. (4.1) or (4.2) with the time-derivatives of the identities (4.7) we obtain the following system of equations for the variations of α, ι and Φ, equivalent to the system (4.1),

dα dt = -ω prec sin Φ sin ι , (4.10a
)

dι dt = -ω prec cos Φ , (4.10b) dΦ dt = ω + ω prec sin Φ tan ι . (4.10c)
On the other hand, using the total angular momentum (3.18) together with the components of the orbital angular momentum given by Eqs. (3.19) -notably the fact that L n and L λ are due to SO terms dominantly of order O(c -3 ), we deduce that sin ι is a small quantity of order O(1/c). From this fact, we conclude by direct integration of the sum of Eq. (4.10a) and Eq. (4.10c) that

Φ + α = φ + O 1 c 4 , (4.11)
in which we have defined the "carrier" phase as

φ = ω dt = ω(t -t 0 ) + φ 0 , (4.12)
with φ 0 the value of the carrier phase at some arbitrary initial time t 0 . We recall that the orbital frequency (3.14) is constant in first approximation for circular motion.

The combination Φ + α being known by Eq. ( 4.11), we can further express the precession angles ι and α in first approximation in terms of the components S n and S λ of the total spin S = S 1 + S 2 . From (4.4) we find (discarding non-linear spin contributions)

sin ι = S 2 n + S 2 λ cL N + O 1 c 3 , (4.13)
where we recall that L N = mνr 2 ω denotes the Newtonian orbital angular momentum. On the other hand, using also Eq. (4.9) and the relation (4.11) we obtain at leading order

e i α = S λ -i S n S 2 n + S 2 λ e i φ + O 1 c 2 . (4.14)
[See also the more precise equations (4.21)-(4.

22).]

It remains now to obtain the explicit time variation of the components of the individual spins S A n , S A λ and S A ℓ . Using (4.13) and (4.14) [and also (4.11)] we shall then be able to obtain the explicit time variation of the precession angles and phase. Combining (3.16) and (4.1) we obtain the precession equations for the three unknowns S A n , S A λ and S A ℓ in the form of the following first-order system (valid at any PN approximation)

dS A n dt = ω -Ω A S A λ , (4.15a) dS
A λ dt = -ω -Ω A S A n -ω prec S A ℓ , (4.15b) dS
A ℓ dt = ω prec S A λ , (4.15c)
where Ω A is the norm of the precession vector of the spin A as given by (3.17), and the precession frequency ω prec is explicitly given by (4.3). Actually the terms involving ω prec in the right-hand sides of (4.15) can be neglected because they are quadratic in the spins. Thus, staying at the linear SO level, we find that the equations (4.15) can be decoupled and integrated as

S A n = S A ⊥ cos ψ A , (4.16a)
S A λ = -S A ⊥ sin ψ A , (4.16b)
S A ℓ = S A . (4.16c)
Here S A ⊥ and S A denote two constants for each spins A, and agree with the projections (which are constant at the linear SO level) of the spins perpendicular and parallel to the direction of the orbital angular momentum ℓ. The phase of each of the spins is given by

ψ A = (ω -Ω A )(t -t 0 ) + ψ 0 A , (4.17)
where ψ 0 A is the constant initial phase at time t 0 . With those results we obtain an explicit solution for the precession angles by substituting Eqs. (4.16) into the results (4.13) and (4.14). We find that ι(t) is given at the 0.5PN level by

sin ι = x 1/2 ν (s 1 ⊥ ) 2 + (s 2 ⊥ ) 2 + 2s 1 ⊥ s 2 ⊥ cos(ψ 1 -ψ 2 ) + O 1 c 3 , (4.18)
where we recall that s A ⊥ = S A ⊥ /(Gm 2 ). Knowing ι(t) we deduce α(t) from sin ι e iα = -i

x 1/2 ν e iφ s 1 ⊥ e -iψ1 + s 2 ⊥ e -iψ2 + O 1 c 3 . (4.
19) The difference of spin phases ψ 12 ≡ ψ 1 -ψ 2 readily follows from Eq. (4.17) and Eq. (3.17) at 1PN order as

ψ 12 = ψ 0 12 + 3 2 ω x δ(t -t 0 ) + O 1 c 4 . (4.20)
Moreover, Eq. (4.9) can be written more explicitly at the 1.5PN level as sin ι e -i Φ = -i

J + L ℓ + O 1 c 4 , (4.21)
where J + ≡ J n + i J λ is given at the 1.5PN order by

J + = S 1 ⊥ c e -iψ1 1 + x - 1 8 - 3 4 + δ 8 (4.22) + e iψ1 x 3 8 + ν 4 - 3δ 8 + 1 ↔ 2 + O 1 c 4 ,
and the 1PN orbital angular momentum L ℓ is known from Eq. (3.19a).

As a check of the previous solution we observe that if we take the time derivative of Eq. (4.4), then evaluate the total angular momentum J given by (3.18) together with the components of the orbital angular momentum L provided in (3.19), and use the solution (4.15) for the evolution of the spin components, we obtain

dι dt = -ω prec S λ S 2 n + S 2 λ + O 1 c 4 , (4.23)
which is consistent with (4.10b) once (4.8) is employed.

Finally we express the triad vectors n(t), λ(t) and ℓ(t) in terms of the precession angles and the carrier phase, and in terms of the initial triad and angles at the initial instant t 0 , modulo terms of order O(c -4 ). To do this we notice that the triad {n, λ, ℓ} at time t is obtained from the inertial triad {x, y, z} by the rotation associated with the three Euler angles α, ι and Φ. Similarly the initial triad {n 0 , λ 0 , ℓ 0 } at time t 0 is obtained by the rotation associated with α 0 , ι 0 and Φ 0 . So, combining those two rotations we readily obtain {n, λ, ℓ} in terms of {n 0 , λ 0 , ℓ 0 }. Using Eq. (4.11) to eliminate the phase Φ in favor of the carrier phase φ -this introduces small remainder terms O(c -4 ) -and neglecting all terms quadratic in the spins, we get

n = cos(φ -φ 0 )n 0 + sin(φ -φ 0 )λ 0 + sin ι sin(φ -α) -sin ι 0 sin(φ -α 0 ) ℓ 0 + O 1 c 4 , (4.24a) λ = -sin(φ -φ 0 )n 0 + cos(φ -φ 0 )λ 0 + sin ι cos(φ -α) -sin ι 0 cos(φ -α 0 ) ℓ 0 + O 1 c 4 , (4.24b) ℓ = ℓ 0 + -sin ι sin(φ 0 -α) + sin ι 0 sin(φ 0 -α 0 ) n 0 + -sin ι cos(φ 0 -α) + sin ι 0 cos(φ 0 -α 0 ) λ 0 + O 1 c 4 . (4.24c)
Since we have found in Eqs. (4.18)-(4.19) an explicit solution for the precession angles ι(t) and α(t), the time dependence of n, λ and ℓ is completely known. Note that in practical computations it is often more convenient to work not with n and λ but with the complex null vector m = (n + i λ)/ √ 2 and its conjugate.

V. COMPUTATION OF THE WAVEFORM

Here we shall compute the SO terms coming from all non-linear (i.e. of formal order G 2 ) contributions associated with tails consistent with the 2.5PN and 3PN orders in the waveform. We shall need only to focus on the tails entering the mass and current quadrupoles U ij and V ij (having ℓ = 2) and on the current octupole V ijk (ℓ = 3).

They indeed contain, when specialized to spinning compact binary systems, the SO contributions we are interested in. The reason is that the leading SO terms start at the 0.5PN order O(1/c) in the current moments, but only at the 1.5PN order O(1/c 3 ) in the mass moments [see Eqs. (3.11) above].

In addition to the tail integrals shown in Eq. (2.4), we shall also compute some terms of order G 2 at 2.5PN or 1.5PN order, but which are non-hereditary, i.e., merely depend on the instantaneous retarded time T R . Those corrections are given in full form by Eqs. (5.4)-(5.6) of Ref. [47], but here we shall need only, for the same reason as before, two terms involving the source current dipole moment or angular momentum J i ; the other terms will not contribute to the SO effect at 3PN order. Furthermore, we have to include some additional corrections of similar nature originating from the differences between the canonical and the source moments in Eqs. (2.6). Those are given in the general case by Eqs. (5.9)-(5.10) of Ref. [47], but we can check that only the contribution in the mass quadrupole U ij a priori matters here.

The relevant contributions in the radiative moments, including only the terms needed for the applications below (see Ref. [47] for more complete expressions) read:

δU ij = I (2) ij + 2Gm c 3 TR -∞ dt ln T R -t 2τ 0 + 11 12
I (4) ij (t) + G c 5 1 3 ε ab i I (4) j a J b + 4 W (2) I ij -W (1) I (1) ij (2)
, (5.1a)

δV ij = J (2) ij + 2Gm c 3 TR -∞ dt ln T R -t 2τ 0 + 7 6 J (4) ij (t) , (5.1b
)

δV ijk = J (3) ijk + 2Gm c 3 TR -∞ dt ln T R -t 2τ 0 + 5 3 J (5) ijk (t) - 2G c 3 J i I (4) jk , (5.1c)
where we have replaced M with m = m 1 + m 2 which is valid at the dominant order. Moreover, we shall find that the terms in the right-hand side of Eq. ( 5.1a) which depend on the moment W vanish at the considered order. Indeed, by inserting the value of σ i from (3.10b) into the potential W defined by (2.8b) and integrating by part, we find that the result is zero.

The corresponding gravitational-waveform, for which all three moments U ij , V ij and V ijk are important, is then given by

δh TT ij = 2G c 4 R P TT ijkl δU kl - 4 3c N a ε ab k δV l b - 1 2c 2 N am ε ab k δV l bm . (5.2)
The main task consists in evaluating the tail integrals (2.10) whose Fourier transforms have already been obtained in Eq. (2.16). The result (2.16) heavily relied on a physical assumption concerning the system in the remote past, namely that it was formed by freely falling incoming particles, see Eq. (2.11). We reviewed in Sec. II B that one can insert in the result (2.16) the binary's current frequency spectrum, i.e. at time T R , modulo small error terms of the order of the adiabatic parameter of the inspiral, or, more precisely, of negligible order O(ln c/c 5 ). In the case of spinning compact binaries this means that we have to include in the spectrum the current orbital frequency ω ≡ ω(T R ), and also the precession frequencies Ω 1 ≡ Ω 1 (T R ) and Ω 2 ≡ Ω 2 (T R ) of the two spins. This follows from the explicit solution of the triad {n, λ, ℓ} and of the precession equations (see Sec. III).[Notice that the precession angles α and ι always appear through the product sin ι e iα given by equations such as (4. 19).] Hence we can take for the Fourier components of the multipole moments

I L (Ω) = 2π n,n1,n2 A n,n1,n2 L δ(Ω -ω n,n1,n2 ) , (5.3)
where the frequency modes are some ω n,n1,n2 = n ω + n 1 Ω 1 + n 2 Ω 2 . The result (2.16) then becomes

U L (T R ) = n,n1,n2 i A n,n1,n2 L (-iω n,n1,n2 ) ℓ+1 e -i ωn,n 1 ,n 2 TR × π 2 s(ω n,n1,n2 ) + i ln(2|ω n,n1,n2 |τ ′ 0 ) + γ E . (5.4)
We recall from Eq. (3.17) that the precession frequencies Ω 1 and Ω 2 are small quantities of order 1PN. This means in particular that because of the explicit factor ω ℓ+1 n,n1,n2 in Eq. (5.4) (which arises from taking the time derivatives of the multipole moment and integrating), the modes for which n = 0 in tail integrals are very small, at least of order 4.5PN, and can be neglected. The SO terms in the radiative tails originate primarily from the spins present in the sources of the integrals, i.e. appropriate derivatives of multipole moments as shown in Eqs. (5.1). The SO contributions in the source moments have been given in Eqs. (3.11). However there are other crucial SO terms which originate from the nonspin parts of the moments. They come from time differentiations of the triad vectors using the evolution equations (4.1), which produce spin terms contained in the precession frequency ω prec , and from time differentiations of spins themselves via the precession equations (3.16). SO contributions may also be generated by the tail integration itself due to the precession of the triad basis {n, λ, ℓ} according to the formula (4.24). On the other hand, other SO terms come from the reduction to circular orbits when we eliminate the orbital separation r in favor of a function of the orbital frequency ω obtained from inverting the relation (3.14). Note that the two latter corrections being of 1.5PN relative order, they cannot come from anywhere but the tail integral of the Newtonian quadrupole moment.

During the practical computation, we make explicit the time dependence of the derivatives of multipole moments, computed in the center-of-mass frame as functions of the relative position, the relative velocity and both spins. For circular orbits, x and v depend only on r and ω, which are approximately constant on dynamical time-scales, and on the unit vectors n and λ. Thus, the whole time dependence arises through that of n and λ (and ℓ = n × λ), and is provided by our explicit solution (4.24), together with the precessing angles α(t) and ι(t) given by Eqs. The complete results for the spin dependent parts of the radiative moments, in which we use the short-hand notation for spins (3.9) and where the basis vectors {n, λ, ℓ} are evaluated at the current time T R , are then

δU ij = 2mνx 4 c 2 1 3 (71s n + 35δσ n ) π n (i ℓ j) -2λ (i ℓ j) ln(4ωτ 0 ) + γ E - 11 12
- 1 3 (29s λ + 17δσ λ ) π λ (i ℓ j) + 2n (i ℓ j) ln(4ωτ 0 ) + γ E - 11 12
+ 4 s ℓ + δ 3 σ ℓ -4n (i λ j) ln(4ωτ 0 ) + γ E - 11 12 + π(n ij -λ ij ) -2 n (i ℓ j) 19 3 s λ + 3δσ λ + λ (i ℓ j) 19 3 s n + 3δσ n - 8 3 s ℓ n (i λ j) , (5.5a
)

δV ij = -3mνx 7/2 c 3 λ i σ j ln(2ωτ 0 ) + γ E - 7 6 - π 2 n i σ j , (5.5b
)

δV ijk = -16mνx 4 c 4 s k (n ij -λ ij ) + 2(s k + δσ k ) (n ij -λ ij ) ln(4ωτ 0 ) + γ E - 5 3 + π n i λ j . (5.5c)
Insertion of the above quantities into Eq. (5.2) yields the non-linearly induced SO contributions at 2.5PN and 3PN orders in the waveform as

δ ĥTT ij = x 5/2 λ k ln(2τ 0 ω) + γ E - 7 6 - π 2 n k (N × σ) i δ jk + N a ε aki σ j + x 3 4 3 s k (n cd -λ cd ) + 2(s k + δσ k ) (n cd -λ cd ) ln(4τ 0 ω) + γ E - 5 3 + πn (c λ d) × × N ak ε aci δ jd + N ca ε aki δ jd + N ad ε aci δ jk - 2 5 N aj ε iac δ kd + 4 s ℓ + δσ ℓ 3 -4n i λ j ln(4τ 0 ω) + γ E - 11 12 + π(n ij -λ ij ) - 8 3 s ℓ n i λ j -2 n i ℓ j 19 3 s λ + 3δ σ λ + λ i ℓ j 19 3 s n + 3δσ n + 1 3 (71s n + 35δσ n ) πn i ℓ j -2λ i ℓ j ln(4ωτ 0 ) + γ E - 11 12
- 1 3 (29s λ + 17δσ λ ) πλ i ℓ j + 2n i ℓ j ln(4ωτ 0 ) + γ E - 11 12 TT , (5.6)
for which we have conveniently introduced the rescaled waveform δ ĥTT ij defined by

δh TT ij = 4Gmνx Rc 2 δ ĥTT ij .
(5.7)

The two gravitational-wave polarizations h + and h × are given in Appendix C. We have checked that the testparticle limit ν → 0 of the -2 spin-weighted spherical modes resulting from the above waveform agrees with the results of Ref. [21] (given explicitly in Ref. [48]) based on black-hole perturbation theory.

VI. ENERGY FLUX AND ORBITAL PHASING

The case of the gravitational energy flux is simpler than for the waveform, notably because we need only the contributions from the mass and current quadrupole moments, i.e.

δF = G c 5 2 5 U (1) ij δU (1) ij + 32 45c 2 V (1) ij δV (1) ij . ( 6.1)
The 3PN SO effects in the energy flux have been computed in two different ways. In the first way, we compute the time derivative of the radiative moments U ij and V ij whose SO-tail contributions are given in Eqs. (5.5), and then square these radiative moments to get the flux (6.1). The second way is completely equivalent, but entirely done by hands. It consists of writing all the separate pieces composing the energy flux (6.1), made of the coupling between some instantaneous moment (evaluated at current instant T R ) times a hereditary tail integral. The SO terms have to be included in either the instantaneous moment in front of the integral, or in the tail integral itself. This gives then several "direct" SO contributions coming from tails at relative 1.5PN order (for the mass quadrupole tail) or 0.5PN order (for the current quadrupole tail) which are then added together.

In addition there is the crucial contribution due to the reduction to circular orbits of the standard (non-spin) tail integral at 1.5PN order, for which the relation between the orbital separation r and the orbital frequency ω [as given by the inverse of Eq. (3.14)] provides a supplementary SO term at relative 1.5PN order, which thus contributes in fine at the same 3PN level as the "direct" SO tail terms. Finally, we obtain the following net result for the SO tail contribution at 3PN order in the total energy flux:

δF = 32 5 c 5 G x 8 ν 2 -16π s ℓ - 31π 6 δ σ ℓ , (6.2)
where we recall that s ℓ = s • ℓ and σ ℓ = σ • ℓ, with the spin variables s and σ being defined by Eqs. (3.7)-(3.9). Let us remark that in the energy flux the 3PN SO term is entirely constituted by the SO tails we have obtained in (6.2). So the complete 3PN SO term in the flux is provided by Eq. (6.2). Contrary to the waveform computed in Sec. V, there are no other SO terms coming from linear source moments at that order. Because the energy flux and the resulting orbital phasing is so important for gravitational-wave observations, we shall now give the complete formula for the total flux, including all non-spin terms and all linear SO terms up to 3PN order (but neglecting non-linear SS interactions). However we shall not write the known non-spin 3.5PN terms in the flux (due to non-spin tails [17]) because some yet uncalculated SO effects should conjointly appear at that order. The 3PN energy flux, complete except for SS interactions, reads then We are consistently using the constant-magnitude spins S A that are related to the original variables SA of Ref. [19,20] by Eq. (3.7); see also Eqs. (7.4) of Ref. [20].

F = 32 5 c 5 G x 5 ν 2 1 + x - 1247
The non-spin terms are given, e.g., in Ref. [1]. We find perfect agreement in the perturbative limit ν → 0 with black hole perturbation calculations reported in Ref. [21].

On the other hand the total conservative energy E of the binary is not affected by the SO terms at the 3PN order (we check this point in Appendix A), hence we have Following Ref. [20] we shall next use the standard energy balance argument to deduce the evolution of the orbital frequency even in the presence of spins. To this end we have to check that the constant-magnitude spins are secularly constant (i.e., constant over a long radiationreaction time scale) up to the right level, 1.5PN order in the present case. In Ref. [20] we have referred to the work [45] for a proof that this is correct up to relative 1PN order, i.e. considering radiation reaction effects up to 3.5PN order. In Appendix B below we extend the argument of Ref. [45] to the relative 1.5PN order, which essentially means adding the tail-induced part of the radiation reaction at 4PN order. This check being done we can thus neglect ds ℓ /dt and dσ ℓ /dt in average over a radiation-reaction time scale. An alternative way to see this is to directly compute the variation of the projection of the spins along the Newtonian orbital angular momentum, i.e. S A ℓ = S A • ℓ, using the precession equations (3.16) appropriate for constant-magnitude spins. We readily find that dS A ℓ /dt = S A •[dℓ/dt+ℓ×Ω A ], which shows that dS A ℓ /dt is at least quadratic in the spins for circular orbits. This readily follows from the facts that ℓ remains constant in the absence of spins, and that, as we have seen in Eq. (3.17), Ω A for circular orbits points in the direction of ℓ modulo spin corrections. Thus we have dS A ℓ /dt = 0 at the linear SO level (neglecting quadratic SS couplings). The argument is in principle valid up to any PN order, but is restricted to circular orbits.

E = - 1 2 m ν c 2 x 1 + x - 3
The conclusion is that the constant-magnitude spin terms can be considered as constant when computing the averaged evolution dE/dt of the energy given by Eq. (6.4). Equating then dE/dt to -F , where F is given by (6.3), we obtain the secular variation of the frequency ω -denoted ω for simplicity -as (neglecting SS contributions) ω ω 2 = 96 5 ν x We recall that to the carrier phase we have also to add the precessional correction, arising from the changing orientation of the orbital plane. We have proved in Eq. (4.11) that at the 1PN order the total phase Φ is given by Φ = φ -α + O(c -4 ). Thus the precessional correction is given by -α and is explicitly provided by the solution (4.18)-(4.19). Alternatively, the precessional correction can be computed numerically [9].

VII. CONCLUSION

So far, the search for gravitational waves with LIGO and Virgo detectors has focused on non-spinning compact binaries [49][50][51][52][53], although in Ref. [54] single-spin templates were employed, for the first time, to search for inspiraling spinning compact objects. It is timely and necessary to develop more accurate templates which include spin effects. Extrapolating results from the nonspinning case, we expect that, for maximally spinning objects, reasonably accurate templates would need to be computed at least through 3.5PN order.

During the last years, motivated by the search for gravitational waves, SO effects have been computed in the two-body equations of motion through 3.5PN order [19,41,43,44,[55][56][57] and energy flux through 2.5PN order [20,[42][43][44][45]. Moreover, SS effects have been calculated through 3PN order in the conservative dynamics [41,42,[58][59][60][61][62][63][64][65][66][67][68] and multipole moments [69].

In this paper, building on our previous work [19,20], we have improved the accuracy of the energy flux and gravitational waveform by computing SO terms induced by tail effects [10][11][12][13][14][15][16][17][18]. Those effects are due to the backscattering of linear waves in the curved space-time geometry around the source. Using the multipolar PN formalism developed in Refs. [13-15, 23, 24], we have identified and computed the radiative multipole moments responsible of tail terms involving SO couplings. More specifically, we have computed those SO tail contributions to the energy flux at 3PN order and to the gravitational waveforms at 2.5PN and 3PN order. Those SO tails constitute the complete coefficient at 3PN order in the energy flux. In particular we find that the energy flux is in complete agreement with the result of black-hole perturbations in the test-particle limit [21]. Our computation is restricted to quasi-circular inspiraling orbits, and uses the two-body precessional dynamics at 1.5PN order.

The computation of SO tail terms in the waveform is summarized in Sec. V, and some building blocks and foundation for calculating tail effects in the PN formalism were reviewed in Sec. II. For the first time, we have computed tail terms when precession effects in the twobody dynamics are also included. The relevant results for the waveform are given in Eq. (5.6) and in Appendix C. The SO tail effects in the energy flux and phasing at 3PN order are given in Sec. VI, see in particular Eqs. (6.2)-(6.3) and (6.6).

Considering the vigorous synergy which is currently taking place between analytical and numerical relativity for building faithful templates [5][6][7][8], we expect that the results developed in this paper will help the construction of more accurate analytical templates describing the entire process of inspiral, merger and ringdown of black holes in presence of spins.

In the near future we plan to complete the knowledge of SO effects in the gravitational waveform at 3PN order, by computing the non-tail (i.e., instantaneous) SO couplings at 2PN and 3PN orders, and the corresponding -2 spin-weighted spherical harmonics (or gravitational modes). This will constitute a step further with respect to Ref. [46] which computed SO effects in the gravitational modes through 1.5PN order. + 4(nS 1 )v i + 4(vS 1 )n i -2 3 (nv)S i 1 .

[Note that the latter expression depends on the specific definition of the spin one is using, i.e. in the present case the Newtonian spin defined for extended bodies by (B1); for the spin variable used in [19,20] the expression is expected to be different, but still in the form of a total time derivative.] The results (B7) or (B8) show that there is no secular evolution for the spin at the dominant 2.5PN order (see also [71] and references therein for related discussions). Note that this conclusion actually applies to any spin variable at dominant order. However it has been shown in Ref. [45] that in the case of the constant-magnitude spin there is also no secular evolution of the spins at the next-to-leading 3.5PN order in the radiation reaction. At next-to-leading order this result [45] applies specifically to the constant-magnitude spins and uses the radiation reaction acceleration up to the 3.5PN order.

We now extend the argument by including the higherorder 4PN correction term (i.e. 1.5PN radiation-reaction order) which is known to be due to gravitational wave tails [14,72]. That extension is actually straightforward since it basically needs only the Newtonian result (B7). The reason is that including the effect of tails in the radiation reaction simply amounts to replacing the source Newtonian quadrupole moment I ij by the tail-corrected expression [72]

I tail ij (t) = I ij (t) + 4GM c 3 t -∞ dt ′ I (2)
ij (t ′ ) ln t -t ′ 2τ 0 , (B9) where τ 0 denotes some arbitrary time scale, for instance the one which appears in Eqs. (2.4). Note the factor 4GM/c 3 in front of the tail integral which is twice the factor 2GM/c 3 in front of the tail integrals in (2.4). This factor ensures the consistency between the work done by the radiation reaction force in the local source and the total energy flux radiated at infinity from the source [15].

Thus the radiation reaction force including the 4PN tails takes (in harmonic coordinates) the same form as in Eq. (B3) but with I ij replaced by I tail ij . This shows that the previous Newtonian argument still holds for the 2.5PN + 4PN radiation reaction force and that the effect on the precession equation is still in the form of some irrelevant total time derivative:

dS i A dt 2.5PN+4PN = d dt - G c 5 (3)
I ij tail (t) S j A .

(B10)

Hence our conclusion that the constant-magnitude spins are secularly constant up to 4PN order corresponding to 1.5PN radiation-reaction order.

Appendix C: Gravitational-wave polarizations

We derive in this Appendix the two gravitational-wave polarizations. They are computed from the projection formulas (2.2), using the expression (5.6) [together with Eq. (5.7)] for δh TT ij . We adopt the convention shown in Fig. 2 for the polarization vectors. To shorten the result, we denote the projections of the polarisation basis {N , P , Q} onto the moving triad {n, λ, ℓ} by e.g. P n , P λ , P ℓ . With this notation, we have

δh + = Gmν c 2 R x 7/2 π(2(P n Q ℓ + P ℓ Q n )σ ℓ -2(P ℓ Q ℓ -P n Q n + P λ Q λ )σ n + 2(P λ Q n + P n Q λ )σ λ ) + 4 -(P λ Q ℓ + P ℓ Q λ )σ ℓ -(P λ Q n + P n Q λ )σ n + (P ℓ Q ℓ + P n Q n -P λ Q λ )σ λ ln(2τ 0 ω) - 7 6 + γ E + x 4 - 16 3 (P n P λ + N n P n Q ℓ -N λ P λ Q ℓ + N n P ℓ Q n + N ℓ P n Q n -N λ P ℓ Q λ -N ℓ P λ Q λ -Q n Q λ )s ℓ - 4 3 (19P ℓ P λ + 12N n P n Q n -4N λ P λ Q n -4N λ P n Q λ -4N n P λ Q λ -19Q ℓ Q λ )s n - 4 3 (19P ℓ P n + 4N λ P n Q n + 4N n P λ Q n -19Q ℓ Q n + 4N n P n Q λ -12N λ P λ Q λ )s λ + δ -12(P ℓ P λ -Q ℓ Q λ )σ n -12(P ℓ P n -Q ℓ Q n )σ λ + π 8(P n 2 -P λ 2 -Q n 2 + Q λ 2 )s ℓ + 142 3 (P ℓ P n -Q ℓ Q n )s n - 58 3 (P ℓ P λ -Q ℓ Q λ )s λ - 16 3 (N λ P n Q ℓ + N n P λ Q ℓ + N λ P ℓ Q n + N ℓ P λ Q n + N n P ℓ Q λ + N ℓ P n Q λ )(s ℓ + δσ ℓ ) - 32 3 (N λ P n Q n + N n P λ Q n + N n P n Q λ )(s n + δσ n ) + δ 8 3 (P n 2 -P λ 2 -Q n 2 + Q λ 2 )σ ℓ

FIG. 2 :

 2 FIG.2: Similar as Fig.1but with the direction of the source N indicated together with a choice of convention for the two polarization vectors P and Q.

  (4.18)-(4.19), or (4.21) and (4.22) with more precision.

Acknowledgments

We thank Etienne Racine for fruitful collaboration during the early stages of this project. A.B. acknowledges support from NSF Grant PHY-0903631 and NASA grant NNX09AI81G. L.B. acknowledges partial support from Programme International de Coopération Scientifique (CNRS-PICS).

Appendix A: 3PN spin terms in the equations of motion

In this Appendix, we check that there are no SO terms at the 3PN order in the total conservative invariant energy of the binary given by Eq. (6.4). Indeed, we find that the 3PN SO terms in the binary's equations of motion (say, in harmonic coordinates) can be gauged away. The result is to be expected because we know that the first SO modification of the radiation reaction damping force arises at the 4PN order rather than 3PN [45].

We compute the near-zone PN metric by solving the Einstein field equations in harmonic coordinates for the stress-energy tensor (3.1). We find by direct PN iteration of the metric, parametrized by means of retarded potentials V , V i , • • • (see Ref. [19] for more details), that the contribution of SO terms at 3PN order in this gauge is given by

where we indicate with cst(t) some irrelevant O(c -8 ) constant term in space, where we keep the SO parts of the acceleration un-replaced, and where 1 ↔ 2 refers to the same expression but for particle 2.

The metric (A1) yields a 3PN contribution in the equations of motion of spinning particles which can be calculated from the Papapetrou [START_REF] Papapetrou | [END_REF]29] equations of motion (see e.g. Sec. III in [19]). The result for the acceleration of particle 1 is

We observe that δa i 1 is symmetric under the exchange of particles 1 and 2. A closer inspection reveals that it is in fact given by the second total time derivative of a certain vector, namely

and we find that δX = νγ 3 (x × σ) in the notation of Sec. III C. This is precisely the effect of the gauge transformation associated with a shift of coordinates x ′i = x i + δX i . We thus conclude that the 3PN SO terms in the equations of motion are pure gauge and cannot affect the binary's invariant energy (6.4).

Appendix B: 4PN spin secular evolution

Here we show that the constant-magnitude spins are secularly constant, i.e. constant over a long radiationreaction time scale, up to the 4PN order corresponding to the 1.5PN relative order. In Ref. [45] this has already been proved up to 1PN relative order; here we extend the argument to 1.5PN order. [In the main text after Eq. ( 6.4) we present an alternative argument valid at any PN order but restricted to circular orbits.] Following Ref. [45] we describe our source by a set of well-separated extended bodies A, supposed to be Newtonian in a first stage. We define the spin S i A of each of the bodies in the usual Newtonian way as an integral extending over the volume of body A,

where ρ * denotes the Newtonian (baryonic) mass density, and x j A is the Newtonian center of mass position of the body A. At Newtonian level the spin (B1) agrees with the definition employed in the present paper. The "baryonic" spin (B1) is only used for the purpose of this Appendix. The equation of evolution of the baryonic spin reads as

The spin precession equation follows from inserting into (B2) an explicit solution for the acceleration in terms of positions and velocities. The resulting equation is then simplified using some virial relations appropriate to the case where the compact body is "stationary", see Ref. [45]. The secular evolution of the spin is then obtained by considering the radiation reaction piece of the acceleration in Eq. (B2).

At the dominant 2.5PN level the radiation reaction acceleration inside an isolated body in harmonic coordinates is given by (see, e.g., Ref. [70])

where I ij is the source's STF quadrupole moment (at Newtonian order), and U ij is the Newtonian potential tensor defined by

(We have U ii = U , the usual Newtonian scalar potential.) It can be shown [45] that the only contribution at 2.5PN order to the spin precession equation comes from the velocity-dependent part of Eq. (B3), i.e.

The other pieces in the 2.5PN acceleration vanish when the size of the body tends to zero (compact-body limit) and may be ignored. Using a virial relation [45] we readily obtain

Because the spin is constant in the lowest approximation the latter result is a total time-derivative:

which can be moved to the left-hand side and absorbed into a negligible redefinition of the spin variable at 2.5PN order. When specialized to two compact bodies the result (B7) becomes

The authors can provide on demand a file containing the results in Mathematica input format.


Vincent Creuze
email: creuze@lirmm.fr

DISTANCE AND ORIENTATION MEASUREMENT OF A FLAT SURFACE BY A SINGLE UNDERWATER ACOUSTIC TRANSDUCER

In this paper, we detail a new method that allows measuring simultaneously distance and orientation of a flat horizontal or vertical surface (e.g. wall, flat bottom) with a single acoustic transducer. A single mono-frequency ping is necessary. Then the analysis of the peaks in the received echo is sufficient to find out both distance and orientation. This method does not require any a priori knowledge of the target backscattering properties. Its advantages are: a very high sampling rate, the smallest and cheapest sensor, and low acoustic interferences with other acoustic devices (single frequency, short ping). We present theoretical aspect, simulation results, and experiments.

INTRODUCTION

Many applications of underwater vehicles (remotely operated or autonomous) require the orientation and distance measurement of a flat surface in order to achieve partially / fully autonomous survey. For instance, video surveys of hydraulic dams by underwater vehicles allow cracks or weaknesses detection. Many interesting sensing approaches have been developed in the literature. Some rely on digital imaging associated with laser triangulation. However, this method suffers from range limitation due to the possible turbidity of the water. Other methods rely on acoustic sensing. Then distance and orientation of the wall are obtained by achieving segmentation of the image acquired by sector scanning sonars [START_REF] Petillot | Underwater Vehicle Obstacle Avoidance Using a Multi-Beam Forward Looking Sonar[END_REF] [START_REF] Leonard | Directed Sonar Sensing for Mobile Robot Navigation[END_REF], or by using directly the beams of a Doppler Velocity Log [START_REF] Damus | A new paradigm for ship hull inspection using a holonomic hovercapable AUV[END_REF]. These approaches are very well suited to the underwater environment but need expensive and quite big sensors (more than 5k€, and sizing at least 10cm x 6cm²). Other techniques use difference of time-offlight with two or more transducers, and have been implemented on underwater vehicles as well as on terrestrial vehicles [START_REF] Yata | A fast and accurate sonar ring sensor for a mobile robot[END_REF]. There also exist techniques based on a single transducer. Some rely on template matching or spectral analysis [START_REF] Yata | Measurement of reflecting direction from the center frequency of the echo in ultrasonic sensing[END_REF], but require a high computational power. An original single transducer based technique uses the measurement of the time difference of the zero-crossing points of the echo of the surface. It allows accuracy of the order of 1° within a 10° range [START_REF] Yata | Wall Following Using Angle Information Measured by a Single Ultrasonic Transducer[END_REF]. The new method that we propose in this paper is based on the directivity pattern of a single transducer and allows orientation detection within a 20° range with the same accuracy. In addition, it requires only a short ping emission at a single frequency. The processing of the received echo allows us to find out both orientation and distance of a flat surface without a priori knowledge of its backscattering properties. In comparison with other acoustic methods, ours offers the following advantages:

 Highest sampling rate (only one ping is required to compute both distance and orientation) and shortest acquisition time (in comparison with sector scanning sonar for instance).  Single frequency emission thus minimizing interferences with other acoustic devices (e.g. modem).  Low power consumption, small computational cost.  The sensor is basic, cheap and small (in comparison with arrays of sensors).  No rotating part (e.g. rotating antenna).

First, we briefly remind the principles of directivity patterns and underwater sound backscattering. Then, we will detail the method and present an algorithm for automatic computation of the received echo. In the last part, simulation and experimental results will be detailed. Before being published, this work has led to the filing of the European Patent application n°EP11305048.8, on January 17th 2011, entitled "Method and device for measuring distance and orientation using a single electro-acoustic transducer".

PROBLEM STATEMENT

Problem statement

We consider the situation of an underwater Remotely Operated Vehicle (ROV) aiming to navigate at constant distance and orientation of a vertical wall (Fig. 1). 𝑑 denotes the distance between the wall and the transducer along the u axis. 𝜓 is the wall's orientation measured between the v axis and the wall. The orientation of the line of sight of the transducer (= axis of the main lobe) regarding the vehicle axis u is denoted 𝛾 and is supposed to be constant (mechanically fixed) and known. The transducer is disc-shaped but could be replaced by any other type of transducer as long as its directivity pattern contains at least two lobes. The pitching angle of the ROV is supposed to be null, i.e. the bodyfixed axes u and v remain horizontal. Fig. 2 -Directivity pattern of a disc-shaped transducer (ra-dius=13.5 mm, f=200 kHz). The sidelobe attenuation is -17.7dB.

Directivity pattern and backscattering

It is well-known that the directivity pattern 𝑏 𝜃 of a disc shaped-transducer is given by:

𝑏 𝜃 = 2𝐽 1 2𝜋𝑓 𝑟 𝑐 𝑐 sin 𝜃 2𝜋𝑓 𝑟 𝑐 𝑐 sin 𝜃 2 (1)
where 𝐽 1 (. ) is the first order Bessel function, 𝑓 is the acoustic wave frequency (200 𝑘𝐻𝑧), 𝑟 𝑐 is the transducer radius, and 𝑐 is sound speed in water (1500 𝑚. 𝑠 -1 ) [START_REF] Urick | Principles of Underwater Sound[END_REF]. A polar representation of the directivity pattern in decibels highlights the main lobe and the sidelobes (Fig. 2). In the following, "sidelobe" will designate only the first sidelobe (17.7dB attenuation). Others are weak enough to be neglected (their attenuation is larger than 24dB). The sidelobe has usually no effect because it is very weak and hits the target (e.g. seabed in the case of echo-sounder) with large incidence angles.

The backscattering strength S s for the surface reverberation of the wall (or of the bottom in the case of bottom tracking) is given in 𝑑𝐵. 𝑚 -2 and is defined as the ratio, in decibel units, of the backscattered intensity produced by a unit area of the wall to the incident intensity. For incidence angle 𝜑 𝑖 beyond 25°, a good approximation of the backscattering strength 𝑆 𝑆 is given by Lambert's law [START_REF] Lurton | An introduction to Underwater Acoustics -Principles and Applications[END_REF]:

𝑆 𝑆 = 𝑆 0 + 10 log cos 2 𝜑 𝑖 (2)
where 𝑆 0 is a constant parameter depending on the wall characteristics, and 𝜑 𝑖 is the incidence angle.

For incidence angles 𝜑 𝑖 less than 25°, the backscattering strength 𝑆 𝑆 can be approximated by (linear approximation):

𝑆 𝑆 = 𝑆 𝑁 + 𝜑 𝑖 𝜑 0 𝑆 0 -𝑆 𝑁 + 10 log cos 2 𝜑 0 (3)
where 𝑆 𝑁 is the backscattering strength for normal incidence and 𝜑 0 is called the transition direction.

Two lobes, two peaks

If the transducer is tilted, the incidence angle of the main lobe is larger than the one of the sidelobe (Fig. 1). Then, the backscattering strength is lower for the main lobe signal than for the sidelobe. Moreover, due to the difference in length of the sound paths, transmission losses 𝑇𝐿 (mainly spherical spreading, as the distance is short) are lower for the sidelobe than for the main lobe. Consequently, in spite of directivity attenuation, the global attenuation level 𝐻 makes the sidelobe echo level significant in comparison with the main lobe echo level. Moreover, as the distance to the wall is smaller for the sidelobe than for the main lobe (Fig. 1), the echo of this latter will be received after the one of the sidelobe. Thus, the backscattered acoustic echo consists of two separated parts corresponding first to the sidelobe and secondly to the main lobe (Fig. 6).

In the following, we will explain how the time of occurrences of the maximum of both echoes can lead to the computation not only of the distance but also of the orientation.

PRINCIPLE OF DISTANCE AND ORIENTATION MEASUREMENT

Intuitive approaches

In order to find out the orientation of the wall, some approaches consist in analyzing the shape of the second part of the echo (template matching). However, the backscattering parameter 𝑆 0 is a priori unknown. It depends on the reflecting surface (concrete wall, metallic hull, seabed…) and strongly affects the shape of the second part of the echo. Thus, without a good knowledge of 𝑆 0 , the shape analysis of the second part of the echo will not lead to a precise estimation of the orientation of the wall. Another approach could be to consider that the maximum of intensity of the first part of the echo corresponds to the axis of the sidelobe, and the second part of the echo corresponds to the axis of the main lobe. Thus, one could easily recover distance and orientation by simple trigonometric equations. This is wrong. In fact, the echo of sidelobe is so narrow that its maximum could be considered to correspond to its axis, but, due to spherical spreading and backscattering angle dependency, the maximum of intensity of the second part of the echo does not correspond to the axis of the main lobe. In the following, we will explain how to determine in which direction of the beam the maximum of sound echo is received. This is also called the "highest echo direction" in Fig. 1.

Expressing the maximum of the echo as a function of orientation 𝝍

We propose to study the echo intensity as a function of 𝜃, where 𝜃 is an angle difference with the axis of the main lobe (Fig. 1). In the following, 𝜓 is the wall orientation angle and 𝑑 is the distance between the transducer and the wall along the u axis of the vehicle. The incidence angle in the 𝜃 direction is denoted 𝜑 𝑖 (𝜃), and the distance between the transducer and the wall along 𝜃 direction is 𝑟 𝜃 . Thus, we have the following relationships:

𝜑 𝑖 = 𝛾 + 𝜃 -𝜓 (4) and 𝑟(𝜃) = 𝑑 cos (𝜓 ) cos ( 𝛾+𝜃 -𝜓 ) . (5)
Let us now determine 𝐻(𝜃, 𝑑, 𝜓), the function defining the echo attenuation in 𝜃 direction. Absorption will be neglected within transmission losses as its effect is very small in comparison with spherical spreading (the measured distances in this application are less than 15 meters). Thus, considering the two-way configuration of our transducer, the two-way spherical spreading losses, and the backscattering (angle of incidence of the main lobe is beyond 25° so Lambert's law applies), we have:

𝐻 𝜃, 𝑑, 𝜓 = 2𝐽 1 2𝜋𝑓 𝑟 𝑐 𝑐 sin 𝜃 2𝜋𝑓 𝑟 𝑐 𝑐 sin 𝜃 4 1 𝑑 cos (𝜓) cos ( 𝛾 +𝜃 -𝜓)
4 𝑖 𝑓0 cos 2 𝛾 + 𝜃 -𝜓 [START_REF] Yata | Wall Following Using Angle Information Measured by a Single Ultrasonic Transducer[END_REF] where 𝑖 𝑓0 = 10 S 0 10 , with S 0 defined in 2.2.

Then, one can see that the derivative d𝐻 𝜃 ,𝑑,𝜓 d𝜃 neither depends on the distance 𝑑, nor on the backscattering parameter S 0 , which means that the angle 𝜽 𝒎𝒂𝒙 corresponding to the maximum of 𝑯(𝜽, 𝒅, 𝝍) depends only on 𝝍 (other parameters 𝑓, 𝑟 𝑐 , 𝑐, and 𝛾 are known constant values). Thus 𝜃 𝑚𝑎𝑥 can be expressed as 𝜃 𝑚𝑎𝑥 = 𝑔 𝜓 , where 𝑔 is a nonlinear function. Of course, d𝐻 𝜃 ,𝑑,𝜓 d𝜃 = 0 has no simple analytical solution, thus 𝑔 cannot be expressed analytically. However, it is very easy to numerically compute 𝑔 and interpolate it by a second order polynomial called 𝑔 (Fig. 3). As it depends only on the measured directivity pattern (it does not depend on 𝑑 or S 0 ), 𝑔 has to be computed only one time (off-line). Once we know 𝑔 , it becomes easy to estimate the wall orientation. The use of the interpolation 𝑔 instead of 𝑔 will allow easier implementation (on DSP) of the algorithm described in the following. It has also to be noticed that the assumption of Lambert's law incidence of the main lobe is true only if the orientation of the wall with respect to the v axis of the ROV is greater than -10° (this mainly depends on the aperture of the main lobe and on the 𝛾 angle). As orientation and distance of the ROV with respect to the wall is close-loop controlled, 𝜓 can be assumed to be maintained close to 0°. Initialisation of this process is out of the scope of this paper, but one can briefly say that this is achieved by a 360° rotation of the ROV. In fact, when 𝜓 is out of the range [-10° +10°], there is either only one echo (superimposition) or no echo (sound is not emitted in the wall direction).

Orientation and distance

Angles shown on fig. 1 can be expressed as follows:

𝐴𝐵𝑂 = 𝜋 2 -𝜓 + 𝛾 -𝜃 𝑠 (7) 𝐴𝐷𝑂 = 𝜋 2 -𝜓 + 𝛾 + 𝜃 𝑚𝑎𝑥 (8) 𝑂𝐴𝐷 = 𝜋 2 + 𝜓 (9)
where 𝜃 𝑠 is the angle between the axis of the sidelobe and the axis of the main lobe (obtained from the directivity pattern of the transducer). Then, we have the following equations: where 𝑡 𝑠 and 𝑡 𝑚 are respectively the times of occurrence of the maxima of the two parts of the received echo (Fig. 4).

As 𝑐 denotes the sound speed, 𝑐. 𝑡 2 denotes the distance between the source and an impact point. Equations ( 10) and (11) lead to:

𝑡 𝑠 cos 𝜓 -𝛾 + 𝜃 𝑠 -𝑡 𝑚 cos 𝛾 -𝜓 + 𝜃 𝑚𝑎𝑥 = 0 (12)
Replacing 𝜃 𝑚𝑎𝑥 by 𝑔 𝜓 in (12) gives us the equation to solve to find out 𝜓.

𝑡 𝑠 cos 𝜓 -𝛾 + 𝜃 𝑠 -𝑡 𝑚 cos 𝛾 -𝜓 + 𝑔 𝜓 = 0 (13

)
Although there is no analytical solution, equation ( 13) can be easily numerically solved. Once orientation 𝜓 has been determined, the distance 𝑑 is computed by applying (10) or (11).

Ping duration

The duration of the acoustic ping to be emitted depends on the distance 𝑑. The closer from the surface we are, the shorter must be the ping in order to keep the ability to distinguish both peaks. In fact the shape of the peaks 𝑡 𝑠 and 𝑡 𝑚 depends on the ping duration, and a too long ping would lead to superimposed echoes of both lobes.

AUTOMATIC SIGNAL PROCESSING

Automatic processing of the received acoustic echo can be implemented by the following algorithm:

Step 1 Low-pass filtering of the echo to extract the envelope of the echo.

Step 2 Detecting the summits locations 𝑡 𝑠 and 𝑡 𝑚 , for instance by computing the first derivative of the echo envelope.

Step 3 Numerically solving (13), for instance by computing all the values of the left part of the equation starting from 𝜓 𝑚𝑖𝑛 (the smallest detectable value within the angular range of this system) until its sign changes (of course more efficient methods can be applied depending of the computing capabilities of the chosen device).

Step 4 Once orientation 𝜓 has been computed, it is easy to compute the distance 𝑑 by using (10) or (11).

SIMULATION RESULTS

Simulation environment

Simulations have been conducted with Matlab TM environment. The reflecting surface is meshed. Each cell of the grid is far smaller than the product 𝑐. 𝜏, where 𝜏 is the ping duration and 𝑐 is the sound speed. All the contributions of the mesh are weighted (depending on their attenuation), time delayed (depending on their position) and then added to compute the received echo.

Comparison of the simulation results with experimental acquisitions shows the ability of the simulator to reproduce the echo shape.

As the echo intensity strongly depends on both sensor and environmental parameters, the echo intensity plotted on simulations has no unit (they should be in Pa or dB). In fact, on the following plots, the intensity scale is relative and should be multiplied by an appropriate constant to get the true intensity. As we only focus on the shape of the echo and the arrival time of its peaks, the depicted intensity will be only relative.

Simulation results

Simulations have shown that 𝜓 can be properly estimated within the range [-10° +10°]. Depending of backscattering properties of the reflecting surface and depending on the directivity pattern of the transducer, this latter range value could be extended. However for autonomous navigation purposes, the vehicle orientation has to be kept constant (with respect to the tracked surface). Once the sensor is well oriented, the measurement range of our system is sufficient for servoing the vehicle position. If a wider measurement range is necessary, the transducer could be actuated by a servomotor.

Figure 4 shows simulation results with the following set of parameters: 𝑓=200 𝑘𝐻𝑧, 𝑟 𝑐 = 1cm, 𝑐 = 1500 𝑚. 𝑠 -1 , 𝛾 = -37.8°, 𝑑 = 10𝑚, 𝜓 = +10° (first plot) or 𝜓 = 0° (second plot). One can see the simulated received echo (in blue solid line) and the times of arrival of the two maxima (in red). Applying the automatic signal processing algorithm described in §4 leads to the following results: 𝑑 = 9.94𝑚, 𝜓 = +9.9° (first plot) and 𝑑 = 10.08𝑚, 𝜓 = 9.7° (second plot).

For values smaller than -10°, echoes of the two lobes occur too closely or are superimposed (Fig. 5).

For values bigger than +10°, the intensity of the sidelobe echo decreases and it becomes harder to detect its maximum. The errors measured during simulations (within the range -10° +10°) are the following: distance error < 2%, orientation error < 1°.

EXPERIMENTS

We are currently conducting experiments of the system (pool tests). First results confirm the simulation results but the number of experiments is not yet sufficient to give a precise experimental estimate of the measurement errors. Peaks can easily be detected and lead to an appropriate estimate of distance and orientation of the bottom of the pool.

CONCLUSION

More experiments will be conducted in the next months in order to verify the experimental accuracy of the system. This work has been presented for the 2D geometry of a vertical wall. The case of non vertical surfaces (3D) is under study. However, as the above-described principles also apply in the air, many applications exist for terrestrial (e.g. wall following) or aerial robotics (floor or ceiling following).

Figure 1 -

 1 Figure 1 -Top view of the underwater vehicle (ROV) navigating along the vertical wall. The transducer is placed in O. (O,u,v) is the body-fixed frame of the ROV.

Figure 3 -

 3 Figure 3 -With the following set of parameters (𝒇 = 𝟐𝟎𝟎 𝒌𝑯𝒛, 𝒓 𝒄 = 𝟏𝐜𝐦, 𝒄 = 𝟏𝟓𝟎𝟎 𝒎. 𝒔 -𝟏 , 𝜸 = -𝟒𝟕. 𝟖°), function 𝜽 𝒎𝒂𝒙 = 𝒈 𝝍 (blue dots) can be interpolated (solid red line) by 𝒈 𝝍 = 𝟎. 𝟎𝟎𝟏𝟕𝝍 𝟐 + 𝟎. 𝟏𝟑𝟒𝝍 + 𝟒. 𝟓𝟎𝟔

Figure 6

 6 shows an acoustic acquisition made under the following conditions: 𝑑 = 2.2𝑚, 𝜓 = +5°. The values obtained by computing the time of arrival of the two peaks are 𝑑 = 2.18𝑚, and 𝜓 = +4.9°.

Figure 4 -

 4 Figure 4 -Simulations of acoustic echos (𝑓=200 𝑘𝐻𝑧, 𝑟 𝑐 = 1cm, 𝑐 = 1500 𝑚. 𝑠 -1 , 𝛾 = -37.8°, 𝑑 = 10𝑚, 𝜓 = 10° 𝑜𝑟 0°).The first small peak corresponds to the echo of the sidelobe, while the second peak corresponds to the echo of the main lobe.

Figure 5 -

 5 Figure 5 -Simulation of an acoustic echo outside the range (𝑓=200 𝑘𝐻𝑧, 𝑟 𝑐 = 1cm, 𝑐 = 1500 𝑚. 𝑠 -1 , 𝛾 = -37.8°, 𝑑 = 10𝑚, 𝜓 = -5°). Both peaks are superimposed and thus cannot be distinguished.

Figure 6 -

 6 Figure 6 -Experimental acquisition of an acoustic echo (in a pool). Peaks can easily be detected and lead to an appropriate estimate of distance and orientation of the bottom of the pool.