Luc Blanchet

Jérome Novak

External field effect of modified Newtonian dynamics in the Solar system

Keywords: planets and satellites:general, galaxies:kinematics and dynamics, methods:numerical, dark matter

The Modified Newtonian Dynamics (MOND) has been formulated as a modification of the Poisson equation for the Newtonian gravitational field. This theory generically predicts a violation of the strong version of the equivalence principle, and as a result the gravitational dynamics of a system depends on the external gravitational field in which the system is embedded. This so-called external field effect has been recently shown to imply the existence of an anomalous quadrupolar correction, along the direction of the external galactic field, in the gravitational potential felt by planets in the Solar System. In this paper we confirm the existence of this effect by a numerical integration of the MOND equation in the presence of an external field, and compute the secular precession of the perihelion of planets induced by this effect. We find that the precession effect is rather large for outer gaseous planets, and in the case of Saturn is comparable to published residuals of precession obtained by Saturn range tracking data. The effect is much smaller for inner planets, but in the case of the Earth it appears to be in conflict for most of the MOND functions µ(y) with the very good constraint on the perihelion precession obtained from Jupiter VLBI data. The MOND functions that are compatible with this constraint appear to have a very rapid transition from the MONDian regime to the Newtonian one.

INTRODUCTION

Motivation

The Modified Newtonian Dynamics (MOND) -Milgrom (1983a,b,c)-is a successful alternative framework for interpreting the galactic rotation curves and the empirical Tully-Fisher relation without relying on dark matter halos (see Sanders & McGaugh 2002 for a review). At the nonrelativistic level, the best modified-gravity formulation of MOND is the modified Poisson equation originally proposed by [START_REF] Bekenstein | [END_REF],

∇ • µ g a0 ∇U = -4πGρ , ( 1
)
where ρ is the density of ordinary (baryonic) matter, U is the gravitational potential, g = ∇U is the gravitational ⋆ E-mail: blanchet@iap.fr † E-mail: jerome.novak@obspm.fr field and g = |g| its ordinary Euclidean norm. The modification of the Poisson equation is encoded in the MOND function µ(y) of the single argument y ≡ g/a0, where a0 = 1.2 × 10 -10 m/s 2 denotes the MOND constant acceleration scale. The MOND function interpolates between the MOND regime corresponding to weak gravitational fields y = g/a0 ≪ 1, for which it behaves as µ(y) = y + o(y), and the Newtonian strong-field regime y ≫ 1, where µ reduces to 1 so that we recover the usual Newtonian gravity.

Relativistic extensions of MOND modifying general relativity have been proposed (Bekenstein 2004;Sanders 2005) and extensively studied (see [START_REF] Bekenstein | the Proc. of the 28th John Hopkins Workshop[END_REF] for a review). Alternatively the MOND equation (1) can be reinterpreted as a modification of dark matter rather than gravity by invoking a mechanism of "gravitational polarisation" -a gravitational analogue of the electric polarisation (Blanchet 2007); this yields to the concept of dipolar dark matter which has been formulated as a relativistic model in standard general relativity by Blanchet & Le Tiec (2008, 2009).

An important consequence of the non-linearity of c 2010 RAS Eq. (1) in the MOND regime, is that the gravitational dynamics of a system is influenced (besides the well-known tidal force) by the external gravitational environment in which the system is embedded. This is known as the external field effect (EFE), which has non-trivial implications for non-isolated gravitating systems. The EFE was conjectured to explain the dynamics of open star clusters in our galaxy (Milgrom 1983a,b,c), since they do not show evidence of dark matter despite the involved weak internal accelerations (i.e. below a0). The EFE effect shows that the dynamics of these systems should actually be Newtonian as a result of their immersion in the gravitational field of the Milky Way. The EFE is a rigorous prediction of the Bekenstein-Milgrom equation (1), and is best exemplified by the asymptotic behaviour of the solution of (1) far from a localised matter distribution (say, the Solar System), in the presence of a constant external gravitational field ge (the field of the Milky Way). At large distances r = |x| → ∞ we have [START_REF] Bekenstein | [END_REF])

U = ge • x + GM/µe r 1 + λe sin 2 θ + O 1 r 2 , ( 2
)
where M is the mass of the localised matter distribution, where θ is the azimuthal angle from the direction of the external field ge (see also Fig. 1), and where we denote µe ≡ µ(ye) and λe ≡ yeµ ′ e /µe, with ye = ge/a0 and µ ′ e = dµ(ye)/dye. In the presence of the external field, the MOND internal potential u ≡ Uge • x shows a Newtonian-like falloff ∼ r -1 at large distances but with an effective gravitational constant G/µe. 1 However, contrary to the Newtonian case, it exhibits a non-spherical deformation along the direction of the external field. The fact that the external field ge does not disappear from the internal dynamics can be interpreted as a violation of the strong version of the equivalence principle. 2 For the reader's convenience we derive the result (2) in the Appendix A.

In a recent paper, Milgrom (2009) has shown that the imprint of the external galactic field ge on the Solar System (due to a violation of the strong equivalence principle) shows up not only asymptotically, but also in the inner regions of the system, where it may have implications for the motion of planets. This is somewhat unexpected because gravity is strong there (we have g ≫ a0) and the dynamics should be Newtonian. However, because of the special properties of the equation (1), the solution will be given by some modified Poisson-type integral, and the dynamics in the strong-field region will be affected by the anomalous behaviour in the asymptotic weak-field region. Thus the results apply only to the nonlinear Poisson equation (1), not to modified inertia formulations of MOND. The anomaly expresses itself as an anomalous quadrupolar contribution to the MONDian internal field u, as compared to the Newtonian potential, given 1 Recall that in the absence of the external field the MOND potential behaves like U ∼ -√ GM a 0 ln r, showing that there is no escape velocity from an isolated system (Famaey et al. 2007). However since no object is truly isolated the asymptotic behaviour of the potential is always given by (2), in the approximation where the external field is constant. 2 However the weak version of the equivalence principle, that all test particles in the MOND gravitational field have universal acceleration a = g, remains satisfied [START_REF] Bekenstein | [END_REF].

by δu = 1 2 x i x j Qij , (3)
where x i is the distance from the Solar System barycentre, and Qij is a trace-free quadrupole moment of the type

Qij = Q2 e i e j - 1 3 δij , (4)
with Q2 being the quadrupole coefficient and e = (e i ) the preferred direction of the galactic centre (i.e. e = ge/ge).

The anomalous term (3) is a harmonic solution of the Laplace equation, i.e. ∆δu = 0.

The radial dependence of this anomaly is ∝ r 2 and can thus be separated from a quadrupolar deformation due to the Sun's oblateness which decreases like ∝ r -3 . This result is valid whenever r is much less than the MOND transition distance for the Solar System, defined by r0 = GM/a0 with M the mass of the Sun and a0 the MOND acceleration scale. This radius corresponds to the transition region where the Newtonian acceleration becomes of the order of the MOND acceleration a0 and therefore, MOND effects become dominant. We have r0 ≃ 7100 AU so the effect (3)-( 4) is valid in a large volume around the Sun including all the planets (recall that Neptune's orbit is at 30 AU).

In addition to (3) we have also the usual MOND effect. In the special case of spherical symmetry, the equation (1) reduces to µ(g/a0)g = gN, with gN the ordinary Newtonian gravitational field. For a MOND function behaving like µ(y) ∼ 1 -ǫy -γ when y → ∞ (i.e. in the region r ≪ r0), the Newtonian field is modified by δgN ≃ ǫa0(r/r0) 2γ-2 . For γ = 1 this corresponds to a constant acceleration similar to the Pioneer anomaly, and for γ 2 the effect is very small and the motion of planets is almost unaffected. This effect is independent of the presence of the external field, and is spherical, so can be distinguished from the quadrupolar deformation (3) induced by the external field; we neglect this effect here.

Summary

In the present paper we shall confirm the existence of the effect (3)-( 4) and shall derive it by an independent numerical integration of the equation (1) using an elliptic solver originally built for numerical relativity purposes. 3 The magnitude of the quadrupole coefficient Q2 depends on the dimensionless ratio between the external field ge and the MOND acceleration a0, say

η = ge a0 , (5)
and on the particular MOND interpolating function µ in use. For the Milky Way field at the level of the Sun we have ge ≃ 1.9 × 10 -10 m/s 2 which happens to be slightly above the MOND scale, i.e. η ≃ 1.6. Our calculation implemented for various standard choices for the µ-function gives

2.1 × 10 -27 s -2 Q2 4.1 × 10 -26 s -2 . ( 6
)
The quadrupole coefficient Q2 is found to be positive corresponding to a prolate elongation along the quadrupolar axis. We shall verify numerically that Q2 is indeed constant over a large portion of the inner Solar System, and starts decreasing when the distance becomes comparable to the MOND transition scale r0, at a few thousands AU from the Sun. We shall also check that Q2 scales approximately like the inverse square root of the central mass, in agreement with dimensional analysis which shows that we should approximately have Q2 = q2 a0/r0, where q2(η) is a dimensionless quadrupole coefficient depending on the ratio η between the external field and the MOND acceleration. We find that for η ≃ 1.6, q2 ranges over 0.02 q2 0.36, in good agreement with the values found by Milgrom (2009).

We then study some consequences of the anomalous quadrupole moment for the motion of planets, notably the secular advances of planetary perihelia. Applying standard perturbation equations of celestial mechanics we obtain the anomalous precession rate

∆2 = dω dt = Q2 √ 1 -e 2 4n 1 + 5 cos(2ω) . (7)
To simplify, we have assumed that the direction of the galactic centre lies in the plane of the ecliptic (this is only approximately correct) and that all planets move in this plane. We defined the origin of the precession angle ω to be the direction of the galactic centre. The orbital frequency of the planet is n = 2π/P (P is the orbital period), and a, e and ω = ω + Ω denote the standard orbital elements. The effect seems to be the most interesting for Saturn for which we find 0.3 mas/cy ∆2 5.8 mas/cy. ( 8)

