Steven Hergt

Jan Steinhoff

Gerhard Schäfer

Reduced Hamiltonian for next-to-leading order Spin-Squared Dynamics of General Compact Binaries

Keywords: numbers: 04.25.Nx, 04.20.Fy, 04.70.Bw, 97.60.Jd

Within the post Newtonian framework the fully reduced Hamiltonian (i.e., with eliminated spin supplementary condition) for the next-to-leading order spin-squared dynamics of general compact binaries is presented. The Hamiltonian is applicable to the spin dynamics of all kinds of binaries with self-gravitating components like black holes and/or neutron stars taking into account spin-induced quadrupolar deformation effects in second post-Newtonian order perturbation theory of Einstein's field equations. The corresponding equations of motion for spin, position and momentum variables are given in terms of canonical Poisson brackets. Comparison with a nonreduced potential calculated within the Effective Field Theory approach is made.

Introduction

A crucial prediction of Einstein's theory of General Relativity is the existence of gravitational waves (GWs), e.g. resulting from the inspiralling and merging process of two compact objects. Up till now those waves are purely theoretical predictions with lack of direct experimental verification, but their direct detection is under preparation by gravitational wave observatories on Earth, e.g., LIGO, VIRGO, GEO, or LISA, a future space mission [START_REF] Rowan | Gravitational wave detection by interferometry (ground and space)[END_REF].

A first indirect evidence for the existence of GWs was the observation of energy loss in the orbital motion of the Hulse-Taylor binary pulsar PSR B1913+16 being in full agreement with the predictions of Einstein's theory using the quadrupole radiation formula. This discovery was awarded the Nobel Prize in 1993. In the meantime another strong indirect evidence occurred with the double pulsar system PSR J0737-3039A and B [START_REF] Kramer | Tests of general relativity from timing the double pulsar[END_REF]. For analysis of the measured GW patterns one has to provide very accurate templates following from theory. This can be achieved by numerical calculations with the matching of functions to the results or directly using analytic tools. In the latter case, waveforms can be obtained only perturbatively due to missing analytic solutions of the Einstein field equations for two or more (spinning) compact objects (black holes, neutron stars). One of the most successful approximation methods is the post-Newtonian one in which the metric keeps close to the flat spacetime relying on the assumption that the typical velocity v in a system divided by the speed of light c is always small, v/c ∼ 1. The deviations from the flat metric can be characterized by the Newtonian potential Φ; so for a binary system, Φ/c 2 ∼ v 2 /c 2 ∼ 2 . In an appropriate limit (as → 0), the post-Newtonian (PN) approximation yields Newton's equations.

The merging process of two compact objects is divided into three time scale sectors: inspiral, merger, and ring-down. Each sector delivers characteristic theoretical wave patterns, that are hoped to be matched against measured signals in the future. The PN approximation provides an excellent analytic handling for the inspiral phase. If the PN calculations are very accurate and thus high in order one can make very sensible predictions when comparing with measured signals. In this article we focus on calculations of the next-to-leading order (NLO) dynamics of spin-induced quadrupolar deformation effects. Surely, the most compact dynamical object is the Hamiltonian, generating the equations of motion, so we calculate in section 2 the NLO spin-squared one including a constant C Q parameterizing spin-induced quadrupolar deformation effects. C Q can be given definite values describing black holes (BHs) or neutron stars (NSs). For neutron stars C Q also depends on the model or equation of state (EoS). Thus our result, as it seems to be necessary to accurately measure C Q within future GW astronomy, could help to find the right EoS. The Hamiltonian in the present paper is calculated within the canonical formalism of Arnowitt, Deser, and Misner (ADM) [START_REF] Arnowitt | The dynamics of general relativity[END_REF]. It should be noted that our Hamiltonian is fully reduced in the sense that the spin supplementary condition (SSC) is eliminated on the level of the Hamiltonian. Further we make a formal counting of the spin as c 0 and do not distinguish between fast and slowly spinning objects (see, e.g., [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF] and also Appendix A of [START_REF] Steinhoff | Canonical formulation of gravitating spinning objects at 3.5 post-Newtonian order[END_REF]).

By now there are a lot of results regarding spin effects at the conservative orders in the PN approximation. The leading order (LO) spin effects are well-known for black holes, see, e.g., [START_REF] Barker | Gravitational two-body problem with arbitrary masses, spins, and quadrupole moments[END_REF][START_REF] Peter | Interaction of two black holes in the slow-motion limit[END_REF][START_REF] Barker | The gravitational interaction: Spin, rotation, and quantum effects-a review[END_REF][START_REF] Thorne | Laws of motion and precession for black holes and other bodies[END_REF][START_REF] Thorne | Multipole expansions of gravitational radiation[END_REF]. The LO C Q -dependence is given in [START_REF] Barker | The gravitational interaction: Spin, rotation, and quantum effects-a review[END_REF][START_REF] Poisson | Gravitational waves from inspiraling compact binaries: The quadrupole-moment term[END_REF]. The NLO spin effects were only tackled recently. The first derivation of the NLO spin-orbit (SO) equations of motion (EoM) is given in [START_REF] Tagoshi | Gravitational field and equations of motion of spinning compact binaries to 2.5 post-Newtonian order[END_REF] which became further developed in [START_REF] Faye | Higher-order spin effects in the dynamics of compact binaries. I. Equations of motion[END_REF],

both in harmonic gauge. Later, within the ADM canonical formalism, a Hamiltonian presentation was achieved [START_REF] Thibault Damour | Hamiltonian of two spinning compact bodies with next-to-leading order gravitational spin-orbit coupling[END_REF] (see also [START_REF] Steinhoff | ADM canonical formalism for gravitating spinning objects[END_REF]). The NLO spin(1)-spin(2) dynamics was found in [START_REF] Steinhoff | Next-to-leading order gravitational spin(1)spin(2) dynamics in Hamiltonian form[END_REF][START_REF] Steinhoff | ADM canonical formalism for gravitating spinning objects[END_REF] and confirmed by [START_REF] Porto | Spin(1)spin(2) effects in the motion of inspiralling compact binaries at third order in the post-Newtonian expansion[END_REF][START_REF] Levi | Next to leading order gravitational spin-spin coupling with Kaluza-Klein reduction[END_REF]. Higher PN orders linear in spin were tackled recently in [START_REF] Steinhoff | Canonical formulation of self-gravitating spinning-object systems[END_REF][START_REF] Barausse | Hamiltonian of a spinning testparticle in curved spacetime[END_REF][START_REF] Steinhoff | Canonical formulation of gravitating spinning objects at 3.5 post-Newtonian order[END_REF]. In particular, Ref. [START_REF] Steinhoff | Canonical formulation of self-gravitating spinning-object systems[END_REF] extended the point-mass ADM formalism to spinning objects, valid to any order linear in spin. Even Hamiltonians of cubic and higher order in spin were obtained for binary black holes (BBHs) [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from source terms of Kerr geometry in approximate ADM coordinates[END_REF][START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF]. Besides quadrupolar deformations induced by proper rotation (spin) and treated in the present paper, tidal deformations induced through the gravitational field of the other object were also treated, see, e.g., [START_REF] Hartle | Tidal shapes and shifts on rotating black holes[END_REF][START_REF] Damour | Effective one body description of tidal effects in inspiralling compact binaries[END_REF][START_REF] Taylor | Nonrotating black hole in a post-Newtonian tidal environment[END_REF].

A nonreduced potential (i.e., with SSC not eliminated on the level of the potential) corresponding to the result of the present paper was already calculated in [START_REF] Porto | Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF], however, a term relevant for the center-of-mass motion was missing and only found recently [START_REF] Porto | Erratum: Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF]. A comparison with the result of the present paper is not trivial, if one wants to avoid comparing all (rather long) equations of motion; instead it is more efficient to stay on the level of the (relatively short) potential. In section 3 we sketch how to transform the potential from [START_REF] Porto | Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF][START_REF] Porto | Erratum: Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF] into a reduced Hamiltonian where we will find full agreement with our result of the present paper. Considering the special case of black holes (or C Q = 1), we already succeeded in calculating the Hamiltonian of the present paper in [START_REF] Steinhoff | Spin-squared Hamiltonian of next-to-leading order gravitational interaction[END_REF][START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF], providing for the first time both the spin and correct center-of-mass dynamics in this case. There we were only able to find agreement with [START_REF] Porto | Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF] in the spin precession equation, see [START_REF] Steinhoff | Comment on two recent papers regarding next-to-leading order spin-spin effects in gravitational interaction[END_REF] (after identifying a sign typo in [START_REF] Porto | Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF]). With the correction in [START_REF] Porto | Erratum: Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF] a full comparison can now be attempted. It will be provided in section 3.

