Galactic kinematics with modified Newtonian dynamics
We look for observational signatures that could discriminate between Newtonian and modified Newtonian (MOND) dynamics in the Milky Way, in view of the advent of large astrometric and spectroscopic surveys. Indeed, a typical signature of MOND is an apparent disk of “phantom” dark matter, which is uniquely correlated with the visible diskdensity distribution. Due to this phantom dark disk, Newtonian models with a spherical halo have different signatures from MOND models close to the Galactic plane. The models can thus be differentiated by measuring dynamically (within Newtonian dynamics) the disk surface density at the solar radius, the radial mass gradient within the disk, or the velocity ellipsoid tilt angle above the Galactic plane. Using the most realistic possible baryonic mass model for the Milky Way, we predict that, if MOND applies, the local surface density measured by a Newtonist will be approximately 78 within 1.1 kpc of the Galactic plane, the dynamically measured disk scalelength will be enhanced by a factor of 1.25 with respect to the visible disk scalelength, and the local vertical tilt of the velocity ellipsoid at 1 kpc above the plane will be approximately 6 degrees. None of these tests can be conclusive for the presentday accuracy of Milky Way data, but they will be of prime interest with the advent of large surveys such as GAIA.
keywords: gravitation – Stars: kinematics – Galaxy: fundamental parameters – Galaxy: kinematics and dynamics – Galaxy: structure
1 Introduction
Our understanding of Galactic stellar populations and kinematics achieves regular progress with the introduction of new large surveys of stars with photometry, distances, radial velocities or proper motions, enabling us to probe the 6dimensional space of stellar positions and velocities (e.g., Perryman et al. hip (); Hog et al. ÷hog00 (); Nordström et al. no04 (); Famaey et al. f05 (); Zwitter et al. zwi08 ()).
For instance, the shape of the Galactic 3dimensional potential has been probed by, e.g., the orbits of the Sagittarius stream (Ibata et al iba01 (); Helmi helmi (); Read & Moore RM05 (); Haghi et al. hag06 ()), the tidal tails of Palomar 5 (Grillmair & Dionatos GD06 (); Grillmair & Johnson GJ06 ()), or the kinematics of halo stars. In the solar neighbourhood, the potential can also be analysed by measuring the Galactic escape speed of high velocity stars (e.g., Smith et al. smi07 (); Famaey, Bruneton & Zhao fbz07 ()), the force perpendicular to the Galactic plane (e.g., Oort oor60 (); Crézé et al cre98 (); Siebert et al. sie03 (); Nipoti et al. nip07 (); Holmberg & Flynn hol00 (), hol04 ()), the coupling between the three components of the velocity distribution in the solar neighbourhood (e.g., Bienaymé bie99 ()), or the orientation of the velocity ellipsoid above the Galactic plane (Ollongren oll62 (); Hori & Liu hor63 (); LyndenBell lyn62 (); Siebert et al. sie08 ()). In the future, much progress is still expected to be made in our understanding of the Galactic potential with the advent of the JASMINE and GAIA missions.
This mapping of the Galactic potential is of great importance because it holds the key to our understanding of the behaviour of dark matter on galactic scales. Nowadays, the dominant paradigm is that dark matter is made of nonbaryonic weaklyinteracting massive particles, called “cold dark matter” (CDM). However, on galactic scales, the observations appear to be at variance with a sizeable list of CDM predictions (see, e.g., Famaey et al. fam07 ()). The observed conspiracy between the mass profiles of baryonic matter and dark matter at all radii in spiral galaxies (e.g., Famaey et al. fam07 ()) rather lends support to modified Newtonian dynamics (MOND, Milgrom mil83 ()), a paradigm postulating that for accelerations below the true gravitational attraction approaches , where is the usual Newtonian gravitational field. Without resorting to CDM, this paradigm is known to reproduce galaxy scaling relations, as well as the rotation curves of individual galaxies ranging over five orders of magnitude in mass (e.g., Sanders & McGaugh San02 ()). In particular, the kinematic analysis of tidal dwarf galaxies by Bournaud et al. (bou07 ()) is difficult to explain within the classical CDM framework, while it is in accordance with MOND (Gentile et al. Gen07 ()). On the other hand, some problems arise with this paradigm in the subgalactic and extragalactic scales (e.g., Zhao Z05 (); Angus et al. ang07 (); see e.g., Zhao Z08 (); Angus ang09 (), and Bruneton et al. br09 () for possible solutions).
The present study examines possible observational signatures of MOND gravity in the Milky Way, which could be detected (or rejected) with the advent of large Galactic surveys. This is in direct continuity with the works of Famaey & Binney (fam05 ()); Nipoti et al. (nip07 ()); Wu et al. (wu08 ()), and McGaugh (mcg08 ()). We build these predictions using the MOND Milky Way model of Wu et al. (wu08 ()), and refer the reader to this paper for a full description of the model^{1}^{1}1We use the model labelled ”MOND ” in Wu et al. (wu08 ()), meaning that the modulus of the external gravitational field acting on the Milky Way is chosen to be .. This model is based on one of the most realistic possible baryonic mass models of the Milky Way, the Besançon model (Robin et al. rob03 ()).
Once the gravitational potential of the model is known, one can apply the Newtonian Poisson equation to recover the density distribution that would have yielded this potential within Newtonian dynamics. In this context, MOND predicts a disk of “phantom” dark matter (Sect. 2), allowing us to differentiate the MOND prediction from these of a Newtonian model with a dark halo by (i) measuring dynamically (within Newtonian dynamics) the disk surface density at the solar radius (Sect. 3), (ii) measuring the radial dynamical mass gradient within the disk (Sect. 4), or (iii) measuring the velocityellipsoid tilt angle above the Galactic plane (Sect. 5). We show that none of these tests are yet conclusive with presentday accuracy, but they are extremely useful to future highprecision astrometric and spectroscopic surveys.
2 The “phantom” disk of MOND
Bekenstein & Milgrom (bek84 ()) proposed a modification of the Newtonian dynamics to reproduce the flat rotation curve of disk galaxies without the introduction of enormous amounts of dark matter on galactic scales. They modified the Poisson equation by relating the gravitational potential to the mass density according to
(1) 
In the case of high acceleration, the function and the equation becomes the usual Poisson equation. The case of low acceleration is explained and described in a long series of papers (see for instance Milgrom mil02 ()). At low accelerations (called the deep MOND regime), when , the systems differ most from Newtonian ones: assuming , the equation becomes:
(2) 
Spherical densitypotential pairs are easily built and some flattened spheroids can be numerically described (see for instance Ciotti et al. cio06 ()). In the Appendix, we develop a new numerical potential solver for this equation.
Within MOND, one finds that any disk or spheroidal density distribution of finite mass produces a spherical potential at large radius (see Appendix). It is the Mondian explanation of the dark mater halo detected by a Newtonist. However, a disky distribution also produces a supplementary disk potential (Milgrom mil83 (), mil01 (); see also Appendix).
Once the MOND gravitational potential is known, one can always apply the Newtonian Poisson equation to it, to recover the density distribution corresponding to the given potential in Newtonian dynamics. Thus, how can we differentiate^{2}^{2}2 MOND is compatible (i.e., not falsified) with recent determination of the Galactic halo flattening, for instance deduced from the kinematics of the Sagittarius stream (Ibata et al iba01 (); Read & Moore RM05 ()). Newtonian from Mondian dynamics, since the difference of density (given by the two theories from the same potential) can always be attributed by a Newtonist to a dark matter component (the socalled “phantom” dark matter of MOND)? The answer lies within the properties of the supplementary disk of “dark matter” that appears in the phantom density .
The properties of this phantom disk are very specific and unique to MOND. Of course, two different kinds of dark matter components (with distinct kinematics) could explain both the round and the thin dark components (see e.g., Read et al. Read08 ()). However, if future observations reveal such a thin darkmatter disk, of the exact charateristics predicted by MOND, they will clearly make MOND dynamics more relevant to explaining Galactic dynamics, while failing to detect it could falsify the Bekenstein & Milgrom (bek84 ()) formulation of MOND. Our goal hereafter is to quantify the predicted effects of this dark matter disk on Galactic kinematics in the Milky Way. This specific dark (phantom) disk coud be easily identified using GAIA data; however, the exact predictions of MOND depend on its ability to describe the visible matter distribution. For that purpose, models such as the Besançon model are necessary because it takes into account important details such as the change in both radial scale length, and the mass distribution between thin, thick disks and other baryonic components.
We note that, hereafter, when we speak about the “dynamical mass distribution in MOND”, we always speak about the mass distribution that would have yielded the MOND gravitational potential from the Newtonian Poisson equation, i.e., the baryonic mass plus the phantom dark mass.
3 The local dynamical surface density in MOND
Milgrom (mil83 ()) already established that the phantom dark disk would enhance the force perpendicular to the Galactic plane. To be compatible with the analysis of Hipparcos data in the solar neighbourhood, this enhancement must not be too strong and must not imply an extremely massive unseen disk. Based on a Galactic model, Nipoti et al. (nip07 ()) determined that at the solar position MOND almost doubles the force produced in Newtonian gravity by the baryonic disk. For this model (based on model 1 of Dehnen and Binney, deh98 ()), the baryonicdisk surfacemass density is =43 M pc. When the disk is embedded in a spherical darkmatter halo, the surface density becomes 65 M pc (see their Figure 5), while for MOND (baryons+phantom) it becomes 85 M pc. This corresponds respectively to an increase of 51 percent and 98 percent compared to the surface density of the baryonic disk.
Using the Milky Way MOND model of Wu et al. (wu08 ()), we found instead an increase of 57, 62, and 66 percent at , 8, and 8.5 kpc respectively (78 at kpc). This difference with Nipoti el al. (nip07 ()) results from different baryonicmass concentrations and different local values between the models. This dynamical surface density of 78 must be compared with the force determination at 1.1 kpc, M pc by Holmberg & Flynn (hol04 ()). Considering the various uncertainties, these measurements are compatible both with the Newton+spherical halo and MOND Besançon models. Indeed, these accurate determinations of the are based on clearly defined, but small, samples of a few hundreds of stars. Due to the small size of these samples, only one parameter can be efficiently recovered: the surface density below some height , the mean distance of the samples from the Galactic plane. This does not break the degeneracy between two effects: a small flatenning of the halo and a small increase of the baryonic disc density. These two effects can compensate without modifying the surface density measured at 1.1 kpc. In contrast, the shape at smaller (and also the Oort limit) will change. In the future, large samples from RAVE and very large ones from GAIA, will make it possible to recover precisely versus z, and to determine dynamically the detailed vertical mass density distribution .
4 The dynamical disk scalelength in MOND
A second and different test concerns the prediction of the force perpendicular to the Galactic plane at various Galactocentric radii, to measure dynamically the scalelength of the (visible+phantom) disk mass distribution and to compare this to the visible disk scalelength. The interesting property of this test is that only the gradients in the density distributions of both the visible and visible+phantom disks have to be measured. These gradients are certainly easier to measure than the true total mass distribution.
In Fig. 1, we plot (dashed line) the radial density distribution of the baryonic matter in the Galactic plane for the Besançon model. By adopting the MOND gravitational potential (Wu et al. wu08 ()), we then determine the baryonic+phantom dynamical density (solid line) that would be inferred by a Newtonist measuring the force perpendicular to the Galactic plane at various Galactocentric radii. A factor of 1.25 is predicted between the two scale lengths, determined between 5 kpc and 10 kpc from the Galactic centre.
Since the Galactic thin disc dominates the disk mass distribution with a radial scale length of 2.5 kpc (Robin et al. rob03 ()), the expected disk dynamical scalelength differs significantly, 3.1 kpc, between measurements by a Newtonist and by a MONDian. This measurement will be a test of prime interest when data from large surveys such as GAIA or JASMINE become available.
5 The tilt angle of the velocity ellipsoid in MOND
A third new and interesting test of MOND relies on the measure of the vertical tilt of the stellar velocity ellipsoid above the Galactic plane (e.g., Siebert et al. sie08 ()). At any position above the Galactic plane, solar neighbourhood stars are distributed along isodensity ellipsoids in the velocityspace (where is the velocity towards the Galactic centre, the velocity in the direction of Galactic rotation, and the vertical veloicty, all with respect to the Sun). At a given position above the Galactic plane, the inclination angle of the ellipsoid minor axis with respect to the axis (see Fig. 6 of Siebert et al. sie08 ()) is intimately linked to the gravitational potential. In the plane, this angle is defined by
(3) 
where , , and are the secondorder velocity distribution moments. Cuddeford & Amendt (cu91 ()) analysed, to the fourth order, the consecutive moments of the Boltzmann equation by developing these moments by expanding a Taylor series. Combining these moment equations, they derived expressions for the velocity moments and obtained an approximate expression for the tilt angle , which depend only on the shape of the gravitational potential (but is exact only for completely separable axisymmetric potentials). More precisely, the ellipsoid inclination close to the Galactic plane is related to the disk scalelength and also to the darkmatter halo flattening (Bienaymé bie09 ()). This dependence of the vertical tilt of the velocity ellipsoid on the shape of the potential has also been checked by numerical simulations (e.g., Siebert et al. sie08 ()).
Here, we determine the vertical tilt angle as a function of (at the Galactocentric radius of the Sun) for the Newtonian Besançon model embedded in a dark halo (Robin et al. rob03 ()), and compare this inclination with that obtained for the MOND model (Wu et al. wu08 ()). For this, we follow Siebert et al. (sie08 ()), and compute the inclination of the velocity ellipsoid by integrating orbits in the meridional plane (we note however that the Besançon model is nonaxisymmetric, so that we average it). We sample the initial conditions from a Shu distribution function (Bienaymé bie99 ()), yielding an orbit library of over 20 million orbits. Each orbit was integrated over 120 rotations, using a 4th order RungeKutta algorithm. We then randomly selected 1 phasespace point per orbit in the last 80 rotations. From there, we computed the inclination angle of the velocity ellipsoid at different positions. Figure 2 displays the resulting vertical tilt angle (calculated with Eq. 3) as a function of at three different Galactocentric locations for both the Newtonian (+dark matter) and MONDian Besançon models. This tilt is also compared with the tilt expectation from a purely spherical potential (dotted lines on Fig. 2): clearly, the flattening of the MOND potential due to the presence of a phantom dark disk makes the inclination angle lower than for a purely spherical potential, and lower than for the Newtonian Besançon model embedded in a spherical halo. The tilt differences are most significant in the inner part of the Galaxy. At the solar position, however, the differences remain small when 1kpc, and become significant only at kpc (a difference of, 3.5 and 2 degrees, respectively at =7.5 and 8.5 kpc, the inclination being smaller within the MOND model). At kpc and kpc (the solar position in the Besançon model), where the strongest observational constraints on the tilt have been obtained, both models predict an inclination of 6 degrees, which agrees with the observational determination of degrees, obtained with a RAVE sample of 580 stars (Siebert et al. sie08 ()). Other observational results on the measure of the vertical tilt remain a bit controversial and conflicting (Fuchs et al. fuchs09 (); Smith et al. smi09 ()), and will require more (and more precise) data to avoid any potential bias, and confirm or disprove the sphericity of the potential.
Besançon MOND  Observations  