Those values are within the range of published residuals for the precession of Saturn as obtained from global fits of the Solar System dynamics in Pitjeva (2005) and [START_REF] Fienga | Relativity in Fundamental Astronomy: Dynamics, Reference Frames and Data Analysis[END_REF]. However we shall find that in the case of the Earth the predicted perihelion precession may be greatly constrained by the best estimates of postfit residuals obtained thanks to the Jupiter VLBI data using the INPOP planetary ephemerides [START_REF] Fienga | Relativity in Fundamental Astronomy: Dynamics, Reference Frames and Data Analysis[END_REF]. This is turn reduces the possible choices for the MOND function µ(y) to those that exhibit a sharp transition from the Newtonian to MONDian regime.

To conclude, we present very accurate numerical computations performed within the well-defined theory (1), and compare the results with recent observations of the motion of planets in the Solar System. We find that the observational constraints are rather strong, and may even conflict with some of the predictions. Thus our results provide motivation for investigating more systematically possible anomalous effects in the Solar System predicted by alternative theories. In particular, the external field effect associated with the preferred direction of the galactic centre could be seen as a typical one arising in generic attempts to modify gravity with motivation coming from dark matter and/or dark energy. Indeed we expect that generic departures from general relativity will violate the strong equivalence principle [START_REF] Will | Theory and experiment in gravitational physics[END_REF].

On the other hand let us cautiously remark that MOND and more sophisticated theories such as TeVeS (Bekenstein 2004), which are intended to describe the weak field regime of gravity (below a0), may not be extrapolated without modification to the strong field of the Solar System. For instance it has been argued by [START_REF] Famaey | [END_REF] that a MOND interpolating function µ which performs well at fitting the rotation curves of galaxies is given by µ1 defined by (34a) below. However this function has a rather slow transition to the Newtonian regime, given by µ1 ∼ 1 -y -1 when y = g/a0 → ∞, which is already excluded by Solar System observations (Sereno & Jetzer 2006). Indeed such slow fall-off -y -1 predicts a constant supplementary acceleration directed toward the Sun δgN = a0 (i.e. a "Pioneer" effect), which is ruled out because not seen from the motion of planets. Thus it could be that the transition between MOND and the Newtonian regime is more complicated than is modelled in Eq. ( 1). This is also true for the modified dark matter model (Blanchet & Le Tiec 2008, 2009) which may only give an effective description valid in the weak field limit and cannot be extrapolated as it stands to the Solar System. While looking at MOND-like effects in the Solar System we should keep the previous proviso in mind. The potential conflict we find here with the Solar System dynamics (notably with the Earth orbital precession) may not necessarily invalidate those theories if they are not "fundamental" theories but rather "phenomenological" models.

The plan of this paper is the following. In Sec. 2, we derive an expression of the solution of the MOND equation near the origin in terms of multipole moments. Section 3 is devoted to numerical techniques and results, with a description of the particular formulation and assumptions used for the numerical integration in Sec. 3.1, of the various tests passed by the code (Sec. 3.2), and the numerical results and values for the first multipoles of the potential in Sec. 3.3. Section 4 details the consequences of this modified gravitational potential on the orbits of Solar-system planets, with first the derivation of the perturbation equations in Sec. 4.1 and numerical values for the planetary precessions (Sec. 4.2). Finally, Section 5 summarises the results and gives some concluding remarks.

MULTIPOLAR EXPANSION OF THE MONDIAN POTENTIAL

In this paper we shall solve the modified Poisson equation (1) with the boundary condition given by the constant external gravitational field ge (defining a preferred spatial direction denoted e = ge/ge), i.e. that the gravitational field g = ∇U asymptotes to limr→∞ g = ge. The external field ge should consistently obey a MOND equation, but here we shall simply need to assume that ge is constant over the entire Solar System. The MONDian physicist measures from the motion of planets relatively to the Sun the internal gravitational potential u defined by

u = U -ge • x , (9)
which is such that limr→∞ u = 0. Contrary to what happens in the Newtonian case, the external field ge does not disappear from the gravitational field equation (1) and we want to investigate numerically its effect. The anomaly detected by a Newtonian physicist with respect to the MONDian physicist is the difference of internal potentials,

δu = u -uN , ( 10
)
where uN denotes the ordinary Newtonian potential generated by the same ordinary matter distribution ρ, and thus solution of the Poisson equation ∆uN = -4πGρ with the boundary condition that limr→∞ uN = 0. Hence uN is given by the standard Poisson integral

uN(x, t) = G d 3 x ′ |x -x ′ | ρ(x ′ , t) . (11)
A short calculation shows that the anomaly obeys the ordinary Poisson equation ∆δu = -4πGρ pdm , where ρ pdm is the density of "phantom dark matter" defined by

ρ pdm = 1 4πG ∇ • (χ∇U ) , (12)
where we denote χ ≡ µ -1. The phantom dark matter represents the mass density that a Newtonian physicist would attribute to dark matter. In the model by Blanchet & Le Tiec (2008, 2009) the phantom dark matter is interpreted as the density of polarisation of some dipolar dark matter medium and the coefficient χ represents the "gravitational susceptibility" of this dark matter medium.

The Poisson equation ∆δu = -4πGρ pdm is to be solved with the boundary condition that limr→∞ δu = 0; hence the solution is given by the Poisson integral

δu(x, t) = G d 3 x ′ |x -x ′ | ρ pdm (x ′ , t) . (13)
We emphasise here that, contrary to the Newtonian (linear) case, the knowledge of the matter density distribution does not allow to obtain any analytic solution for the potential; however, we can still infer the structure of the multipolar expansion near the origin, and the moments will be computed numerically. We can check that the phantom dark matter behaves like r -3 when r → ∞, so the integral ( 13) is perfectly convergent. We want to discuss the influence of the external galactic field in the inner part of the Solar System where the gravitational field is strong (g ≫ a0); thus µ tends to one there, so χ tends to zero. For the discussion in this section we adopt the extreme case where χ is exactly zero in a neighbourhood of the origin, say for r ε, so that there is no phantom dark matter for r ε; in later sections devoted to the full numerical integration we shall still make this assumption by posing χ = 0 inside the Sun (in particular we shall always neglect the small MOND effect at the centre of the Sun where gravity is vanishingly small). If ρ pdm = 0 when r ε we can directly obtain the multipolar expansion of the anomalous term (13) about the origin by Taylor expanding the integrand when r = |x| → 0. In this way we obtain

δu = +∞ l=0 (-) l l! x L QL , (14)
where the multipole moments near the origin are given by 4

QL = G r>ε d 3 x ρ pdm ∂L 1 r , (15)
4 Our notation is as follows:

L = i 1 • • • i l denotes a multi-index composed of l multipolar spatial indices i 1 , • • • , i l (ranging from 1 to 3); ∂ L = ∂ i 1 • • • ∂ i l is the product of l partial derivatives ∂ i ≡ ∂/∂x i ; x L = x i 1 • • • x i l is the product of l spatial positions x i ; similarly n L = n i 1 • • • n i l = x L /r l is the product of l unit vectors n i = x i /r; the symmetric-trace-free (STF) projection is indicated with a hat, for instance xL ≡ STF[x L ],
and similarly for nL and

Because the integration in ( 15) is limited to the domain r > ε and ∂L(1/r) is symmetric-trace-free (STF) there [indeed ∆(1/r) = 0], we deduce that the multipole moments QL themselves are STF. This can also be immediately inferred from the fact that ∆δu = 0 when r ε, hence the multipole expansion ( 14) must be a homogeneous solution of the Laplace equation which is regular at the origin, and is therefore necessarily made solely of STF tensors of type xL . Hence we can replace x L in ( 14) by its STF projection xL . It is clear from the formula (15) that the MONDian gravitational field (for r r0) can influence the near-zone expansion of the field when r → 0.

With the expression ( 12) for the phantom dark matter and the MOND equation ( 1), we can further transform QL as

QL = - 1 4π r>ε d 3 x 4πGρ + ∆U ∂L 1 r . ( 16
)
Approximating the central matter distribution as being spherically symmetric (i.e. ignoring the "back-reaction" of the non-spherical anomalous field δu on the matter which generates the field), we see that the first term is non-zero only in the monopolar case l = 0 where it reduces in the limit ε → 0 to minus the Newtonian potential evaluated at the origin. On the other hand the second term in ( 16) can be transformed as a surface integral. We find

QL = -uN(0)δ l,0 + 1 4π r=∞ r=ε d 2 Si U ∂iL 1 r -∂iU ∂L 1 r .
(17) Our notation means that the surface integral is composed of two contributions, an inner one at r = ε (denoted Q ε L ) and an outer one at infinity r = ∞ (say Q ∞ L ). In Eq. ( 17) we implicitly assume that the limit ε → 0 is to be taken after evaluating the inner surface integral.

Inserting U = u + ge • x into the surface integral at infinity Q ∞

L we find that it contributes only to the dipolar term and reduces to the external galactic field; thus

Q ∞ i = g i e
with all other Q ∞ L 's being zero. On the other hand the inner surface integral Q ε L is found to have a well-defined limit when ε → 0 given by Q ε L = (-) l ( ∂LU )(0), where we recall that ∂L denotes the STF part of ∂L = ∂i 1 • • • ∂i l . We then find that the galactic field in U = u + ge • x cancels the dipole term in the surface integral at infinity, so that the result is

QL = -uN(0)δ l,0 + (-) l ( ∂Lu)(0) . ( 18
)
The internal MONDian potential u admits therefore the following STF multipolar expansion (equivalent to a near-zone expansion when r → 0)

u = uN -uN(0) + +∞ l=0 1 l! x L ( ∂Lu)(0) . (19)
In Appendix B we derive an alternative proof of this result.

At this stage we have elucidated the structure of the multipole expansion of the anomaly δu near the origin, but still we need to resort to a numerical integration of the non-linear MOND equation in order to obtain quantitative values for the multipole moments. Section 3 will be devoted to this task.

Finally we show that the dipole moment Qi (with l = 1) is actually zero (Milgrom 2009). This follows from the weak equivalence principle satisfied by the Bekenstein-Milgrom theory. As a consequence, the total acceleration of the centre of mass of the Solar System in the galactic external field ge does not depend on the Solar System's internal dynamics and is simply given by ge (for a detailed proof see [START_REF] Bekenstein | [END_REF]. The centre of mass of the Solar System does not differ much from that of the Sun, so we deduce that the "self-acceleration" of the Sun, i.e. the total acceleration due to the internal potential u when integrated over the whole volume of the Sun, must be (approximately) zero:

d 3 x ρ ∂iu = 0 . ( 20
)
We have obviously the same result for the Newtonian potential uN in Newtonian gravity so the same will be true for the anomalous field δu = u -uN. In other words the phantom dark matter which is the source for the anomaly exerts no net force on the Sun and we get, for a spherically symmetric star,

Qi = - 1 M d 3 x ρ ∂iδu = 0 . (21)
In Sec. 3.3 we shall numerically verify that the dipole moment Qi is indeed zero within our numerical error bars.

NUMERICAL INTEGRATION OF THE MOND EQUATION

Formulation and implementation

Using spherical spatial coordinates {r, θ, ϕ}, we consider a star represented by a spherically symmetric distribution of matter ρ(r) obtained from a hydrostatic equilibrium model in Newtonian theory (polytrope). These models are obtained by integrating the hydrostatic spherical equilibrium equation (relating the pressure p and the gravitational field) and the Newtonian equation for the gravitational field (relating the gravitational field and the density distribution ρ), together with a polytropic equation of state (EOS) of the form p = Kρ γ , where K and γ are two constants. As it shall be shown later in Sec. 3.2, the results do not depend on the specific EOS used to obtain this hydrostatic equilibrium, so we do not discuss the particular EOS used for obtaining the density distribution. It also means that we neglect the effect of MOND theory on the structure of the star. This is not a severe restriction, since inside the star the gravitational field is much higher in amplitude than the constant a0, making µ ∼ 1 and Newtonian theory is recovered up to very high accuracy. We then solve the MOND equation ( 1), with the boundary conditions given by the constant galactic gravitational field ge (see Fig. 1), i.e. limr→∞ g = ge = gee.

In order to be closer to a Poisson-like form of the partial differential equation (PDE), we are solving in terms of the internal potential u = Uge • x such that limr→∞ u = 0. Defining h = ∇u = g-ge, the MOND equation ( 1) becomes the following PDE:

∆u = 1 µ(g/a0) -4πGρ - µ ′ (g/a0) a0g g i g j ∂ihj , (22)
with µ ′ ≡ dµ/dy being the derivative of the function µ with respect to its argument y = g/a0. This PDE is solved numerically, using the library lorene which implements spectral methods (for a review in the case of numerical relativity see Grandclément & Novak 2009) in spherical coordinates. In our case (see Fig. 1), the fields do not depend on the azimuthal angle ϕ and the problem is axisymmetric. Since Eq. ( 22) is a non-linear elliptic PDE, the algorithm used is the fixed point iteration method in which one starts from an initial guess u0(r, θ), here the solution of the Newtonian Poisson equation ∆u0 = -4πGρ. Knowing un(r, θ) at a given iteration step n, one computes the non-linear source terms in the right-hand-side of Eq. ( 22), say σ(u, ρ), to obtain a new value of the potential un+1, solving a linear Poisson equation

∆ũn+1 = σ(un, ρ) , (23)
and using a relaxation technique

un+1 = λũn+1 + (1 -λ)un , (24)
with λ ∈ (0, 1] being the relaxation parameter (often taken to be 0.5). This iteration is stopped when the relative variation of the potential u between two successive steps becomes lower than a given threshold (usually 10 -12 ). The linear equation ( 23) is solved decomposing each field on a truncated base of spherical harmonics Y m l (θ, ϕ) (with m = 0 because the problem is axisymmetric) in the θ-direction, and Chebyshev polynomials Tn(r) in the r-direction. For this last coordinate, a multi-domain technique is used, with linear mappings of the r coordinate to the interval [-1, 1], except for the last domain where a mapping of the type s = 1/r allows to treat spatial infinity in our computational domain and to impose the right boundary conditions at r → ∞ (see Grandclément & Novak 2009 for details). Moreover, this multi-domain technique imposes that the potential u, as well as its first radial derivative, be continuous across the stellar surface.

Finally, the only modification made to the exact equation ( 22) in order to be numerically integrated is the setting µ(g/a0) = 1 inside the star. Indeed, there are two problems here. First, the MOND gravitational field g(r) does not admit a regular Taylor expansion near the origin, as does the density in Eq. ( B6) and the associated Newtonian gravitational field which follows as

gN = 4πG +∞ k=0 r 2k+1 (2k + 3)(2k + 1)! (∆ k ρ)(0) . ( 25
)
This can be immediately seen from the fact that in spherical symmetry the MOND equation ( 1) reduces to µ(g/a0)g = gN. Using (25) we see that the expansion of g when r → 0 starts with a term proportional to the square root of r and is therefore not regular; more precisely we have

g = 4πGρ0a0r 3 1 + O(r 2 ) . (26)
With our polynomial decomposition of fields in the radial direction, where we decompose in Chebyshev polynomials Tn(r), we therefore get an incompatibility when trying to solve the Poisson equation ( 23) in the vicinity of the origin at the centre of the Sun. The second problem is that, at the very centre of the star g → 0 and therefore µ → 0, which makes the division present in the right-hand-side of Eq. ( 22) difficult to handle numerically.

Here we have chosen to avoid both problems by imposing µ = 1 in the star so it is entirely Newtonian. This approximation may produce a noticeable change in the value of µ(y) only within a sphere of radius r 3a 0 4πGρ 0 ∼ 10 -15 R⊙, where ρ0 is the central density and R⊙ is the radius of the Sun. It is thus supposed to have a completely negligible effect on the results.

The result which is finally used to study the influence of the modification of the Newtonian gravity on the orbits of planets is the value of the trace-free multipole moments QL defined in Sec. 2. In the case where all the multipole moments are induced by the presence of the external field ge in the preferred direction e, the situation is axisymmetric and all the moments QL will have their axis pointing in that direction e, and we can define the multipole coefficients Q l by

QL = Q l êL , (27)
where êL denotes the STF part of the product of l unit vectors e L = e i 1 • • • e i l . Then the multipole expansion ( 14)

reads as 5 δu(r, θ) = +∞ l=0 (-) l (2l -1)!! r l Q l (r) P l (cos θ) , (28)
where P l is the usual Legendre polynomial and θ is defined in Fig. 1. Although from the considerations of Sec. 2 the multipole moments Q l should be approximately constant within the MOND transition radius r0, here we compute them directly from the numerical solution of (1) and shall obtain their dependence on r; we shall check numerically that Q2(r) and Q3(r) are indeed almost constant in a large sphere surrounding the Solar system. With our definition the monopole, quadrupole and octupole pieces in the internal field are given by

δu1 = -r Q1 cos θ , ( 29a
) δu2 = 1 2 r 2 Q2 cos 2 θ - 1 3 , ( 29b
) δu3 = - 1 6 r 3 Q3 cos θ cos 3 θ - 3 5 cos θ . ( 29c
)
From dimensional arguments, it is possible to see that Q l should scale as Q l ∼ a0/r l-1 0 [see for instance Eqs. ( 33) and ( 37)] and therefore

6 δu ∼ l Q l r l ∼ a0r0 l r r0 l . (30)
This last expression shows that higher-order multipoles should have lower influence on Solar-system planets, for which r ≪ r0. This shall be exemplified with the influence of the octupole, as compared to that of the quadrupole, in Sec. 4.2.

Tests of the code

A first test of our code is the standard comparison between the Laplace operator applied to the potential u, and the right-hand-side of Eq. ( 22). In all results displayed here, the maximum of this error has always remained lower than 10 -11 and shall not be considered as an interesting error indicator.

Next, an indication of the error comes from the use of the Gauss theorem:

µ g a0 g • d 2 S = -4πGM, ( 31
)
where M is the mass of the star, computed from initial data at a very high precision. Then, another check (which is however not a priori independent) is the asymptotic behaviour of the potential u(r, θ) when r → ∞, which is 34b)] and the standard value a 0 = 1.2 × 10 -10 m.s -2 . The abscissa gives the ratio (5) between the norm of the asymptotic gravitational field ge and the parameter a 0 . Note that both curves do coincide.

u = GM rµe 1 + λe sin 2 θ + O 1 r 2 , (32)
where µe = µ(ye) and λe = yeµ ′ e /µe with ye = ge/a0 (see the proof in Appendix A). The two tests are not really independent because the Gauss theorem is obtained by integrating the asymptotic field over a sphere at infinity. The main interest of these tests is that while the field is computed asymptotically on a sphere at infinity, the mass M is obtained "locally" by integrating numerically the density over the volume of the star. We emphasise that in our code the correct asymptotic behaviour (32) comes out directly and needs not to be imposed by hands at the beginning of the calculation.

From the use of the mapping s = 1/r (see Sec. 3.1), it is numerically possible to check the angular dependence of the quantity r × u when r → ∞ and thus to estimate independently the error on the potential u. Both tests -Gauss theorem (31) and asymptotic behaviour (32) -have been performed and have given accuracy levels of ∼ 10 -3 -10 -4 . Some examples are given in Fig. 2, for a sequence of different values for the asymptotic gravitational field ge. Both tests give the same numbers, up to five digits. In all results shown hereafter, the accuracy level defined by these two tests is better than 10 -2 .

Among the parameters entering the numerical model, we have checked that the details of the density profile ρ(r) do not contribute to the value of the quadrupole moment Q2(r). As an illustration, we give here the quadrupole obtained with two different density profiles, displayed in Fig. 3, both giving a star of exactly one solar mass. With the same other parameters [choices of µ(y), a0 and ge], the relative difference in the induced quadrupole is 4 × 10 -3 , while the error indicators give a relative numerical uncertainty of ∼ 2 × 10 -3 . We therefore consider that the particular form of the density profile used to model the star (i.e., the choice of EOS) has little influence on the results. Both profiles yield a one-solar mass star in hydrostatic equilibrium. The two profiles yield very little (same level as the overall numerical accuracy ∼ 10 -3 ) difference on our results.

Results on the induced quadrupole and octupole

In what follows, unless specified, the mass of the star is that of the Sun. As a first result, we show in Fig. 4 the profile of the quadrupole induced by the MOND theory through the function Q2(r) defined in Eq. (29b). We find that this quadrupole is decreasing from the star's neighbourhood to zero, on a typical scale of 10000 astronomical units (AU). However, this quadrupole can be considered as almost constant closer to the Sun, as it has a relative variation lower than 10 -4 within 30 AU (see the zoomed region in Fig. 4). We shall therefore refer to the quadrupole as a simple number, noted Q2(0) or simply Q2, when evaluating its influence on the orbits of Solar-system planets. Our numerical results for the quadrupole are given in Table 1.