More work needs to be done for an application of the result of the present paper to GW astronomy. In order to obtain the NLO radiation field (for the SO case see [START_REF] Blanchet | Higher-order spin effects in the dynamics of compact binaries. II. Radiation field[END_REF][START_REF] Blanchet | Erratum: Higher-order spin effects in the dynamics of compact binaries. II. Radiation field[END_REF]) the stress-energy tensor has to include spin-squared corrections. This stressenergy tensor arises from the one with a general quadrupole [START_REF] Steinhoff | Multipolar equations of motion for extended test bodies in general relativity[END_REF] by a spin-squared ansatz for the mass-quadrupole, see [START_REF] Steinhoff | Spin-squared Hamiltonian of next-to-leading order gravitational interaction[END_REF]. Moreover, the NLO spin contribution should be of importance for data analysis. In a recent publication [START_REF] Reisswig | Gravitational-wave detectability of equal-mass black-hole binaries with aligned spins[END_REF] it has been shown that for maximal spins (aligned with the total orbital angular momentum), the event rates are roughly thirty times larger than of those matter systems with anti-aligned spins to orbital angular momentum and eight times as large as for non-spinning binaries. So especially considering such sources the event rate will increase with the inclusion of spin effects. Further, for the creation of templates, it is useful to find a parametrization of the orbits by solving the EoM. It is common to describe the conservative dynamics in terms of certain orbital elements. Spin precession and dissipative effects can then be formulated as secular EoM of the orbital elements. For explicit solutions including spin at LO SO see, e.g., [START_REF] Königsdörffer | Post-Newtonian accurate parametric solution to the dynamics of spinning compact binaries in eccentric orbits: The leading order spin-orbit interaction[END_REF][START_REF] Tessmer | Gravitational waveforms from unequal-mass binaries with arbitrary spins under leading order spin-orbit coupling[END_REF].

The NLO Spin-Squared Hamiltonian

We start with giving a short overview concerning the calculation of the Hamiltonian in question. This calculation is done within the ADM canonical formalism [START_REF] Arnowitt | The dynamics of general relativity[END_REF]. We use units in which 16πG = c = 1, where G is the Newtonian gravitational constant and c the velocity of light. Greek indices will run over 0, 1, 2, 3, Latin over 1, 2, 3. For the signature of spacetime we choose +2. We employ the following notations: x = x i (i = 1, 2, 3) denotes a point in the 3-dimensional Euclidean space R 3 endowed with a standard Euclidean metric and a scalar product (denoted by a dot). Letters a and b are body labels (usually they are set to 1 or 2), so x a ∈ R 3 denotes the position of the ath point mass. We also define r a := xx a , r a := |r a |, n a := r a /r a ; and for a = b, r ab := x ax b , r ab := |r ab |, n ab := r ab /r ab ; | • | stands here for the length of a vector. The linear momentum vector of the ath body is denoted by p a = (p ai ), and m a denotes its mass parameter. The usual flat space spin vector of the ath body (in local coordinates) is denoted by S a = S a(i) in correspondence with our paper [START_REF] Steinhoff | Spin-squared Hamiltonian of next-to-leading order gravitational interaction[END_REF] and its associated antisymmetric tensor by S a(i)(j) = ijk S a(k) with the total antisymmetric -symbol defined as 123 = 1. We abbreviate δ (xx a ) by δ a . The partial differentiation with respect to x i is denoted by ∂ i or by a comma, i.e., ∂ i φ ≡ φ ,i ; the partial differentiation with respect to x i a we denote by ∂ ai . Following the ADM canonical formalism, the independent degress of freedom of the gravitational field are described by h T T ij , the transverse-traceless part of

h ij = g ij -δ ij (h T T ii = 0, h T T ij,j = 0)
, and by conjugate momenta π ij T T . The needed energy and linear momentum density expressions are given by

γ 1 2 T µν n µ n ν = H m(atter) , (1)
-γ

1 2 T µ i n µ = H m(atter) i , (2)
where γ = det(g ij ), γ ij is inverse to g ij , n ν is a unit timelike normal to hypersurface x 0 = const, and T µν is the energy-momentum tensor of the matter system. Hereafter, we call its constituents the "particles", but they may well represent neutron stars or black holes. This is substantiated by "general relativity's adherence to the strong equivalence principle": black holes and other compact bodies, to some approximation, obey the same laws of motion as test bodies; see, e.g., [START_REF] Peter | Black holes: Gravitational interactions[END_REF]. Also, the analysis of the initial-value solutions for black holes shows that as in electromagnetism, where image charges are described by delta functions, black holes in full general relativity can be represented by "image masses" with delta functions support [START_REF] Jaranowski | Bare masses in time-symmetric initial-value solutions for two black holes[END_REF]. It is convenient to choose the following four coordinate conditions

π ii = 0 , g ij = ψ 4 δ ij + h T T ij , ψ = 1 + 1 8 φ . (3)
The standard ADM Hamiltonian (cf. [START_REF] Arnowitt | The dynamics of general relativity[END_REF])

H = dS i (g ij,i -g jj,i ), (4)
then becomes, using the Gauss theorem,

H = -d 3 x∆φ . (5)
The integrand ∆φ = ∂ i ∂ i φ can be expressed in terms of x a , p a , S a , h T T ij and π ij T T using the constraint equations. By expansions of the field equations in powers of G and after adopting suitable regularization procedures of integrals involved (see, e.g., Ref. [START_REF] Jaranowski | Technicalities in the calculation of the 3rd post-Newtonian dynamics[END_REF] and the Appendix in [START_REF] Jaranowski | Third post-Newtonian higher order ADM Hamilton dynamics for two-body point-mass systems[END_REF]), one can determine the Hamiltonian.

Following the procedure outlined in our previous papers [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF][START_REF] Steinhoff | Spin-squared Hamiltonian of next-to-leading order gravitational interaction[END_REF], the Hamiltonian and the other generators are constructed as to fulfill the Poincaré algebra up to 2PN order depending on standard canonical variables

{x i a , p bj } = δ i j δ ab , (6)
{S (i) a , S (j) a } = ijk S (k) a , (7)
with all other brackets being zero. The coefficient equations resulting from this procedure will change due to a modified Hamiltonian and CoM vector entering the crucial relation {G i , H} = P i see Equation (2.4) in [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF]. The modification leading to spin quadrupolar deformation effects of a general compact object has to be made in the leading order spin-squared Hamiltonian labeled as H S 2 1 in equation (2.8) in [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF], which now has to include a general spin-quadrupole constant C Q see [START_REF] Poisson | Gravitational waves from inspiraling compact binaries: The quadrupole-moment term[END_REF]; additionally an appropriate S 2 1 -CoM vector has to be found. Both can be accomplished by adopting the static source expression for H m from our paper [START_REF] Steinhoff | Spin-squared Hamiltonian of next-to-leading order gravitational interaction[END_REF] equation ( 4) and incorporating the C Q constant reading

H matter S 2 1 ,static = c 1 m 1 I ij 1 δ 1 ;ij + 1 8m 1 g mn γ pj γ ql γ mi ,p γ nk ,q Ŝ1ij Ŝ1kl δ 1 + 1 4m 1 γ ij γ mn γ kl ,m Ŝ1ln Ŝ1jk δ 1 ,i , (8)
c 1 = - 1 2 C Q , (9)
meaning C Q = 1 for BH. Symbolic abbreviations in this formula are taken unaltered from the original paper thus denoting the same mathematical objects. This means that S (i) being given in an Euclidean basis can be related to a spin tensor Ŝij in a coordinate basis with the help of a triad (dreibein) e i(j) by Ŝij = e i(k) e j(l) klm S (m) .

The dreibein as a function of the metric is just e i(j) = ψ 2 δ ij because the metric can be taken as conformally flat, g ij = ψ 4 δ ij , in our approximation. The mass-quadrupole tensor of object 1, I ij 1 , is given by

I ij 1 ≡ γ ik γ jl γ mn Ŝ1km Ŝ1nl + 2 3 S 2 1 γ ij , (10)
2S 2 1 = γ ik γ jl Ŝ1ij Ŝ1kl = const . (11)
The relation to