78  74 6 (Holmberg & Flynn hol04 ())  
Tilt at kpc  6 degrees  7.3 1.8 degrees (Siebert et al. 2008) 

A measurement of the necessary accuracy will possibly be achieved at the solar Galactic radius with large RAVE samples in the near future, and with GAIA data in a few more years. It will be important to account for secondary effects (see, e.g., Famaey et al. fam05 ()), such as streams, resonances, and the bar influence, and to correctly take into account their impact on the vertical tilt of the ellipsoid. However, we note, that, in the immediate solar neighbourhood, only the Hyades are clearly identified in the distribution that looks extremely well phase mixed, in contrast to the distribution where many substructures (from resonances or nonstationnarity) have been identified. This is also true for the distributions of distant samples at 1 kpc from the Galactic plane (Soubiran et al. sou08 () and Siebert et al. sie08 ()). The effect of the bar and specific vertical resonances needs to be understood in greater detail since they can modify significantly the verticaltilt orientation (Ollé & Pfenniger oll98 ()). However, this family of orbits should also be visible as substructures in the plane.
6 Conclusions
In addition to the predictions made by Wu et al. (wu08 ()) about the rotation curve and the escape speed, MOND also makes very specific predictions about Galactic kinematics related to the presence of a disk of phantom dark matter. These predictions will be tested with future large surveys:
(i) By measuring the force: at the solar radius, MOND predicts a 60 percent enhancement in the dynamical surface density at 1.1 kpc relative to the baryonic surface density, a value not excluded by current data. This enhancement would become more apparent at large Galactic radii, where the stellardisk massdensity becomes negligible.
(ii) By determining dynamically the scalelength of the diskmass densitydistribution. This scale length is a factor of 1.25 larger than the scalelength of the visible stellar disk if MOND applies. This test could be applied with existing RAVE data (Steinmetz et al. ste06 (); Zwitter et al. zwi08 ()), but the accuracy of available proper motions still limits the possibility of exploring the gravitational forces too far from the solar neighbourhood.
(iii) By measuring the velocity ellipsoid tilt angle within the meridional Galactic plane. This tilt is different for the two dynamics in the inner part of the Galactic disk. However, the tilt of about 6 degrees at z=1 kpc at the solar radius agrees with the determination of degrees obtained by Siebert et al. (sie08 ()). The difference between MOND and a Newtonian model with a spherical halo becomes significant at z=2 kpc.
These tests of gravity could be applied with future GAIA or JASMINE data that will allow us to reconstruct the 3dimensional gravitational potential of the Galaxy. To assess the values of the current local constraints, the predictions of the Besançon MOND model are compared with the relevant observations in Table 1.
Acknowledgements.
The authors gratefully acknowledge Annie Robin for her help. BF acknowledges financial support of the FNRS.References
 ()
 (2009) Angus, G.W. 2009, MNRAS, 394, 527
 (2007) Angus, G.W., Shan H.Y., Zhao H.S., Famaey, B. 2007, ApJ, 654, L13
 (1984) Bekenstein, J., Milgrom, M. 1984, ApJ, 286, 7
 (1999) Bienaymé, O. 1999, A&A, 341, 86
 (2009) Bienaymé, O. 2009, A&A, in press
 (2007) Bournaud, F., Duc, P.A., Brinks, E. et al. 2007, Science, 316, 1166
 (1995) Brada, R., Milgrom, M. 1995, MNRAS, 276, 453
 (2009) Bruneton, J.P., Liberati, S., Sindoni, L., Famaey, B. 2009, JCAP, 3, 31
 (2006) Ciotti, L., Londrillo, P., Nipoti, C. 2006 ApJ, 640, 741
 (1998) Crézé, M., Chereul, E., Bienaymé, O., Pichon, C. 1998, A&A, 329, 920
 (1991) Cuddeford, P., Amendt, P. 1991, MNRAS, 253, 427
 (1998) Dehnen, W., Binney, J. 1998, MNRAS, 294, 429
 (2005) Famaey, B., Binney, J. 2005, MNRAS, 363, 603
 (2005) Famaey, B., Jorissen, A., Luri, X., et al. 2005, A&A, 430, 165
 (2007) Famaey, B., Bruneton, J.P., Zhao, H.S. 2007, MNRAS, 377, L79
 (2007) Famaey, B., Gentile, G., Bruneton, J.P., Zhao, H.S. 2007, Phys Rev D 75, 063002
 (2009) Fuchs, B., Dettbarn, C., Rix, H.W. et al. 2009, AJ, 137, 4149
 (2007) Gentile, G., Famaey, B., Combes, F., Kroupa, P., Zhao, H.S., Tiret, O. 2007, A&A, 472, L25
 (2006) Grillmair, C.J., Dionatos, O. 2006, ApJ, 641, L37
 (2006) Grillmair, C.J., Jonhson, R. 2006, ApJ, 639, L17
 (2006) Haghi, H., Rahvar, S., HasaniZonooz, A. 2006, ApJ, 652, 354
 (2008) Hecht, F., Pironneau, O., Le Hyaric, A., Ohtsuka, K. 2008 http://www.freefem.org/ff++
 (2004) Helmi, A. 2004, MNRAS, 351, 643
 (2000) Hog, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27
 (2000) Holmberg, J., Flynn, C. 2000, MNRAS,313, 209
 (2004) Holmberg, J., Flynn, C. 2004, MNRAS,352, 440
 (1963) Hori, L., Liu, T.P. 1963, PASJ, 15, 100
 (2001) Ibata, R., Lewis, G.F., Geraint F., Irwin, M., Totten, E., Quinn, T. 2001, ApJ, 551, 294
 (1962) LyndenBell, D. 1962, MNRAS, 124, 95
 (2008) McGaugh, S. 2008, ApJ, 683, 137
 (1983) Milgrom, M. 1983, ApJ, 270, 371
 (2001) Milgrom, M. 2001 ApJ, 326, 1261
 (2002) Milgrom, M. 2002 New Astronomy Review, 46, 741
 (2007) Nipoti, C., Londrillo, P., Zhao, H. S., Ciotti, L. 2007, MNRAS, 379, 597
 (2004) Nordström, B., Mayor, M., Andersen, J., et al., 2004, A&A, 418, 989
 (1998) Ollé, M, Pfenniger, D 1998, A&A, 334, 829
 (1962) Ollongren, A. 1962, Bull. Astron. Inst. Netherlands, Vol. 16, 241
 (1960) Oort, J.H. 1960, Bull. Astron. Inst. Netherlands, Vol. 15, p.45
 (1997) Perryman, M.A.C., et al. 1997, A&A, 323, L49
 (2002) Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,2002, Numerical Recipes, Cambridge: University Press
 (2005) Read, J.I., Moore, B. 2005, MNRAS, 361, 971
 (2008) Read, J.I., Lake, G., Agertz, O., Debattista, V. 2008, MNRAS, 389, 1041
 (2003) Robin, A. C., Reylé, C., Derrière, S., Picaud, S. 2003, A&A, 409, 523
 (2002) Sanders, R.H., McGaugh, S. 2002, ARA&A, 40, 263
 (2003) Siebert, A., Bienaymé, O., Soubiran, C. 2003, A&A, 399, 531
 (2008) Siebert, A., Bienaymé, O., Binney, J., BlandHawthorn, J. et al. 2008, MNRAS, 391, 793
 (2007) Smith, M.C., Ruchti, G., Helmi, A. et al. 2007, MNRAS, 379, 755
 (2009) Smith, M. C., Evans, N. W., An, J. 2009, arXiv:0902.2709
 (2008) Soubiran, C., Bienaymé, O, Mishenina, T.V., Kovtyukh, V.V. 2008 A&A, 480, 91
 (2006) Steinmetz, M., Zwitter, T., Siebert, A., Watson, F.G., Freeman, K. C., Munari, U. et al. 2006, AJ, 132, 1645
 (2007) Tiret, O., Combes, F. 2007 A&A, 464, 517
 (2008) Wu, X., Famaey, B., Gentile, G., Perets, H., Zhao, H. 2008, MNRAS, 386, 2199
 (2005) Zhao, H.S. 2005, A&A, 444, L25
 (2008) Zhao, H.S. 2008, J. Phys.: Conf. Ser., 140, 012002
 (2008) Zwitter, T., Siebert, A., Munari, U., Freeman, K. C.; Siviero, A. et al. 2008 AJ136, 421