We now briefly study other multipole moments, namely for l = 1 and l = 3. The profiles for the dipole Q1(r) and octupole Q3(r), defined by (29a) and (29c), are displayed in Fig. 5. We find that near the Sun and the solar planets the dipole can be considered as zero, since it is three orders of magnitude lower than the typical value for an acceleration in the problem (i.e. a0 or ge). Indeed we expect on dimensional analysis that Q1 should scale with the MOND acceleration a0 [see (30)],

Q1 = a0 q1(η) , (33)
where the dimensionless coefficient q1 depends on the ratio η between ge and a0 as defined by Eq. ( 5). Our values given in Table 1 show that q1 is actually very small. It means that Q1 is lower than the numerical error and confirms the analytical proof given in Sec. 2 that the dipole moment is approximately zero. On the other hand, we find that the octupole tends to a non-zero value, because the corresponding q3, which has also an influence on the orbits of planets and ge = 1.9 × 10 -10 m.s -2 . The MOND transition radius r 0 = GM/a 0 is shown by a dash-dotted line at r ≃ 7100 AU. Right panel: zoom of the central region (r 50 AU), where the quadrupole is almost constant (the spikes for r → 0 are due to the high value of the monopole term, from which it is numerically difficult to extract the quadrupole). Note the difference in the y-axis scales.

Figure 5. Profiles of dipole Q 1 (r) (left) and octupole Q 3 (r) (right) in the solar system, with same settings as those of Fig. 4. The MOND transition radius r 0 = GM/a 0 is shown by a dash-dotted line at r ≃ 7100 AU. Numerically we get |Q 1 (0)| ∼ 9 × 10 -14 m.s -2 ≪ a 0 which shows that the dipole is in fact zero.

(see the next section), is found to be of the order of one. Numerical values of the octupole are also given in Table 1.

The dependence of the quadrupole moment Q2 upon the value of the galactic gravitational field is displayed in Fig. 6, for four different coupling functions µ(y). Here we consider four cases widely used in the literature:

µ1(y) = y 1 + y , ( 34a
)
µ2(y) = y 1 + y 2 , ( 34b
)
µexp(y) = 1 -e -y , (34c)
µTeVeS(y) = √ 1 + 4y -1 √ 1 + 4y + 1 . ( 34d
)
The function µ1 has been shown to yield good fits of galactic rotation curves [START_REF] Famaey | [END_REF]. However because of its slow transition to the Newtonian regime it is a priori incompatible with Solar System observations. The function µ2 is generally called the "standard" choice and was used in fits such as those of [START_REF] Begeman | [END_REF]. We here use the general notation for any positive integer n:

µn(y) = y n √ 1 + y n . ( 35
)
We include also the function µexp having an exponentially fast transition to the Newtonian regime. The fourth choice µTeVeS is motivated by the theory TeVeS (Bekenstein 2004).

One should note that none of these functions, except maybe the fourth one, derives from a fundamental physical principle.

One observes that, up to the standard value of ge = 1.9 × 10 -10 m.s -2 the four curves are quite close, giving an interval of values for Q2:

2.2 × 10 -26 s -2 Q2 4.1 × 10 -26 s -2 . ( 36
)
However, for other choices of the MOND function µ(y), Q2 can take lower values down to 2.1×10 -27 s -2 . In all cases, we find that Q2 is positive, which means a prolate deformation of the field toward the galactic centre. Two profiles of µ(y) seem to give a maximum for Q2 near the standard value of ge, namely the standard choice µ2 and the exponential choice µexp. For the two other choices of µ(y), we have not been able to increase the value of η = ge/a0 further than 10, because the error indicators would become too large and the results could not be trusted any longer. Note that in Fig. 6, a0 has been kept fixed, while we have varied ge. We have further explored the dependence of the EFE induced quadrupole Q2 on the type of MOND function.

We have used several different functions of type µn, as defined in Eq. ( 35), and the results for the variation of Q2 depending on n are displayed in Fig. 7. From this figure, one can notice that the value of Q2 decreases with n, that is with a faster transition from the weak-field regime (where µ(y) ∼ y) to the strong field regime (µ(y) ∼ 1). We have been unable to study higher values of n and to determine a possible limit for Q2 as n goes to infinity. In Table 1, we give values of the dipole, quadrupole and octupole near the Sun (r 50 AU), where they are observed to be constant, for the four different expressions (34) of the interpolating function.

We have also explored different values of a0 at fixed ge, with the results that Q2 was increasing for small a0 as a power-law, reaching a maximum value for a0 ≃ 10 -10 -10 -9 m • s -2 and then slowly decreasing. Finally, we have varied the mass of the central star. Results are displayed in Fig. 8 where, in particular, the M -1/2 dependence of the quadrupole moment is recovered. This is not a surprise since on general grounds one expects that the quadrupole moment (and similarly the octupole moment) should scale like (Milgrom 2009)

7 Q2 = a0 r0 q2(η) , (37a)
Q3 = a0 r 2 0 q3(η) , (37b)
where r0 = GM/a0 is the MOND transition radius and q2, q3 denote some dimensionless coefficients depending on the ratio η = ge/a0 and on the choice of the interpolating function µ. Having obtained the expected dependence on the mass in Fig. 8 is another check that our code behaves correctly, since we are able to recover the analytic results known for this system. The values of the dimensionless multipole coefficients q2 and q3, which are expectedly close to unity, are reported in Table 1. 7 To compare with the results of Milgrom (2009) for the quadrupole, one should note the different conventions in use, namely q Milgrom ij = -Q BN ij and q Milgrom (η) = -2 3 q BN 2 (η). Our results for the quadrupole are in good agreement with values given by Milgrom (2009).

Table 1. Numerical values of the dipole Q 1 , the quadrupole Q 2 and the octupole Q 3 together with the associated dimensionless quantities q 1 , q 2 and q 3 defined by Eqs. ( 33) and ( 37). All values are given near the Sun. We use different choices of the function µ(y) defined in Eqs. ( 34) and ( 35). These values have been obtained for a 0 = 1.2 × 10 -10 m • s -2 and ge = 1.9 × 10 -10 m • s -2 (and M = 1 M ⊙ ).

MOND function

µ 1 (y) µ 2 (y) µ 5 (y) µ 20 (y) µexp(y) µ TeVeS (y) Q 1 [m s -2 ]
-8.9 × 10 -14 -9.8 × 10 -14 -1.1 × 10 -13 -1.2 × 10 -13 -9.7 × 10 -14 -3.5 × 10 -14 q 1 -7.4 × 10 -4 -8.2 × 10 -4 -9.2 × 10 -4 -10 -3 -8.1 × 10 -4 -2.9 × 10 -4

Q 2 [s -2 ] 3.8 × 10 -26 2.2 × 10 -26 7.4 × 10 -27
2.1 × 10 -27 3.0 × 10 -26 4.1 × 10 -26 q 2 0.33 0.19 6.5 × 10 -2 1.8 × 10 -2 0.26 0.36

Q 3 [m -1 s -2 ] -1.2 × 10 -40 -9.3 × 10 -41 -4.9 × 10 -42 -2.3 × 10 -42 -1.2 × 10 -40 -1.1 × 10 -40 q 3 -1.1 -0.87 -4.6 × 10 -2 -2.1 × 10 -2 -1.1 -1

EFFECT ON THE DYNAMICS OF SOLAR SYSTEM PLANETS

Perturbation equations

In this section, we investigate the consequence for the dynamics of inner planets of the Solar System of the presence of an abnormal multipole moment QL oriented toward the direction e of the galactic centre, in the sense of Eq. ( 27).

Recall that the domain of validity of this anomaly is expected to enclose all the inner Solar System (for distances r r0 ≈ 7100 AU), with the quadrupole coefficient being constant up to say 50 AU (see Fig. 4). As we have seen, the anomaly induces a perturbation on the Newtonian gravitational potential, namely u = uN + δu, where the perturbation function δu is given for various multipole moments by Eqs. (29). Following the standard practice of celestial mechanics we denote the perturbation by R ≡ δu. The perturbation function R is such that the perturbing force (or, rather, acceleration) acting on the Newtonian motion is

F = ∇R.
The unperturbed Keplerian orbit of a planet around the Sun is described by six orbital elements (say cA with A = 1, • • • , 6). For these we adopt the semi-major axis a, the eccentricity e, the inclination I of the orbital plane, the mean anomaly ℓ defined by ℓ = n(t -T ) where n = 2π/P (n is the mean motion, P is the orbital period and T is the instant of passage at the perihelion), the argument of the perihelion ω (or angular distance from ascending node to perihelion), and the longitude of the ascending node Ω. We also use the longitude of the perihelion defined by ω = ω+Ω.