Q ij 1 and a 2 1 is I (i)(j) 1 = m 2 1 Q ij 1 and S 2 1 = m 2 1 a 2 1
, so in leading order the related quadrupole-moment tensor Q ij 1 is just given by

Q ij 1 = a (i) 1 a (j) 1 - 1 3 a 2 1 δ ij (12)
This static source alone is also enough to determine all the G 2 terms (static, free of linear momenta) of the Hamiltonian in question. The LO spin-squared Hamiltonian and the S 2 1 -CoM vector are calculated via the formulae H =d 3 x∆φ and G i =d 3 x x i ∆φ, respectively, with a post-Newtonian perturbatively expanded ∆φ and H matter S 2 1 ,static according to equations (4.14) -(4.16) in [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF]. The results are

H CQ S 2 1 = G 2 m 1 m 2 r 3 12 C Q 3 (S 1 • n 12 ) 2 m 2 1 - S 2 1 m 2 1 , (13)
and

G S 2 1 = G m 2 m 1 ν 1 (S 1 • n 12 ) S 1 r 2 12 + (S 1 • n 12 ) 2 r 3 12 (ν 2 x 1 + ν 3 x 2 ) + S 2 1 r 3 12 (ν 4 x 1 + ν 5 x 2 ) ( 14
)
with coefficients

ν 1 = - 1 2 - 3 2 C Q , ν 2 = 3 4 C Q , ν 3 = 3 4 C Q , ν 4 = 1 2 + 1 4 C Q , ν 5 = - 1 2 - 3 4 C Q . ( 15
)
The non-static (momenta based) terms of the Hamiltonian will be determined via the same ansatzes for the source terms in H m and H m i as in [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF] with the LO quadrupole moment [START_REF] Tagoshi | Gravitational field and equations of motion of spinning compact binaries to 2.5 post-Newtonian order[END_REF] reading

H m = 2 b=1 - m b 2 C Q Q ij b ∂ i ∂ j - 1 2 p b • (a b × ∂) + γ ij p bi p bj + m 2 b 1/2 + λ 1 p 2 b 2m b Q ij b ∂ i ∂ j + λ 2 m b (p b • ∂)Q ij b p bi ∂ j + λ 3 m b a 2 b (p b • ∂) 2 -λ 8 p b • (a b × ∂) Q ij b ∂ i ∂ j δ b , (16)
H m i = -2 2 b=1 Q kl b λ 5 p bk ∂ l ∂ i + λ 6 p bi ∂ k ∂ l + λ 7 (p b • ∂) δ li ∂ k + λ 4 a 2 b (p b • ∂)∂ i + m b 4 (a b × ∂) i 1 - 1 6 Q kl b ∂ k ∂ l - 1 2 p bi δ b . ( 17
)
Notice that the static term in H m involves the C Q constant which follows from the expansion of ( 8) being the only modification of Eq. (4.11) in [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF]. These sources allow the calculation of the NLO spin-squared Hamiltonian (with yet undetermined coefficients) leading to the same coefficient equations (4.50) -(4.62) in [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF] except for the test particle terms β 2 and β 4 which are just multiplied by the quadrupole constant C Q . These equations have to be matched to the ones resulting from the requirement of fulfilling the Poincaré algebra which now include the C Q constant and for that reason will slighty differ from the equations (3.8) - (3.22) in [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF]. The matching procedure then fixes all the coefficients left attributing to the source term coefficients the values

λ 1 = 7 4 - 3 2 C Q , λ 2 = - 5 4 + 3 2 C Q , λ 3 = - 1 24 , λ 4 = 0, (18)
λ 5 = 1 12 - C Q 6 , λ 6 = - 1 8 + C Q 4 , λ 7 = 1 8 , (19)
which agrees with our results obtained in [START_REF] Hergt | Higher-order-in-spin interaction Hamiltonians for binary black holes from Poincaré invariance[END_REF] when setting C Q = 1. In view of section 3 we label from now on standard canonical variables with a 'hat' specifying its affiliation to the Newton-Wigner (NW) SSC except for the momentum which is chosen to be the same for NW and covariant SSC. The 'hated' variables are then called NW variables in the sense that they are standard canonical, meaning

{x i a , pbj } = δ i j δ ab , (20)
{ Ŝ(i) a , Ŝ(j) a } = ijk Ŝ(k) a , (21) pai
= p ai , (22)
all other brackets being zero. Subtleties arising from that definition of NW variables are discussed in detail in our Comment [START_REF] Steinhoff | Comment on two recent papers regarding next-to-leading order spin-spin effects in gravitational interaction[END_REF]. The resulting NLO spin-squared Hamiltonian for general compact binaries reads

H ADMcan NLO S 2 1 = G r3 12 m 2 m 3 1 - 21 8 + 9 4 C Q p 2 1 ( Ŝ1 • n12 ) 2 + 15 4 - 9 2 C Q (p 1 • n12 )( Ŝ1 • n12 )( Ŝ1 • p 1 ) + - 5 4 + 3 2 C Q ( Ŝ1 • p 1 ) 2 + - 9 8 + 3 2 C Q (p 1 • n12 ) 2 Ŝ2 1 + 5 4 - 5 4 C Q p 2 1 Ŝ2 1 + 1 m 2 1 - 15 4 C Q (p 1 • n12 )(p 2 • n12 )( Ŝ1 • n12 ) 2 + 3 - 21 4 C Q (p 1 • p 2 )( Ŝ1 • n12 ) 2 + - 3 2 + 9 2 C Q (p 2 • n12 )( Ŝ1 • n12 )( Ŝ1 • p 1 ) + -3 + 3 2 C Q (p 1 • n12 )( Ŝ1 • n12 )( Ŝ1 • p 2 ) + 3 2 - 3 2 C Q ( Ŝ1 • p 1 )( Ŝ1 • p 2 ) + 3 2 - 3 4 C Q (p 1 • n12 )(p 2 • n12 ) Ŝ2 1 + - 3 2 + 9 4 C Q (p 1 • p 2 ) Ŝ2 1 + C Q m 1 m 2 9 4 p 2 2 ( Ŝ1 • n12 ) 2 - 3 4 p 2 2 Ŝ2 1 + G 2 m 2 r4 12 2 + 1 2 C Q + m 2 m 1 1 + 2C Q Ŝ2 1 + -3 - 3 2 C Q - m 2 m 1 1 + 6C Q ( Ŝ1 • n12 ) 2 , (23)
for C Q = 1 being in full agreement with the result for BH presented for the first time in [START_REF] Steinhoff | Spin-squared Hamiltonian of next-to-leading order gravitational interaction[END_REF].

Comparison with the NLO spin(1)spin(1) potential

In order to transform the Routhian R obtained within the Effective Field Theory (EFT) approach [START_REF] Porto | Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF][START_REF] Porto | Erratum: Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF] to a nonreduced Hamiltonian H, we first have to eliminate the acceleration term [START_REF] Porto | Erratum: Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF] with the help of the Newtonian equations of motion (corresponding to a redefinition of the position variables, see [START_REF] Schäfer | Acceleration-dependent Lagrangians in general relativity[END_REF]). This generates corrections terms to the order G 2 in the NLO spin-squared potential. Next, one must replace velocities by canonical momenta p a = ∂R ∂va to get the Hamiltonian by a Legendre transformation, i.e.,

H = v 1 • p 1 + v 2 • p 2 -R . (24)
The canonical momentum p 1 necessary to cover all NLO spin effects explicitly reads

p 1 = 1 + 1 2 v 2 1 m 1 v 1 + Gm 1 m 2 2r 12 [6v 1 -7v 2 -(n 12 • v 2 )n 12 ] + G r 2 12 [m 2 (n 12 × S 1 ) + 2m 1 (n 12 × S 2 )] , (25)
and similarly for particle 2. The Poisson brackets at this stage are

{x i a , p bj } = δ i j δ ab , (26)
{S (i) a , S (j) a } = ijk S (k) a , (27)
{S (i) a , S (0)(j) a } = ijk S (0)(k) a , (28)
{S (0)(i) a , S (0)(j) a } = -ijk S (k) a , (29)
and zero otherwise. ‡ Notice that these are not yet the reduced or standard canonical brackets as S (0)(i) a is still an independent degree of freedom and was not eliminated using the covariant SSC S µν a u ν = 0. It is well known that one has to proceed to Dirac brackets (DBs) if S (0)(i) a is going to be eliminated from the Hamiltonian H using a SSC, see, e.g., [START_REF] Paul | Lectures on Quantum Mechanics[END_REF][START_REF] Hanson | The relativistic spherical top[END_REF]. However, it is possible to find new variables xi a , pbj and Ŝ(j) a for which the Dirac brackets take on the standard form,