Appendix A Full multigrid MOND potential solver
We developed a new full multigrid (FMG) solver in the same spirit as Brada & Milgrom (bra95 ()) and Tiret & Combes (tir07 ()). It solves Eq. (1) by means of relaxation performed on the density field sampled at different resolutions. Largescale feature convergence is obtained from a coarse sampling of the density, while highfrequency features result from fine grids (see e.g., Press et al. NR02 ()). Relaxation is acheived with a 4 colours GaussSeidel, with 5 pre/post iterations.The solver is materially accelerated by means of Graphical Processor Units (GPUs) using the CUDA API developed by NVidia for its hardware. The Poisson equation is solved in second over a grid using a GeForce 8800 GTX device. Typical relative residuals obtained for the MOND Poisson equation are .
We tested our potential solver with disks of finite thickness, namely an exponential disk (scale length =1, thickness ), and a thickened Kuzmin disk (scale length =1, and thickness =0.2):
(4)  
(5) 
The FMG results for the effective (visible+phantom) density are plotted in Figs. 3 and 4. We note that the phantom disk clearly appears in the outskirts of the Galaxy. The FMG results were previously checked by solving Eq. (2) using the FreeFem++ software, a partial differential equation solver (Hecht et al. hecht ()). For spherical potentials the accuracy of the computed potential with FreeFem++ is excellent. With very flattened systems (axis ratio 1/10), the achieved accuracy remains good thanks to an adaptive grid within FreeFem++. For instance, in the case of the exponential disk, our mean relative accuracy is 1.5 per cent at radii of between 0.4 and 6 scale lengths, and lower than 5 per cent for radii between 0.2 and 10 scale lengths. Excellent agreement is found between the FMG and FreeFem++, with relative differences of smaller than within 10 scale lengths.