The perturbation is a function R(cA) of the orbital elements of the unperturbed Keplerian ellipse. With our choice for the orbital elements {cA} = {a, e, I, ℓ, ω, Ω} the perturbation equations of celestial mechanics read as (see e.g. [START_REF] Brouwer | Methods of celestial mechanics[END_REF]

da dt = 2 an ∂R ∂ℓ , (38a)
de dt = √ 1 -e 2 ea 2 n 1 -e 2 ∂R ∂ℓ - ∂R ∂ω , (38b)
dI dt = 1 a 2 n √ 1 -e 2 sin I cos I ∂R ∂ω - ∂R ∂Ω , (38c)
dℓ dt = n - 1 a 2 n 2a ∂R ∂a + 1 -e 2 e ∂R ∂e , (38d)
dω dt = 1 a 2 n √ 1 -e 2 e ∂R ∂e - cos I sin I √ 1 -e 2 ∂R ∂I , (38e)
dΩ dt = 1 a 2 n sin I √ 1 -e 2 ∂R ∂I . (38f)
We shall be interested only in secular effects, so we average these equations over one orbital period P ; denoting the time average by brackets we have

∂R ∂cA = 1 P P 0 dt ∂R ∂cA = 1 2π 2π 0 dℓ ∂R ∂cA . (39)
In particular this implies that we shall always have ∂R/∂ℓ = 0 hence the perturbation equation (38a) gives da/dt = 0 (at first order in perturbation theory). Indeed, a perturbative force of the type F = ∇R is conservative, so the energy of the orbit is conserved in average and there is no secular change in the semi-major axis.

Let us now apply the perturbation equations ( 38) to the specific case of the perturbation function corresponding to the quadrupole anomaly, namely

R = δu2 = 1 2 r 2 Q2 cos 2 θ - 1 3 . ( 40
)
Here cos θ = e • n, with n = x/r being the unit direction of the planet and r 2 = x 2 + y 2 + z 2 , where (x, y, z) are the coordinates of the planet in an absolute Galilean frame centred on the focus of the unperturbed Keplerian ellipse, and with respect to which the orbital elements {a, e, I, ℓ, ω, Ω} are defined. In the following, to simplify the presentation, we shall choose the x-direction of this Galilean frame to be the direction of the galactic centre e = ge/ge. That is, we assume that the origin of the longitude of the ascending node Ω lies in the direction of the galactic centre. This means that cos θ = x/r so that

R = Q2 6 2x 2 -y 2 -z 2 . ( 41
)
We express the planet's absolute coordinates (x, y, z) in terms of the orbital elements {a, e, I, ℓ, ω, Ω} by performing as usual three successive frame rotations with angles Ω, I and ω, to arrive at the frame (u, v, w) associated with the motion, where (u, v) is in the orbital plane, with u in the direction of the perihelion and v oriented in the sense of motion at perihelion. The unperturbed coordinates of the planet in this frame are

u = a (cos U -e) , (42a)
v = a 1 -e 2 sin U , (42b)
w = 0 , (42c)
where U denotes the eccentric anomaly, related to ℓ by the

These equations are general and give in particular the precession of the node dΩ/dt and the perihelion precession dω/dt or, rather, dω/dt where ω = ω + Ω is the longitude of the perihelion.

In order to make some estimate of the magnitude of the effect, let us approximate the direction of the galactic centre (which is only 5.5 degrees off the plane of the ecliptic)8 as being located in the plane of the orbit; consequently we choose I = 0. In this case ω = ω + Ω is the relevant angle for the argument of the perihelion. The non-zero evolution equations then become

de dt = 5Q2e √ 1 -e 2 4n sin(2ω) , (44a)
dℓ dt = n - Q2 12n 7 + 3e 2 + 15(1 + e 2 ) cos(2ω) , (44b) dω dt = Q2 √ 1 -e 2 4n 1 + 5 cos(2ω) . (44c)
We recall that with our notation ω is the azimuthal angle between the direction of the perihelion and that of the galactic centre (approximated to lie in the orbital plane). Of particular interest is the secular precession of the perihelion dω/dt due to the quadrupole effect which we henceforth denote by9

∆2 = Q2 √ 1 -e 2 4n 1 + 5 cos(2ω) . ( 45
)
The precession is non-spherical, in the sense that it depends on the orientation of the orbit relative to the galactic centre through its dependence upon the perihelion's longitude ω.

The effect scales with the inverse of the orbital frequency n = 2π/P and therefore becomes more important for outer planets like Saturn than for inner planets like Mercury. This is in agreement with the fact that the quadrupole effect we are considering increases with the distance to the Sun (but of course will fall down when r becomes appreciably comparable to r0, see Fig. 4).

We have also computed the secular planetary precession induced by the octupole moment Q3, for which the perturbation function reads with the same conventions

R = - 1 6 r 3 Q3 cos 3 θ - 3 5 cos θ = - Q3 30 x 2x 2 -3y 2 -3z 2 .
(46) Redoing the perturbation analysis we find the octupolar precession in the case I = 0 as

∆3 = Q3a √ 1 -e 2 32en
2 -13e 2 + 35e 2 cos(2ω) cos ω . (47)

Numerical evaluation of the planetary precession

We now compute numerically this effect for the various planets of the Solar System; the relevant orbital elements for the planets we used in this calculation are provided in Table 2.

Our numerical values for the quadrupole and octupole anomalies ∆2 and ∆3 are reported in Tables 3 and4 respectively. As we see from Table 3 the quadrupolar precession ∆2 is in the range of the milli-arc-second per century which is not negligible. In particular it becomes interestingly large for the outer gaseous planets of the Solar System, essentially Saturn, Uranus and Neptune. The dependence on the choice of the MOND function µ is noticeable only for functions µn(y) defined by ( 35) with large values of n, where the effect decreases by a factor ∼ 10 between n = 2 and n = 20. However, the octupolar precession ∆3 in Table 4 is much smaller, being rather in the range of the micro-arc-second per century.

We give for comparison in Table 5 the best published postfit residuals for any possible supplementary precession of planetary orbits (after the relativistic precession has been duly taken into account), which have been obtained from global fits of the Solar System dynamics by Pitjeva (2005) and [START_REF] Fienga | Relativity in Fundamental Astronomy: Dynamics, Reference Frames and Data Analysis[END_REF]. In particular the INPOP planetary ephemerides by [START_REF] Fienga | Relativity in Fundamental Astronomy: Dynamics, Reference Frames and Data Analysis[END_REF] use information from the combination of very accurate tracking data of spacecrafts orbiting different planets.

From Table 5 we see that our results for ∆2 are smaller Table 2. Orbital elements of planets used in the computation of the abnormal precession rates ∆ 2 and ∆ 3 . The longitude of the perihelion ω is defined here with respect to the direction of the external galactic field (assuming that the galactic centre lies within the plane of the ecliptic). Thus we pose ω = ωAllen -λ, where λ = -93.2 • is the longitude of the galactic centre and ωAllen is given in [START_REF] Allen | Astrophysical Quantities[END_REF]. 2009) publish an offset with respect to the general relativistic value for the precession: namely they quote -10 ± 8 mas/cy as obtained from the Cassini tracking data and 200 ± 160 mas/cy from the VEX data. The values we find for ∆2 are smaller but grossly within the range of these postfit residuals. However we find that ∆2 is positive for Saturn while the offset reported in [START_REF] Fienga | Relativity in Fundamental Astronomy: Dynamics, Reference Frames and Data Analysis[END_REF] from Cassini data is negative. But note that this may not be relevant because the INPOP ephemerides are used by [START_REF] Fienga | Relativity in Fundamental Astronomy: Dynamics, Reference Frames and Data Analysis[END_REF] to detect the presence of an eventual abnormal precession, not to adjust the value of that precession.

Mercury

In the case of the Earth the INPOP ephemerides find a tight constraint 0 ± 0.016 coming from the Jupiter VLBI data [START_REF] Fienga | Relativity in Fundamental Astronomy: Dynamics, Reference Frames and Data Analysis[END_REF]. This constraint seems already to exclude most of our obtained values for ∆2, except for MOND functions of type µn (35) with rather large values of n. In particular, one needs to take n > 8 in order to have an EFE precession compatible with this constraint on the Earth orbit. On the other hand we note that for all the planets the octupolar precession rate ∆3 is very small and clearly completely negligible in current fits of the Solar System dynamics.

Thus it seems that in the case of the Earth the constraint from the Jupiter VLBI data is already quite tight, and it is difficult to accommodate our anomalous quadrupolar precession rate ∆2 for most choices of MOND function µ(y). However let us remark that the postfit residuals in Table 5 are obtained by adding by hands an excess of precession for the planets and looking for the tolerance of the data on this excess. But in order to really test the anomalous quadrupolar precession rate ∆2, one should consistently work in a MOND picture, i.e. consider also the other effects predicted by this theory, like the precession of the nodes, the variation of the eccentricity and the inclination, and so on. Then one should perform a global fit of all these effects to the data; it is likely that in this way the quantitative conclusions would be different.

On the other hand, as we have commented in the Introduction, it is non-trivial to know if testing MOND-like theories (including TeVeS and the model proposed by Blanchet & Le Tiec 2008, 2009) in the Solar System could invalidate those theories, if they represent only some phenomenological models describing the weak field regime of gravity and as such should not be extrapolated as they are in the strong field of the Solar System.

CONCLUSION

In MONDian gravity, the influence of the external galactic gravitational field onto the orbits of Solar-system planets appears mainly through the presence of a quadrupolar perturbation to the gravitational potential, as it has been shown through a multipolar analysis of the solution of the MOND equation. This contribution has been computed using a spectral, highly accurate, numerical Poisson solver to obtain a solution of the equation (1). This numerical solution has been extensively tested and compared to known analytical formulas, and all the expected properties have been recovered. Using the most commonly used MOND coupling functions µ(y) and value of the MOND acceleration a0, we have obtained quadrupole contributions that are compatible with those published by Milgrom (2009). We have also extended this calculation to include octupole contributions. With the derivation of detailed perturbation equations, the influence of the quadrupole and octupole on the perihelion precession of planets have been quantified and compared to two sets of observations, in particular the INPOP planetary ephemerides by [START_REF] Fienga | Relativity in Fundamental Astronomy: Dynamics, Reference Frames and Data Analysis[END_REF]. At first sight, these observational results seem to exclude the presence of a MOND-EFE induced quadrupole, for some choices of the MOND function µ and scale a0. This is particularly true for the Earth orbit, where constraints from the Jupiter VLBI observations are very strong. However, for other choices of the MOND function the observational constraints on the precession of the Earth orbit are still compatible with the computed effect. These compatible MOND functions share the property of having a very sharp transition from the weak-field (modified) to the strong field (Newtonian) regime. Moreover one should be very cautious and keep in mind that:

1. The observational constraints on the Solar-system orbits have been obtained by a global fit, using a fully relativistic first-post-Newtonian model. Taking only into account a particular MOND effect like precession without considering a fully-MONDian model is in principle incoherent since all the other effects should also be considered;

2. The MOND theory is a phenomenological approach trying to describe gravity in the weak field regime (where evidence for dark matter arises). It is therefore unclear if such theory can be applied to stronger gravitational fields and, in particular, in the Solar System. Conversely, the constraints obtained here in strong field may not be relevant for the weak-field case unless some more fundamental theory, valid in both regimes, is known.