{x i a , pbj } DB = δ i j δ ab , (30)
{ Ŝ(i) a , Ŝ(j) a } DB = ijk Ŝ(k) a , (31)
and zero otherwise. These new variables can only be unique up to canonical transformations. This freedom allows us to choose p 1i = p1i , as for the flat space case [START_REF] Hanson | The relativistic spherical top[END_REF]. A possible transition to xi a and Ŝj a then reads,

S 1(i)(j) = Ŝ1(i)(j) - p 1[i Ŝ1(j)](k) p 1k m 2 1 1 - 3p 2 1 4m 2 1 - 2Gm 2 m 2 1 r12 p 1[i Ŝ1(j)](k) p 1k + 3G m 1 r12 p 1[i Ŝ1(j)](k) p 1k + G m 1 r12 p 1[i Ŝ1(j)](k) nk 12 (n 12 • p 2 ) + 2Gm 2 m 2 1 r2 12 p 1[i Ŝ1(j)](l) Ŝ1(k)(l) nk 12 + 2G m 1 r2 12 p 1[i Ŝ1(j)](l) Ŝ2(k)(l) nk 12 , (32)
x 1i = x1i - 1 2m 2 1 p 1k Ŝ1(i)(k) 1 - p 2 1 4m 2 1 -G m 2 m 2 1 p 1k Ŝ1(i)(k) r12 + 3 2 G p 2k Ŝ1(i)(k) m 1 r12 + G 2 nk 12 (n 12 • p 2 ) Ŝ1(i)(k) m 1 r12 + G m 2 m 2 1 Ŝ1(k)(l) Ŝ1(i)(l) nk 12 r2 12 + G nk 12 Ŝ1(i)(l) Ŝ2(k)(l) m 1 r2 12 , (33)
where the antisymmetrization of indices pertaining to a tensor A ij is defined as

A [ij] = 1/2(A ij -A ji ).
The rather complicated form of these variable transformations reflects the complicated structure of the DBs for self-interacting spinning objects in the covariant SSC. We will elaborate on its specific calculation in another paper. Notice that these results are applicable to all NLO spin effects (for the spin-orbit contributions the corrected form of the tetrad from [START_REF] Steinhoff | Comment on two recent papers regarding next-to-leading order spin-spin effects in gravitational interaction[END_REF] has to be inserted into the SSC). To best of our knowledge this is the first time that DBs are applied to gravitationally selfinteracting spinning objects. In [START_REF] Barausse | Hamiltonian of a spinning testparticle in curved spacetime[END_REF] test-spinning objects are considered and [START_REF] Hanson | The relativistic spherical top[END_REF] covers the flat space case only. The Routhian from [START_REF] Porto | Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF] now leads us to the reduced NLO spin-squared Hamiltonian in the form [START_REF] Porto | Spin(1)spin(2) effects in the motion of inspiralling compact binaries at third order in the post-Newtonian expansion[END_REF][START_REF] Porto | Next to leading order spin(1)spin(1) effects in the motion of inspiralling compact binaries[END_REF] a different sign convention was used for the Poisson brackets of the spin. As usual, we are not showing the canonical conjugate of the spin here; see [START_REF] Hanson | The relativistic spherical top[END_REF].

H EFTcan NLOS 2 1 = G r3 12 m 2 m 3 1 - 21 8 + 9 4 C Q p 2 1 ( Ŝ1 • n12 ) 2 ‡ Notice that in
+ 21 4 - 9 2 C Q (p 1 • n 12 )( Ŝ1 • n12 )( Ŝ1 • p 1 ) + - 7 4 + 3 2 C Q ( Ŝ1 • p 1 ) 2 + - 21 8 + 9 2 C Q (p 1 • n12 ) 2 Ŝ2 1 + 7 4 - 9 4 C Q p 2 1 Ŝ2 1 + 1 m 2 1 - 15 4 C Q (p 1 • n12 )(p 2 • n12 )( Ŝ1 • n12 ) 2 + 3 - 21 4 C Q (p 1 • p 2 )( Ŝ1 • n12 ) 2 + -3 + 9 2 C Q (p 2 • n12 )( Ŝ1 • n12 )( Ŝ1 • p 1 ) + -3 + 3 2 C Q (p 1 • n12 )( Ŝ1 • n12 )( Ŝ1 • p 2 ) + 2 - 3 2 C Q ( Ŝ1 • p 1 )( Ŝ1 • p 2 ) + 3 - 15 4 C Q (p 1 • n12 )(p 2 • n12 ) Ŝ2 1 + -2 + 13 4 C Q (p 1 • p 2 ) Ŝ2 1 + C Q m 1 m 2 9 4 p 2 2 ( Ŝ1 • n12 ) 2 - 3 4 p 2 2 Ŝ2 1 + G 2 m 2 r4 12 2 + 1 2 C Q + m 2 m 1 1 2 + 3C Q Ŝ2 1 + -3 - 3 2 C Q - m 2 m 1 1 2 + 6C Q ( Ŝ1 • n12 ) 2 . ( 34
)
This Hamiltonian and the one calculated in the last section should differ only up to a canonical transformation. It should thus be possible to generate the difference ∆H NLOS 2 1 = H EFTcan

with H N being the Newtonian Hamiltonian of a two-body system and g being an appropriate generator. It turns out that with the generator

g can NLOS 2 1 = G r2 12 m 2 m 2 1 - 1 2 + C Q (p 1 • n12 ) Ŝ2 1 + 1 2 ( Ŝ1 • n12 )( Ŝ1 • p 1 ) , (36)
equation ( 35) can be fulfilled and so agreement is achieved. This means that the Hamiltonian H ADMcan NLOS 2 1 calculated with the aid of the ADM method, in terms of invariant physical quantities, agrees with the Routhian from above, hence there is great confidence that H ADMcan NLOS 2 1 correctly describes a binary consisting of BHs and/or NSs or other kinds of compact objects in post-Newtonian Einsteinian theory. The new Hamiltonian may find immediate application in the problem of motion of orbiting binaries as investigated and solved in e.g., [START_REF] Tessmer | Gravitational waveforms from unequal-mass binaries with arbitrary spins under leading order spin-orbit coupling[END_REF][START_REF] Memmesheimer | Third post-Newtonian accurate generalized quasi-Keplerian parametrization for compact binaries in eccentric orbits[END_REF].

Acknowledgments

The authors wish to thank M. Tessmer for useful discussions. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) through SFB/TR7 "Gravitational Wave Astronomy".


Ian Stavness
email: stavness@ece.ubc.ca

John E Lloyd

Yohan Payan

Sidney Fels

Towards Speech Articulation Simulation with a Dynamic Coupled Face-Jaw-Tongue Model

Introduction

Understanding the motor control strategies underlying human speech production is an active area of research. In order to separate features of observed speech movements that arise from central motor commands from those that are due to the mechanics of the peripheral musculoskeletal system, a number of researchers have developed detailed biomechanical models jaw, tongue, larynx, lips, and face. The finite-element method (FEM) is commonly used to characterize the non-linear mechanical properties of the face, lips, tongue and other soft-tissue structures of the upper airway. A number of models have been proposed, but typically focus on one single articulator. 3D FEM models of the face include [START_REF] Sifakis | Simulating speech with a physics-based facial muscle model[END_REF] and [START_REF] Nazari | Shaping by stiffening: a modeling study for lips[END_REF]. Coupled jaw-tongue-hyoid models include an early 2D model by [START_REF] Sanguineti | A dynamic biomechanical model for neural control of speech production[END_REF], a 3D mass-spring tongue and rigidbody jaw model by [START_REF] Fang | A model-based investigation of activations of the tongue muscles in vowel production[END_REF], and a 3D dynamic FEM tongue by [START_REF] Buchaillard | A biomechanical model of cardinal vowel production: Muscle activations and the impact of gravity on tongue positioning[END_REF] coupled with a rigid-body jaw-hyoid model by [START_REF] Stavness | Coupled hard-soft tissue simulation with contact and constraints applied to jaw-tongue-hyoid dynamics[END_REF].

