{'fraction_non_alphanumeric': 0.03167191677278583, 'fraction_numerical': 0.015392864323853355, 'mean_word_length': 4.496234446627374, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Post-Newtonian Lagrangian of an N-body System with Arbitrary Mass and Spin Multipoles

The present paper derives the post-Newtonian Lagrangian of translational motion of N arbitrary-structured bodies with all mass and spin multipoles in a scalar-tensor theory of gravity. The multipoles depend on time and evolve in accordance with their own dynamic equations of motion. The Lagrangian is retrieved from the post-Newtonian equations of motion by solving the inverse problem of the Lagrangian mechanics and generalizes a well-known Lagrangian of pole-dipole-quadrupole massive particles to the particles of higher multipolarity. Analytic treatment of the higher-order multipole contributions is important for more rigorous computation of gravitational waveform of inspiralling compact binaries at the latest stage of their orbital evolution before merger when tidal and rotational deformations of stars are no longer small and rapidly change in time. The Lagrangian of an N-body system with arbitrary mass and spin multipoles is instrumental for formulation of the post-Newtonian conservation laws of energy, momenta and the integrals of the center of mass.

Modern astrophysics has reached a significant progress in understanding the nature of ultra-relativistic objects like neutron stars and black holes through precise and accurate measurements of compact binary systems conducted in different bandwidths of electromagnetic spectrum spreading out from radio to gamma-ray wavelength [1][2][3][4]. Gravitational wave astronomy had opened a new window to study dynamics of the binary systems with compact stars having fast intrinsic rotation and ultra-short orbital periods [5]. Such binary systems consist either of two neutron stars or a neutron star (NS) and a black hole (BH) or two black holes which are significantly deformed by centrifugal and tidal forces. Moreover, the gravitational radiation-reaction force causes a rapid secular decrease of the orbital period and separation between the two bodies while the intrinsic angular momentum (spin) of each body is subject to wide precession which changes the rotational and tidal deformations of the bodies in the course of time. Precise measurement of these deformations by means of gravitational wave detectors is a key to more profound understanding of the interior structure of compact astrophysical objects and the equation of state (EoS) of hadronic matter at super-critical density and pressure of the degenerate Bose condensate of baryons [6] and/or quark-gluon plasma [7] just before the merger. It also provides a deeper insight into non-linear geometrodynamics of binary black hole [8,9].
Information about the equation of state and other fundamental physical parameters of matter of the tidallydisturbed bodies in the binary system (like the Love numbers, viscosity, shear stress, etc.) is mapped through the gravity field equations to a set of tidally-induced multipole moments and their time derivatives [10,11]. There are two, physically-meaningful types of the multipole moments in general relativity called mass and spin multipoles [12]. Scalartensor theory of gravity contains one more type of multipoles called scalar-field multipoles [13,14]. The multipoles parameterize the equations of orbital evolution of the inspiralling binary systems. The mass multipole moments appear already in the Newtonian approximation while the spin and scalar field multipoles come up for the first time at the first post-Newtonian (PN) approximation. Translational equations of motion of extended, arbitrary-structured bodies in an N -body system with the full account for the intrinsic mass and spin multipoles of the bodies have been derived in general relativity [15,16] and in scalar-tensor theory of gravity [14,17,18]. These equations extend the domain of analytic applicability of equations of motion of pole-dipole massive particles in higher PN approximations from inspiral to coalescence stage of the orbital evolution of compact binary systems.
General physical characteristics of a gravitational N -body system can be, in principle, understood without explicit integration of equations of motion if the first integrals of the equations are known. The first integrals are functions of dynamic variables and their time derivatives in the configuration space of the N -body system. The dynamic variables are coordinates of the center of mass and multipole moments of the bodies which are functions of time. The first integrals can be obtained by applying the Nöther theorem to the Lagrangian of an N -body system [19]. The knowledge of the analytic form of the first integrals is important for formulating the conservation laws of energy, angular and linear momenta, and the integrals of the center of mass of the N -body system. Knowledge of their analytic form is critical for minimization of the errors of numerical integration of the equations of motion of inspiralling compact binaries when computing the tidal and tidal-resonant effects [20] as well as for refining the statistical estimates of the parameters of gravitational waveforms [21].
In principle, the Lagrangian of an N -body system can be derived directly from the Hilbert action of general relativity (or any other, viable alternative gravity theory) as proposed by Infeld [22] who extended similar idea of Fokker [23] from electromagnetism to general relativity. Fokker-Infeld method was applied by Brumberg [24] who derived the post-Newtonian Lagrangian of an N -body system consisting of rigidly rotating and oblated bodies possessing the intrinsic spin-dipole and mass quadrupole moments. The Fokker-Infeld Lagrangian of an N -body system was also derived for structureless (non-spinning) point-particles 1 in the second PN approximation of general relativity with different techniques [25][26][27][28][29][30][31]. Similar method was used to derive the Arnovitt-Deser-Misner (ADM) Hamiltonian of point-like particles of an N -body system in the third [32,33] and forth PN approximations [34][35][36]. Corresponding ADM Hamiltonian of an N -body problem for spinning point-like particles has been derived by Jaranowski and Schäfer [37], Schäfer and Jaranowski [38] who employed the regularization technique to deal with the formally-divergent integrals from distributions which are used to model the point-like bodies of the N -body system [39].
An independent method to get the Lagrangian L of an N -body system is based on solving the inverse problem of the Lagrangian mechanics by writing down a system of self-adjoint differential equations in partial derivatives for, yet unknown, function L such that its solution, if it exists, yields the variational derivative corresponding to the dynamical equations of motion of the N -body system. The Lagrangian is defined not uniquely but up to a total time derivative from an arbitrary scalar function [19]. Mathematical foundations justifying the existence and uniqueness of the solution of the inverse problem for a system of ordinary differential equations of the second order have been worked out by Douglas [40], significantly rectified and improved by Santilli [41] and further refined and enhanced by other researchers -see, for instance, [42,43] and reference therein. Solution of the system of the self-adjoint differential equations relies upon a rigorous mathematical background but is hardly applicable for practical computation of the Lagrangian of a gravitational N -body problem in the post-Newtonian approximations due to the complexity of the equations. Technically, it is much simpler to look for the Lagrangian in the form of a linear combination of some scalar functions with yet undetermined numerical coefficients, then, to take the variational derivative from this combination and equate it to the original equations of motion of the N -body problem. Identification and comparison of similar terms yield a system of algebraic equations for the above-mentioned numerical coefficients which can be solved by applying a standard linear algebra technique. If such solution exists it immediately yields the post-Newtonian Lagrangian of the N -body system. This "linear algebra" approach to solving the inverse problem of the Lagrangian mechanics was applied in general relativity to build the Lagrangian of an N -body system made up of pole-dipole particles [30,33,44,45].
We notice that in spite of the enormous success in developing the post-Newtonian approximations for solving twobody problem of point-like masses [46] the post-Newtonian Lagrangian of a gravitational N -body system consisting of arbitrary-structured, rotating and dynamically-evolving extended bodies has not been formulated so far. The problem is that such a Lagrangian should contain contributions from all mass and spin multipoles of the bodies and match analytically with the Lagrangian of the point-like pole-dipole particles. Some ideas of how to fulfill this task have been proposed at the Newtonian level in [47,48]. Brumberg [24] had employed the Fokker-Infeld approach combined together with Fock's method [49] in order to get the post-Newtonian Lagrangian for particles endowed with dipole and quadrupole moments defined in the global coordinates of an N -body system. These multipoles are not directly measurable and must be converted to the body-adapted local coordinates in order to get physically meaningful results [31,[49][50][51]. Definite progress was made in deriving the Lagrangian for particles with spin [39,[52][53][54][55] and quadrupole moments [56][57][58][59]. However, we are not aware of any other attempt in scientific literature aimed to build the post-Newtonian Lagrangian of an isolated gravitating system comprised of N arbitrary-structured extended bodies possessing all time-dependent mass and spin multipole moments beyond the quadrupole order. Therefore, the present paper focuses on solution of this problem in the first PN approximation of scalar-tensor theory of gravity.
We find out the Lagrangian by solving the inverse problem of the Lagrangian mechanics with the post-Newtonian equations of motion of N bodies derived in our previous paper [18]. Section II is a brief introduction to the post-Newtonian formalism of the local and global coordinate charts used for description of the dynamics of an astronomical N -body system. We also define in this section the post-Newtonian multipole moments as well as the mass and center of mass of an extended body from the N -body system. Section III explains the configuration space of an N -body problem, describes the post-Newtonian equations of motion of arbitrary structured bodies and discusses the method of solving variational equations of the inverse problem of the Lagrangian mechanics. Section IV describes preparatory transformations of the gravitational force from its original form given in section III to that which is required for retrieving the Lagrangian by solving the inverse problem. Section V brings the Newtonian terms in the equations of motion to the Lagrangian form. The post-Newtonian gravitational force is transformed to the Lagrangian form in section VI. A general formula for the Lagrangian of an N -body system is derived in section VII. Section VIII specifies the generic Lagrangian of the previous section to the case of a binary system and analyzes its structure in monopole, dipole, and quadrupole approximations. It also compares our Lagrangian in the quadrupole-particle approximation with the results of previous researchers. Conclusions are given in section IX. Appendix A explains the virial theorems which are used for interconnecting the external and internal dynamic variables in the configuration space of an N -body system in the process of derivation of the Lagrangian.
x α = {t, x i } are the global (single chart) coordinates covering the entire spacetime in a single chart.
-∂ i = ∂/∂x i is a partial derivative with respect to a spatial coordinate x i , -the multi-index partial derivatives with respect to coordinates x α are denoted as -the round parentheses embracing a group of tensor indices denote full symmetrization, for example, T (αβγ) = 1 3! (T αβγ + T βγα + T γαβ + T βαγ + T αγβ + T γβα ), etc. , -the square parentheses around a pair of tensor indices denote anti-symmetrization, for example, T [αβ]γ = 1 2 T αβγ − T βαγ ,etc.
-the angular brackets around spatial tensor indices denote a symmetric trace-free (STF) projection. For example, the STF projection T <L> of tensor T L is constructed from its symmetric part S L = T (L) = T (i1i2...i l ) by subtracting all the permissible traces [12,60] where [l/2] is the largest integer less than or equal to l/2.
Other notations will be introduced and explained in the main text of the paper as they appear. In particular, in most equations of the present paper we use a geometic system of units in which the fundamental speed, c, and the universal gravitational constant, G, are set equal to unity.

II. COORDINATE SYSTEMS
We consider an isolated astronomical system consisting of N extended, arbitrary-structured bodies which are bound by attractive gravitational forces. We label the bodies with the capital letters A,B,C each of which takes values from the set of integers {1,2,...,N}. Each body A is characterized by yet unspecified but physically-admissible distribution of mass density ρ, internal thermodynamic energy Π and stress tensor σ ij . The bodies are assumed to be well-separated in space and consist of matter occupied a finite volume with a distinctive boundary. The present paper does not admit physical singularities of spacetime like black holes and assume that the bodies have weak self-gravity and move sufficiently slow in space. This is the realm of the PN approximations which are applied to solve the gravity field equations in relativistic celestial mechanics of the solar system [50,51], binary pulsars [1] and for calculating templates of gravitational waves emitted by inspiralling compact binaries [46,61].
An adequate mathematical description of dynamics of an N -body system requires introduction of N+1 coordinate systems [62][63][64]. The primary coordinate chart, x α = (t, x i ), is global and covers the entire spacetime manifold while a local coordinate chart is introduced in the neighborhood of world line of each body A∈ {1, 2, ..., N}. From the point of view of the Lagrangian mechanics the global coordinates are necessary to parameterize the orbital motion of each body A in the N -body system in terms of coordinates x i A , velocity v i A =ẋ i A and acceleration a i A of its center of mass which are treated as external dynamic variables in the configuration space. The local coordinates, w α A = (u A , w i A ), are introduced for each body A in order to describe the temporal evolution of its mass and spin multipole moments considered as internal dynamic variables. The global and local coordinates of body A overlap in a neighborhood of the body and can be transformed one to another [50,51,65] so that the local coordinates can be represented as function of the global coordinates, w α A = w α A (t, x), and vice versa, In what follows, we shall drop the body label A attached to the local coordinates of the body and identify w α ≡ w α A if there is no confusion with the local coordinates of another body B =A. In all other cases, the local coordinates are tagged with the label of the corresponding body.
The global coordinates specify the orbital motion of the center of mass of each body. The local coordinates characterize the changes in the orientation of the body and its rotational and tidal deformations which cause temporal changes of the multipole moments of the body as the body moves along its orbit. Conversely, the intrinsic changes of the body's multipole moments affect the body's orbital motion. It is desirable to separate the body's internal and external degrees of freedom as much as possible. It can be achieved through a suitable construction of the local coordinates adapted to each body and their subsequent splicing with the global coordinates by the method of matched asymptotic expansions.
We utilize a scalar-tensor theory of gravity where the scalar field Φ and metric tensor g αβ are solutions of the corresponding field equations of the theory -for more details see [18]. We assume that the bodies of an N -body system move slowly, and gravitational and scalar fields are weak everywhere. The field equations are solved iteratively by the method of post-Newtonian (PN) approximations [50] with a small parameter ǫ ∼ v/c ∼ (U/c 2 ) 1/2 where v is a characteristic velocity of matter, U is a typical value of gravitational potential, and c is the ultimate speed of gravity which is equal to the speed of light in vacuum [66,67]. The scalar field Φ and the metric tensor g αβ are expanded in Taylor series around their background values with respect to the small parameter ǫ. In particular, the post-Newtonian expansion of scalar field Φ is where Φ 0 is the background value of the scalar field and φ is a dimensionless perturbation of the field. The scalar field is assumed to be self-interacting and the strength of the self-interaction is given by function ω(Φ). This function is also expanded in a Taylor series with respect to the perturbation of the field, where ω 0 ≡ ω(Φ 0 ) is the background value of the self-interaction, and ω ′ 0 ≡ (dω/d ln Φ) Φ0 is its first derivative. The metric tensor g αβ is expanded in the post-Newtonian series around the background Minkowski metric, η αβ , where h αβ is the perturbation of the metric tensor. The perturbations φ and h αβ are of the order of ǫ 2 . Exact form of the perturbations φ and h αβ is given below in next sections. The residual terms in (2)-(4) can be found in [13,18,50] but they are not required in the present paper. Specific functional form of the perturbations, φ and h αβ , depends on the type of the coordinates introduced in spacetime of an N -body system. We distinguish the global coordinates covering the entire spacetime manifold from the local coordinates adapted to each body and covering a close neighborhood of the body all along the world line of its center of mass.

A. The global coordinates
In the global coordinates x α = (t, x i ) the post-Newtonian perturbation of the metric tensor reads where is the parameter of the scalar-tensor theory, and here and everywhere else the symbol O(n) denotes the post-Newtonian terms of the order of ǫ n . Parameter γ is formally similar but not physically-equivalent to the parameter γ PPN of the parametrized post-Newtonian (PPN) formalism [68] which is introduced in the PPN metric tensor formally to characterize the curvature of space-like hypersurface of constant time t. On the other hand, the covariant parameter γ defined in the scalar-tensor theory of gravity by (8), is a measure of the strength of self-coupling interaction of the long-range scalar field Φ given in terms of the background value ω 0 of function ω(Φ).
The other functions appearing in (5)-(7) include scalar (Newtonian-type) and vector gravitational potentials of the N -body system, which are linear superpositions of scalar and vector gravitational potentials of each body of the N -body system, and Here V A denotes the spatial volume occupied by body A, ρ * ≡ √ −gu 0 ρ is the invariant mass density [49], u 0 is time component of 4-velocity of matter u α , g = det|g αβ | is a determinant of the metric tensor, ρ is the local mass density entering the energy-momentum tensor, and v i = dx i /dt = u i /u 0 is 3-dimensional velocity of matter in the global coordinates. The invariant mass density ρ * is more convenient for calculations as its product with the volume element is invariant, ρ * d 3 x = inv, with respect to coordinate transformations as well as to the Lie transport along matter worldlines [49,50]. Scalar field perturbation is mathematically identical to the Newtonian gravitational potential [18], Nonetheless, the reader should keep in mind that the laws of transformation of the scalar field and the metric tensor are different. It requires to carefully distinguish transformations of the scalar potential in (9) from that in (11) which is a part of the metric tensor component g 00 , as they carry out information about geometrically-different types of the field. Subsequent derivation requires to single out one of the bodies, let say a body A, and split the gravitational potentials in two parts -internal and external, where U A and U i A denote the internal gravitational potentials produced by body A alone, denote the external gravitational potentials of all other (external) bodies of the N -body system but the body A, and denote the gravitational potentials of the external body B. We do not impose any specific limitations either on the distribution of mass and velocities inside each body of the N -body system or on the shape of the volume of integration. They remain arbitrary to the extent admitted by the gravitational field equations of the PN approximations. In particular, the mass density and velocity must obey the microscopic equations of motion of matter explained in Appendix A 1. The boundary of the volume integration V A of each body A is determined in accordance with the solution of the internal problem for this body in the local coordinates -see, for example, [69,70].

B. The local coordinates
There are N local coordinate charts each of which is adapted to a particular body in the N -body system [50,51]. By default, we denote the local coordinates of body A as w α = (u, w i ). If one needs to distinguish the local coordinates adapted to body A from those adapted to body B, we use a corresponding label, for example, w α A = (u A , w i A ) are the local coordinates adapted to body A, w α B = (u B , w i B are the local coordinates adapted to body B, and so on. In the local coordinates of body A∈ {1, 2, ..., N} the metric tensor perturbation takes on the following form [18], where Q L and C L are the external gravitational multipoles generated by the bodies being external to body A, P L are the external multipoles moments of the scalar field, and U A and U i A are the internal gravitational potentials generated solely by the matter of body A. The multipoles Q L and C L are also called gravitoelectric and gravitomagnetic multipoles as they are connected with the corresponding electric-type and magnetic-type components of the curvature tensor [71]. We did not show explicitly in the local metric (15)-(17) the residual post-Newtonian terms though the terms of the order of O(4) in h 00 are requited for calculation of the post-Newtonian equations of motion [14,15,18,64,72] and for defining the post-Newtonian mass multipole moments of the bodies -see section II C below. These terms can be found in [18].
The internal gravitational potentials are defined as integrals over volume V A occupied by the body, where ν i = dw i /du is velocity of matter of body A with respect to the local coordinates adapted to this body. Definitions of the internal potentials in (18) look similar to those given in (10) in the global coordinates. However, since the local and global coordinates are interrelated by the post-Newtonian transformation, the numerical values of the internal potentials U A (u, w), U i A (u, w) coincide with the numerical values of U A (t, x), U i A (t, x) only in the Newtonian approximation -see more details about the correspondence between the gravitational potentials in [50, §5.2.3].
The internal potentials are expanded in the multipolar series outside the body A as explained in detail in [73,74] where M L A are the mass multipoles, and S L A are the spin multipoles of body A, r = |w| is the radial distance from the origin of the local coordinates of body A to the field point w i , and the overdot denotes the total time derivative with respect to the local time u. The mass and spin multipoles are defined in terms of the matter variables below in section II C.
The metric tensor perturbation h αβ (u, w) in the local coordinates is a solution of the field equations of the scalartensor theory of gravity. Therefore, it must be consistent with the metric tensor perturbation h αβ (t, x) in the global coordinates as they represent exactly the same gravitational field written down in terms of different variables. This fact allows us to establish a transformation from the global to local coordinates of body A (and its inverse) by applying the technique of matched asymptotic expansions [75] to the metric perturbations. This technique also allows to find out the explicit form of the external multipoles Q L , C L in terms of partial derivatives from the gravitational potentials U (t, x),Ū i (t, x) of the external (with respect to body A) bodies of the N -body system [50].
In case of l ≥ 2 the gravitoelectric multipoles Q L represent physically a tidal gravitational field in the neighborhood of body A [18] where we have used multi-index notation for the STF partial derivative, ∂ L ≡ ∂ i1i2...i l . For l = 1, the gravitoelectric dipole moment [18] is the difference between the Newtonian gravity force and acceleration of the center of mass of body A, a i A =ẍ i A , where the overdot is time derivative with respect to time t. There is no monopole gravitoelectric multipole, Q ≡ 0, as it is eliminated by re-scaling of local coordinates.
The gravitomagnetic multipoles C L are [18] where we have used notation for the STF partial derivative, the overdot denotes the time derivative with respect to time t, and the square parentheses around a pair of indices denote anti-symmetrization. The external gravitomagentic dipole C i has a physical meaning of the dynamic angular velocity of rotation of spatial axes of the local coordinates [62]. We assume that the local coordinates of body A∈ {1, 2, ..., N} do not rotate dynamically, that is C i ≡ 0. This assumption is equivalent to the statement that the local triad of spatial axes of the local coordinates undergoes the Fermi-Walker transport along the world line of the center of mass of body A. The external multipoles of the scalar field Φ look similar to the gravitoelectric multipole, However, it should be emphasized that P L = Q L + O(2), if and only if, the index l ≥ 2. Monopole of the scalar field defines the value of the gravitational potential of external bodies on the world line of the center of mass of body A, and dipole of the scalar field, formally coincides with the Newtonian gravity force computed on the world line of the center of mass of body A. Dipoles P i and Q i are linked one to another by equation (22) that is, Post-Newtonian corrections of the order O(2) are required only in (21) for derivation of the post-Newtonian equations of motion. They can can be found in [14,18] but we don't need their explicit form in the present paper C. The internal mass and spin multipoles of a single body from an N -body system Like in general relativity, the metric tensor in scalar-tensor theory of gravity depends on two types of canonical STF multipoles -the mass multipoles M L A and the spin multipoles S L A . However, in contrast to general relativity there are two possible definitions of the multipoles affiliated respectively to the Einstein and Jordan frames of the scalar-tensor theory. It turns out that the Jordan-frame multipoles define the force of gravitational interaction between massive bodies in the equations of motion of an N -body system [14,18]. On the other hand, the principle of the effacing of internal structure of the bodies [76] from the equations of motion suggests that it is the Einstein-frame monopole and dipole moments which are to be used for definition of the inertial mass and the center of mass of each body. The mass multipoles in the Jordan frame are called active multipoles and those defined in the Einstein frame are termed conformal multipoles [50].
The active mass multipoles referred to the local coordinates of body A, read [18] where the active post-Newtonian mass and current density distributions [14] σ(u, w) = ρ * (u, w) 1 depend on the invariant density ρ * , velocity ν i , kinetic and potential energy density of matter as well as on trace σ kk = δ ik σ ik of the stress tensor σ ik of matter of body A. The stress tensor and the internal structure of bodies in the N -body system are constrained by equation of state which is not required to be specified for our calculations. Thus, the bodies can be either solid or liquid or made of some other kind of ordinary or condensed matter. Definition (28) of the mass multipoles includes post-Newtonian terms which come from the post-Newtonian terms O(4) in the metric (15) -see more details in [50,Chapter 4].
The active mass multipoles depend on two parameters of the scalar-tensor theoryγ and β. Parameter γ has been explained above in (8). Parameter β is a measure of the dependence of self-interaction of the scalar field Φ on its strength [18] where ω ′ 0 = (dω/d ln Φ) Φ=Φ0 is a derivative of the self-coupling function ω(Φ) with respect to the dimensional scalar field perturbation φ = Φ/Φ 0 − 1, where Φ 0 is the background value of the scalar field [18]. In the present paper β is formally similar to but not fully identical to the parameter β PPN of the PPN formalism [68] where it was introduced to the metric tensor phenomenologically to parameterize the possible deviation of the strength of self-interaction of the metric tensor perturbations from its general-relativistic value, β PPN = 1. Notice that since Φ 0 depends on time due to the expansion of universe the parameters β and γ are functions of time as well, and their secular evolution can be measured, at least in principle, with precise astronomical techniques like pulsar timing, lunar laser ranging, etc. [77].
The mass multipole (28) also depends on the gravitoelectric external multipoles Q L and the scalar field multipoles P L which are expressed in terms of the partial derivatives from the gravitational potential of external bodies. The reason for inclusion these terms to the definition of the multipoles is that it fully eliminates non-canonical multipoles (35) from equations of motion and validates the effacing principle [76]. It is also required by the law of conservation of energy in tidal interaction [78]. Making use of equations (21), (22) and (24) we can bring the mass multipoles (28) to the form that is more suitable for derivation of the Lagrangian of the N -body problem. It reads where are the Blanchet-Damour mass multipoles of body A [73] in scalar-tensor theory of gravity that are functionals of the active mass density distributioñ and non-canonical multipoles The symbol denotes a symmetric mass multipole of body A in the body-adapted local coordinates which is required only in the post-Newtonian terms. It is worth noticing some useful relationships between various multipoles, The mass multipoles of body A are functions of the local coordinate time u which can be also interpreted as the proper time of a fictitious observer placed at the origin of the local coordinates adapted to body A. It is useful to notice that the multipoles M L A and I L A differ by terms of the post-Newtonian order of magnitude and can be used interchangeably in the post-Newtonian terms. We prefer to use in the post-Newtonian terms the multipoles M L A like in (32). It should be emphasized that the mass multipoles (32) are not simply the integrals of the internal dynamic variables which appear in the definition of the Blanchet-Damour mass multipoles (33) and the non-canonical multipoles (35), (36). The active mass multipoles also depend on the external dynamic variables -coordinates x i A of the centers of mass of the bodies and their accelerations a i A . This dependence complicates derivation of the Lagrangian of the N -body problem since the partial derivatives from the mass multipoles with respect to the external variables do not vanish -∂M L A /∂x i A = 0 and ∂M L A /∂a i A = 0. Hence, the variational derivative from the active mass multipoles do not vanish either -see section V B.
The active multipoles of body B are defined similarly to those of body A. Namely, we have where are the Blanchet-Damour active mass multipoles of body B. The mass multipoles of body B are functions of the coordinate time of the local coordinates of body B. The local coordinate time of body A does not coincide with that of body B as the bodies move with respect to each other and are exposed to different gravitational potentials. We shall discuss the correspondence between the local times of the bodies in section IV A. In addition to the internal mass multipoles, M L A of body A, there exists a set of internal STF spin multipoles S L where ε kpq is the fully antisymmetric Levi-Civita symbol. The internal spin multipoles depend on the local coordinate time of the corresponding body. Definition (40) is sufficient for deriving the post-Newtonian translational motion of extended bodies in the N -body system. Derivation of the post-Newtonian equations of rotational motion of the bodies requires extension of (40) to the post-Newtonian domain. We don't discuss it in the present paper but details can be found in [18]. Both the mass and spin multipole moments of each body are subject to tidal and rotational deformations which progressively change as the bodies of the N -body system move in the course of time. The time evolution of multipole moments is determined by corresponding equations of motion for these moments. We assume that these equations of motion are known and can be solved, at least, in principle. In the present paper we will need three equations of motion which govern the time evolution of mass monopole, mass quadrupole and spin dipole. These equations are given in equations (A9), (A22) of Appendix A and in (70) respectively. D. Definition of the inertial mass and the center of mass of a single body from an N -body system There are two distinctive world lines inside the world tube covered by the local coordinates adopted to body A. They are representing a world line of the origin of the local coordinates and that of the body's center of mass. Naturally, the origin of the local coordinates adapted to body A has spatial coordinates w i = 0 which are identified with the point x i A in the global coordinates. The center of mass of body A can move along a different world line and, as a rule, the correspondence between the two world lines is established by imposing an additional vector-like supplementary condition [79].
We shall denote the local coordinates of the center of mass of body A as w i A . Definition of the center of mass in the post-Newtonian approximation depends on our choice of the mass density function describing the post-Newtonian distribution of matter inside the body. We want to find out the most optimal definition of the center of mass which eliminates as much spurious effects in the equations of motion of the body as possible. As mentioned above, in scalartensor theory of gravity there are two mass density distribution functions corresponding to the conformal (Einstein's frame) and active (Jordan's frame) multipole moments. Our study [14,17,18] shows that in the scalar-tensor theory of gravity the center of mass of body A is defined the most optimally if one uses the conformal mass density distribution inside the body, which depends on the invariant mass density ρ * = √ −gu 0 ρ, internal velocity of matter ν i = dw i /du, its internal energy Π, and gravitational potential U A of the body. The conformal mass density ̺ does not depend on the PPN parameters β, γ It is instructive to introduce a general-relativistic post-Newtonian mass of body A by formula where the general-relativistic mass density in the post-Newtonian approximation is [49,68] σ(u, w) = ρ * (u, w) 1 + 1 2 The general-relativistic mass (42) depends exclusively on the internal dynamic variables of body A and does not depend on the external variables like coordinates and velocity of the centers of mass of bodies in the N -body system. The post-Newtonian inertial mass of body A is defined as a conformal monopole moment of the body by the following equation [18,Eq. 164] where P =Ū (x A ) is the external monopole moment of the scalar field and Q L are the gravitoelectric multipole moments (21). Notice that the inertial mass (44) depends on the external dynamic variables through P and Q L . Position of the center of mass w i A of body A in the local coordinates adapted to the body is defined by the conformal dipole moment of the body [18], where (see [18,Eq. 170]) and a vector function I i A = I i A (t) characterizes an additional post-Newtonian displacement of the center of mass with respect to the origin of the local frame of body A which can depend on time. The quantity m A in (46) is the inertial mass of body A (44) and the STF non-canonical multipole, N L A , is given in (35a). The displacement vector function I i A should be specified by imposing a particular supplementary condition which choice is dictated by computational circumstances but, in principle, is arbitrary. Notice that in accordance to the computational scheme of derivation of the equations of motion, the point on spacetime manifold having spatial global coordinate x i A is always identified with the origin of the local coordinates adapted to body A (the point with w i = 0) independently of the value of the residual displacement I i A of the center of mass with respect to the origin of the local coordinates. We postulate that at any instant of time the conformal dipole moment vanishes at the origin of the local coordinates adapted to body A. In other words, the condition J i A ≡ 0 defines the point with the local coordinates w i = 0. The center of mass of body A, w i A , coincides with the origin of the local coordinates, w i = 0, if the residual displacement, I i A = 0. It is evident that the condition, I i A = 0, can be always fulfilled, at least, at some fixed instant of time. However, this condition should be propagated in time as derivation of the equations of motion requires to know behavior of I i A along the entire worldline of the origin of the local coordinates. It turns out that the condition I i A (t) = 0 can indeed be satisfied continuously. The proof was furnished in the Newtonian approximation in [63,65] and in the first PN approximation by Kopeikin and Vlasov [14], Damour et al. [72] in the local coordinates. Extension of this result to the global coordinates involves more comprehensive computations which have been conducted by Racine and Flanagan [15], Racine et al. [16] in general relativity and, independently, by Kopeikin [18] in scalar-tensor theory of gravity. These computations yield translational post-Newtonian equations of motion of extended bodies in the global coordinates with all multipoles taken into account and with the center of mass of each body lying at the origin of the local coordinates. We have discovered (see discussion in section IX.A.6 of [18]) that the constrain, I i A (t) = 0, may be not the most optimal one as a sufficiently large number of terms with second time derivatives entering equations of translational motion can be eliminated by choosing another condition, [18,Eq. 289]. This simplifies the equations of motion. For the time being we leave the vector function I i A unconstrained and discuss its specific choice in section V B more explicitly. The fact of the matter is that derivation of the Lagrangian from the equations of motion produces terms with the second time derivatives which are not reduced to the form of the Lagrangian (variational) derivative and must be associated with the specific choice of function I i A -see equation (148). It is worthwhile to notice that the conformal dipole moment J i A of body A is not equivalent to the active dipole moment M i A of the body. Thus, M i A = 0 even if J i A = 0 which means that extended bodies in scalar-tensor theory of gravity possess a non-vanishing active mass dipole moment as contrasted to general relativity [14,80]. This is why the lowest order of gravitational waves emitted by self-gravitating systems in scalar-tensor theory of gravity is dipole with multipole index l = 1 [68] while in general relativity it is quadrupole with l = 2 [81]. For the reader's convenience we provide a formula for the active dipole moment of body A in case when the conformal moment J i A = 0. This formula is derived by making use of a vector virial theorem (A15) in the definition of the active dipole moment and reads [18,Eq. 171] where η = 4β − g − 3 is called the Nordtvedt parameter [68], P i is the external dipole moment of the scalar field defined in (26), the local non-geodesic acceleration Q i is defined in (22), the external gravitoelectric multipoles Q L are shown in (21) and the non-canonical multipole moments are given in (35). The reader can see that the active dipole M i A vanishes in general relativity if the center of mass is chosen at the origin of the local frame of body A. In scalar-tensor theory the active dipole of a single body differs from the conformal dipole even for a single, non-spherical and/or non-stationary body by terms entering the first line of equation (47).
Equation (47) is also useful for discussing the dipolar tidal effects in scalar-tensor theories which cause a peer interest for their measurement may be possible in coalescing binary systems with an advanced third generation of gravitational wave detectors [80,82]. The tidal active dipole moment of each body in the binary is caused by the external dipole moment P i of the scalar field induced by its companion. Assuming that tides are stationary and neglecting higher-order multipoles in (47) we get the tidal part of the active dipole of body A where N A is the moment of inertia of body A. Expression (48) tells us that in general relativity there is no dipole effects as the active dipole vanishes in the limit when β → 1. Formula (48) describes one of the numerous finite-size effects in scalar-tensor theory of gravity and agrees with the expectation that these effects are of the first post-Newtonian order of magnitude [14,83]. In case, of black holes and/or neutron stars the moment of inertia N A ≃ G 2 M 2 A /c 4 , where M A is relativistic mass of the object. It makes clear that for the compact relativistic objects the tidally-induced active dipole moment D i A ≃ G 3 c −6 that affects motion of the finite-size bodies at the third post-Newtonian approximation [14] unlike general relativity where the finite-size tidal effects appear only in the fifth (∼ c −10 ) post-Newtonian order [31,84]. We also notice that the coefficient in the right hand side of equation (48) can be interpreted as the scalar tidal deformability parameter similar to that having been introduced in paper by Bernard [80,Eq. 2] λ (s) The scalar deformability parameter defines a dimensionless scalar tidal Love number (sTLN) k which extends the concept of the TLN in general relativity where tides are caused solely by the quadrupole and higher harmonics of tidal field [10,11,85]. Of course, the sTLN are absent in general relativity.

A. Configuration Space of Dynamic Variables
We consider post-Newtonian dynamics of N extended, arbitrary-structured bodies with a continuous distribution of matter inside each body. The distribution of matter is fully characterized by the energy-momentum tensor T αβ depending on the internal dynamic variables like the density of matter, its velocity, internal energy and stress/pressure. The internal variables of the matter's continuum obey a set of equations of motion which govern the dynamic evolution of the variables (see Appendix A for the case of the Newtonian approximation). The N -body system is considered in the present paper as being completely isolated from other astrophysical systems in the universe and we assume that the background spacetime manifold is asymptotically-flat with the Minkowski metric η αβ at infinity 2 . We postulate that each body of the N -body system is electrically neutral and does not exchange mass with the other bodies of the system. Thus, the only force of interaction between the bodies is the force of gravity carried out by tensor and scalar fields.
Self-consistent description of celestial dynamics of an N -body system calls for N+1 coordinate charts which are used for solving the internal and external problems of motion. A single global coordinate chart x α = (t, x i ) covers the entire spacetime manifold. It is introduced in order to characterize the relative motion of the centers of mass of the bodies with respect to each other. The other N local coordinate charts are introduced in a neighborhood of each body A∈ {1, 2, ..., N}. The local coordinates w α A = (u A , w i A ) adapted to the body A and are carried out along with the body. The local coordinates are used for the description of dynamics of the internal motions of matter inside the body and evolution of the multipole moments of the body.
The independent dynamic variables of the external problem are coordinates

and accelerations a i
A =ẍ i A of the center of mass of body A in the global coordinates. We have included the acceleration as it explicitly appears in the equations of motion of the N -body problem possessing all multipoles [15,16,18]. The accelerationdependent Lagrangians are fairly common in relativistic mechanics. For example, in case of description of the dynamic evolution of spin [52,53] and in the Lagrangians of point-like particles in the post-Newtonian approximations of higherorder [29,30,44,45]. The configuration space of the external variables x i A has dimension 3N . Choosing a set of the internal dynamic variables in the local coordinates of body A depends on the approach to the solution of the internal problems. In case of an extended bodies the most natural choice would be a set of hydrodynamics variables like density of matter, pressure, internal energy, etc. interconnected through the equations of state. This approach is clearly unsuitable in case of black holes which internal dynamic variables consists of a set of the mass and spin multipole moments. The multipole moments of the extended bodies can be also used as a set of basic dynamic variables of the internal problem of an extended body A because the multipole moments are defined as volume integrals from the distribution of hydrodynamic variables which can be considered as a mapping from one set of the dynamic variables to another. In any case the configuration space of the internal variables is infinite.
In this paper we use the multipole moments of the body (both mass and spin type) as the primary internal dynamic variables. They are not supposed to change under a virtual variation of the external variables and should be defined in such a way that make them split distinctively from the external degrees of freedom. The mass and spin multipoles M L A and S L A which enter equations of motion of the center of mass (see section III B below), have been defined by equations (32) and (40) respectively. The spin multipoles enters only the post-Newtonian terms of the equations of motion and its Newtonian-like definition definition (40) is fully satisfactory and sufficient to consider it as a primary dynamic variable in the post-Newtonian Lagrangian of translational motion. This is not true for the mass multipole M L A of body A enters the Newtonian force while its definition (32) includes explicitly post-Newtonian terms depending on the coordinates and accelerations of the centers of mass of external bodies. These external variables should be excluded so that only the first term I L A in the right hand side of (32) can be considered as a primary internal dynamic variable. There are also secondary internal dynamic variables of body A which include other characteristics of matter distribution inside the body like the non-canonical multipole moments N L A and R L A defined in (35) as well as various integrals defining the amount of the total internal energy of matter and/or gravitational field of the body, etc., discussed in Appendix A. It is worth pointing out over here that the mass and spin multipoles of body A are residing on the world line of the center of mass of the body and can be treated as purely spatial Cartesian tensors in tangent spacetime being orthogonal to 4-velocity of the center of mass of the body. Time evolution of the multipoles is governed by the local equations of motion of the internal dynamic variables, e.g., the equation for time derivative of spin (70) or that for mass (A9), etc.
Dynamical behavior of an N -body system can be analyzed with the help of a Lagrangian which serves as an indispensable tool for deriving conserved quantities that are obtained by employing the property of relativistic invariance of the Lagrangian with respect to local diffeomorphisms in the configuration space of dynamic variables in accordance with Nöther's theorem [19]. The Lagrangian of the N -body problem can be formulated directly starting from the full action of gravity theory [24,29,35,36,87,88]. It is known as the Fokker Lagrangian since the method of its derivation has been first formulated by Fokker [23]. Especially successful for the derivation of the Fokker Lagrangian became the effective field theory approach applied to point-like massive particles [89,90] and to covariant description of dynamics of extended bodies with spin and quadrupole moments [58,59].
The Lagrangian of an N -body system can be also derived independently by solving the inverse problem of the Lagrangian mechanics which comes down to the integration of a system of partial differential equations being adjoint to the equations of motion of the N -body problem [41]. Of course, the Fokker Lagrangian and the Lagrangian of the inverse problem are to be identical up to the terms which do not affect the dynamics of the bodies. These are the total time derivatives and/or functions of the "double zero" [53]. Correspondence between the various presentations of one and the same Lagrangian can be established by the method of redefinition of position variables [91] which is based on application of contact transformations of dynamic variables in the jet prolongation of extended configuration space. In the present paper we derive the Lagrangian by solving the inverse problem for the post-Newtonian equations of motion of the N -body system taken from [18] and given explicitly in next section.

B. Post-Newtonian Equations of Motion
Translational equations of motion for each body A∈ {1, 2, ..., N} from an N -body system in the framework of a scalar-tensor theory of gravity have the following form [18,Eq. 290] where m A is the inertial (conformal) mass of body A (44), a i A = d 2 x i A /dt 2 =ẍ i A is the coordinate acceleration of the center of mass of body A, and the net gravitational force in the right hand side of (50) is split in three parts -the Newtonian gravitational force F i N , the post-Newtonian gravitational force F i pN , and the post-Newtonian term I i A depends on the post-Newtonian shift of the origin of the local coordinates of body A from its center of mass as explained above in section II D. Vector-function is not subject to dynamical equations as it characterizes the residual coordinate freedom in choosing position of the center of mass of the body in the local coordinates. The original derivation of the post-Newtonian equations of motion (50) of the N -body problem has been given in [18] under condition that the center of mass of body A is located at the origin of the local frame at any instant of time which yields I i A ≡ 0. Below, we demonstrate that solution of the inverse problem of the Lagrangian mechanics does not allow to keep on this condition and requires to displace the body's center of mass from the origin of the local frame yielding I i A = 0. The inertial mass m A is not constant but depends on time. This is because body A is affected by the tidal forces from other N-1 bodies which interact with the multipole moments of body A causing post-Newtonian variation of the mass as the bodies move [14,15,64,72]. Depending on the equation of state and viscosity the tidal deformation of the body causes dissipation of the internal energy of the body that produces a secular decrease of the body's mass. If the dissipation of the internal energy can be neglected the mass of the body is subject merely to periodic orbital variations. In any case, the time derivative of the body mass,ṁ A = 0. This effect should be computed analytically and taken into account in the subsequent derivation of the Lagrangian.
The force components are represented in the form of a linear superposition of the STF partial derivatives from functions R AB and R −1 AB representing respectively the coordinate distance between bodies A and B and its inverse. More specifically, are the components of the active STF multipoles of external body B referred to the local frame of body B, is the coordinate distance between the centers of mass of body A and B, and repeated indices mean the Einstein summation rule. The mass and spin multipole moments of body A that appear in the force of gravity (51), (52), are given by equations (32), (40) and similarly for body B. The multipoles of body A are spatial tensors defined on the world line W A of the center of mass of the body. They are functions of the coordinate time u A of the local coordinates adapted to body A. We use in (51), (52) the following notations: The multipoles of any other body B = A are spatial tensors defined on the world line of the center of mass of body B and depending on the coordinate time u B of the local coordinates adapted to the body B. We have denoted (51), (52). Particular functional dependence of the multipoles on time is determined by solving equations of the time evolution of the multipoles of the corresponding body.
The spatial partial derivatives from the coordinate distance are computed in accordance with definitions which are taken at point x A ≡ (x i A ) -the center of mass of body A, and we understand that R B ≡ |x − x B |. In other words, the partial derivatives from R −1 AB and R AB are to be always understood as the derivatives taken with respect to the first argument of R AB or R −1 AB . This convention means, for example, that the l-order partial derivative and similar equation is valid for R BA . It is also worth noticing that At the same time, △R B = 0 and, hence, ∂ <L> R AB = ∂ L R AB . Exact correspondence between ∂ <L> R AB and ∂ L R AB is given in [18,Eq. (268)] The coefficients of the post-Newtonian force (52) are given below exactly as they appear in [18,Eqs. (302)-(313)] with some typos corrected in Eqs. (304), (306), (309) and (310). We have The PN coefficients (57)-(68) depend on the active mass and spin multipoles of the bodies of the N -body system and their time derivatives taken in the local frame adapted to the corresponding body. They also depend on velocities of the centers of mass of the bodies and their accelerations defined with respect to the global coordinates. Coefficient (68) describes dependence of the force on the matrix F ij A of relativistic precession which obeys an ordinary differential equation [18,Eq. (151)] The matrix of the relativistic precession consists of the contribution of the Lense-Thirring, de Sitter and Thomas precessions which correspond to the first, second and third terms in the right hand side of (69) respectively. The mass and spin multipoles of each body appearing in the equations of motion (50)-(52) describe the internal degrees of freedom of motion of matter of the corresponding body. They are treated in the equations of motion as purely spatial tensors residing in the tangent spacetime manifold being orthogonal to the world line of the origin of the local coordinates adapted to the body. The time dependence of the multipoles is not constrained by the condition of rigid rotation. We assume the most general case that the multipoles evolves over time in accordance with their own equations of motion derived from the microscopic equations of motion of matter inside the corresponding body. In the present paper we need three equations of motion of the multipoles in the Newtonian approximation. One of them is the equation of time evolution of the mass quadrupole moment which is derived in Appendix A. Second equation is the equation of evolution of the spin dipole of body A. It reads [18,72,92] where ε ipq is fully-antisymmetric symbol of Levi-Civita. The third equation is the equation of time evolution of the inertial mass (136).
Post-Newtonian equations of motion (50)- (52) of the N -body problem with all multipoles taken into account represent a system of ordinary differential equations of the forth order. Indeed, the reader can observe that the last line in equation (52) explicitly contains the first and second time derivatives of the acceleration of the center of mass of body A coupled with the spin-dipole and mass quadrupole of body A. This fact has been noticed by Damour et al. [72]. Equation of motion depending on the higher-order time derivatives of the center of mass of moving bodies may lead to the appearance of spurious, self-accelerating solutions [93] which make no physical sense. A standard way to eliminate the self-accelerated modes is to replace the terms with higher-order time derivatives with the equations of motion of lower approximations which reduces the order of the time derivatives [81,94]. This procedure is applied for calculating the gravitational force in inspiralling compact binaries in the second, third and forth PN approximations [33,38,61,95,96] where the terms with the higher order time derivatives appear due to the finite speed of gravity and retardation effects [31,44,45,94]. It should be noticed, however, that the nature of the terms depending on the higher-order time derivatives in the post-Newtonian force (52) is not associated with the retardation of gravity but rather with the accelerated motion of the center of mass of body A and non-linear (quadratic) dependence of the transformations between global and local coordinates adapted to the body, on spatial coordinates.
Finally, we notice that the post-Newtonian equations of translational motion derived by Racine and Flanagan [15], Racine et al. [16] in general relativity are in complete agreement 3 with our equations (50)- (70) in the limiting case of β = γ = 1. We draw attention of the reader to the fact that it is important to know the dependence of body's multipoles on the external variables for derivation of the Lagrangian of the N -body problem. Therefore, we have computed the multipole moments as functions of the internal and external variables and shown their explicit form in section II C of the present paper. This part of the computation is omitted in the paper by Racine and Flanagan [15].

C. Lagrangian and Variational Derivative
The central task of the present paper is to solve the inverse problem of the Lagrangian mechanics [41][42][43] that is to find out the Lagrangian L corresponding to the translational equations of motion (50)-(52) with the PN coefficients given in (57)- (68). The post-Newtonian Lagrangian L of translational motion of an N -body system with arbitrary structured bodies possessing the mass and spin multipoles depends on the external dynamic variables which are coordinates and velocities of the centers of mass of the bodies. It can also depend on accelerations of the bodies as was shown in case of pole-dipole particles with spin by Damour [52] and Barker and O'Connell [53,54] (see also [97,98]). The Lagrangian of point-like massive particles depends on accelerations of the bodies non-linearly starting from the second PN approximations [30,44,45]. We have observed that the post-Newtonian force (52) depends on the second order time derivatives of the accelerations. Hence, it suggests that the post-Newtonian Lagrangian of N bodies having quadrupole and higher-order multipole moments should depend on accelerations as well.
We consider the principle of the least action for the external problem of motion of N bodies in the configuration space of 3N external dynamic variables which are the global coordinates of the center of mass of body A∈ {1, 2, ..., N}. The internal dynamic variable is any integral quantity taken over volume of body A from the interior distribution of mass density, velocity, etc. Such integrals do not depend explicitly on the external variables. The internal dynamic variables are identified with the Blanchet-Damour mass multipoles (33), the non-canonical multipoles (35), the spin multipoles (40) and a number of other integrals characterizing the internal structure of bodies like those entering the virial theorems -see appendix A.
The primary internal dynamic variables of the N -body system are the mass and spin multipoles -I L A and S L A . We notice that the multipole I L A makes up the basic (Newtonian) part of the mass multipole moment M L A of body A given in (32) which tells us that M L A is a composite function of both internal and external variables. Nonetheless, M L A is convenient for more economic formulation of the equations of motion (50)- (52). Definition (40) of spin multipoles S L A is sufficient in description of translational equations of motion in the post-Newtonian approximation and does not depend on the external variables explicitly.
The internal dynamic variables are related to the external dynamic variables through the differential equations of motion of the internal problem. These equations are required for finding the Lagrangian of translational motion of the N -body system. In this paper we need the Newtonian equations of motion for the mass-monopole, spin-dipole and mass-quadrupole of body A. The equations of motion for mass-monopole and mass-quadrupole are derived in Appendix A. The equation for spin-dipole is given in (70). These equations are employed in section VI D for derivation of the Lagrangian of translational motion.
We denote the Lagrangian of the external problem as L ≡ L ext . We postulate that in the global coordinates the Lagrangian is given by a scalar function are coordinates, velocity and acceleration of the center of mass of the bodies, Q A = {M L A , S L A } denotes a set of the internal dynamic variables which are the mass M L A and spin S L A multipoles, and Λ A is a set of the internal integral variables like those entering the virial theorems explained in appendix A. We denote the components of the multipoles in the global coordinates as M L A and S L A to distinguish from their components M L A and S L A computed in the local frame of body A. Both external and internal dynamic variables are functions of global coordinate time t but the internal variables enter the Lagrangian of the external problem as time-dependent parameters. Equations of motion for the internal variables are derived not from L ≡ L ext but the Lagrangian of the internal problem L int which we don't discuss in the present paper.
Translational motion of bodies of an N -body system is determined by the action Variation of the action leads to the variation of the Lagrangian where the variation of coordinates, , is taken on a hypersurface of constant coordinate time t, and denotes δ-variational derivative of the Lagrangian in the extended configuration space of the external dynamic variables, [19,30] taken under condition that all partial derivatives are computed with keeping the internal variables The least action principle states that the corresponding variation of the action (71) vanishes on-shell, δS = 0, that yields the Euler-Lagrange equations of motion of the centers of mass of bodies in terms of the variational derivatives, The δ-variational derivative represents a covariant variational derivative expressed in a single global coordinate chart. The internal dynamic variables Q A ,Q A , Λ A in this representation are expressed in the same coordinates as the external dynamic variables. It is the δ-variational derivative which is used in covariant calculations of equations of motion in the effective field theories [59,89,99]. Nonetheless, physically measurable are the components of multipoles expressed in the local frame adapted to each body and these, locally-frame multipoles appear in translational equations of motion (50)-(52) of extended bodies in the global coordinates. Therefore, it is conceivable to look for the Lagrangian of the external problem considered as a scalar function L ≡ L x A , v A , a A ; q A ,q A , λ A with the internal dynamic variables q A = {M L A , S L A } and λ A computed in the local frame of body A. To make the Lagrangians L and L consistent, the variables {q A , λ A } must be related to {Q A , Λ} by the frame and/or gauge transformations of the internal variables, let say -see sections IV A and VII B for specific details. It is clear that the concept of the variational derivative must be modified to accommodate the new set of variables.
Let us introduce a new symbol d for the variation of the action S and corresponding dynamic variables which respects the variables q A , λ A defined in the local frame of body A, The d-variation of the external variables in (75) is now performed with keeping constant values of the internal variables where, by the new convention, the partial derivatives from the Lagrangian L are computed under the condition that the internal variables {q A ,q A , λ A } in the Lagrangian L are kept constant. The least action principle, dS = 0, applied to the Lagrangian L yields translational equations of motion which should be equivalent to the equations of motion (74) after accounting for the law of transformation of the internal variables.

D. The inverse problem of the Lagrangian mechanics
The primary goal of the present paper is to find a Lagrangian such that its variational derivative would be equal to the equations of motion (50). Two possible approaches are conceivable based on two types of the Lagrangians and the internal dynamic variables explained in the previous section. We have or The problem of derivation of the Lagrangian by solving either (78) for L of (79) for L is known as the inverse problem of the Lagrangian mechanics. It has been thoroughly analyzed by a number of previous researchers and the most notable results are summarized in [41][42][43]. A general method developed for solving the inverse problem of the Lagrangian mechanics relies upon mapping of equations of motion to an adjoint system of partial differential equations for a scalar function associated with the Lagrangian if solution of the adjoint system exists. Finding solution of the adjoint system is not an easy problem of the theory of differential equations in partial derivatives and will not be implemented in the present paper due to the complexity of post-Newtonian equations of motion with all multipoles. In fact, there is no need to solve equations (78) or (79) separately as they suggest that the two variational derivatives must be equal, which means that the Lagrangians L and L are to be interrelated by some transformation which can be found by solving (80). Assuming that the transformation between the internal variables is of the post-Newtonian order and the dependence of the Lagrangians onq A and λ A variables is only at the post-Newtonian terms, equation (80) is reduced to which suggests that L can be obtained from L by doing some transformation of the internal variables which is not obvious in the most general case of the bodies depending on all multipoles. In what follows, we shall focus on solving equation (78) and finding the Lagrangian L. Its counterpart, the Lagrangian L, which is a solution of (79) will be found in section VII B by applying a series of transformations of the internal variables and the virial theorems from appendix A which is equivalent to solving equation (81). The reason for this approach to the problem is simply due to the fact that we could solve the inverse problem by integrating (78) and to find the Lagrangian L while solution of the inverse problem by integrating (79) turned out to be less obvious to tackle. The procedure that shall be used in the present paper to find out the Lagrangian L is based on the observation that the right hand side of (78) is represented in the form of the differential operator from coordinates of the center of mass of the bodies and depend on velocities and accelerations algebraically. The Lagrangian must have a similar algebraic-differential structure and can be found by a series of transformations bringing the terms with the second time derivatives of the multipole moments and accelerations of the bodies in the post-Newtonian force (52), to the form of the variational derivative (73) accompanied with additional terms preserving the identity of the transformation. Some of the additional terms are, then, canceled out with similar terms in the post-Newtonian force. Repeating this procedure several times reveals that all terms in the right hand side of (78) are unambiguously reduced to the form of the variational derivative (73) taken from a linear combination of scalar functions which, thus, represents the major part of the Lagrangian that we are looking for. Nonetheless, there remain a few troublesome terms in the right hand side of (78) which cannot be directly reduced to the variational derivative by applying the above procedure. These terms are associated with the implicit coupling of the internal and external degrees of freedom, and their reduction to the variational derivative requires utilization of the equations of motion of the internal problem for spin-dipole and mass-quadrupole moments. These equations map the troublesome terms in the translational equations of motion to a number of the integrals and their time derivatives which are the intrinsic variables. Rather remarkably, these intrinsic variables couple to the velocity of body A in such a way that they can be reduced to the variational derivative. It completes derivation of the Lagrangian.
It should be clarified that solution of the inverse problem of the Lagrangian mechanics is not unique. Indeed, the δvariational derivative (73) taken from a total time derivative of an arbitrary function F (t, x A , v A ; Q A ,Q A , Λ A ) depending on coordinates and velocities of the bodies vanishes identically. Therefore, the Lagrangian L = L(t, are equivalent in the sense that they yield identical equations of motion. However, it may be necessary to take into account that adding a total time derivative to the Lagrangian must be accompanied with a corresponding change of the initial and/or boundary conditions for the N -body system. Naturally, the total time derivative added to the Lagrangian also contributes to the laws of conservation of the N -body system [19]. We shall discuss these subtleties somewhere else. Besides the Lagrangian transformations mentioned above, there are three auxiliary transformations which we have to perform before looking for the solution of the inverse problem. First, equations of motion (50)- (52) have been derived and are given in terms of the STF partial derivatives taken from the coordinate distance R AB = |x A − x B | between each couple of bodies in the N -body system. However, the variational derivative (73) has a free spatial index which is clearly split from the STF projection. It points out that we are to convert all STF partial derivatives to a symmetric form to peel off the index of the variational derivative. Second, the procedure of the matched asymptotic expansions used for derivation of the translational equations of motion makes the multipole moment of each body B = A of the N -body system a function of the local coordinate time u B which is simultaneous with the time argument u A of multipoles of body A. However, the least action principle applied to the action (71) The coordinate basis vectors E α are neither mutually orthogonal in the tangent spacetime nor they are unit vectors.
Transformation between the basis vectors on a tangent manifold affiliated with body A∈ {1, 2, ..., N} reads where the components of the matrix of transformation Λ α β are given by the law of transformation of the coordinate bases of the local and global coordinates taken at the origin of the local coordinates [18,.
and F ij A is the skew-symmetric matrix of relativistic precession of the spatial axes of the L-frame with respect to the G-frame defined by solution of equation (69).
We have used calligraphic letters to denote four-dimensional tensor components of the mass, M α1..  (83) demonstrates that the variation of the basis vectors of the L-frame will not vanish because the components of the matrix Λ β α depend on the coordinates and velocity of the center of mass of body A. Therefore, the L-frame components of multipoles M L A and S L A are functions of the external variables through the components of the matrix of transformation Λ β α . It means that the variational derivative from the L-frame components of the multipoles M α1...α l A and S α1...α l A do not vanish and must be computed explicitly. For this purpose it is convenient to transform the L-frame components of multipoles to the G-frame to make the dependence of the L-frame components on the external variables explicit.
To this end, we recollect that, by definition, the multipoles of body A are tensors which have only spatial components in the L-frame of the body. It means that the multipoles of body A are orthogonal to four-velocity of the world line of its center of mass. We denote the four-velocity as U α in the L-frame, and as U α in the G-frame. The condition of orthogonality of the multipole moments of body A to the four-velocity of its center of mass reads, In the L-frame of body A we have, U 0 = 1 and U i = 0, so that the condition (88) agrees with the condition that in the L-frame multipoles have only the spatial components M L A and S L A . This property does not hold in the G-frame that is the time-like components of M α1...α l−1 β A and S α1...α l−1 β A do not vanish. However, all time components of the G-frame multipoles like M 0i1..i l , etc., can be expressed solely in terms of their spatial components by making use of (88).
The components of the multipoles in the L-frame are defined by volume integrals (28), (40), where the local coordinates w α = (u, w i ) should be treated as the affine coordinates of the L-frame. The affine coordinates of the G-frame in the tangent manifold are denoted as X α = (T, X i ) where T coincide with the global coordinate time, Transformation between the affine coordinates of the tangent spacetime is equivalent to the transformation of 1-forms on the tangent manifold which is inverse to the transformation (83) of the basis vectors of the two frames, The frame transformation of the mass multipoles of body A should be performed with the post-Newtonian accuracy as the mass multipoles enter the Newtonian force in equations of motion (50)- (52). At the same time, the spin multipoles enter only the post-Newtonian force (52) and their frame transformation is sufficient with the Newtonian accuracy. Thus, for the spin multipoles we simply have The L-frame tensor components, M L A , of the mass multipoles are transformed to the G-frame components, M L A , with the help of the coordinate transformation (89). It yields This transformation includes the time component, M α1...α l−1 0 A , of the mass multipole in the G-frame but it can be found from the condition (88) of orthogonality of the multipoles to the four-velocity of body A, Four-velocity of body A in the G-frame is given in the Newtonian approximation by U 0 = −1 and that allows us to compute all time components of the multipoles in the G-frame from the corresponding spatial components of the multipole. Hence, the frame transformation (91) of the mass multipole components of body A can be written down solely in terms of the spatial components, and similarly for body B, where the G-frame multipoles M L A , M N B are explained below in (97), (101), and are the G-frame components of the non-canonical multipoles defined earlier in the L-frames of the corresponding bodies as N L A and N N B . Time T A in the G-frame of body A corresponds to the coordinate time u A in the L-frame of this body in accordance with (89).
The active mass multipoles M L A of body A take in the G-frame on the following form where the Blanchet-Damour mass multipoles and are the non-canonical multipoles computed in the G-frame, and the densityσ(T A , X) is exactly the same as defined above in (34) after replacing coordinates u → T A , w → X. Integration in (98), (99) is performed over the hyperplane of constant time T A of body A with making use of the post-Newtonian law of transformation the invariant mass density ρ * = ρ √ −gu 0 (see, e.g., [68,Eq.4.81]) Similar expression defines the active mass multipoles of body B where and Integration in (102), (103) is performed over the hyperplane of a constant time T B of body B. Time T B in the G-frame of body B corresponds to the coordinate time u B in the L-frame of body B in accordance with (89). Correspondence between the times T A and T B will be explained in the next section IV B. Post-Newtonian frame transformations (94) and (95) do not affect the multipoles in the post-Newtonian force (52) but bring about additional post-Newtonian corrections to the Newtonian force (51) which now takes on the following form, where we have explicitly shown the time dependence of the mass multipoles in the first, purely Newtonian term in the right hand side of (104). Times u A and u B are simultaneous in the L-frame of body B but they are not simultaneous in the G-frame of this body [18]. Since times T A and T B are related to the local times u A and u B by linear transformation (89) the times T A and T B are not simultaneous either. We consider the time T A as a reference time for all bodies of the N -body system and take T A = t. The time delay between T A and T B = T A + ∆ AB makes the multipoles of body B dependent on the external variables like the coordinates and velocities of the centers of mass of the bodies from the N -body system through the time delay which depends on coordinates and velocities of the center of mass of body B, ∆ AB = ∆(v B , R AB ). This dependence of the time arguments of multipoles on the external variables should be made explicit in order to compute the variational derivative from the multipoles. This is done in the next section.

B. Time transformation of the Newtonian force
Coordinates of the center of mass of each body of the N -body system are computed in the equations of motion (50)-(52) at the same instant of global coordinate time t, that is On the other hand, the mass multipole moments in the Newtonian force (104) are taken at different times T A and T B = T A . The time difference is a consequence of the application of the matching procedure to establish a correspondence between the local and global coordinates which yields [18, Section IX.B] where t is the value of the global coordinate time. The mass multipoles of body A are computed at the time t, that is M L A (T A ) = M L A (t). The mass multipoles M N B (T B ) are expanded in the Taylor series around the instant of time t by making use of (105), From now on, we adopt the following notations for the G-frame components of the mass and spin multipoles taken on the hypersurface of constant time t, Our next step is to substitute (106) to the first term in the right hand side of the Newtonian force (104) and discard all residual terms of the order of O(4) which are redundant. It yields, which can be recast to STF form with the help of Eq. (B15) from paper [18], The term with the STF derivative in (109) should be transformed to a symmetric derivative by applying the peeling off decomposition of indices [18, Eq. (A1)], It allows us to bring (108) to the following form, which should be substituted to the very first (Newtonian) term in the right hand side of (104).

C. Transformation of the STF partial derivatives
Computation of the Lagrangian by solving the inverse problem equation (78), is significantly simplified if all STF partial derivatives in the post-Newtonian gravity force F i pN are replaced with the symmetric partial derivatives entering definition of the STF projection (1). To this end, we notice that the STF derivatives from R −1 AB are already equivalent to the symmetric derivatives from R −1 AB , in the sense that This property is a direct consequence of the Laplace equation AB . However, the post-Newtonian force F i pN contains STF derivatives from R AB which are not equivalent to the symmetric derivatives, that is ∂ <L> R AB = ∂ L R AB . Terms with the STF derivatives from R AB appear in the third line of (52) and we are to bring all of them back to the symmetric form.
For converting the term depending on the second time derivative from multipole M N B we employ [18, Eq. (271)] which yields, where we have also utilized the frame transformation of the multipoles which is reduced in the post-Newtonian terms to a trivial replacement of all multipoles defined in the L-frame of a body to its G-frame.
To convert the term depending on the acceleration a p B of the center of mass of body B in the third line of (52) we resort to [18,Eq. (273)]. It results in where again we have replaced in the right side of (114) the mass multipoles defined in the L-frame to those defined in the G-frame. For transforming the term being quadratic with respect to velocities of the center of mass of the bodies in the third line of (52) we employ [18,Eq. (272)] along with (110) of the present paper. It brings forth the following identity, This completes the requested transformation of the STF partial derivatives to a standard symmetric derivatives.

D. Reduced form of the equations of motion
Now, we summarize the results of the preceding transformations from sections IV A-IV B by substituting them to equations (50)-(52) and contracting similar terms. It yields the equations of motion (50)-(52) in a reduced form that is suitable for construction of the Lagrangian by solving equation (78) of the inverse problem. It is convenient to split various terms in the reduced equations of motion in several distinctive groups depending separately on the mass and spin multipoles. as they can be treated independently. More explicitly, the reduced form of the equations of motion (50)-(52) reads where m A is the inertial mass of body A defined above in (44), F i n is the Newtonian force, F i m and F i s are two components of the post-Newtonian force depending respectively on the mass and spin multipole components, We emphasize that all spatial derivatives in (116)-(119) are now the symmetric partial derivatives, and the differentiation with respect to spatial coordinates of the bodies is understood in the sense of equations (54), (55). The PN coefficients of the differential operator in (118), (119) depend on time and are divided in two clusters depending on the mass and spin multipoles. These coefficients are explained below in full detail.

Mass multipole coefficients
The mass multipole coeffcients form three groups. The first group of mass multipole coefficients includes those which depend solely on the mass multipoles, their time derivatives and/or velocities of the bodies. They are: The second group of the mass multipole coefficients depend on the mass multipoles and accelerations of the center of mass of the bodies: The third group consists of just a single coefficient describing gravitational coupling of the mass multipoles with the gravitational potential of the bodies being external to body B in the N -body system.

Spin multipole coefficients
The spin multipole coefficients are split in two groups. The first group depends on the products of spin and mass multipoles and either their time derivatives or velocities of the bodies. They are: α pLN The second group includes only one coefficient describing mutual gravitational coupling of the spin multipoles of bodies from the N -body system: Equations of motion (116)-(119) replace the right hand side of (78) and will be used below for solving the inverse problem of the Lagrangian mechanics in the form We shall find the Lagrangian L by consecutively picking up the coefficients in the right hand side of (134) depending on the second time derivatives and bringing them to the form of the variational derivative (73) with simultaneous reduction and simplification of the remaining terms. After converting all those terms to the variational derivatives and mutual cancellation of all similar terms we shall be left with a few residual terms which are not directly reducible to the variational derivative. Nonetheless, the residual terms can be converted to the variational derivative after transforming them to the time derivatives from the internal dynamic variables by applying the equations of motion (70) for spin dipole S i A and the virial theorem (A24) for the mass symmetric moment ‫ג‬ ij A .

V. THE LAGRANGIAN FORM OF THE NEWTONIAN FORCE
We start solution of the inverse problem of translational motion of N bodies from transformation of the first two terms in the right hand side of (134) which constitute the Newtonian second law of classic mechanics. The inertial term, m A a i A , corresponds to the kinetic energy of translational motion of the bodies and the Newtonian force, F i n , complies with the potential energy of their Newtonian gravitational attraction.

A. Kinetic energy
The kinetic energy of body A depends on its inertial mass that was defined in (44). The inertial mass contains the intrinsic part m A which does not depend explicitly on the external variables of the configuration space of the N -body system. However, it also contains the external multipoles moments Q L and P which are the explicit functions of the coordinates of the centers of mass of the external bodies in accordance with definitions (21) and (25). Due to the fact that active dipole moments of body A have a post-Newtonian order of magnitude, M i A = O(2), its exact structure in the post-Newtonian terms is not important as the difference goes over to the second post-Newtonian order. To make equations look more homogeneous we replace the dipole term in the last (post-Newtonian) term of (44) Q i M i → P i M i . We can also neglect in this term the difference between the multipole moments defined in the local and global frames and equate, M L A = M L A . We also use in the post-Newtonian terms notation M A for gravitational mass. Under these conditions the inertial mass of body A (44) can be written down in the following form The inertial mass m A is not constant. The rate of change of the inertial mass has been computed in [14,72,92] and can be determined directly by taking a time derivative from both sides of (135) with accounting for the time derivative of the mass m A -see [18,Eq. 163] and (A9). It yieldṡ Let us define the Newtonian kinetic energy K N of translational motion of all bodies in the N -body system by a standard formula of classic mechanics The variational derivative (73) from the kinetic energy K N is computed with taking into account that the inertial mass m A of body A depends explicitly on the coordinates of the external bodies of the N -body system. We have where in the first term of the very last line the partial derivatives are taken with respect to coordinates of body B, that is Some terms in the right hand side of (138) can be grouped together to form a variational derivative which can be taken to the left hand side of (138) to combine with the variational derivative from K N . More specifically, we notice that Moreover, it is straightforward to show that in the first PN approximation the following identity is valid, This term represents a post-Newtonian correction to the kinetic energy (137) and it is reasonable to add it to the variational derivative (138) of the kinetic energy along with the term (139). Subtracting (139) from (138) and adding (140), we notice that some similar terms in the so-obtained expression cancel out and the inertial, acceleration-dependent term in the right hand side of equation (134) takes on the form of the variational derivative plus some other terms, where we have denoted We shall use (141) in section VI D to combine it with other terms in the right hand side of (134) to derive the Lagrangian of translational motion of the N -body problem.

B. Potential energy
The Newtonian potential energy of gravitational interactions in the N -body system is associated with the gravitational force F i n presented by (117). It depends on the active mass multipoles (97) (97), (101) in the Newtonian gravitational force F i n , we get a fairly long formula, (98), (102), and taken at the same instant of time t as M L A and M N B -see (107). We have also introduced new notation for the symmetric moments in the G-frame, It is worth noticing that in the case l = 0 the monopole symmetric moment coincides with the mass of body A, that is ‫ג‬ A = M A . We are also allowed to consider the dipole ‫ג‬ i A = M i A = 0 in any post-Newtonian term. Almost each term in the decomposition (143) of the Newtonian force can be separately converted to the variational derivative either exactly or up to a term being equal to the ordinary, second-order time derivative from a linear combination of some functions which will be specified below. Indeed, it is straightforward to establish that the very first term in the right hand side of (143) is equivalent to the variational derivative where we have taken into account that the variational derivatives from the internal variables vanish: The terms with the accelerations in the third and forth lines of the right hand side of (143) are converted to The terms with the second-order time derivatives in the right hand side of (146), (147) can be removed from the equations of motion by appropriately choosing the shift function I i A in (134). It corresponds to making translation of the coordinates of the center of mass of body A from the origin of its local coordinates, point w i = 0, to another point w i A defined by equation (45) under the condition of vanishing conformal dipole moment J i A = 0. More exactly, we choose which means that the local coordinates of the center of mass of body A are now It should be emphasized that the active mass and spin multipole moments M L A and S L A with l ≥ 1 are still computed with respect to the origin of the local coordinates of the corresponding body, that is with respect to the point w i = 0 which is defined by the condition J i A = 0. The shifted values of the mass multipoles,M L A are related to the mass multipoles defined with respect to the origin of the local frame, by transformationŝ where we have taken into account that the translation w i A has a post-Newtonian order of magnitude according to (148), (149). Because the equations of motion depend on the mass and spin multipoles, M L A and S L A , defined with respect to the origin of the local coordinates, we continue to use these multipoles in the rest of the present paper.
The terms in the very last line of (143) represent the most laborious part in the computation of the variational derivative of the gravitational potential energy of N bodies. These terms describe two-point and three-point non-linear interactions of the mass multipoles. We follow Fock [49, §77] who did similar computations but only for monopole terms (the indices l = k = 0). Extending Fock's calculation for all other indices l ≥ 0, k ≥ 0, we have where are defined in terms of fully-symmetric (with respect to the summation indices A,B,C) functions Ψ AB = Ψ (AB) , Ψ ABC = Ψ (ABC) . Function describes the two-point gravitational interaction, and function describes the three-point interaction. Notice that the non-canonical symmetric multipole moments ‫ג‬ LN A that appear in (153), (154), have been introduced in (144). They can be replaced with their STF projections I IJ by making use of formula (A22b) from paper by Blanchet and Damour [73]. For example, where min(l, n) denotes the minimum of two numbers l and n, the indices I ≡ L − P , J ≡ N − P , and we have taken into account that the symmetric derivative from R −1 AB is equivalent to its STF derivative 5 . Similar STF-decomposition equations can be written for other terms in (153) and (154). Accounting for the results obtained in this section, we can express the Newtonian force (143) in the form of the variational derivative supplemented with terms depending on the external gravitational potentialsŪ (x A ),Ū (x B ) for body A and B respectively, and a second time derivative of vector function (148) which cancel out the corresponding shift function in the right hand side of the inverse-problem equation (134), where is the Newtonian-like potential energy of gravitational field where the mass multipole M L A and M N B belong to the G-frame of the bodies and have been defined in (98), (102).
Equation (157) can be further simplified to where we have used the Newtonian equations of motion, M B a i A = F i n , to reduce the term depending on the gravitational potentialŪ (x A ) of the bodies being external to body A. This replacement elucidates that the term depending on U (x A ) in (159) is exactly equal to a similar term in the right hand side of (141) so that the two terms will mutually cancel out in the final reduction of similar terms in the right hand side of (134).

VI. THE LAGRANGIAN FORM OF THE POST-NEWTONIAN FORCE
In this section we consider transformation of the post-Newtonian terms in the right hand side of (134) which contain time derivatives of the first, second and higher order, to the variational derivative. The terms with the time derivatives consist of four groups. The first group includes terms being proportional to the products of the partial derivatives from the coordinate distance R AB between the bodies and the time derivatives of the mass multipoles or acceleration of the bodies. The second group consists of terms being proportional to the products of the partial derivatives from the inverse of the coordinate distance R −1 AB with time derivatives of the mass multipoles or accelerations of the bodies which appear in force F i m . The third group consists of terms being proportional to the products of the partial derivatives from the inverse of the coordinate distance R −1 AB and the first time derivatives of the spin multipoles which constitute the force F i s . The forth group consists of the residual terms depending on the higher-order time derivatives of the multipoles and the external variables like velocity and acceleration of body A and their time derivatives. Below we consider these groups of terms one by one.
A. Transformation of the higher-order derivative terms coupled to distance RAB Terms being proportional to the symmetric partial derivatives from R AB appear in the second line of equation (134). First, we transform to the variational derivative the term which is proportional to the second time derivative, Next term is proportional to the acceleration. It is transformed as follows, Making use of the above formulas we get that the combination of terms in the second line of (134) is completely reduced to a variational derivative without residual terms, where B. Transformation of the higher-order derivative terms in the force F i m The force F i m is given in (118) with the PN coefficients listed in section IV D 1. The higher-order derivative terms include the second time derivatives of the mass multipole moments and accelerations of the center of mass of the bodies of the N -body system.

Terms with the second time derivatives of the mass multipoles
The second time derivatives of mass multipole moments appear only in the mass multipole coefficient α iLN m shown in (120). The transformation equations for the four different terms with the second time derivatives, which appear in this coefficient, are as follows: Transformations (164)-(167) convert all terms with the second time derivatives of the mass multipoles in the force F i m to the variational derivative and a number of other terms which bring additional contributions to other mass multipole coefficients entering (118). After reduction of similar terms the post-Newtonian force (118) is modified and takes on the following form and the coefficients in the force F i m , which have been modified by the above transformations, are marked with hat. They arê The remaining coefficients in (168) (without the hat) are not affected by the transformations (164)-(167) and remain the same as shown in (124)- (129). Now, we are to single out the variational derivatives from acceleration-dependent terms in (168).

Terms with the accelerations
Equation for the force (168) includes four terms depending explicitly on accelerations of the center of mass of the bodies of the N -body system. These terms appear in the PN coefficients β iLN m and β ipLN m . We transform these terms to the form of the variational derivative and combine the rest of the terms, obtained after these transformation, with similar terms in the remaining PN coefficients of the force (168).
Two terms depend on the acceleration of body A. Their corresponding transformations are: There are also two terms which depend on the acceleration of body B. Their transformations read: After replacing the above transformations of the accelerations in (168), the force F i m takes on the following form: where L 4 has been derived previously in (169). A new contribution to the variational derivative reads and the PN coefficients which have been changed in F i m after the reduction of similar terms, are marked with tilde and read, The coefficients α pqLN m and γ LN m are not affected by the above transformations and remain the same as in (125) and (129) respectively.
C. Transformation of the spin multipole terms in the force F i s In this section we transform the spin multipole-dependent force F i s given in (119) to the Lagrangian form. To this end, we compute the following variational derivative, This formula, after having been compared with expression (119) for the spin multipole-dependent force F i s , elucidates that the force can be exactly reduced to the variational derivative, where D. Reduction of the residual terms. Application of the virial theorem.
The post-Newtonian Lagrangian is obtained by replacing (141), (159), (162), (178), and (184) in the corresponding terms of equation (134) of the inverse problem of the Lagrangian mechanics. This makes the right hand side of (134) consisting of two groups of terms. The first group is represented in the form of the variational derivative from a linear combination of scalars L 1 , L 2 ,...,L 6 . The second group includes a few residual terms which have not yet been represented in the form of the variational derivative. The goal of this subsection is to demonstrate that the residual terms are also reduced to the variational derivative after some additional transformations.
First of all, we notice after inspection of structure of the residual terms that some of them can be directly converted to the Lagrangian form. Indeed, the terms in the last line of (134) depending on the STF quadrupole moment M ij A of the body A and its higher-order time derivatives can be written as with Alternatively, the terms depending on the STF quadrupole moment M ij A of the body A and its higher-order time derivatives can be removed from the equations of motion by making a post-Newtonian translation of the position of the center of mass of body A with respect to the origin of the local coordinates. This is what we did in our paper [18] to circumvent the appearance of possible spurious solutions of the equations of motion with higher-order time derivatives [93]. However, since the acceleration-dependent terms (186) are fully reduced to the variational derivative from the scalar L 7 , we prefer in the present paper to include them to the total Lagrangian of the N -body problem instead of making a complementary shift in the position of the center of mass of body A in addition to that shown in (149).
After putting together different pieces of the equations of motion given by equations (141), (159), (162), (178), (184) and (186) and canceling similar terms, the equation of the inverse problem (134) acquires the following form Despite the fact that we have moved quite far in reducing equation (134) of the inverse problem to the integrable form, the equation (188) still contains in the right hand side some residual terms that are not looking like the variational derivative. Our goal is to find out a next chain of transformations which will bring these terms to the explicit form of the variational derivative, thus, making equation (188) fully-integrable.
To this end, we notice that the residual terms in the second line of (188) can be reassembled as follows, where the last line was obtained by making use of the rotational equations of motion (70) for the spin dipole of body A written in the reverse order, Such a replacement of some terms in the post-Newtonian approximation with the help of the Newtonian equations of motion for dynamic variables is fully legitimate as this is done on-shell in the configuration space of the dynamic variables before derivation of the Lagrangian [19]. Now, we algebraically split the second term in the last line of (189) in its symmetric and anti-symmetric parts, and make use of (190) once again to replace the anti-symmetric piece with the time derivative from the spin S i A of body A. The symmetric piece in (191) is transformed on-shell to a set of intrinsic dynamic variables by applying a tensor virial theorem (A24) which is explained in Appendix A in more detail. Adding the results of these on-shell transformations all together yields, This reveals that all residual terms in the right hand side of (188) are fully reduced to the Lagrangian form, It is instructive to recast L 8 by making use of tensor virial theorem (A24). It yields where we have introduced a shorthand notation for a virial tensor function of a single body A Notice that function K pq A vanishes on-shell due to the virial theorem (A24) but it must be kept in the Lagrangian as the variational derivative from K pq A is not nil. Results of the previous section demonstrate that the right hand side of equation (78) of the inverse problem of the Lagrangian mechanics can be put to the Lagrangian form which means that (78) is reduced to a fully integrable form, where we have denoted and L ab , L abc are given in (152)-(154). We can now integrate equation (197) and obtain the post-Newtonian Lagrangian L of the translational equations of motion of bodies with arbitrary mass and spin multipole moments expressed in the global G-frame. The Lagrangian is a linear combination of several components, where is an arbitrary function of coordinates and velocity of the bodies describing the gauge freedom in definition of the Lagrangian [19,102]. This function is omitted from subsequent discussion as it does not contribute to the equations of motion. Nonetheless, the reader should keep in mind that different choices of the gauge function F yield dissimilar forms of the Lagrangian which can make more difficult the comparison with the results of other authors . The explicit form of L a is where m A being the post-Newtonian (Tolman) mass in general relativity defined in (42), the multipole moments M L A and M L B are defined by (98), (102) in the G-frames on the world lines of the corresponding bodies. They are Blanchet-Damour multipoles [73] in scalar-tensor theory for all l ≥ 0 including the mass monopole.
There are several points to comment on.
1. The Lagrangian L is identical to the Fokker Lagrangian which would have been obtained by integrating the action of the scalar-tensor field theory. All external and internal dynamic variables of the Lagrangian L are unequivocally computed in a single global G-frame. The multipole moments entering the Lagrangian are the Blanchet-Damour multipole moments [73] extended to tensor-scalar theory [13,14].
2. The Lagrangian L depends on accelerations of the bodies. The linear acceleration terms associated with spin of body A and its quadrupole moment have been known [52,53,59]. We have extended these results to the multipoles of higher order.
3. The Lagrangian L also depends on terms being quadratic with respect to accelerations of the bodies which couple with the quadrupole moments of the bodies. The quadratic-acceleration term means that the Hessian of L computed with respect to the accelerations is non-degenerate. This brings about the Ostrogradsky instability [103,104] that may be problematic in some applications. It is tempting to replace the accelerations of bodies in the Lagrangian with the equations of motion from the lower-order approximation. This is a legitimate procedure in action-based gravitational theories but it leads to an implicit (and, hence, uncontrolled) change in the original coordinates used for derivation of equations of motion and the laws of conservation [30,45,105]. Correct mathematical procedure of elimination accelerations and higher-order time derivatives from the Lagrangians relies upon the method redefinition of position variables [91] which is a generalization of the method of the "double zero" proposed by Barker and O'Connell [53]. We notice that the quadratic acceleration terms can be also removed from the Lagrangian L by redefinition of the center of mass of body A -see [18,section IX.6].
4. The Lagrangian L explicitly depends on the non-canonical multipoles like ℵ L A which are multipole extension of the second moment of inertia ℵ A of body A depending on the finite radius of the body. This fact has been noticed by Nordtvedt [106]. The moment of inertia for a compact body ℵ A ≃ G 2 m 2 A /c 4 , where m A is the Tolman mass of the object. Since the moment of inertia ℵ A appears in the post-Newtonian Lagrangian explicitly, the finite-size effects of compact bodies enter dynamics at the 3PN (∼ c −6 ) order. This conclusion is valid, however, only in the scalar-tensor theory with the parameter β = 1. In general relativity β = 1 so that all terms with the moment of inertia ℵ A vanish in the equations of motion (see section VIII E for the details) that postpones the appearance of the finite-size effects in general relativity to the 5PN (∼ c −10 ) order. Previous results of some other researchers, for instance, [24,[107][108][109], suggested that 1PN equations of motion in general relativity do depend on the moment of inertia ℵ A of the bodies (not associated with spin). These results are based on underdeveloped mathematical technique and are unjustified as the authors of the papers [24,[107][108][109] did not take into account that physical moment of inertia ℵ A must be calculated in the local coordinates adapted to body A but not in the global coordinates [14,50,110]. 5. The Lagrangian L a depends on tensor function K pq A given in (196) which explicitly includes the integral characteristics of the bodies depending on their specific internal structure -the second symmetric moments ‫ג‬ ij A , the kinetic energy of the internal motion of matter T ij A , the potential energy of gravitational field U ij A , and the integral of stresses S ij A -see their definitions in Appendix A. It can make an impression that the Lagrangian L can be applied only to the bodies with the known internal structure. However, this is not true since, after taking the variational derivative, we can take function K pq A = 0 on-shell due to the virial theorem and, thus, it does not enter explicitly the results of all possible applications like computation of the equations of motion or conserved quantities. Furthermore, the explicit dependence on the integrals characterizing the internal structure of the bodies can be completely removed from the Lagrangian by switching from the internal dynamic variables (multipoles) expressed in the global G-frame to physical multipoles computed in the local L-frame along with applying the virial theorems (A24), (A25) for transformation of gravitational masses M A of the bodies. This procedure is worked out in next section.
B. Lagrangian L with the multipoles defined in the local frame

Gauge transformation of the internal dynamic variables
This section solves the inverse problem in the form of equation (79) with the Lagrangian L expressed in terms of the physical multipoles defined in the local L-frame of each body. First of all, we transform the monopole which is the gravitational mass M A given by equation (98) for index l = 0. It reads where, m A , is the Tolman mass, a constant coefficient η = 4β − γ − 3 is known as the Nordtvedt parameter [68], the quantities ℵ A and ℜ A are scalar non-canonical monopoles defined in (96) and (99), the integrals T A , U A , ℵ A , S A characterize the internal structure of the body A as explained in equations (A6), (A7), (A20) and (A21) respectively. Accounting for the relationship,‫א‬ A = 2Ṙ A , and applying the virial theorem (A25) we bring (201) to the following form where the scalar function K A is the trace of tensor function K pq A , Function K A vanishes on-shell but must be kept explicitly in the Lagrangian as its variational derivative does not disappear.
The terms in the square brackets of (202) are of the post-Newtonian order of magnitude. Hence, the dipole (l = 1) term in the sum could, in principle, be neglected as it is of the post-Newtonian order of magnitude itself. Expression (202) helps us to understand that in case of general relativity, when parameters γ = β = 1, η = 0, the inertial, m A , and gravitational, M A , masses are identical for spherically-symmetric bodies in complete agreement with the strong principle of equivalence (SEP). In scalar-tensor theory of gravity the two masses are not identical -mainly due to the self-gravity term ηU A in (202) -and the difference between them causes violation of the SEP that was first noted by Nordtvedt [111,112,113]. Testing this violation in the Earth-Moon system and binary pulsars is one of the most effective experimental probes for detecting the presence of a hypothetical long-range scalar fields in nature [113][114][115][116][117][118][119][120][121][122]. Our calculations reveal that there is a non-static term contributing to the SEP violation. Indeed, according to (202) the gravitational mass M A of body A disagrees with its inertial mass m A not only due to the static Nordtvedt effect (term with η) but also because the body's moment of inertia ℵ A (defined in (A20)) changes in the course of time (for example, radial pulsations of star) and/or due to a tidal gravitational interaction. Let us introduce new notation for the G-frame multipoleŝ which allows to decompose the Newtonian gravitational energy as follows is the Newtonian potential of all bodies being external to body A. Now, we consider a linear combination of three terms which appears in the Lagrangian L a along with function K A emerging after transfromation (202). We take the variational δ-derivative from W in the sense of equation (73). It yields where we have accounted for the fact that after taking the variational derivative functions K pq A = K A = 0 and K pq B = K B = 0 on-shell due to the virial theorem but their partial derivatives do not vanish. After developing the partial derivatives in (208) by using definitions (196), (203), we get However, the right hand side of (209) is equivalent (up to the relativistic precession term) to the transformation (104) of the Newtonian gravity force due to the conversion of the multipoles from the global to local frame. Making use of the inverse transformations (94), (95) of the G-frame multipolesM L A to L-frame multipoles where I L A are given by (33), we find out that all terms in the right hand side of (209) cancel out with the corresponding terms from the multipole transformation formulas and it reduces to a simple Newtonian-like form where the matrix F kp A of the relativistic precession is a solution of the precessional equation (69), and we replaced the G-frame multipoles M L A with their L-frame counterparts M L A . Now, we define function that is similar but not identical to function W in which the relativistic frame transformation of multipoles (94) from the global to local frame has been performed. The frame-transformed function W differs from W by terms depending on the internal structure of bodies which are not present explicitly in W. Applying the variational d-derivative to W in the sense of equation (76) (that is the moments M L A do not vary) we find out that This equation allows us to transform the Lagrangian L with the mass multipoles M L A defined in the global frame to the Lagrangian L with the mass multipoles M L A defined in the local frame. More specifically, we make use of the transformations of the Newtonian potential energy described above and use them to replace the corresponding terms in the Lagrangian L from previous section. It yields the Lagrangian L with all physical multipoles expressed in the local frame of the bodies.

Lagrangian
The Lagrangian of the N -body problem with all multipoles expressed in the local frame of each body is given by equation where where m A being the post-Newtonian (Tolman) mass of body A in general relativity defined in (42), the mass multipole moments M L A and the spin multipoles S L A are respectively defined by (32) and (40) in the L-frames of the body A, and the non-canonical multipoles N L A are given in (35). Notice that the multipoles M L A are the scalar-tensor variant of the Blanchet-Damour multipoles [73] only for l ≥ 1. The monopole M A (l = 0) is not the Blanchet-Damour mass but a generalization of the Tolman mass in scalar-tensor theory of gravity found by Nordtvedt [112].
The terms L ab and L abc in (214) are given by combining equation (152) along with the very last term in the right hand side of the transformation equation (205). We have where functions , and the symmetric (not STF) moments I L A have been introduced in (36). We have extracted in (216), (217) the terms being proportional to the factor (1 − β) since these terms originate from the presence of self-interaction of the scalar field and characterize its strength in scalar-tensor theories. It is these terms which lead to the appearance of dipole tidal effects in the scalar-tensor theories of gravity as has been recently pointed our by [80]. Susceptibility of each body A to the dipolar tidal field is characterized by the deformability parameter λ (s) A introduced in (49). We emphasize that the variational derivative that must be applied to the Lagrangian L in (214) is d-derivative introduced in (76). The two Lagrangians L in (199) and L in (214) are equivalent in the sense of the identity (80) that is they both yield the same equations of translational motion (50)- (68). It is interesting to notice that the Lagrangian L can be obtained from the Lagrangian L by formally replacing the G-frame multipolesM L A to the L-frame multipoles M L A and using the virial theorems to nullify functions K pq A and K A . Indeed, the Lagrangian L is a function of the Gframe multipoles, L = L M L A while the Lagrangian L is a function of the L-frame multipoles, L = L M L A . Previous section tells us that 6 where functions K pq A and K A vanish on-shell due to the virial theorem, which is the equation of motion for the internal variables. Hence, it is tempting to take K pq A = K A = 0 in (220), and render a formal replacementM L A → M L A that makes the impression that L = L. However, this procedure is not a right way to establish the correspondence between the two Lagrangians. In fact, it is illegitimate since the virial theorem is valid solely on-shell and the equality K pq A = K A = 0 is not allowed in the Lagrangian before taking a variational derivative from the Lagrangian [91]. Moreover, replacing the multipoles from one frame to another is not a matter of changing their labels but a corresponding mathematical transformation. The correct procedure of establishing the correspondence between the two Lagrangians must rely upon the legitimate transformations and application of the variational derivatives in the way explained in previous section VII B 1. The Lagrangians are equivalent in the sense of the equivalence of two variational derivatives as have been discussed in section III D. The role of the virial theorem in derivation of the Lagrangian of the N -body problem with all multipoles taken into account and its importance for making a correct gauge transformation of the multipolar Lagrangian from one frame to another has not been pointed out previously in works of other researchers. The Lagrangian L with the physical multipoles M L A defined in the local frame of each body is more preferable in astrophysical applications as it contains the least number of terms depending on the internal structure of the bodies. The most important for gravitational wave astronomy are binary systems consisting of two objects of comparable masses. We discuss the Lagrangian L of such two-body system in next section.

VIII. LAGRANGIAN OF A BINARY SYSTEM
Two-body system is probably the most interesting case for gravitational wave astronomy. The inspiralling compact binaries are main sources of gravitational wave signal for present-day and near-future generations of gravitational-wave observatories [5,123]. Current model of inspiralling compact binary treats the bodies of the system as pole-dipole particles with constant mass and precessing spin. This model is tractable for analytic computations up to the 4PN approximation and was extensively used to construct a bank of templates of waveforms [124,125]. At the latest stages of coalescence tidal and rotational deformations of the bodies in the binary become so significant that the pole-dipole model is getting no longer accurate and should be extended to include higher-order multipoles in the Newtonian and post-Newtonian approximations [10,11,58,85]. This section presents the Lagrangian L for description of the orbital evolution of two-body system with all mass and spin multipole moments of the bodies taken into account.

A. General Formula
In case of two-body system, A={1,2}, the post-Newtonian three-body interaction Lagrangian L abc vanishes, and the overall Lagrangian is simplified to a linear sum of two terms The first term in (222) reads where the multipoles M L 1 , M L 2 are given by (210), R ≡ R 12 = |x 1 − x 2 |, all partial derivatives from R are taken with respect to x 1 like ∂ i R = ∂R 12 /∂x i 1 , a dot between two spatial vectors denote their Euclidean dot product while a cross between two vectors does the Euclidean cross product. The second piece of the Lagrangian (222) is where we have singled out the piece of the Lagrangian which is generated exclusively by the self-interaction of the scalar field having external multipole moments P 1L and P 2L produced by the first and second body respectively. Definition of the external scalar multipoles has been given in (24) and reads more explicitly as Formally, the structure of the term with the coefficient 1 − β in (224) is identical with the second term in the right hand side of this equation but the physical meaning of the two terms is drastically different because the second term in (224) describes exclusively the effects of general relativity.
In what follows, we consider a series expansion of the Lagrangian (222) with respect to the degrees of multipole order more explicitly, and write this series as where L m contains only monopole terms, L d includes dipole terms, L q takes in quadrupole terms while the terms of higher multipole order are denoted with a lower dots symbol.

B. Monopole approximation
The monopole part L m of the Lagrangian corresponds to the case of massive bodies treated as structureless point-like particles. It consists of all terms in (222)-(224) constrained to the value of indices l = n = 0. It reads, where the gravitational mass M A (A=1,2) is defined in (210) for the value of index l = 0. Notice that the term with the parameter β comes from the Lagrangian L ab describing the nonlinear interaction of the external scalar field. Because M A is not equal to the Tolman mass m A the Newtonian dynamics is modified as pointed out by Nordtvedt [111,112]. We have neglected the difference between M A and m A in the post-Newtonian terms of (227) as it leads to the residual terms of the 2PN order. Notice that our treatment of the SEP accounts for the effects of non-stationarity of orbiting bodies 7 which can be important in the treatment of dynamical tides in coalescing binary systems beyond the adiabatic tidal approximation in scalar-tensor theory [80,82].
The Lagrangian (227) for monopole massive particles has been derived independently many times by various researchers who have also used different mathematical techniques. Among these derivations the most notable works are by Lorentz and Droste [126] and Fikhtengolts [127] who had derived the Lagrangian L m from EIH equations of motion (see also [81, §106]). PPN parameter γ has been introduced to the Lagrangian L m by Estabrook [128]. Later on, PPN parameter β was added to the Lagrangian of point-like bodies by Will [68] and by Barker and O'Connell [129]. We notice that the β − γ Lagrangian written by Will [68,Eq. (6.80)] treats the inertial and gravitational masses entering the Lagrangian as being equal while Barker and O'Connell [53] clearly distinguish the inertial and gravitational masses to take into account the Nordtvedt effect. where is the mass dipole part of the Lagrangian, and is its spin dipole part. Here, the unit vector n = (n i ), the dot between two spatial vectors denote their Euclidean dot product and the cross "×" between vectors is the Euclidean cross product. It is worth noticing that the part of L dm with the scalar dipole moments D i is taken from the first term in (224) which is proportional to 1 − β.
In general relativity the particle dipole approximation refers merely to spins of the bodies and the dipole Lagrangian of general relativity includes only L ds . The mass dipole moments of the bodies in general relativity depend merely on the displacement between the origin of body's local coordinates and its center of mass. Thus, the mass dipole is not physical and vanishes in general relativity if the origin of the local coordinates is chosen at the center of mass of the body. In scalar-tensor theory of gravity the definition of the center of mass of each body A is given in terms of the conformal dipole moment J i A , introduced above in (45), while the gravitational force in the equations of motion is governed by the active multipole moments of the bodies. The active dipole moment M i A of body A (A=1,2) is taken from (210) and (33) for l = 1, and reads where the general-relativistic mass densityσ(u, w) is defined in (43), and vector integrals T i A , U i A , S i A depend on the internal distribution of matter as explained in Appendix A 3. The dipole M i A does not vanish even if the conformal dipole moment J i A = 0. Indeed, the very first term in the right hand side of (231) can be computed from the condition J i A = 0 by taking into account definitions (45) and (148), where the external gravitoelectric moments Q L have been defined in (21). Replacing (232) in (231) and making use of vector virial theorem (A15) we get which demonstrates the presence of the Nordtvedt effect in the dipolar approximation of the scalar-tensor theory of gravity. Notice that M i A = 0 for spherically-symmetric bodies because all terms in the right hand side of (233) are integrals which angular part vanishes in case of the spherical symmetry except of N A . However, this term couples with the local non-geodesic acceleration Q i which is nil for spherically symmetric bodies. This consideration reveals that the active mass dipole M i A is sufficiently small in practical situations and can be neglected in most cases. Nevertheless, the scalar tidal dipole D i does not vanish and gives contribution to L dm which stems from the dipolar part of the term being proportional to 1−β in (224). According to definition (48), the scalar tidal dipole is proportional to the moment of inertia N A coupled with the external scalar dipole P i and reads where λ It brings about the tidal dipolar effects which measurement can be used in order to put additional constraints on the parameters of the scalar-tensor theory by gravitational-wave detectors [80]. Equations of motion with the scalar tidal dipole effects included, are obtained from the Lagrangian L dm after taking the variational derivative. It yields The relative acceleration a i = a i 1 − a i 2 of the two bodies is where m = m 1 + m 2 is the total gravitational mass of the system and the parameter α = M 1 M 2 /m 1 m 2 characterizes the Nordtvedt effect. Equation (238) coincides in the post-Newtonian approximation with similar formula derived by Bernard [80, eq. 7] after accounting for the definitions of the PPN parameters and the sensibility of the bodies. The spin dipole part of the Lagrangian, L ds , has been derived by Brumberg [24] in general relativity who used the harmonic gauge but replaced the body's accelerations with the Newtonian equations of motion in the Lagrangian. It was Damour [52] who included the acceleration-dependent part to the Lagrangian L ds in general relativity while Barker and O'Connell [53] extended this Lagrangian (230) to the case of scalar-tensor theory of gravity and developed the method of the double zero to eliminate the acceleration-dependent terms from the Lagrangian without imposing an implicit gauge transformation. Barker and O'Connell [54] have also analyzed the structure of the spin-dipole Lagrangian under different spin complementary conditions and in different post-Newtonian coordinates.

D. Quadrupole approximation
The quadrupole approximation includes all terms depending on mass and spin quadrupole moments of the bodies. The quadrupole Lagrangian consists of three pieces L q = L qm + L qs + L qa .
Here, the first term in (239) n p n q R 2 + ( 1 2 − β) (3m 1 + 4m 2 ) m 2 M pk 1 + (3m 2 + 4m 1 ) m 1 M pk 2 n p n q R 4 − 3γ m 2 1 M pk 2 + m 2 2 M pk 1 n p n q R 4 − The second term in the right hand side of (239) is the spin quadrupole Lagrangian that coincides exactly with the corresponding Lagrangian derived recently by Banihashemi and Vines [130,Eq. 29]. The third term in the right hand side of (239) describes the acceleration-dependent part of the mass quadrupole Lagrangian which can be reduced to a "canonical" form depending only on the coordinates and the first time derivatives of the dynamic variables. More specifically, the terms depending quadratically on accelerations in the Lagrangian (242) can be eliminated by applying the method of the double zero having been worked out by Barker and O'Connell [53]. Let us introduce two functions where we have used the monopole approximation for gravitational force as it is sufficient in the quadrupole approximation. These functions vanish on-shell but they are non-vanishing if the equations of motion are not implied. Equations (243) allows us to recast the Lagrangian (242) to the following form n q R 2 .
The products Z p 1 Z q 1 and Z p 2 Z q 2 are the double-zero terms [45,53] with the property that each factor in the product vanishes on-shell independently one from another. Therefore, the double-zero terms do not contribute to the equations of motion after taking the variational derivative and can be discarded from the Lagrangian.
Terms which are linear with respect to accelerations in the Lagrangian can be eliminated with the help of the Leibniz chain rule by interchanging time derivatives in each term of (244) which depends on acceleration linearly. For example, and similar procedure is applied to other terms in (244). After employing (245) the total time derivative appears in the Lagrangian (244) but it has no effect on the equations of motion of the bodies and, as such, can be discarded. The remaining terms do not depend on accelerations and can be combined with similar terms in L qm . The elimination of the acceleration-dependent terms from the Lagrangian may be useful under certain circumstances and for comparison with the Lagrangians derived by other authors.

E. Comparison of the quadrupole Lagrangians of different authors
First attempt to derive the post-Newtonian Lagrangian of an N -body system with accounting for the quadrupole moments of the bodies had been undertaken by Brumberg [24] who conducted calculations in global coordinates. Brumberg used the Newtonian definition of the center of mass of the bodies and the multipoles defined in the global coordinates. For these reason, Brumberg's Lagrangian is expressed in terms of not physically-measurable quantities and has a lot of coordinate-dependent effects which are formidably difficult to track down. In particular, reduction of Brumberg's Lagrangian to the case of spherically-symmetric bodies lead to the appearance of the finite-size effects depending on the moments of inertia of the bodies but they can be eliminated after transformation of body's multipoles to the local coordinates as pointed out by Vincent [110] and detailed in [50, section 6.3.3]. The unphysical definitions of the center of mass and body's multipoles inhibit the range of practical applications of Brumberg's Lagrangian in astrophysics and relativistic celestial mechanics.
More recently, Xu et al. [131], Wu et al. [132] have striven to derive the post-Newtonian Lagrangian in quadrupole approximation by applying the DSX-formalism [64,72] of the local coordinates in the N -body problem. However, they obtained equations of motion that significantly diverge from two subsequent and completely independent analytic derivations by Racine and Flanagan [15], Racine et al. [16] and Kopeikin [18] who came up to exactly the same result. We could not track down the origin of the disagreement and do not comment on the Lagrangian corresponding to the equations derived by Xu et al. [131], Wu et al. [132].
The post-Newtonian Lagrangian L VF qm of two-body problem in the mass quadrupole approximation has been worked out by Vines and Flanagan [56] from the equations of motion having been derived previously in [15,16] for N bodies possessing all mass and spin multipoles. Later on, the quadrupole Lagrangian L S qm has been derived by Steinhoff et al. [58] in the framework of the effective field theory [89]. It was shown [58] that the two Lagrangians differ by a total time derivative where the function F is shown in [58,Eq. 3.13]. Since the difference between the two Lagrangians is a total time derivative, the Lagrangians L VF and L S are effectively the same. For computational reasons the Lagrangian L S is more convenient for the comparison with our Lagrangian L which will be done under simplifying assumption that the first body is characterized by both mass and quadrupole moment while the second body has only mass. We shall also reduce our expressions to general relativity so that the parameters β = γ = 1 yielding η = 0. The Lagrangian L S reads where L pm is the point particle (monopole) Lagrangian, the quadrupole potential V EQ is given in [58,Eq. 3.11], and we neglected for the sake of simplification the spin and quadrupole-quadrupole (tidal) terms. Following the paper [58] we denote masses of the first and second body, m 1 and m 2 , and the quadrupole moment of the first body Q ij 8 . The second body is assumed to be point-like mass m 2 having no internal structure. Then, the monopole Lagrangian − 3 4 m 2 Q pq n p n q R 3 3v 2 1 − 7 (v 1 · v 2 ) + 3v 2 2 − 5 (n · v 1 ) (n · v 2 ) + m 2 Q pq n p n q R 4 where notations for coordinates, velocities and accelerations are the same as in our equation (240). Our Lagrangian has the monopole term where T A is the total kinetic energy of the internal motions of matter of body A, U A is the total energy of gravitational self-interaction of body A, and P A is the free thermodynamic energy of matter of body A. Making use of microscopic equations of motion (A1) of matter of body A, it is straightforward to prove a scalar virial theorem relating the above three quantities, by taking time derivatives from each of them and adding the derivatives all together. It yields the law of the time evolution of the general-relativistic mass of body A defined in (42). After accounting for (A5), we havė where the multipole moments M L A are defined in (32). This equation is employed in section V A for computing the time derivative from the inertial mass of body A.

Vector Virial Theorem
We consider time derivatives from two non-canonical dipoles Let's introduce the following notations for vector integrals Vector virial theorem states [18,Eq. 169] that This equation is used in computation of an active dipole moment in section VIII C.

Tensor Virial Theorem
Let us now introduce tensor integrals Their scalar counterparts obtained by contraction of the two indices are and the remaining two scalars are given above in (A6), (A7). Now, we take a second time derivative from I ij A and use (A1). We get a tensor virial theorem Contracting over two indices yields a scalar virial theorem for the central moment of inertia After substituting the gravitoelectric potential Q L from (21), and reshuffling terms in (A22) and (A23) we obtain more explicit tensor and scalar virial theorem which connects the internal and external dynamic variables. The scalar and tensor virial theorems are used for calculation of the Lagrangian in section VI D and in section VIII E for comparison of the Lagrangians derived by different authors.