In any case, further studies are to be done if one wants to obtain more stringent conclusions about constraints imposed by Solar-system observations onto MOND-like theories. More precise observations could also give valuable informations about an eventual EFE due to MOND theory and also restrict the number of possible MOND functions that are compatible with the observations. of r, so we write accordingly uN(r) = f (r 2 ). Acting on a function of r 2 the STF derivative operator gives ∂Lf (r 2 ) = nL (2r) l f (l) (r 2 ) where f (l) denotes the l-th derivative of this function (see Blanchet & Damour 1986). Hence we find that ∂LuN is a perfectly regular function and we conclude that ( ∂LuN)(0) = 0 except for the monopole case l = 0. We can therefore replace ( ∂Lδu)(0) in (B4) by ( ∂Lu)(0) for any l 1 and we have finally established the result δu = -uN(0) + +∞ l=0 1 l!

x L ( ∂Lu)(0) . (B7) which is the same as in Eq. ( 19).

Figure 1 .

 1 Figure1. Setting of the spherically symmetric star and the asymptotic galactic field ge, using spherical coordinates {r, θ, ϕ}.

Figure 2 .

 2 Figure 2. Relative error of the code on the two tests described in Sec. 3.2, for a 1M ⊙ spherical mass distribution, for the standard MOND function µ 2 (y) [Eq. (34b)] and the standard value a 0 = 1.2 × 10 -10 m.s -2 . The abscissa gives the ratio (5) between the norm of the asymptotic gravitational field ge and the parameter a 0 . Note that both curves do coincide.

Figure 3 .

 3 Figure3. Two density profiles for testing the result dependence on the equation of state. Both profiles yield a one-solar mass star in hydrostatic equilibrium. The two profiles yield very little (same level as the overall numerical accuracy ∼ 10 -3 ) difference on our results.

Figure 4 .

 4 Figure4. Left panel: profile of Q 2 (r) in the solar system, for a standard choice of function µ 1 (y) [see Eq. (34a)], a 0 = 1.2 × 10 -10 m.s -2 and ge = 1.9 × 10 -10 m.s -2 . The MOND transition radius r 0 = GM/a 0 is shown by a dash-dotted line at r ≃ 7100 AU. Right panel: zoom of the central region (r 50 AU), where the quadrupole is almost constant (the spikes for r → 0 are due to the high value of the monopole term, from which it is numerically difficult to extract the quadrupole). Note the difference in the y-axis scales.

Figure 6 .Figure 7 .

 67 Figure 6. Quadrupole in the solar system Q 2 ≡ Q 2 (0) as a function of the external galactic gravitational field ge, for four different expressions of the coupling function µ(y) [see Eqs. (34) for a description]. The standard value of ge = 1.9 × 10 -10 m.s -2 is shown by a thin dash-double-dotted line.