We are developing a biomechanics simulation toolkit ArtiSynth (www.artisynth.org) in order to simulate coupled hard and soft tissue structures. Here we present the first biomechanical model that dynamically couples the face and lips with a jaw-tongue-hyoid model, as pictured in Figure 1.

Modeling requirements

Physiologically-based simulation of speech articulation presents a number of unique challenges as compared to traditional musculoskeletal modeling of limb movements. Firstly, the upper airway involves a combination of hard tissues (skull and jaw bones) and softtissues (face, lips, tongue, and soft-palate). While non-linear FEM models can be used to represent biological soft-tissue mechanics, the technique is computationally expensive and using detailed FEM models to represent the predominantly rigid movement of skeletal structures is particularly inefficient. For this reason, we have focused on developing constraint mechanisms in ArtiSynth to handle mechanical coupling between rigid and deformable bodies. Speech movements are typically observed only in the mid-sagittal plane (with electromagnetic articulography, dynamic MRI, ultrasound, or video fluoroscopic recordings). However, the 3D shape of the vocal tract is important factor in the modulation of acoustics during speech production. Also, the mid-sagittal shape of the tongue is determined by its 3D volume (lateral expansion and contraction of the tongue body) as the tongue is a nearly-incompressible tissue. For these reasons, 3D FEM modeling is important as well as incompressibility. ArtiSynth is a fully 3D simulation environment with tools for viewing and manipulating the model in 3D. ArtiSynth includes hyper-elastic FEM materials that are nearly-incompressible and includes the option to enforce FEM incompressibility as a constraint.

Another crucial component of speech production is contact. Contact is observed between different articulators and includes deformable to deformable body contact (between the upper and lower lip) as well as deformable to rigid body contact (between the lips and teeth and the tongue and palate). ArtiSynth supports mesh-based collision detection and contact handling using dynamic constraints. Sub-regions of an FEM model can be defined for purposes of handling contacts between two regions of an FEM model, e.g. between the upper and lower lips of the face mesh, and we are working on other techniques for robust self-collision detection.

Face-Jaw-Tongue model

Our face-jaw-tongue-hyoid model couples two previous reference models that were adapted and registered to Computed Tomography (CT) data for a specific subject, including a dynamic FEM face model [START_REF] Bucki | Finite element speaker-specific face model generation for the study of speech production[END_REF] and a jaw-tongue-hyoid model [START_REF] Stavness | Coupled hard-soft tissue simulation with contact and constraints applied to jaw-tongue-hyoid dynamics[END_REF]. The volumetric face mesh and the surface jaw and skull meshes had been registered to the subject's CT data; however, the inner-surface of the face mesh did not exactly conform to the jaw and skull surfaces due to inaccuracies in the original adaptation. To improve the fit, we fixed the outer surface of the face mesh and used contact detection to conform the inner face surface to the jaw/skull surfaces. Once the contact-based morphing was complete, the outer-face-surface nodes were set free and a number of inner-face-surface nodes were attached to the jaw and skull, with the nodes around the lips and cheeks left free. Face muscles were defined by serial line segments in the original face model and were adapted to the subject along with the FEM mesh [START_REF] Bucki | Finite element speaker-specific face model generation for the study of speech production[END_REF]. The adapted muscle paths were then manually smoothed to correct for registration artifacts.

Speech simulations

We have identified a few particular speech articulations that we are investigating with our new modeling platform. Example articulation simulations are shown in Figure 2. Contact between the tongue and the lower teeth occurs in the /a/ simulation. Also, contact between the upper tongue surface and the hard palate occurs during the alveolar stop /t/ and the French syllable /ta/. In our previous work we have simulated maximum voluntary tonguepalate pressure and shown it to be within a physiologically valid range. Muscle forces required to create the physiological level of pressure observed in speech articulations is of interest.

Labial phonemes are also of interest, including articulations that include lip-toteeth contact, such as /fa/ and /va/, and bilabial lip-lip contact, such as /ba/ and /pa/. Lip rounding and protrusion have been simulated with the face model alone [START_REF] Nazari | Shaping by stiffening: a modeling study for lips[END_REF] and we are interested in modulating the degree of jaw opening by activating jaw opening muscles in the coupled face-jaw-tongue model.

We will report our progress on this speech articulation simulations at the ISSP 2011 meeting. Also, the ArtiSynth modeling platform is open source and freely available for use by the ISSP community. Therefore, we hope to engage ISSP 2011 attendees to discuss the application of biomechanical modeling to the analysis of speech production mechanics and motor control.

Figure 1 .

 1 Figure 1. Oblique, oblique cut-away and sagittal cut-away views of the face-jaw-tongue model at rest.

Figure 2 .

 2 Figure 2. Oblique and sagittal cut-away views of the face-jaw-tongue model during muscle activated vowel articulations: /a/ (upper panels) and /u/ (lower panels).

Acknowledgements

We thank Pierre Badin at GIPSA-Lab, Grenoble for providing the CT data used in the subject-specific model. We also thank Mohammad A. Nazari and Pascal Perrier at GIPSA-Lab, Grenoble for their assistance with simulations of lip protrusion.


Florent Duguet

George Drettakis

Flexible Point-Based Rendering on Mobile Devices

Point-based rendering is a compact and efficient means of displaying complex geometry. Our goal is to enable flexible point-based rendering, permitting local image refinement, required for example when zooming into very complex scenes, and efficient shadow computations. We use hierarchical packed point representations, based on recursive grid data structures. Such compact structures are particularly well adapted to devices with limited memory and display resolution, such as PDA's. To achieve flexible rendering we store intermediate attributes such as normals and colors at internal nodes of the hierarchy. We examine the memory and computation tradeoffs involved in the type of structure used, and find that tri-grids (3x3x3 hierarchical grids) are a suitable compromise for many cases. We show implementations of our method on PC and on a PDA, for octrees and tri-grids. The PDA version can render objects sampled by 1.3 million points at 2.1 frames per second.

Introduction

Recent years have seen the proliferation of complex virtual environments with the deployment of the more ubiquitous computer, which we now use at work, in the household, for leisure etc. With the increasing demand for detailed and high quality geometric models, typical scene size, containing for example scanned 3D objects, easily reaches millions of geometric primitives. At the same time, display devices are becoming more diverse, including mobile or embedded platforms such as personal digital assistants (PDAs) or mobile phones, thus significantly increasing the user base for complex virtual environments.

The main motivation of our work is to allow interactive display of complex scenes on different kinds of platforms with heterogeneous rendering capabilities. The traditional 3D object representation with vertices and polygons (faces), coupled with the traditional rendering pipeline are not always well suited in such a context. The gain of standard rasterization coherence is lost when many polygons project onto a single pixel of the screen. The visualization of such complex scenes quickly becomes challenging: models do not fit in main memory even on high performance graphics workstations. This reasoning is even more relevant for mobile platforms such as PDAs where screen resolutions around a hundred thousand pixels are typical (320x240), memory is quite limited, and computational capabilities are limited to integer arithmetic with simple instructions.

To accommodate all the above constraints, we choose a packed hierarchical point-based representation for rendering. Point-based rendering has gained popularity in recent years, offering a simple-to-use level of detail mechanism, since the number of points rendered can be adapted to the screen size of the underlying object. Initial work in point-based rendering concentrated on the generation [START_REF] Stamminger | Interactive sampling and rendering for complex and procedural geometry[END_REF][START_REF] Pfister | Surfels: surface elements as rendering primitives[END_REF] and structuring [START_REF] Rusinkiewicz | Qsplat: a multiresolution point rendering system for large meshes[END_REF] of points and rendering quality [START_REF] Zwicker | Surface splatting[END_REF]. More recent work has introduced hierarchical structures for points using coding techniques, which allow very compact storage of massive models [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF].

Our work strives to achieve flexible rendering, that is the ability to render only parts of the hierarchical data structure as necessary, including interior hierarchy nodes. In particular, we need to avoid the traversal of the entire hierarchy, and the reconstruction of normals/colors for interior nodes, since both can be prohibitively expensive operations. To this end, we store model attributes such as normals and color information at intermediate nodes of the hierarchy, allowing very efficient rendering of large models, since only a fraction of the 1 total number of samples is rendered. This approach also allows us to efficiently perform view-frustum culling and level of detail selection without the need to traverse the entire hierarchy. Flexible rendering also allows the traversal of the hierarchy in a specific order, resulting in a fast, one-pass shadow mapping algorithm. The data structures we use are recursive grids. In what follows, we refer to a ρxρxρ recursive grid as a ρ-grid. Using this terminology, an octree is 2-grid, used for example by Botsch et al. [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF]. The data structures used can be stored as a packed data bit stream encoding the hierarchy, both with and without the normal/color attributes. For rendering, additional pointers need to be stored in-core allowing random access to the hierarchy.