Figure 8 .

 8 Figure8. Quadrupole moment Q 2 as a function of the stellar mass M . The M -1/2 -law is recovered as expected from Eq. (37a).

  Kepler equation ℓ = U -e sin U . Inserting the resulting expression of R into the perturbation equations (38), and performing the time average (in practice, an average over the eccentric anomaly U , after the appropriate change of vari-

	able), then yields
	da dt de dt	= 0 , = 5Q2e	√ 4n 1 -e 2	[cos I cos(2ω) sin(2Ω)	(43a)
	dI dt dℓ dt	(43b) 2 + 3e 2 + 5e 2 cos(2ω) cos Ω + sin(2ω)(cos 2 Ω -cos 2 I sin 2 Ω) , Q2 sin I sin Ω 4n √ 1 -e 2 (43c) -5e 2 cos I sin(2ω) sin Ω , = n + = Q2 96n -15(1 + e 2 ) cos(2ω) [2 -2 cos(2I)
			+ cos(2I -2Ω) + 6 cos(2Ω) + cos(2I + 2Ω)] +12 -8 + 3e 2 cos 2 Ω -12 1 -6e 2 + (7 + 3e 2 ) cos(2I) sin 2 Ω +20 2 -3e 2 +6(1 + e 2 ) cos I sin(2ω) sin(2Ω) , (43d)
	dω dt	= -	Q2 √ 4n 1 -e 2 +(-1 + e 2 ) (3 + 5 cos(2ω)) cos 2 Ω 2 -2e 2
				+	5 2	e 2 cos I sin(2ω) sin(2Ω)
	dΩ dt	=	-10 cos 2 I sin 2 ω sin 2 Ω +5(1 -e 2 ) cos I sin(2ω) sin(2Ω) , Q2 sin Ω 4n √ 5e 2 cos Ω sin(2ω) 1 -e 2 + cos I -2 -3e 2 + 5e 2 cos(2ω) sin Ω . (43e)

c 2010 RAS, MNRAS 000, 1-15

Table 3 .

 3 Results for the precession rates of planets ∆ 2 due to the quadrupole coefficient Q 2 . We use the values for Q 2 corresponding to various MOND functions defined in Eqs. (34)-(35) (see Table1). All results are given in milli-arc-seconds per century (mas/cy).

			Venus	Earth	Mars	Jupiter	Saturn	Uranus	Neptune
	a [AU]	0.397	0.723	1.00	1.523	5.203	9.537	19.229	30.069
	e	0.206	0.007	0.017	0.093	0.048	0.054	0.047	0.009
	P [yr]	0.24	0.62	1.00	1.88	11.86	29.46	84.01	164.8
	ω [deg]	170.7	224.7	196.1	69.2	108.0	185.6	264.2	138.2
			Quadrupolar precession rate ∆ 2 in mas/cy
	MOND function	Mercury	Venus	Earth	Mars	Jupiter	Saturn	Uranus	Neptune
	µ 1 (y) µ 2 (y) µ 5 (y) µ 20 (y) µexp(y) µ TeVeS (y)	0.04 0.02 7 × 10 -3 2 × 10 -3 0.03 0.05	0.02 0.01 3 × 10 -3 10 -3 0.02 0.02	0.16 0.09 0.03 9 × 10 -3 0.13 0.17	-0.16 -0.09 -0.03 -9 × 10 -3 -0.13 -0.17	-1.12 -0.65 -0.22 -0.06 -0.88 -1.21	5.39 3.12 1.05 0.3 4.25 5.81	-10.14 -5.87 -1.97 -0.56 -8.01 -10.94	7.93 4.59 1.54 0.44 6.26 8.56
	or much smaller than the published residuals except for
	the planets Earth, Mars and Saturn. The case of Saturn
	is interesting because the INPOP planetary ephemerides by
	Fienga et al. (

Table 4 .

 4 Results for the precession rates of planets ∆ 3 due to the octupole coefficient Q 3 . Beware of the results given here in microarc-seconds per century (µas/cy).Octupolar precession rate ∆ 3 in µas/cy

	MOND function	Mercury	Venus	Earth	Mars	Jupiter	Saturn	Uranus	Neptune
	µ 1 (y) µ 2 (y) µ 5 (y) µ 20 (y) µexp(y) µ TeVeS (y)	2 × 10 -3 2 × 10 -3 10 -4 5 × 10 -5 2 × 10 -3 2 × 10 -3	0.2 0.1 5 × 10 -3 2 × 10 -3 0.2 0.2	0.2 0.2 0.01 5 × 10 -3 0.2 0.2	-0.03 -0.03 -2 × 10 -3 -7 × 10 -4 -0.03 -0.03	1.4 1.1 0.06 0.03 1.4 1.3	19.5 15.1 0.8 0.4 19.5 17.9	12.1 9.4 0.5 0.2 12.1 11.1	1499 1162 61.2 28.7 1499 1374

Table 5 .

 5 We reproduce some published residuals for any supplementary orbital planetary precession. All results are given in milli-arcseconds per century (mas/cy).

			Postfit residuals for ∆ = dω/dt in mas/cy
	Origin	Mercury	Venus	Earth	Mars	Jupiter	Saturn	Uranus	Neptune
	Pitjeva (2005) Fienga et al. (2009)	-3.6 ± 5 -10 ± 30	-0.4 ± 0.5 -4 ± 6	-0.2 ± 0.4 0 ± 0.016	0.1 ± 0.5 0 ± 0.2	-142 ± 156	-6 ± 2 -10 ± 8	-0 ± 2 × 10 4	-0 ± 2 × 10 4

The code is based on the library lorene, available from the website http://www.lorene.obspm.fr. c 2010 RAS, MNRAS 000, 1-15

c 2010 RAS, MNRAS 000, 1-15

Here, the expansion is defined for any range in radius r, contrary to Sec. 2 and Appendix B, where the multipoles were only numbers defined for r → 0 and not functions of r.

 6 The radial dependence of the contribution of the l-th multipole moment is the same as that of the usual MOND spherical effect corresponding to a MOND function behaving like µ ∼ 1 -ǫy -γ when y → ∞, with γ = (l + 1)/2 (see the discussion at the end of Sec. 1.1). But of course the effect we are considering here is non-spherical. c 2010 RAS, MNRAS 000, 1-15

The latitude β and longitude λ of the galactic centre in the standard ecliptic coordinate system are β = -5.5 • and λ = -93.2 • (see e.g.[START_REF] Allen | Astrophysical Quantities[END_REF]).

The result found byMilgrom (2009) for this effect [his Eqs. (37)-(38)] looks more complicated, but can be simplified and is seen to be equivalent to ours. c 2010 RAS, MNRAS 000, 1-15

ACKNOWLEDGMENTS

It is a pleasure to thank Christophe Le Poncin-Lafitte, Alexandre Le Tiec and Mordehai Milgrom for a careful reading of the manuscript and very useful comments.

APPENDIX A: ASYMPTOTIC BEHAVIOUR OF THE FIELD

In this Appendix we prove that the asymptotic behaviour of the solution of the Bekenstein-Milgrom equation (1) at large distances from an isolated matter source in the presence of an external field ge, is given by Eq. ( 32). The proof has already been given by [START_REF] Bekenstein | [END_REF]; here we present a slightly different derivation. For simplicity we suppose that the source is spherically symmetric, so that the asymptotic form of the MOND potential is axisymmetric around the direction of the external field and of the type

where f (θ) is a function of the azimuthal angle θ measured from the direction e of the external field. We naturally assume that the function f (θ) is regular. With the ansatz (A1) we obtain successively the leading order terms in the expansion of various quantities as

where we denote f θ ≡ ∂f /∂θ, µe ≡ µ(ye), λe ≡ yeµ ′ e /µe (with ye = ge/a0 and µ ′ e = dµ(ye)/dye). Next, we impose that ∂i(µg i ) = 0 holds asymptotically far from a localised matter distribution. This, together with the fact that f is a regular function of the azimuthal angle θ, yields the firstorder differential equation

whose solution is given in terms of a constant K by

Finally we determine this constant by using the Gauss theorem (31), i.e.

where we choose a coordinate sphere since the result does not depend on the chosen surface at infinity. Inserting into (A5) the asymptotic expansion of the field given by (A2c), together with the solution (A4) found for f , we obtain

which establishes the looked-for formula (32).

APPENDIX B: EXPRESSION OF THE MULTIPOLE EXPANSION

The result (19) can also be proved directly, without resorting to the computation of surface integrals of Sec. 2. We first expand δu = u -uN when r → 0 using Taylor's formula as

Here the derivative operator ∂L is a priori non STF but we can decompose it into a sum of STF components using the formula (see the Appendix A of Blanchet & Damour 1986 for a compendium of useful formulas)

where k ranges from 0 to the integer part of l/2 denoted [l/2], where δ2K denotes a product of k Kronecker symbols carrying 2k indices chosen among those of L = i1 • • • i l , where the parenthesis indicate full symmetrisation over the l indices i1 • • • i l , where ∆ ≡ ∂ii denotes the ordinary Laplacian, and where the coefficients depending on l and k are given by

Using now the fact that ∆δu = 0 in a neighbourhood of the origin (r ε) we see that only the term k = 0 will survive in the latter sum, and since a l k=0 = 1 we conclude that we can replace in fact the derivative operator in (∂Lδu)(0) by its STF projection ∂L, hence

Next, we prove that the contribution from the Newtonian potential uN in (B4) is purely monopolar. This comes from our assumption that the density of ordinary matter ρ in the Sun is both regular (i.e. C ∞ ) and spherically symmetric, i.e. depending only on r (neglecting the back-reaction of the field on its source). As we shall now prove, this means that the Taylor expansion of ρ(r) in spherical symmetry must contain only even powers of r. The density being regular its expansion when r → 0 is given by Taylor's formula as

Next we again use the formula (B2) to decompose (B5) into a sum of STF pieces. Because ρ is spherically symmetric all these STF pieces must be zero except the monopolar contribution; hence only the terms for which l = 2k survive in the latter decomposition. Using the value of the coefficient a 2k k = 1/(2k + 1) we are then left with ρ = +∞ k=0 r 2k (2k + 1)! (∆ k ρ)(0) .

(B6)

As we see ρ(r) admits an expansion containing only even powers of r. Now the spherically-symmetric Newtonian potential uN will also admit an expansion only in even powers


C Aschan-Leygonie

& S Baudet-Michel

Building comparable synthetic health-related indicators of air quality in cities

Keywords: air pollution, spatial variation, indicator, cities, environmental epidemiology

Working on the relationship between respiratory health and atmospheric pollution in French cities led us to adopt a perspective that considers the risk for a single city of showing high levels of respiratory health problems. The risk level in a city is dependent on population vulnerability, the city context and pollutant levels. This presentation focuses on the observation of the unequal pollution hazards in French cities on a macro-scale. The study field covers 34 French cities of more than 100,000 inhabitants. The aim is to build comparable synthetic indicators of atmospheric pollution that consider the temporal and spatial diversity within a city, and to compare the different results. The focus is on the results for nitrogen dioxide (NO 2 ) over the 2007 and 2008 winter season. The paper presents the database and methodology used. The findings show that the spatial structure of air pollution between the cities remains globally the same whatever the indicator chosen. However, the results also show significant changes in the relative position of certain cities. When comparing the more conventional indicator calculated from mean daily concentrations in urban monitoring sites with an indicator also taking into account the concentrations observed in roadside stations, relative changes appear for more than ten cites. These findings highlight the importance of the criteria used in the choice of the indicators, and the benefit of using a set of complementary indicators in epidemiological studies.

Introduction

This paper is part of a research program that aims at investigating the relationships between respiratory health and urban environment, at both local and regional level, for a large set of 56 French urban areas, from 100,000 inhabitants to 10 million. Following the recommendations of Huang and Batterman (1), who highlight the fact that exposure assessment practices in epidemiological studies are highly dependent on the geographical exposure indicators, the aim is to compare different indicators. The comparison is based on the results obtained from different synthetic indicators regarding their levels of air pollution. Two questions are addressed: How can one build comparable pollution indicators, at the global level of a city, that take into account the strong spatial variety in the pollution level within a city? How can one build comparable pollution indicators that take into consideration the daily pollution cycle with alternating peaks and lows? The paper focuses on cities of more than 100,000 inhabitants, defined as built-up areas, over the 2008 winter season (December 2007to February 2008), as this is a high season for respiratory disease with no confounder such as pollens. The methodology is applied to nitrogen dioxide (NO 2 ) as this is reported to be one of the pollutants most closely related to respiratory disease and is the best-surveyed pollutant within urban areas in France.

There are two main types of epidemiological study that consider the link between air pollution and health in urban environments. The first type is based on an exploration of the relationship over a time period at city level, and although the study often includes different cities, the aim is not to compare the urban situations.

In the early 1990s, the APHEIS studies in Europe (2,3) and PSAS (4) studies in France focused on the short-term relationship between morbidity and mortality levels for cardiovascular and respiratory diseases and levels of pollutants such as sulphur dioxide, nitrogen dioxide, ozone and particles of more than 10 or 2.5 microns. The relationship between the pollutant time-series and health timeseries was calculated using Poisson regressions. The health risk was calculated in each city usually for a particle size elevation of 10 µg/m 3 : those risk levels rarely exceeded 1.5% (4).

The second type, which has been developed more recently (5), investigates health inequalities within cities and their link to air pollution. The local air pollution is estimated with a dispersion model with input data such as emission inventories, meteorological data, and concentrations of background pollution, and sometimes pollution registered by proximity monitors.

The present research aims to assess a city's global level of pollution, in order to compare it with other comparable French cities.

Using the BDQA database to build comparable urban pollution indicators

In France, air pollution monitoring has been obligatory since 1998 for cities with more than 100,000 inhabitants. ADEME (Agence De l'Environnement et de la Maîtrise de l'Energie) collects the pollutant concentration levels measured through a national network of independent agencies, the AASQA (Associations Agréées de la Surveillance de la Qualité de l'Air

Estimating ambient air pollution: a method integrating spatiotemporal variation

The construction of indicators allowing the comparison of NO 2 levels between cities should integrate both the spatial and temporal variation in pollutant concentrations in order to better approximate global exposure to pollution.

Dealing with daily cycles: making use of four parameters

As shown for Lyon (fig 2), all cities feature two peaks with high concentration of NO 2 and two lows. The time of day at which the peaks and lows occur varies slightly between cities, but the cycle observed in Lyon (peaks around 7-10 am and 5-8 pm and lows around 2 to 5 am and 1 to 3 pm) is fairly characteristic of French urban areas (6). Therefore, parameters such as mean and median are not accurate synthesisers of daily pollution and other parameters have to be used. We believe comparison of air pollution in cities might be improved if observed according to different parameters. A set of four different parameters were created in order to estimate the global NO 2 concentration level for each city: high and low percentiles and daily peaks and lows. Firstly, daily average values for each city were calculated from peak hours alone over winter 2008 for urban and roadside sites. Secondly, another set of daily mean values, this time from the lowest concentrations, was calculated for each of the 91 days of the study period, for urban and roadside sites. Thirdly, for each day of winter 2008, the 90 th percentile of the hourly NO 2 concentrations was determined. Finally, the 10 th percentile of the hourly average was determined per day over the three month period.

Addressing the spatial variability of NO 2 concentration

As shown for Lyon, nitrogen dioxide levels do not only differ from one hour to the next, but also naturally between different monitor types (figure 3). The averages observed in urban background sites are much lower in general and present a more homogenous profile than the traffic-related sites. The two urban monitors (Saint Just and Gerland) located within Lyon city centre present a mean value of the peak hours during winter 2008 of around 65 µg/m 3 , whereas the traffic-related sites register much higher means over the period (> 80 µg/m 3 ). The roadside sites show greater differences between the sites as they are directly affected by the intensity of the passing traffic (Garibaldi and the site near the A7). Surprisingly, the variability within the traffic-related sites over this period is generally lower for the traffic sites (relative standard deviation) of around 0.4 in urban sites and between 0.2 and 0.3 in roadside sites, except for the Grand Clément and Eastern ring road where variations are high).

A hierarchical ascendant classification (HAC) was used to evaluate the global differences in NO 2 levels measured in urban background sites and roadside sites.

The classification was carried out on the daily average of the eight-hour peak (7-10 am and 5-8 pm) during winter 2008. The dendrogram using Ward's criterion showed a perfect dichotomy between the NO 2 concentrations observed in background monitoring sites (LYCe, ST J to COTI) on the one hand, and in traffic-related sites (BERT to GARI) on the other hand. These findings and the fact that the spatial variability of NO 2 may mainly be allocated to two categories of land use (main roads and densely built-up areas) suggest that it is relevant to estimate the global air pollution of a city including the two monitor types.

Our work is based on the hypothesis that measurements in one type of monitoring site are representative of all those locations with similar features: urban background monitors are representative of all areas with a population density above 3,000 inhabitants per square kilometre, more than 100 meters away from a road with 4,000 to 7,000 vehicles per day and at least 200 meters away from a road with more than 7,000 vehicles per day. In the same way, ambient air concentrations measured in a roadside site are considered representative of all zones near to a road with more than 4,000 vehicles per day. Suburban background monitors and monitors located close to industrial plants have not been included in the study for two reasons: the location criteria are defined in such a way that these monitors survey specific pollutions (for example photochemical pollution in suburban sites) and few cities monitor NO 2 in suburban or industrial areas.

In order to estimate the share of the two types of area within a city we used GIS. Two categories of land use defined by CORINE Land Cover (250 m grid-based database) were used to estimate the share, within the city limits, of continuous urban built-up areas (combined with a population density above 3,000 inhabitants per square kilometre) and road networks and associated land. A 50meter buffer zone (7) from the roadside of the larger roads was applied, and for the urban canyon streets a smaller buffer zone was used (25 meters). Tested on Lyon, it resulted in a respective share of 70% for the continuous urban built-up areas and 30% of the roadside zones. Thus, to build up each global city indicator of NO 2 concentration levels, the levels in both types of monitoring station were weighted; a coefficient of 0.7 was applied for urban monitors and 0.3 for roadside monitors.

Results

Six parameters using median values were calculated in order to compare the cities' relative positions with regard to NO 2 concentrations. All were calculated for the 2008 winter. The parameters relying on low and peak hours and 10 th and 90 th percentile were calculated and used in order to produce a synthetic global indicator accounting for the daily variability in NO 2 concentration, based on the use of traffic and urban monitors. Two other synthetic parameters, similar to those use in epidemiological studies such as APHEIS, were calculated: one using daily means for urban monitors only, the other using daily means for traffic and urban monitors. The first observation when comparing the parameters (table 1) is the obvious fact that the mean levels differ from one parameter to another. The introduction of traffic monitors into the calculations leads to an increase in all levels observed for the median of the daily means: the minimum rises from 15.7 ug/m 3 to 23 ug/m 3 , the maximum from 55.2 to 63.5 ug/m 3 and the mean from 33.6 to 39.3 ug/m 3 . Other parameters calculated on traffic and urban monitors show the highest level of concentration for the 90 th percentile (61.9 ug/m 3 ) and the peak hours (51.3 ug/m 3 ) and the lowest level for the 10 th percentile (18.6 ug/m 3 ) and the low hours 29 ug/m 3 . Very high linear correlations (>0.84) are shown between all parameters (except with P10 R²=0.7).

In order to assess the pollutant levels within the cities a principal component analysis was conducted on the six parameters that render the daily variability of NO 2 concentration levels. The PCA gives rise to a very strong first factor explaining 92% of interurban variations in NO 2 levels. This very high variance shows that the results obtained by the different indicators are very similar. This result refutes our hypothesis that introducing spatial variability and not using mean values would radically change the results. Its the second axis that explains 8% of the inter-urban variations that shows in what way the choice of the indicators is important. This axis shows that the main difference is to be found between the indicators that are based on peak hours and those calculated from the daily lows. The factor scores place opposite each other cities with high peaks and low levels of NO 2 during the daily lows (at night and early afternoon) like Toulon, Nîmes and Clermont-Ferrand, to cities with low levels in peak hours and high concentrations during the daily lows (Strasbourg, Dijon, Le Havre, Rouen).

The first component can be used to characterise the cities regarding their global pollutant level (fig. 3). Lyon, Valence, Paris, Marseille, Grenoble, Nice and Montpellier stand out as the most polluted cities (high scores for all indicators), opposite St-Nazaire, Brest, Rennes, Dunkerque, Orléans, Tours, Bayonne for the latter (low scores for all indicators). Around the gravity centre eleven cities with medium concentrations, like Bordeaux, Toulouse, La Rochelle. A regionalisation of the urban system appears that produces three classes of cities: the southeastern ones, the north-western ones, the south-western ones and finally Paris, Lyon and Valence.

Figure 3 : Cities position on the first axis of the PCA (4 parameters to describe their pollution level : peak and low hours, P10 and P90)

A comparison between the spatial organisation obtained from the parameters using only urban monitors and using both traffic and urban monitors (weighted according to land use) showed significant changes of ranking for two groups of cities. The parameter used here to calculate the two synthetic indicators was in both cases the daily mean concentration, as is usual for epidemiological studies. The mean level for the sample of cities rises from 33.6 (urban monitors only) to 39.3 µg/m 3 (urban and traffic). The minimum and maximum values show similar elevations. The comparison between these two indicators shows that cities such as Valence, Paris, Nîmes, Avignon, Rouen, Montpellier and Nice not only show higher NO 2 concentrations, but their relative position in the sample of cities also changes, appearing more polluted in relation to the other cities than previously when using only urban monitors. Meanwhile Toulon, Amiens, Dunkerque, Dijon, Lille and Montbéliard appear to have NO 2 concentration that, although higher, rank them better in relation to the other cities. To us these results show that the use of urban, or urban and roadside monitors significantly changes the geography of polluted cities. It broadens differences between western (least polluted), south-eastern (most polluted) and south-western (averagely polluted)

cities. Therefore, we believe that the traffic monitors should not be systematically put aside in epidemiological analysis.

Discussion

The global indicators proposed in this paper have been built using a set of monitors characterising different location types in order to respond to the fact that the use of one sole background monitor may not reflect the situation of large cities and may lead to accepting the zero hypothesis concerning the relationship between pollution level and morbidity.

The introduction of traffic monitors into the analysis, the use of most city monitors, responds to the remark made by Katsouyanni et al. (8) that a reliable representation of human exposure to atmospheric pollution would involve at least three monitors (regardless of the type of site). This allows a more accurate assessment of the diversity of inner city pollution. As shown in 3.2., traffic monitors and urban monitors capture different pollution situations that both reflect the diversity of urban spaces. As assumed by a number of authors (2,3,4), this paper demonstrate the roadside and urban monitor's global indicator shows a rise (compared to urban monitors alone) in the NO 2 concentration level of all cities. However, this rise may be a better estimation of the global exposure situation of the city than the estimation from background monitors alone. Boudet, Zmirou (9) compared PM pollution exposure using mobile monitors and the AASQA stationary monitors. Mobile monitors were implemented by a set of individuals living and working in different parts of Grenoble city. The results showed that including proxy monitors in the analysis enabled an estimation of pollution levels closer to that of the mobile monitors than the observation of background monitors alone.

The interurban spatial structure remains globally the same whatever the indicator chosen. However, the results also show significant changes in the relative position of several cities. Comparing the more conventional indicator calculated from daily mean concentrations in urban monitoring sites with an indicator taking into account the concentrations registered by roadside monitors as well as urban ones introduces changes for almost a third of the cites. These findings highlight the importance of the criteria used in the choice of the indicators, and the benefit of using a set of complementary indicators, like the one proposed in the PCA. The use of daily means also appears to be less accurate than the use of for example percentiles.

Limits to our methodology are due to the comparative ambition that lead to the examination of a large set of cities and therefore do not allow the use of detailed information such as NO 2 cadastre and mobile measurement series or results from dispersion models.

Firstly, the hypothesis that areas with similar location features have a similar level of NO 2 pollution can be discussed. Even though the HAC showed a clear similarity between the pollutant levels observed in urban background locations on the one hand, and in traffic-related sites on the other hand, clearly there are numerous zones within a city for which the hypothesis of homogeneous pollution level derived from either traffic, either urban background sites is not defendable : as an example pollution level for traffic sites are different wether measured in large or small roads, in urban residential or urban working places... However, as it is impossible to make use of detailed information while looking at such a large set of cities, introducing variability relating to the identification of a monitor's location type seems to be a fair compromise. The weighting of NO 2 concentrations according to the type of monitors could be improved and adjusted for each city using Corine Land Cover database.

Secondly, and these are the main limits, there are great inequalities between cities regarding the quality of the surveillance of NO 2 (and other pollutants), as 50% host less than four monitor sites. Amongst the 56 largest cities in France, 19 do not survey traffic-related NO 2 . Although theoretically most monitor types should be found in every city, in reality, most cities only survey pollutants in two kinds of location. The AASQA have clearly favoured the background location, used in most health studies, while traffic monitors are too scarce. Eventhough pollution surveillance is now compulsory, the AASQA are responsible, in each city, for atmospheric pollution surveillance. The AASQA benefit from ADEME's guidelines and expertise in implementing their stationary network (10), but each ASQA decides for itself where to install permanent measurement stations, which pollutant to observe, when to expand their measurement network, and when to start observation of an emission campaign. They are independent agencies, but placed under the supervision of a management board where civil society, politicians and manufacturing companies operating in the city are represented. Therefore, strong interests are in conflict in AASQAs and the geography of each city's stationary network may result more in a compromise between political forces than in a scientific objectivity. Therefore, stationary networks can vary in density, seniority, variety of pollutant surveyed and monitoring sites. Thus, there are monitors that do not accurately reflect the air pollution of the type of area that they are supposed to represent.

Conclusion

The purpose of this work was to propose an alternative method for building indicators allowing registration of different urban exposure situations, which may be used in epidemiological and risk assessment studies. Suggested improvements to the conventional use of mean concentrations, are based on the integration of spatial variability by weighting the measurements according to the share of different types of land use and to avoid smoothing the daily variations of NO 2 . Depending on the types of monitoring sites used in the indicator and the way of dealing with the daily variations in the pollutant concentration, the relative position of the cities with respect to the global pollutant level varies.

These findings suggest that further investigations should be carried out on the relevance of certain indicators used in health studies.

Figure 2 :

 2 Figure 2: Daily variation of hourly median concentrations of NO 2 (Lyon, winter season 2007-2008) (data source: ADEME, 2009)

Table 1 :

 1 Parameters to describe NO 2 concentration in cities

	Variable name	Mean	Standard deviation	Minimum	Maximum
	Médiane_low urb/traf pond 6h	29.297	8.432	17.000	51.000
	Médiane_peak urb/traf pond 6h	51.365	11.103	31.800	73.400
	MédianesP10	18.568	7.428	8.100	37.600

Acknowledgements

This research has been supported by a contract from PIR Ville Environnement, an interdisciplinary research programme financed by CNRS (Centre National de la Recherche Scientifique) and the French ministry of Ecology and Sustainable Development (MEEDDM).

(1) Huang, Y. L. and S. Batterman (2000). "Selection and Evaluation of Air Pollution Exposure Indicators Based on Geographic Areas." The Science of the Total Environment(253): 127-144.

(2) APHEIS (2000)(2001) Atmpos. Environ 1998 ;32 :1921-1930. (8) Katsouyanni K., Zmirou D., Spix C., Sunyer J., Schouten J.P., Pönka A., Anderson H.R., Le Moullec Y., Wojtyniak B., Vigotti M.A., Bacharova L., Short-term effects of air pollution on health : a European approach using epidemiological time-series data. European Respiratory Journal (8), pp. 1030-1038, 1995. (9) Boudet, C., D. Zmirou, et al. (2001)