We examine the storage and rendering time tradeoffs of these structures. This study indicates that for our flexible rendering approach tri-grids (3-grids), are an appropriate compromise, since they permit memory savings of about 33% for the in-core data structure used for rendering, compared to octrees. Our conclusions are backed up with experimental statistics on a set of typical models.

In the next section, we briefly present previous related work. We discuss the method of Botsch et al. [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF] in some detail, since we have been inspired by some of the choices they made. We next describe our storage scheme, the various tradeoffs depending on the data structure used, and why we choose to use the tri-grid. We then discuss our flexible rendering algorithm and our fast shadow approach. Compact and efficient structures such as the tri-grid become particularly important when the geometric model does not fit in main memory, especially for mobile platforms. We have implemented our structure and rendering algorithm both on a standard PC, and on a 200Mhz Compaq iPAQ H3850 running PocketPC 2002. We present results and statistics of our method, justifying our choices, and conclude.

Previous Work

For very complex scenes, the traditional graphics pipeline can be wasteful, with much effort spent transforming and rasterizing large numbers of geometric primitives which may cover less than a pixel. Point based rendering, introduced as early as 1985 by Levoy and Whitted [START_REF] Levoy | The use of points as display primitives[END_REF], addresses this issue by point-sampling complex geometry, and rendering an appropriate number of points depending on the screen size of the complex object.

Several papers have been published in this domain over the last few years. Grossman and Dally [START_REF] Grossman | Point sample rendering[END_REF] generated point representations of objects in a preprocess, and presented efficient rendering algorithms of these point sets.

One approach is to generate unstructured point sets, by generating point-samples stochastically as a preprocess or procedurally on the fly [START_REF] Stamminger | Interactive sampling and rendering for complex and procedural geometry[END_REF]. Another approach is to generate structured point sets, for example by creating a hierarchy. The QSplat algorithm is such an approach permitting the visualization of very complex models [START_REF] Rusinkiewicz | Qsplat: a multiresolution point rendering system for large meshes[END_REF], in particular those that do not fit in main memory. A hierarchy of bounding volumes is created, and flexible rendering is achieved by halting the descent into the hierarchy depending on rendering speed requirements and screen size. Intermediate normal and color attributes are stored in this hierarchy. This approach was later extended to handle network transmission [START_REF] Rusinkiewicz | Streaming qsplat: a viewer for networked visualization of large[END_REF].

Pfister et al. [START_REF] Pfister | Surfels: surface elements as rendering primitives[END_REF] create an octree representation in their Surfels method. Surfels [START_REF] Pfister | Surfels: surface elements as rendering primitives[END_REF], surface splatting [START_REF] Zwicker | Surface splatting[END_REF] and point-set surfaces [START_REF] Alexa | Point set surfaces[END_REF] also address, among other aspects, high-quality rendering issues.

Inspired by techniques related to geometric coding [START_REF] Devillers | Geometric compression for interactive transmission[END_REF], a more compact representation is presented by Botsch et al. [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF], in which an octree is used to code point positions implicitly. This representation is based on an octree cell. If the cell hits the surface, then it is marked. All marked cells are then split up to some level. The maximal depth of the octree has to be known in advance. For each non-leaf cell of the octree, an eight bit code is stored: the childhood code. Each bit of this code corresponds to a leaf cell. The bit is set if the cell is marked. Using this technique, marked cells are stored as eight bits codes and empty cells come for free. It is important to note that sample positions are coded implicitly using this bit-code mechanism.

As a result, this approach is extremely memory efficient, and the authors present a very rapid rendering algorithm. In particular, the rendering algorithm factorizes transformations during hierarchy traversal, thus achieving very fast point projection for display. This work has inspired our approach and our choices for both storage and rendering.

The flexible rendering approach we present here combines the advantages of hierarchical and flexible rendering of QSplat [START_REF] Rusinkiewicz | Qsplat: a multiresolution point rendering system for large meshes[END_REF][START_REF] Rusinkiewicz | Streaming qsplat: a viewer for networked visualization of large[END_REF] and the compact and efficient storage of Botsch et al. [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF]. In particular, the compact nature of the data structure allows us to store and render very complex models, enabling high-quality rendering with limited resources, and the flexibility of the approach facilitates efficient rendering.

Point Sampling and Storage

We next explain how we create our packed hierarchical structure, and its cost in memory. As mentioned previously, we allow flexible rendering which requires the storage of intermediate sample attributes, such as normals and colors, at interior nodes of the hierarchy. We will examine the memory tradeoffs which result from this choice. This section is dedicated to a general study of storage cost for any value of ρ and σ. These values will be fixed in section 4.

Structure Creation

3D objects are often given in their traditional representation as vertices and polygons (faces). Since we need to generate our hierarchical structure from this input face set, we sample the geometry. Our sampling strategy is simple and general: for each cell of the hierarchy where geometry is present, whether interior or leaf, we sample the object at the point of the object closest to the center of the cell. If the cell hits no geometry, the cell is flagged as empty. The position of the sample point is not stored since it is implicitly given by the cell position in the hierarchy, and attributes such as normal and material index are coded using a fixed number of bits, noted σ.

In what follows we use a 16-bit code for a 13-bit quantized normal index (as in [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF]) and for indexing 8 materials. Evidently, we need more bits if we have more materials or colors.

Storage Requirement Analysis Number of non-leaf cells in a ρ-grid

To study the cost of our hierarchical structures, we separate the structure into two parts: the hierarchy without the leaf cells and the leaf cells attributes. Depending on the chosen value of ρ, we need to compute the number of interior cells, i.e., cells that are not leaf cells. The set of leaf cells are stored in an array of sample attributes and will be the same whatever the hierarchical structure chosen.

The hierarchical representation is especially well suited for grids with a small number of non-empty cells. In the typical example of a plane in space, only a small subset of the cells are not empty. For the octree, at any level, the number of non-empty subcells for an axis aligned plane is 4 = 2 2 , while the number of subcells in an octree cell is 8 = 2 3 . For the general case of a ρ-grid, the number of non-empty subcells for a plane is ρ 2 .

In the following, we suppose that, on average, we have ρ 2 non-empty subcells per non-leaf cell. As shown in the results and statistics (Section 5), this assumption is valid for many surfaces at a fine enough level.

Let n be the number of leaf samples, and k the depth (i.e., the number of subdivision levels) of the structure, then we shall have ρ 2k = n. The number N of non-empty interior cells is thus given by:

N = k-1 ∑ i=0 ρ 2i = n -1 ρ 2 -1 (1)
Note that this number of intermediate samples in the structure depends on the number of leaf samples n and on 1/ρ 2 . Thus, we expect to have fewer intermediate samples for the tri-grid than for the octree.

Storage Requirement for single level rendering The efficient encoding scheme introduced by Botsch et al. [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF] can be extended to general ρ-grids, using ρ 3 bits for each non leaf cell childhood code. The memory cost C of this representation in bits, without the cost of intermediate samples attributes (such as normals or colors), is the following:

C = N * ρ 3 = n -1 ρ 2 -1 ρ 3 (2)
The cost c 0 per leaf sample is thus given as (Eq. ( 2), ( 1)):

c 0 = C n = n -1 n * ρ 3 ρ 2 -1 ≈ ρ 3 ρ 2 -1 (3)
Note that the cost C accounts for the hierarchy of cells. The total cost of the structure also includes the cost of the leaf points attributes (normals/colors) which is the same whatever the hierarchical representation.

Additional storage for intermediate samples If no intermediate samples are stored in the interior nodes of the hierarchy, the entire hierarchy needs to be traversed since the leaf cell attributes need to be accessed. This can be prohibitively expensive, since traversing the hierarchy can be much more expensive than projecting and displaying a sample representing an intermediate node.

We thus store attributes (e.g., normals or colors) for intermediate samples at each interior hierarchy cell. These intermediate sample attributes are used to efficiently render a whole branch of the tree as a single sample. The cost per non-leaf cell is thus increased by the cost of an attribute: σ. The cost of an internal node is now ρ 3 + σ. The storage cost per leaf sample c σ is thus: The bold numbers are the most efficient structure for a given σ. σ=0 corresponds to no sample attribute, as in [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF].

c σ = ρ 3 + σ ρ 2 -1 (4)
The structure can be stored as a bitstream: for each cell, we store the attribute, the code of child cells (zero for empty cell, and one for non-empty cell) using ρ 3 bits, then we store the non-empty child cells. Note that leaf sample attributes can be interleaved in the structure providing cache coherence when accessing the bitstream.

Runtime Storage for Flexible rendering

To achieve flexible rendering, random access to the hierarchical structure is required. In particular, each cell needs to be accessed in any order.

We use the following structure to provide such flexibility: for each interior cell we store the intermediate attributes (normals and colors), and a pointer to a table of child cells. The size of the child cells table and random access to it can be obtained using the childhood code of the cell (see Section 2). This pointer to the table requires an additional 32 bits of storage per non-leaf cell. Using this approach, the traversal is completely flexible, a whole sub-branch of the tree can be skipped and child cells can be traversed in any order.

In this configuration, the memory requirement for such a structure is given Table 1 on the bottom rows. For an attribute size of 16 bits, which is the coding we use here, σ is 48 (32 for the pointer and 16 for the attributes).

All the tables of child cells are contiguously stored in a large table so that there is no additional cost of storage for the tables themselves. The pointer to a table can be seen as an index to the large table.

Flexible Rendering

We extend the rendering algorithm of Botsch et al. [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF] to enable flexible rendering. We reuse the same concept of transformation factorization, but introduce a test for each cell whether it should be rendered at this level as a splat or a single point (splat condition), or whether sub-cells should be rendered. This extension can also be seen as the adaptation of the QSplat [START_REF] Rusinkiewicz | Qsplat: a multiresolution point rendering system for large meshes[END_REF] rendering algorithm to this new data structure. This is possible because we can access the sub-cells independently from one another and for any level of the hierarchy. The pseudo-code is as follows:

Render (cell,

Intermediate Node Rendering and Splats

To render intermediate nodes, we compute a conservative approximation of the projection of a cell (a cell being an axis aligned cube), and render a rectangular splat. The user selects a splat size s in pixels; in all examples presented here s is equal to 1. A larger s results in faster frame rates. We compute the screen-space bounding rectangle of the cell. This is done by computing a bounds table in a manner similar to cell centers displacements [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF]. The screen-space bounds are given as (min x , min y , max x , max y ); we also compute min z/w , which is the minimum homogeneous depth of the bounding box. See Figure 2 for different values of the splat size.

For each level, if the extent of the projected cell bounding box is less than the splat size (or a pixel), we draw the intermediate sample. As with previous point-based approaches we render rectangular splats to represent these samples. Denote dx = (max x -min x ) the extent in x and s the splat size.

The test to determine whether we stop at this level is performed in homogeneous coordinates for each dimension (dx and dy), which costs two multiplications and two tests, as shown next:

dx w ≤ s, w > 0 (5)
which is equivalent to dx ≤ ws, w > 0 [START_REF] Pfister | Surfels: surface elements as rendering primitives[END_REF] This flexible rendering also permits view-frustum culling. If the conservative projection of the cell is out of the screen bounds, or the cell is completely behind the viewer, then the cell is not rendered at all.

Shadows

Our approach is well adapted to efficient shadow map computation. If we use the standard shadow-map algorithm, for each pixel of the image (with depth), we transform to the light source coordinate system using a 3x3 matrix. This operation requires 9 multiplications, 6 additions, and a test in the shadow map per pixel, without shading. For a typical iPAQ screen, this step would require 690,000 multiplies, and 540,000 additions or tests.

For directional light sources, we can modify our flexible rendering algorithm to compute shadows in a single pass. Given the light direction, we can use a strict ordering in the rendering of the hierarchical structure. As shown in Figure 3, we can define an ordering on subcells so that cells are rendered in front to back order. This idea has been proposed for the depth buffer by Botsch et al. [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF], but could not be implemented using their rendering algorithm, as the authors mention.

This ordering is precomputed once for the ρ 3 sub-cells for a given light source direction by projecting subcell centers along this direction. By traversing the hierarchy in this order, we guarantee that lit samples are rendered first, and samples in shadow are rendered afterwards. Thus, we compute the shadow map and the view from the camera at the same time, and perform the depth comparisons as we render. Note that computing both positions requires an additional displacement table, and an additional temporary position. With this technique, we avoid the expensive matrix multiplication for each pixel of the screen and interleave the rendering from the light source and rendering from the viewpoint into a single pass. This also results in better cache coherence. The performances of this one-pass approach are given in Table 3.

It is important to note that these algorithms (multi-level, splatting and shadows) are given for general ρgrids and could be applied to any specific type of grid, the octree being one of them. Our choice of the tri-grid for both storage and rendering is justified by Tables 1 and3. The tri-grid is more efficient than the octree in memory consumption when storing intermediate attributes. From Table 1, 4 or 5-grids may seem more appropriate. The drawback of these structures however is the fact that the "jump" between levels becomes problematic since an interior cell is replaced by 64 or 125 children resulting in a large amount of overdraw, and in very visible popping artifacts if the splat size is more than one pixel. Different grid types could be used within the same runtime environment with a minimal overhead. The value ρ = 3 has given good results both in terms of runtime storage (σ = 48), and in terms of rendering quality/performances.

Implementation and Results

In this section we present the results of our data structure construction and point-based rendering. The trigrids and octrees are computed on a workstation, and transferred to the iPAQ. The iPAQ we use has a 200Mhz processor with 64 Megabytes of main memory. In practice however we have about 8-9 Mbytes available for our data structure, with the graphics mode enabled.

We have used six models, four from the Stanford database (Buddha, Dragon, Blade and Lucy), and two tree models. All statistics are for 16 bit attribute samples, that is 13 bits for normals and 3 bits for a material index (but this layout is not restrictive). The disk storage requirements are close to the ones expected (σ = 16), as well as the run-time storage requirements (σ = 16 + 32 = 48). We also implemented the algorithm on a PC using floating point arithmetic and fixed-point arithmetic. All implementations are in C++.

We have included an electronic version of a movie showing our system in use on the PDA1 .

Implementation issues

Most mobile platforms do not have floating point units. As a result direct porting of graphics code will have significantly reduced performance. In addition, certain complex instructions may not be implemented; for example a division can be more expensive than a look-up in a precomputed table. This was confirmed on the platform we use. We implemented fixed point arithmetic for the purpose of our renderer, i.e., multiplications and additions. We also implemented an approximation of the inverse up to a given precision using a look-up table. Since we use 32-bit fixed point numbers, a global look-up table could not be used. We use a shift on numbers instead, with table size depending on expected precision.

Note that even on platforms equipped with FPUs, rendering is more efficient using fixed-point arithmetic rather than floating-point arithmetic. This is due to the final conversion from screen coordinates in floats to array indices as integers. This final step is a simple shift with fixed-point values, but is more expensive in floating-point.

Basic Statistics

In Table 2 we present a set of statistics on memory usage of our data structures. The basic statistics for the three main data structures are shown. We show statistics for two levels of both the octree and the tri-grid, for a set of models. The structures are evidently different; however we have chosen subdivision levels which are more or less equivalent for each. For example, 4 levels of subdivision for the tri-grid of the dragon model results in 143,077 samples, and 7 levels of the octree in 158,737. These subdivision levels gives structures which can fit in the memory of the iPAQ we used in our experiments. Subdivision level 10 for the octree and 6 for the tri-grid are used in the same manner for visualization on a PC. We show both memory usage in megabytes, and bits per sample in parentheses. We first show statistics for the recursive grid structure without intermediate attributes (NIS), which cannot be used for flexible rendering. In this case, the octree (which is the same as [START_REF] Botsch | Efficient high quality rendering of point sampled geometry[END_REF]) is on average 4% more efficient than the trigrid in terms of bits used per sample. However, once intermediate attributes are introduced (IS), permitting flexible rendering, the tri-grid is on average 12% more efficient in terms of memory. These numbers confirm the predictions presented in Table 1. A final compression pass using general purpose compression schemes (such as gzip) can be used. However, the compression results were inconsistent from one model to the other and from one level to the next.

Model

More interestingly, we also examine the memory usage of the structure when we include the additional pointers required at runtime for flexible rendering; the octree requires 33% more memory per sample, which is significant.

This packed hierarchical representation also permits visualization of huge geometries with a reasonable number of samples; such as the Lucy model from the Stanford Michelangelo project. This model originally contains more than 28 million triangles that do not fit in main memory, and is stored in a tri-grid at level 7 with more than 53 million point samples and attributes, for around 215 MB of RAM. This multi-level data structure permits an efficient storage of the geometry with a small approximation error (1.5x10 -4 ) of the entire model in memory. The data structure is constructed level-by-level in an out-of-core preprocessing step.

Rendering Performance and Shadows

In Table 3 we show statistics for our rendering algorithm using the tri-grid structure. The statistics are average frames per second on the iPAQ. An example of the kind of view chosen is shown in the top row ("far view") of Fig. 5. As expected, overall rendering all of the points is more expensive than the multi-level display, especially Undersampling problems are evident at level 4 subdivision when using points only. At level 5, the increase in frame rate is significant. All fps on the iPAQ. at maximum subdivision level 5. Nonetheless, for level 4 subdivision, the expense of computing the splats can outweigh the gain from multi-level rendering (Buddha(4) and Blade(4) in Table 3). Note, however, (Fig. 5 upper left) that Buddha at level 4 has visible undersampling artifacts with pure point rendering.

In Fig. 5, we show the problem of undersampling for a close view of a model, which is solved using splats. In addition, we show the effect on frame rate of multi-level display. Recall that multi-level display is possible because we store intermediate samples.

In the case of our largest model (Lucy), we run two versions of the rendering algorithm: flexible and not flexible for level 7 of the tri-grid with more than 53 million samples, on a window size of 512x512.

Model

Points Splats Shadows Buddha ( 4 3: Statistics for our rendering algorithm (in frames per second). "Points" corresponds to direct display of points, "Splats" refers to the flexible rendering algorithm with splatting, and "Shadows" to the one-pass shadow algorithm.

Conclusion

In this paper we have presented a general framework for point-based sampling of complex models, which we call ρ-grids. We have performed analysis of the storage requirements of these structures, and we have shown that the tri-grid is the most memory efficient when we store sample attributes for intermediate nodes.

We have used the tri-grid structure for efficient rendering of very complex models on mobile devices with limited computation and memory capabilities and a small display. We have introduced a flexible multi-level rendering algorithm which allows efficient rendering when zooming into the model, and an efficient one-pass shadow mapping algorithm. Our implementation on a 200 Mhz Compaq iPAQ PDA allows the display of 1.3M point-based representation of a 4.7 million polygon model at 2.1 frames per second.

In future work, we will be examining better reconstruction algorithms, and a progressive transmission approach which will take into account network parameters. We will also investigate hybrid rendering combining different representations such as polygons, lines and points, taking into account the specific constraints of the different platforms we consider. Also, in the case of manifold objects, adaptive sampling based on surface curvature could be investigated. This might result in better compression and rendering for such objects.

Figure 1 :

 1 Figure 1: (Left) Our structure used for hierarchical point-based rendering of a 4.7M polygon model at 2.1 fps on a 200Mhz IPAQ, sampled at 1.3 M points total. The multi-level approach restricts the number of points rendered depending on the view. Notice the shadows in the scene. (Right) The dragon model at 2.3 fps.

Figure 2 :

 2 Figure 2: Using splats for rendering. Left, two pixels splats, Right, five pixels splats.

Figure 3 :

 3 Figure 3: Ordering for shadow map rendering.

Figure 4 :

 4 Figure 4: (Left) The "Big" (see Table 3) scene rendered on the iPAQ. (Right) The Lucy model rendered as a 2048x2048 image with level 7 of the tri-grid on a PC (53 million points).

Figure 5 :

 5 Figure 5: Above, a far view of the Buddha model shown (left) with points at level 4 (4.15 (frames per second (fps)) (middle) with points at level 5 and (0.63 fps) (right) with splats at level 5 (2.85 fps). Below, a close view of the Buddha model shown (left) with points at level 4 (4.15 fps) (middle) with points at level 5 and (0.63 fps) (right) with splats at level 5 (1.42 fps). Undersampling problems are evident at level 4 subdivision when using points only. At level 5, the increase in frame rate is significant. All fps on the iPAQ.

  The nonflexible version renders an image in around 4 seconds, and the flexible version for a splat size of one pixel and shadows renders 7.9 frames per second, on a Pentium IV Xeon 2.0 Ghz PC under Windows 2000. 0.55 fps 2.28 fps

Figure 6 :

 6 Figure 6: Quality and speed of shadows. View of the Dragon model shown (left) with points at level 5 with shadows and (right) multi-level display with one-pass shadows multi-level. Notice that shadows generated with multi-level rendering are practically indistinguishable from shadows generated with full point resolution.

Table 1

 1 summarizes the different costs per leaf sample for various values of ρ and σ. Note that with no intermediate attributes (σ = 0), the octree is the most memory-efficient structure, but for intermediate attributes with size greater than 4 bits, the octree is no longer the best structure. The figures given in Table1are based on equation 4. Experimental results (see Table2) confirm these figures.

	σ ρ = 2 ρ = 3 ρ = 4 ρ = 5
	0	2.66	3.38	4.26	5.21
	4	4	3.88	4.53 5.375
	8	5.33	4.38	4.8	5.54
	12 6.67	4.88	5.06	5.71
	16	8	5.38	5.33	5.88
	32 13.33 7.38	6.4	6.54
	48 18.67 9.38	7.47	7.21

Table 1 :

 1 Variation of ρ-grid storage cost per leaf sample with ρ and the number of bits per sample attribute (σ).

Table 2 :

 2 Basic statistics. The results are presented for the octree at level 7 (oct/7) and level 10 (oct/10), and for the trigrid at level 4 (tri/4) and 6 (tri/6). Samples: the number of leaf samples; NIS: bytes on disk without intermediate attributes; IS: bytes on disk with intermediate samples; Mem: size of in-core structure including pointers and intermediate samples. For all columns, bits per leaf sample are shown in parenthesis.

		samples	NIS (bps)	IS (bps)	mem (bps)
	Lucy (oct/7)	91,207	206 kB (18.5)	263 kB (23.6)	400kB (36.37)
	Lucy (tri/4)	82,078	191 kB (19.1)	210 kB (21.0)	270kB (27.16)
	Lucy (oct/10)	5,708,675 12.7 MB (18.7) 16.4 MB (24.2) 25.6MB (37.73)
	Lucy (tri/6)	6,495,339 15.0 MB (19.5) 16.6 MB (21.5) 21.9MB (28.28)
	Buddha (oct/7)	125,283	284 kB (18.6)	365 kB (23.9)	550kB (37.00)
	Buddha (tri/4)	112,924	265 kB (19.2)	290 kB (21.1)	370Kb (27.40)
	Buddha (oct/10) 7,286,757 16.4 MB (18.8) 21.3 MB (24.5) 33.6MB (38.76)
	Buddha (tri/6)	8,242,567 19.3 MB (19.7) 21.4 MB (21.8) 28.5MB (29.04)
	Blade (oct/7)	248,066	560 kB (18.5)	711 kB (23.5)	1.06MB (35.95)
	Blade (tri/4)	223,123	517 kB (19.0)	565 kB (20.8)	710kB (26.62)
	Blade (oct/10)	14,925,097 33.4 MB (18.8) 43.3 MB (24.4) 68.1MB (38.29)
	Blade (tri/6)	16,929,466 39.5 MB (19.6) 43.8 MB (21.7) 57.9MB (28.69)
	Dragon (oct/7)	158,737	361 kB (18.7)	464 kB (24.0)	700kB (37.25)
	Dragon (tri/4)	143,077	338 kB (19.3)	372 kB (21.3)	480kB (27.90)
	Dragon (oct/10)	9,272,277 20.5 MB (18.6) 26.8 MB (24.2)	43MB (38.93)
	Dragon (tri/6)	10,314,218 24.2 MB (19.7) 26.9 MB (21.9) 35.8MB (29.17)

This is a DIVX 5.0.2 avi movie with a soundtrack.

Acknowledgements

The authors wish to thank Alexandre Olivier-Magnon for modelling the tree, using Maya, thanks to a donation of Alias|Wavefront. The Buddha, Bunny, Blade, Dragon and Lucy models are from the the Stanford 3D Scanning Repository http://graphics.stanford.edu/data/3Dscanrep. Thanks to Marc Stamminger and Frédo Durand for their insightful comments on an early draft.