On the rotation curve of disk galaxies in General Relativity

Recently, it has been suggested that the phenomenology of flat rotation curves observed at large radii in the equatorial plane of disk galaxies can be explained as a manifestation of General Relativity instead of the effect of Dark Matter halos. In this paper, by using the well known weak field, low velocity gravitomagnetic formulation of GR, the expected rotation curves in GR are rigorously obtained for purely baryonic disk models with realistic density profiles, and compared with the predictions of newtonian gravity for the same disks in absence of Dark Matter. As expected, the resulting rotation curves are indistinguishable, with GR corrections at all radii of the order of $v^2/c^2\approx 10^{-6}$. Next, the gravitomagnetic Jeans equations for two-integral stellar systems are derived, and then solved for the Miyamoto-Nagai disk model, showing that finite-thickness effects do not change the previous conclusions. Therefore, the observed phenomenology of galactic rotation curves at large radii requires Dark Matter in GR exactly as in newtonian gravity, unless the cases here explored are reconsidered in the full GR framework with substantially different results (with the surprising consequence that the weak field approximation of GR cannot be applied to the study of rotating systems in the weak field regime). In the paper, the mathematical framework is described in detail, so that the present study can be extended to other disk models, or to elliptical galaxies (where Dark Matter is also required in newtonian gravity, but their rotational support can be much less than in disk galaxies).


INTRODUCTION
Recently, following the original suggestion of Cooperstock and Tieu (2007, see also Balasin and Grumiller 2008), several papers addressed the interesting possibility that the observed phenomenology of rotation curves in disk galaxies can be explained by General Relativity effects (hereafter, GR) peculiar of rotating systems, without the need to invoke the presence of Dark Matter (hereafter DM) halos in order to produce the flat behavior at large galactocentric distances.
Unfortunately, no definite consensus about the importance of GR effects on the rotation curve of disk galaxies seems to be reached, with widely different conclusions ranging from support to the hypothesis, to the identification of possible mathematical issues affecting the disk models used to compute the GR solutions (for a representative, but almost certainly incomplete list of papers representing the different positions, see, e.g., Korzynski 2005, Vogt and Letelier 2005a,b, Cross 2006, Fuchs and Phleps 2006, Crosta et al. 2020, Carrick and Cooperstock 2012, Deledicque 2019, Ludwig 2021, 2022, Toth 2021, Astesiano and Ruggiero 2022. If confirmed, the suggestion above would be a most surprsing result, with consequences extending well beyond the problem of the interpretation of galaxy rotation curves, and perhaps even beyond the problem of the existence of DM. A list of some of these consequences (not necessarily in order of importance) is the following: 1) usually, it is expected that lowest order corrections of GR to newtonian dynamics are of the order of v 2 /c 2 , where v is a characteristic velocity associated with the newtonian gravitational potential of the system, and c is the speed of light. In astronomical systems where the presence of DM is required by newtonian gravity, v/c ≈ 10 −3 or less, therefore such systems are empirically in the GR weak field regime. If the flat region of rotation curves is a GR effect, then we face the problem of explaining in physical terms how an expected effect of the order of ≈ 10 −6 actually becomes more important than the zeroth-order newtonian term: we notice that mathematically such a behavior characterizes singular perturbation theory. 2) In case the previous point is satisfactorily answered, we should then explain why other effects of GR in the weak field approximation (not directly related to rotation) remain at the level of small perturbations, for example by requiring the presence of DM in order to reproduce gravitational lensing, and correctly predicting the small amount of planetary precession left unexplaied once newtonian precession is considered 1 . 3) DM is clearly required by newtonian gravity not only in rotating disk galaxies, but also in velocity dispersion (pressure) supported astronomical systems with low ordered rotation, such as elliptical galaxies and cluster of galaxies (with converging predictions about the structural properties of the inferred DM halos obtained from different diagnostics, such as dynamical analysis, the modeling of X-ray emitting halos, gravitational lensing, see e.g. BT08, Bertin 2014).
Different strategies can be imagined to test the possibility of a GR origin of the flat rotation curve in disk galaxies. In the first, one just try to reproduce the expected GR rotation curve of some specific galaxy starting from its observed baryonic (e.g., stars and gas) density profile. This attempt suffers for some shortcoming. In fact, the mathematical modeling in GR (also considering effects such as sparse data points, error bars, and so on) is more complicated than in newtonian gravity, and great care is needed to interprete the results. Moreover, it should be recalled that the newtonian rotation curve of a purely baryonic razor-thin exponential disk (the common stellar density profile oberved in disk galaxies) of total mass M d and scale-lenght R d , is not Keplerian over a large fraction of the stellar disk (see e.g. Binney andTremaine 2008, Ciotti 2021, hereafter BT08 andC21, respectively), increasing from the center, reaching a quite shallow maximum at ≈ 2R d , and remaining almost flat up to ≈ 3R d (a radius already encircling ≈ 0.8M d , see Section 4.1). It follows that in newtonian gravity observed stellar rotation curves in disk galaxies hardly require any DM over a large fraction of the optical disk, while DM is certainly required by the rotation curves at larger galactocentric distances measured by radio observations in HI gas (e.g., Kalnajs 1983, Kent 1986, van Albada et al. 1985, van Albada and Sancisi 1986, and in particular Chapter 20 in Bertin 2014 for a complete account of the situation), and theoretically by stability arguments (Ostriker and Peebles 1973). Therefore, in this first approach the prediction of a GR rotation curve not declining with R over a large part of the stellar disk would just be what expected in case of small GR corrections to the newtonian rotation curve.
In a second approach (that we follow in this paper), one avoid direct comparison with observational data, but instead construct the newtonian and the (weak field) GR rotation curves produced in the equatorial plane by a rotating stellar disk without DM, with total barionic mass, scale lenght, and density profile similar to those of observed in real disks/elliptical galaxies. The obvious advantage of this approach is that, whatever the solution is, we learn something. In fact, let us assume that the obtained rotation curve in GR and in newtonian gravity are essentially the same (the common expectation): after excluding the case of mathematical/physical errors in the modeling, only three conclusions are possible. Conclusion 1: if we still pretend that GR can explain the flat profile of the rotation curves of disk galaxies at large radii without invoking the presence of DM halos, then we must conclude that the adopted weak-field approximation of GR cannot be used to describe the dynamics of disk galaxies, even though these systems are (empirically) in the weak-field regime. Conclusion 2: the weak field approximation of GR can be used to describe the weak field regime of disk galaxies, but in the explored cases we fail to reproduce the flat region of rotation curvers at large radii (i.e., well beyond the geometrical outskirts of the optical stellar disk) because we assumed a too simple/idealized orbital structure for the stars producing the velocity-dependent component of the GR force in the equatorial plane. Conclusion 3: the weak-field approximation can be used to describe the weak-field regime in disk galaxies, and DM halos are required in GR as they are in newtonian gravity.
Conclusion 1 would be quite formidable, and should be proved convincingly showing that higher order effects not considered in the weak field expansion used are even more important than those considered or, better, by presenting a case obtained by numerically solving the full (non-linear) GR equations for a realistic barionic disk in the weak field regime, with a predicted rotation curve significantly different with respect to the newtonian rotation curve for the same barionic disk, i.e. with GR "corrections" well above 300% and increasing as the square root of the galactocentric radius at larger and larger distances from the center. Conclusion 2 would imply that an universal property of rotation curves of disk galaxies is a GR effect produced by peculiar (actually, never observed) significant streaming motions of the stellar populations along the vertical and radial directions (see Section 4 for details). As we will see, even if with the aid of simple models (however, the most realistic used so far in this problem), in this paper we present quite strong evidence that a commonly adopted weak field approximation of GR predicts rotation curves indistinguishable from the newtonian ones in absence of DM, therefore strongly supporting Conclusion 3 above. In order to illustrate the results in the most transparent way, the mathematical setting is rigorously presented, and all the assumptions made explicitely stated, so that interested reader is in position to repeat and extend the study. We use cylindrical coordinates (R, ϕ, z), and vectors are in bold-face.
The paper is organized as follows. In Section 2 we recall the expression for the gravitomagnetic equations in case of low velocities for the sources, and the equation of motions for a low-velocity test mass, togheter with the general integral expressions of the gravitomagnetic fields in terms of the mass current for axisymmetric systems, in the two alternative fomulations of elliptic integrals and Bessel functions, also discussing their convergence properties. In Section 3 we consider the case of the rotation curve of generic razor-thin disks supported by circular orbits, and provide a regular series expansion of the rotational velocity in terms of the small expansion parameter = v 2 0 /c 2 (the square of the ratio between a charactertistic newtonian velocity of the disk and the speed of light c). In Section 4 we consider the explicit case of baryonic exponential disk, and we obtain essentially identical newtonian and GR rotation curves. The result is then confirmed with the aid of the Kuzmin disk, where the first order GR correction can be computed analytically. In Section 5, we relax the assumption of razor-thin disks, and we derive the gravitomagnetic Jeans equations for collisionless axisymmetric systems supported by two-integrals phase-space distribution function. This allows to investigate finite-thicknes GR effects on the rotational velocity in the equatorial plane. The case of the Miyamoto-Nagai disk is studied, again fully confirming the results for razor-thin disks. Section 6 concludes, while in the Appendix mathematical details are provided.

THE GRAVITOMAGNETIC EQUATIONS
In this paper, following previous works we adopt the gravitomagnetic formulation of GR in the weak-field limit, for low velocities and for steady motions of the sources of the gravitational field, where the equations, truncated at the first order (inclusive) of the source velocity in unit of the speed of light c, reduce to (see, e.g., Landau andLifshitz 1971, andPoisson and, for a detailed description of higher-order post-newtonian expansions; see also Rindler 1997, Lynden-Bell and Nouri-Zonoz 1998, Clark and Tucker 2000, Costa and Natário 2021, Mashhoon et al. 1999, Mashhoon 2008, Ruggiero 2021, Ruggiero and Tartaglia 2002. We indicate with ∧ the vector (cross) product, ∇ is the usual nabla operator, and j = ρv is the gravitational current density, being ρ and v the mass density and velocity field of the sources. Notice that at this expansion order the gravitoelectric field E = −∇φ is just the newtonian field produced by ρ, so that in Equation (1) the GR effects arise only from the gravitomagnetic Ampere law. The equation of motion of a test star at position x * is given by the analogous of the low-velocity limit of the Lorentz force (e.g., see Feynman 1977, Jackson 1998, where B(x * ) is the gravitomagnetic field produced at x * by the total gravitational current density distribution: notice the minus sign in front of the gravitomagnetic force, instead of the plus sign appearing in the electromagnetic case. For future use (see Section 4) it is useful to distinguish between the velocity v * of the test star, and the velocity v(x * ) of the mass current of the sources at x * , even though we anticipate that in the case of razor-thin disks made by circular orbits (see Section 3), necessarily v * = v(x * ).
Equations (1) are formally coincident with the (stationary) Maxwell equations, so the mathematical treatment is the same as in electrodynamics: however, due to some mathematical subtlety of the present astronomical problem, it is useful to list the most important properties of the field B in the general case. As well known (e.g., see J98), for a current density j well behaved at infinity where the first expression is the Biot-Savart law, . . . is the standard euclidean norm, and A is the gravitomagnetic potential vector in the Coulomb gauge, particularly appropriate for the magnetostatic case of steady currents. The Biot-Savart law can be proved by carring the ∇ operator acting on x (hereafter ∇ x ) under integral sign in the second expression above, and then using the general identity ∇ , where in the present case the meaning of the functions g(x, y) and f (y) is obvious, and so ∇ x ∧ (f g) = (∇ x g) ∧ f . A useful (and less known) equivalent expression of the Biot-Savart law can be obtained from the last identity by recognizing that ∇ where the last expression follows from the identity ∇ y ∧ (f g) = g ∇ y ∧ f + (∇ y g) ∧ f . As the volume integral of ∇ y ∧ (f g) vanishes for well-behaved fields at infinity, we finally obtain the alternative expression that we will use later on. Before addressing our specific problem, it is important to recall some convergence property of the fields A and B: notice that the only troublesome points x are those inside the current distribution, when the denominators in the integrands vanish for y = x. The following results can be easily established. 1) For well-behaved, genuinely three dimensional currents, the integrals in Equation (3) converge absolutely, i.e., the volume integral of the norm of the integrand converges, as can be seen by changing variable y = r + x at fixed x, and using spherical coordinates for r. It follows that each component of the A and B fields is also absolutely convergent. The Fubini-Tonelli theorem then assures that the integrals can be evaluated as repeated integrals, and that the order of integration does not matter, even though integrable singularities over sets of null measure can appear: from the physical point of view, such singularities have no consequences, but attention is to be payed in numerical studies.
2) For well-behaved razor-thin currents in the z = 0 plane j(y) = χ(y 1 , y 2 )δ(y 3 ), where δ is the Dirac delta-function and χ is the surface current density, the A and B integrals in Equation (3) again converge absolutely for all points x outside the current plane, and A also for points in the z = 0 plane, as can be proved with the change of variables y = r + x, and expressing r in cylindrical coordinates. The absolute convergence of B in the z = 0 plane cannot instead be established just from boundedness of χ , as the corresponding integral in Equation (3) diverges logarithmically for r → 0. However, from specialization of Equation (4) to razor-thin currents where e z = (0, 0, 1), and proceeding as in point (2), it follows that absolute convergence of B for points in the disk (where the second integral above vanishes from well known properties of δ , and the fact that χ is independent of y 3 ) is guaranteed under the additional request that ∇ y ∧ χ is well-behaved, a condition obeyed by the disks considered in this paper.
Summarizing, the general results above assure that not only for three dimensional currents, but also for points inside razor-thin disks, the gravitomagnetic field cannot develop nasty singularities if the surface density of the disk is sufficiently well behaved, and B and A are given by absolutely converging integrals.

The axisymmetric case
We now restrict to the case of axisymmetric systems, the subject of this paper. In cylindrical coordinates x = Re R + ze z , where e R = (cos ϕ, sin ϕ, 0) and e z = (0, 0, 1), so that ρ = ρ(R, z); we will not discuss how to obtain the associated newtonian gravitational potential φ(R, z), a problem fully addressed in the literature (see, e.g. BT08, C21). For the moment, purely circular orbits are considered for the sources, with a velocity field v = v(R, z)e ϕ , where e ϕ = (− cos ϕ, sin ϕ, 0), so that j = j(R, z)e ϕ , j = ρ(R, z)v(R, z).
The condition of circular orbits will be relaxed in Section 4, however notice that this is a quite natural idealization when considering the rotation curve produced by razor-thin disks in their equatorial plane (see Section 3). Obviously, in a three dimensional system, a purely circular current field outside the equatorial plane requires a vertical pressure gradient (or a vertical velocity dispersion field, as the case in Section 4). In order to evaluate the integrals in Equation (3), we introduce the source coordinates y = ξe R + z e z , with e R = (cos ϕ , sin ϕ , 0), and e ϕ = (− cos ϕ , sin ϕ , 0). It is a simple exercise to show that for axisymmetric currents the fields B and A are rotationally invariant (i.e., independent of ϕ), so that they can be evaluated without loss of generality at ϕ = 0, where x = (R, 0, z). For the potential vector A, it is trivial to prove that for sufficiently regular currents (also in the razor-thin case), not only A z = A · e z = 0 everywhere, in accordance with Point 3) above, but also A R = A · e R = 0 everywhere, while the only non-zero component is where d 3 y = ξdξdϕ dz . An analogous treatment of the B field in Equation (3) shows quite easily that B ϕ = B · e ϕ = 0 everywhere 2 , while therefore B R (R, 0) = 0 for currents with reflection symmetry about the equatorial plane, j(R, z) = j(R, −z), and in particular for razor-thin currents, again in accordance with Point 3) above. Notice that from Equation (A1) it follows that for the axisymmetric currents in Equation (7) ∇ as a sanity check, it is possible to verify by direct evaluation that Equations (8)-(10) are in fact equivalent to Equation (9). Notice also that the analogous of Equation (10) can be applied to the current in Equation (7), to compute the numerator of the integrand in the alternative formulation of the Biot-Savart law in Equation (4). Summarizing, the gravitomagnetic fields in the plane z = 0, produced by the axisymmetric currents in Equation (7) with reflection symmetry plane z = 0 (and in particular in the razor-thin case), can be written in full generality as A = A ϕ (R)e ϕ and B = B z (R)e z , where the functions A ϕ (R) and B z (R) are obtained from Equations (8)-(9) setting z = 0.

The field in terms of complete elliptic integrals
In potential theory it is customary to integrate Equations (8)-(9) first on ϕ , as the knowledge of the specific form of j(R, z) is not required. From the general results mentioned in Points 1) and 2) after Equation (4), the Fubini-Tonelli theorem assures that this is legitimate, and at worst only integrable singularities over zero-measure sets occur. In fact, from Equations (8)

and (A5)
where the function F 1 (R, ξ, ∆z) can be expressed in terms of complete elliptic integrals of first and second kind, K(k) and E(k), and k is given in Equation (A6). Notice that according to Equation (A8) F 1 presents an integrable singularity for k = 1, i.e. over the ring ξ = R at z = z, of easy treatment in numerical applications. The components of the B field can be obtained from differentiation of A ϕ (R, z) by using Equations (10), (11), and (A3): however, when using the fornulation in terms of complete elliptic integrals, it is convenient to avoid explicit differentiation, and instead work with the Biot-Savart law in Equation (9), integrating first over ϕ , and then proceeding as follows. For the radial component B R the resulting kernel in the integrand is a perfect differential of F 1 with respect to both z or z : the first possibility just corresponds to evaluate −∂A ϕ /∂z as required by Equation (10). In the second option one avoid differentiation of the elliptic integrals performing integration by parts with respect to z (provided the current j is well-behaved for |z | → ∞, the usual situation), and noticing that ∂F 1 /∂z = −∂F 1 /∂z . For the vertical component B z we also recognize the kernel as an exact differential with respect to ξ of F 0 in Equation (A4), and again we can avoid differentiation of the elliptic integrals integrating by parts (provided j is well-behaved for ξ → ∞). One finally obtains notice that the first expression for B z can be also obtained from Equation (10), by using the non-trivial identity (A7). Of course, the two second identities above are not unexpected: they are just how the Biot-Savart law in Equation (4) reduces in axisymmetric systems. These expressions are particularly useful when working with regular currents, as the integrable singularity in F 0 and F 1 can be explicitely taken into account. In particular the second expression for B z not only shows again that the field is wellbehaved inside a regular razor-thin current, confirming the conclusion in Point 2) above, but also shows that for a surface density current with an abrupt radial truncation, the field diverges at the edge, in analogy with the well known feature of the rotation curve produced by truncated disks (e.g., see Casertano 1983, see also exercises 5.4-5.6 in C21).
2.3. The field in terms of Bessel functions As well known, the integrals appearing in potential theory can also be expressed in an alternative formulations, for example as Fourier-Bessel series, obtained by using the apparatus of Green functions and Hankel transforms (e.g., Toomre 1963, see also J98, BT08, C21). This approach is very elegant from the analytical point of view, but the numerical implementation is not straightforward, due to the oscillatory nature of Bessel functions and the slow convergence of their integrals. In practice, we substitute Equation (A9) in then second of Equation (3) and we perform a first integration over ϕ . The current in Equation (7) is a vector quantity, and so in principle two integrals should be performed over the components of e ϕ ; however, it is possible to reduce the computation to a single integration, observing that the two components of e ϕ are just the real and imaginary parts of the complex number ie iϕ . We therefore introduce the complex current j * = ie iϕ j(ξ, z ), we determine the complex vector potential A * with a Fourier-Bessel expansion, and finally we switch back to the vectorial representation by separating the real and imaginary parts. The integration of j * over ϕ is elementary, and only the m = 1 component survives, producing 2πiδ m1 j(ξ, z ). Therefore, the final expression is in terms of Hankel transform of index 1 for the current (see, e.g., J98 for the case of the potential vector of a circular spire). The last step is the evaluation of the real and imaginary parts of the resulting expression, and it is immediate to show that the result is just proportional to e ϕ , i.e., we prove again that the potential vector is just A = A ϕ e ϕ . In particular, where is the Hankel transform of order 1 of the current density: reassuringly, by inverting order of integration in the triple integral in Equation (13), and performing first the integration over λ, eq. (2.12.38.1) in Prudnikov et al. (1986, hereafter P86) proves that the resulting expression coincides 3 with Equation (11). From Equation (10) we then obtain the two components of the B field where the expression for B z derives from the identity d[xJ 1 (x)]/dx = xJ 0 (x). The proof of equivalence of Equations (12)-(15) is important but quite laborious, and requires some comment. For B R we move λ in front of J 1 (λR) inside the integral over z , we recognize a derivative with respect to z of the exponential factor and integrate by part, we then exchange order of integration and finally integrate in λ by using eq. (2.12.38.1) in P86. For B z , the approach is to move the λ in front of J 0 (λR) inside the Hankel transform 1 , use the identity λJ 1 (λξ) = −dJ 0 (λξ)/dξ and integrate by parts over ξ. Finally, we invert order of integration and use again eq. (2.12.38.1) in P86. With this first approach we proved the equivalence of Equation (15) with the second expressions for B R and B z in Equation (12). However, there is a second approach that proves the equivalence of Equation (15) with the first identities in Equation (12), and that also help to clarify an important convergence issue. In both the triple integrals in Equation (15) we exchange order of integration, and we evaluate first the integrals over λ, that belong to the family ∞ 0 xe −px J µ (ax)J ν (bx)dx, with the aid of eq. (2.12.38.2) in P86, where in particular p = |∆z| and x = λ. It can be proved (for example by asymptotic expansion of the Bessel functions for large values of their argument) that for ∆z = 0 the two integrals over λ diverge: however, if the limit for ∆z → 0 is evaluated after integration, then the first identities in Equation (12) are recoverd 4 after expressing ∂F 1 /∂z and −∂F 0 /∂ξ in explicit form.

SERIES SOLUTION FOR THE GRAVITOMAGNETIC EQUATIONS: THE RAZOR-THIN DISK WITH CIRCULAR ORBITS
Having established the general setting of the problem, we are now in position to discuss the radial dependence of the circular velocity of a test mass (a star or a gas cloud) in the plane of a razor-thin disk made of field stars in circular orbits, so that where Σ(R) and v(R) are respectively the mass surface density, and the circular velocity of the disk, and is the radial profile of the two-dimensional current; according to the orientation of the coordinate system, a positive v means counterclockwise rotation. From the assumption of circular velocities, it follows necessarily that the modulus of the circular velocity of the test stars and of the rotating disk coincide, i.e. v * (R) = |v(R)|, and quite naturally we also assume that v * (R) = v(R)e ϕ , i.e. that the test stars rotate as the stars of the disk. From Equation (2), where the field E is just the newtonian gravitationl field produced by Σ(R), and from the fact proved in Section 2.1 that in the z = 0 plane B = B z (R)e z , we obtain the scalar equation for the circular velocity v at each radius R where v 0 is the circular velocity of the disk in the newtonian case and where, from Equations (12)-(15), the gravitomagnetic field in the disk plane reduces to the equivalent expressions The general discussion in Sections 2.2-2.3 shows that the first expression contains a logarithmic integrable singularity at ξ = R, while in the second the order of integration cannot be exchanged with the Hankel transformχ 1 given by Equation (14) applied to χ(R). Notice that with B z [χ] we indicate the linear operator acting on the mass current profile (i.e., on the v(R) profile, if the surface density Σ(R) profile is assigned); notice also that Equation (17)  In practice, by solving Equation (17) for assigned Σ(R), we obtain (at the order of the gravitomagnetic equations) the "selfconsistent" circular velocity of the disk produced by the combined effects of the newtonian field and of the gravitomagnetic potential produced by the rotation curve itself.
For a razor-thin disk of total mass M d and scale-lenght R d , it is natural to normalize lenghts to R d , surface densities to M d /R 2 d , and velocities to GM d /R d : a tilde over a quantity indicates normalization to its associated scale. Accordingly, Equation (17) is recast in dimensionless form asṽ where is a dimensionless parameter arising naturally from the normalization of B. In fact, from Equation (18) we obtain where B is the dimensionless gravitomagnetic field in the disk. Notice that for a disk the parameter is nothing else that the natural proxy for the quantity v 2 /c 2 mentioned in the Introduction, where v is a charaterstic velocity associated to the newtonian gravitational field (see the solid lines in Figure 3). Therefore, in real galaxies is very small: for example 10 −6 , in a disk galaxy with a characteristic rotational velocity of 300 km/s. Accordingly, we represent v as a regular asymptotic series in powers of where the exponent α is to be fixed by order balance, and of course for → 0 we reobtain the newtonian case. Inserting the expansion above in Equation (17), from the linearity of B z [χ] on velocity, it follows necessarily that α = 1, and at points wherẽ v 0 = 0ṽ whereχ 0 is the (normalized) surface density current profile corresponding to the newtonian rotation curveṽ 0 : from O(B z [χ 0 ]) = , it follows that the termṽ 1 depends on the gravitomagnetic field produced by the newtonian current only. At the origin, the . The field is computed for the newtonian current χ 0 , as required by Equation (21), and it looks similar to the field that would be produced by the current of a circular spire; the maximum of the newtonian current is located at R 0 0.574R d for the exponential disk, and at R 0 0.535R d for the Kuzmin disk. The maximum of Aϕ is located at R A 1.12R d for the exponential disk and at R A 0.85R d for the Kuzmin disk. Bottom panels: the corresponding B field, in units of GM d /R 3 d . At z = 0, Bz is positive in the inner regions of the disk, and negative outside; the radius R B at which Bz = 0 is determined by the critical point of RAϕ(R), and so does not coincide with the position of the maximum of Aϕ(R), as apparent from the figures. Numerically, R B 2.2R d for the exponential disk, and R B 1.68 for the Kuzmin disk. only point 5 at finite R whereṽ 0 can vanish, the perturbation analysis shows that again α = 1, butṽ 1 =RB[χ 0 ]: however, as R = 0, no difference arises with Equation (21), that therefore can be used uniformly over the whole radial range. Moreover, Equation (21) shows that the GR corrections due to a counterclockwise current lead to a decrease of the rotational speed where B z [χ 0 ] is negative, and an increase where B z [χ 0 ] is positive: of course, it is immediate to verify that the effects on the circular speed are the same for a global inversion of the rotational velocity v 0 → −v 0 , due to the associated change of sign of B z [χ 0 ]. It is important to remark that Equation (21) allows for an alternative interpretation: one could just consider the field B z [χ 0 ] at the r.h.s of Equation (17), and then solve in closed form the resulting quadratic equation for v(R). After selecting the sign so that for → 0 the newtonian profile is reobtained, an expansion of the solution for → 0, and truncation at the first order, gives again Equation (21), proving again the the linear gravitomagnetic formulation of the problem leads to a regular perturbation problem (e.g., see Bender and Orszag 1978, Chapter 7). Of course, even if the problem admits a regular perturbation approach, with a solution reducing to the newtonian one in the limit = 0, the functionṽ 1 could attain very large (but finite) values, compensating the small (but non-zero) value of , and producing a perturbation of the same order of magnitude 6 or even larger than the newtonian term: this would be a strong support to the possibility of a GR origin of the flat rotation curve of disk galaxies. In the next examples, based on realistic disk density profiles, we show however that this is not the case, andṽ 1 remains small, with absolute values well below the unity.
3.1. The exponential disk In our first application we consider the razor-thin exponential disk of total mass M d and scale-lenght R d , the standard model used to describe the stellar density distribution of disk galaxies (e.g., BT08, Bertin 2014). The surface density is given by and in newtonian gravity the gravitational potential in the equatorial plane is so that the associated circular velocity is where I m and K m are the modified Bessel functions of order m (e.g., BT08, C21). As already remarked in the Introduction, the maximum of v 0 is reached at R 2.15R d , and in the range 1.5 <R < 3 the curve is almost flat even in absence of DM: notice that inside 3R d the disk already contains 0.8M d . The newtonian current surface density χ 0 for a disk made by purely circular orbits is then obtained from Equations (22)-(24): for a global counterclockwise rotation, it reaches the maximum at R 0 0.574R d .
For the exponential disk it turns out that the most efficient way to compute the gravitomagnetic field is to use the formulation in terms of elliptic integrals. In particular, as the current density is a regular function of R, the general discussion about convergence leads to use Equation (11) for the evaluation of A ϕ , and the first and second integrals in Equation (12) for the evaluation of B R Figure 3. Radial trend of the rotational velocity components in Equation (20) normalized to GM d /R d , for the exponential (left) and Kuzmin-Toomre (right) disks. The solid line is the newtonian rotation curve v 0 , and the dotted line is the v 1 perturbative term, that should then be added to v 0 after multiplication by 10 −6 . Notice how, in principle, the effect of the GR gravitomagnetic field at large radii is a decrease of the circular velocity with respect to the newtonian case. and B z , respectively: of course, in the disk plane the latter expression reduces to the first case in Equation (18). Numerically, the integrable divergence at ξ = R is easily treated by splitting the integral from the origin to (1 − η)R, and from (1 + η)R to infinity, and reducing η until acceptable convergence is reached (actually, thanks to absolute convergence, the two η's not need to be the same). In the following experiment, convergence was already reached for η = 10 −3 , with stable results at 5 significative digits according to Mathematica NIntegrate function. Importantly, the numerical agreement with the alternative formulation in terms of Bessel functions has been also verified over all the disk.
In the top left panel of Figure 1, the only nonvanishing component A ϕ of the potential vector (in units of GM d /R d ), is shown in the meridional plane, where the overall striking similarity with the field produced by an "effective" circular current loop of radius R A 1.12R d is apparent, where R A is the position of the maximum of A ϕ (R); quite obviously, R A does not coincide with the position R 0 of the maximum of the current. In the bottom left panel, the associated (normalized) B field is shown. As expected, the field is similar to that of a circular (counterclockwise circulating) current loop, with B z positively directed in the inner regions of the disk (where the Biot-Savart fields of each current ring composing χ 0 reinforces), and negatively directed in the outer regions, due to the cooperative (negative) contribution of the current in the inner regions of the disk, while the counteracting positive contributions of the current in the outer regions of the disk are less and less important, due to vanishing of the current at large radii. The radial trend of B z [χ 0 ] in the disk plane is represented by the solid line in Figure 2. Notice that the radius at which B z (R) = 0 is R B 2.2R d , not coincident with R A , as from Equation (10) R B is given by the critical point of RA ϕ (R), while R A is the critical point of A ϕ (R). Having determined the field B z for the exponential disk, from Equation (21) we finally compute v 1 . In the left panel panel of Figure 3 the solid line shows the (normalized) profile of the newtonian rotation curve v 0 , while the dotted line v 1 . After multiplication by 10 −6 , it is clear the the effects of GR on the curve of the exponential disk are well below the possibility of any practical detection, with differences with respect to the newtonian rotation curve well below 1 m/s: at the description level of gravitomagnetism, GR does not produce any effect on the rotation curve of the disk. In any case, it is interesting to notice that for the reasons explained above, the sign of the gravitomagnetic field in the Lorentz equation actually produces a decrease of the rotational speed at large galactocentric distances!

The Kuzmin-Toomre disk
Having determined the solution for the exponential disk, in order to confirm the obtained results we move to explore a different disk model. We focus on the Kuzmin-Toomre razor-thin disk (Kuzmin 1956, Toomre 1963, see also BT08, C21), an idealized model widely used stellar dynamics for its mathematical simplicity. The surface density-potential pair of the disk of total mass M d and scale-lenght R d is given by where the potential is restricted to the disk plane, and the newtonian circular velocity and the surface mass density current are given by The maximum of v 0 is located at R = √ 2R d , and that of χ 0 at R 0 = 2/7R d 0.535R d , a value curiously similar to that of the exponential disk. We now show that v 1 can be obtained in closed form by using the Fourier-Bessel identity in Equation (18). In fact, from eq. (6.565.4) of GR07, or eq. (2.12.4.28) of P86, we obtain the Hankel transform of order 1 for the newtonian mass current asχ where Γ is the complete gamma function, and we used the fact that for modified Bessel functions K −ν = K ν ; notice thatλ is a dimensionless quantity. The last integral in Equation (18) can be expressed analytically thanks to eq. (6.576.3) of GR07 or eq. (2.16.21.1) of P86, and finally from Equation (21) where 2 F 1 is the standard hypergeometric function. Notice that equation above can be recast in terms of elliptic integrals, resulting in perfect agreement over the whole radial range with the values of v 1 (R) obtained from Equation (21) and numerical integration of the last expression of B z in Equation (12), in a reassuring validity check. As a mathematical curiosity, we also notice that for the newtonian current χ 0 of the Kuzmin-Toomre disk, not only the field B z in the disk plane can be obtained explicitely, but also Equations (15) can be solved analytically over all the space by using Fox H function,s after expressing Equation (27) in terms of these functions (e.g., Mathai et al. 2010, and in particular eq. (2.25.3.2) in Prudnikov et al. 1990): the resulting expression is however of no practical use, and so not reported here. In the right panels of Figure 1 we plot the normalized potential vector of the disk in the meridional plane, and the associated gravitomagnetic field: the maximum of the gravitomagnetic potential vector is reached at R A 0.85R d . The similarity with the case of the exponential disk, and also the range of values spanned by the normalized fields is remarkably similar, even if the radial density profiles of the two disks are quite different, especially at large radii; notice also the similar behavior of B z in Figure 2, with positive values in the inner regions of the disk (i.e., for R < R B 1.68), and negative values outside. In Figure  3 we plot the normalized profiles of v 0 and v 1 . It is again apparent how v 1 remains limited over all the radial range, confirming that the contribute to the rotational velocity of the first order perturbative term in the velocity expansion is fully controlled by the smallness of , and again the expected GR corrections to the newtonian rotational velocity of the galaxy are well below the detection limit, with differences between the newtonian and the gravitomagnetic GR curve well below 1 m/s. Finally, in accordance with the change of sign of B z along the equatorial plane, the gravitomagnetic effects produce a decrease 7 of the rotational speed at large distances from the galaxy center.
Therefore, from the analysis conducted so far, it seems a quite robust conclusion that GR, at the level of the gravitomagnetic weak field approximation, does not affect the rotation curve of self-gravitating barionic disks with realistic density profiles, and made mostly by circular orbits. In the next Section we relax the assumption of a two-dimensional density distribution, and we also consider non-circular orbits for the stars producing the field in the system.

THE GRAVITOMAGNETIC JEANS EQUATIONS
After the discussion of razor-thin disks, we address the problem of the expected effects of GR on the internal dynamics of genuinely tridimensional (collisionless) astronomical systems, such as finite-thickness disk galaxies, or axisymmetric elliptical galaxies. We first rigorously derive the gravitomagnetic modification of the Jeans equations, and then we restrict to stationary axisymmetric systems. As usual, we indicate with v the phase-space velocity, so that the density ρ and the streaming velocity field v of the system are given by where f (x, v, t) is the phase-space distribution function (hereafter DF, see BT08, Bertin 2014, C21), and a bar over a quantity indicates its average over the velocity space. In the following we will write where sums over repeated indices hold in Cartesian coordinates. Suppose the considered stellar system is self-consistent, i.e., the motion of each star is determined by the field produced by the combined effects of all the other stars. As we are in the low velocity limit, when the Biot-Savart law can be interpreted as the sum of the Lorentz fields produced by each moving charge (an interpretation in general not true for currents produced by particles of arbitrary large velocities, see e.g. J98, Chapter 5; Griffiths 1999, Chapter 5;Feynman 1977, Chapters 13 and 21;Panofsky and Phillips 1962, Chapter 7), it follows that the acceleration experienced by a star at x with velocity v can be written by summing the r.h.s. member of Equation (2) over the DF, where φ(x) is the newtonian gravitational potential of the system and B[j], by virtue of the linearity of B on the current, is the gravitomagnetic field produced at x by the total streaming current of the system. The hierarchy of the Jeans equations of increasing order is then obtained by taking moments over the velocity space of the collisionless Boltzamnn equation (e.g., see BT08, C21) ∂f ∂t here written in Cartesian coordinates. After multiplication of Equation (31) by 1, v i , v i v j , etc., and integration, the three second order moments equations (the gravitomagnetic modification of the usual Jeans equations of stellar dynamics) are given by where in particular, σ 2 ij are the components of the second-order velocity dispersion tensor. In Appendix B, the general gravitomagnetic time-dependent Jeans equations are also reported in cylindrical coordinates.
We now to restrict to the axysimmetric stationary case. In the usual stellar dynamical case, from the Jeans theorem one would assume a phase-space DF depending on the two classical integrals of the motion (per unit mass) E and J z , i.e. the orbital energy, and the z-component of the angular momentum of stellar orbits. As well known, under this assumption (e.g., BT08, C21), important constraints on the allowed velocity moments follow, namely 1) the only possible streaming motions are in the azimuthal direction, i.e., v = v ϕ (R, z)e ϕ , and 2) the velocity dispersion tensor is in diagonal form at each point in the system, with σ 2 R = σ 2 z ≡ σ 2 , while only v 2 ϕ = σ 2 ϕ + v 2 ϕ remains determined; when σ 2 ϕ = σ 2 everywhere, the system is said isotropic. Therefore, for stationary stellar dynamical two-integral axisymmetric systems with f (E, J z ), the current would be of the "circular type" considered in Equation (7), with j = je ϕ = ρv ϕ e ϕ : notice that, even if the current associated with the streaming velocity field is circular, the orbits of the individual stars contributing to the resulting gravitomagnetic field are in general not circular, and so we are considering an orbital structure more complicated than that of the razor-thin disk in Section 3.
However, in the case of stationary collisionless stellar systems with gravitomagnetic forces, the situation is not so simple, because while E is still an integral of motion, J z is not conserved along orbits and so it cannot be used as an isolating integral in the DF. Fortunately, the problem is solved as follows. As well known, the lagrangian (per unit mass) associated with Equation (30) is where at variance with the EM case (e.g., see Landau and Lifshitz 1971) the + sign in the generalized potential derives from the -sign in the gravitomagnetic Lorentz force. Suppose now that v = v ϕ (R, z)e ϕ , so that A = A ϕ (R, z)e ϕ from the results in Section 2.1. The Euler-Lagrange equations show immediately that J z is not conserved, but a second integral of motion exists, I 2 = J z − RA ϕ . The Jeans theorem then in principle allows for a two-integrals DF f (E, I 2 ), and the interesting question arises wether such a DF is consistent with the assumption of streaming motions with v = v ϕ (R, z)e ϕ used to prove the existence of I 2 .
The answer is in the affirmative, and it is easy to prove that the properties (1) and (2) mentioned above for systems supported by a f (E, J z ) are preserved by the generalized f (E, I 2 ), in particular the fact that v R = v z = 0. Therefore, the vertical and radial Jeans equations (B2) for a stationary, axisymmetric two-integrals system with gravitomagnetic forces, supported by a DF of the family f (E, where ∆ = v 2 ϕ − σ 2 , and in the isotropic rotator case ∆ = v 2 ϕ . Finally, the azimuthal Jeans equation vanishes identically also in the gravitomagnetic case, as obvious from the last of Equations (B2). Notice how at variance with the classical newtonian case, now the vertical and radial velocity dispersions depend on the ordered rotational field v ϕ (R, z), yet to be determined.
Equations (35) can be formally solved in terms of σ 2 and ∆ as in the usual stellar dynamical case, obtaining where σ 2 0 is the newtonian solution, and U (R, z) accounts for the gravitomagnetic effects. Moreover, some algebra shows that where the commutator C[ρ, φ] is the standard newtonian solution (e.g., Hunter 1977, Barnabè et al. 2006, see also C21), while the term contains the gravitomagnetic effects. The triple integral in the second expression above (obtained from the first by exploiting the solenoidal nature of B) is particularly well suited for numerical computations when the current is analitically available (see next Section), while the last expression is obtained from the first taking into account that B = ∇ ∧ A, Equation (10), and finally integrating by parts with jA ϕ = 0 for z → ∞.
As noticed, the difficulty in the solution of the gravitomagnetic Jeans equations is the fact that the force field depends on the azimuthal streaming velocity distribution, a quantity that it is not known in advance. In principle, to solve the gravitomagnetic Jeans equations one should 1) impose a decomposition on the still unknown azimuthal velocity field, for example by using the widely adopted Satoh (1980) ansatz v 2 ϕ = k 2 ∆ (e.g., see C21), 2) solve the resulting nonlinear integrodifferential Equation (37) for v ϕ and obtain B R [j], 3) integrate Equation (36) and obtain σ 2 , 4) conclude by determining ∆ and the tangential velocity dispersion σ 2 ϕ = σ 2 + (1 − k 2 )∆ = σ 2 + (1/k 2 − 1)v 2 ϕ . Quite obviously, a considerable simplification is achieved by exploiting the fact that in Equation (35) the fields B R and B z are proportional to the small parameter = GM * /(r * c 2 ), where M * is the total mass of the system, and r * its scale-lenght. As in the razor-thin disks, it follows immediately that the first order perturbative terms in the natural expansions σ 2 = σ 2 0 + σ 2 1 . . . , ∆ = ∆ 0 + ∆ 1 . . . , and v ϕ = v ϕ0 + v ϕ1 . . . , are related to the functions U and V evaluated over the newtonian current j 0 = ρv ϕ0 , where v ϕ0 (R, z) is the streaming rotation velocity obtained from the classical Jeans equations.
We are now in position to evaluate the effects of GR corrections on the circular velocity v( , R) = v 0 (R) + v 1 (R) + . . . of a test star in the equatorial plane of an axisymmetric collisionless system described by Equations (35). The only difference with the razor-thin disks in Section 3 is that now the current is provided by the azimuthal streaming velocity field of the system. From Equation (30), the analogous of Equations (17) and (19) still apply, and B z = GM * /r 3 * B is computed in the disk plane, and finally j 0 = GM 3 * /r 7 * 0 . The expression for the first order correction term to the streaming velocity (restricting for simplicity to the isotropic rotator case for the stellar distribution of the galaxy) is instead obtained from Equation (37), considering that at the lowest order V = O( ), and so from v ϕ = v ϕ0 + v ϕ1 . . . , after order balance and some simplification, one gets where V = (GM 2 * /r 5 * )V is computed from Equation (38).
4.1. The Miyamoto-Nagai disk We now solve the gravitomagnetic Jeans equations, and evaluate Equations (39)-(39) for the widely used Miyamoto-Nagai (1975, Nagai andMiyamoto 1976, see also BT08, C21) disk of total mass M * and scale-lengths a and b, with newtonian densitypotential pair where ζ = √z 2 + s 2 ,R = R/a,z = z/a, and the newtonian circular velocity in the equatorial plane is given by the dimensionless parameter s = b/a measures the disk flattening, and for b = 0 the model reduces to the Kuzmin disk of total mass M * and scale-lenght a, so we are in the condition to explore the GR effects of a non-zero thickness on the rotation curve, by comparison with the case in Section 3.2. The solid lines are the newtonian components, and the dotted lines are the first perturbative term, that should be added to the corresponding newtonian term after multiplication by 10 −6 . Note how the circular velocity is slightly larger than the streaming velocity, the well known phenomenon of asymmetric drift produced by the disk vertical velocity dispersion support.
We restrict for simplicity to the isotropic rotator, with a conterclockwise streaming velocity field, so that the newtonian azimuthal streaming velocity field and the associated current density are (e.g., Ciotti andPellegrini 1996, Smet et al. 2015). Notice that in the equatorial plane the location of the maximum of j 0 (R, 0) can be evaluated analytically as a biquadratic equation, and for s → 0 the position coincides with that of the maximum for the Kuzmin disk.
In the left panel of Figure 4 we show the results for the circular velocity, to be compared with the analogous Figure 3 for the two razor thin disks made by circular orbits. The solid line represents the newtonian circular velocityṽ 0 in Equation (42), while the dotted line gives the correction term in Equation (39), where the gravitomagnetic field is now produced by stars that in general are not moving in circular orbits. Again, for realistic mass distributions, the perturbation term shoud be multiplied by factor ≈ 10 −6 before addition to the newtonian term: the GR effects appear again to be completely negligible and, in analogy with the case of razor thin disks, positive in the inner regions of the galaxy, and negative at large radii, reducing there the rotational speed of a test mass. In the right panel we show instead the streaming velocity of stars in the disk equatorial plane, both the newtonian component (solid line), and the lowest order GR term (dotted line). Again, it is important to remark that the orbits of the individual stars producing the v ϕ in the equatorial plane are in general not circular, and not restricted to the equatorial plane. Notice how the newtonian streaming velocity is slightly below the newtonian circular velocity, the well known phenomenon of asymmetric drift, due to effect of the vertical component of the velocity dispersion tensor: (e.g., see BT08). Notice also how the asymmetric drift (a small but well measured phenomenon in real galaxies) is several order of magnitude larger than the GR corrections.
The results for the Miyamoto-Nagai disk reinforces the previous conclusions, i.e., also in three dimensional galaxies, with stars moving on non-circular orbits, the GR effects predicted by linear gravitomagnetism on the rotational velocities are a factor of ≈ 10 −6 smaller than the newtonian predictions. If newtonian gravity requires DM, unless "strong" weak-field GR effects (beyond the linear gravitomagnetic description) actually dominates the dynamics of the galaxies, exactly the same amount and distribution of DM is required in GR and in newtonian gravity, in order to reproduce the observed rotational curves.

DISCUSSION AND CONCLUSIONS
Recently, it has been suggested that the observed flat radial profile of the rotation curves of disk galaxies at large galactocentric radii could be a peculiar effect of GR in rotating systems, instead of the dynamical signature of the presence of DM halos. The suggestion is somewhat puzzling, as observed galaxies are empirically in a weak-field regime, with v/c ≈ 10 −3 , and so GR corrections are expected to be very small: however, GR is a non-linear theory, and so important effects cannot be excluded in principle.
In this paper, instead of attempting a GR modelization of the rotation curve of some observed galaxy, we followed a different approach, and we built the rotation curve for purely barionic disks, structurally resembling real galaxies, and with total mass and scale lenght in the observed range. The rotation curves are computed both in newtonian gravity and in the (gravitomagnetic) weak field formulation of GR, and compared. After a detailed discussion of the mathematics, we considered first the case of razor-thin disks made by purely circular orbits. The orbital structure is highly idealized, but the surface density of the disk is the standard exponential profile adopted to model the stellar distribution in real disk galaxies; we also considered another disk model, the Kuzmin-Toomre model, to gain confidence with the results obtained for the exponential disk. Consistently with the gravitomagnetic approximation adopted, we impose "self-consistency" on the rotational velocity, i.e., we assume that the rotational velocity of the test particle (star, or gas cloud), and the rotational velocity of the stellar population producing the gravitomagnetic field, are the same at each radius, and we require that the rotational velocity is the solution of the Lorentz-like equation of motion. The resulting problem is shown to be a regular perturbation problem, with an intrinsinc order parameter ≈ 10 −6 , and the gravitomagnetic field is obtained from two different (but equivalent) methods, i.e., by using elliptic integrals and Bessel functions. It is found that the GR (normalized) perturbative term is of the order of the unity, with no indications of anomalously large emerging effects, and the resulting physical modification of the newtonian rotation curve is therefore of the order of ≈ 10 −6 , well below the possibility to detection in observed rotation curves, in agreement with the order-of-magnitude estimate obtained by Toth (2021). Curiously, in the external regions of the disk the gravitomagnetic GR perturbative term tend to decrease the rotational speed, as can be easily understood by considering the gravitomagnetic field as given by the sum of the fields produced by each current ring of the disk, and by the fact that the disk current in realistic disk galaxies decreases sufficiently fast at increasing radius, so that in the outer regions of the disk the gravitomagnetic field is dominated by the inner current distribution. An important warning follows: if one assumes a disk surface density distribution with a sufficiently slow radial decline at increasing R, it is quite easy to show that the sum of the fields of external rings can be very large (or even diverge) leading to the wrong conclusion of a significant GR effect in real galaxies: in theoretical/numerical studies of the present problem, the use of realistic density profiles for the disks it is of fundamental importance.
As doubts have been casted that disk thickness effects, or motions more complicated than purely circular orbits for the stars producing the gravitomagnetic fields, could be important, we also studied the case of the gravitomagnetic rotation curve in genuinely tridimensional stellar systmes. We first derived the gravitomagnetic Jeans equations starting directly form the collisionless Boltzamnn equation, and the associated Jeans theorem. In this case the equations contain naturally a vertical "pressure" term due to the vertical velocity dispersion. We showed how in axisymmetric systems the gravitomagnetic field depends on the azimuthal streaming velocity field of the system (even though the orbits of the single stars are in general not axisymmetric). In such a system, we can identity two different rotational velocities: the circular velocity of a tracer in the equatiorial plane, and the (circular) streaming velocity of the stellar population producing the gravitomagnetic field. For both cases we found that the GR perturbative term, after multiplication by the expansion parameter ≈ 10 −6 , is completely negligible over the newtonian term, so that also for collisionless axysimmetric systems the GR effects due to rotation appear to be well below detectability, and unable to reproduce a flat rotation curve at large radii any better than a newtonian model in absence of DM.
In conclusion, the study conducted in the paper excludes the possibility that gravitomagnetic GR effects can compensate by any detectable amount the keplerian fall of the rotational velocity that would characterize disk galaxies at large distances in absence of DM halos, and produce the observed flat profiles: from the observational point of view, in rotating disk galaxies DM is required by GR exactly as in newtonian gravity. Of course, if a full GR simulation for the same barionic disks used in this papers (with disk parameters corresponding to the observed weak field regime of real galaxies) convincingly proves the opposite, then among other things we should conclude that gravitomagnetic weak field approximation of GR cannot be used to describe the weak field regime in rotating galactic disks. I thank Giuseppe Bertin, Antonio Mancino, Bahram Mashhoon, Jerry Ostriker, Silvia Pellegrini, Francesco Pegoraro, Renzo Sancisi, Massimo Stiavelli, and an anonymous Referee, for important comments. Alberto Parmeggiani (Department of Mathematics of Bologna University), is especially thanked for discussions about some mathematical aspect of this work.
For the considerations in Sections 2.1 and 2.2, we recall the important identity that can be proved by differentiation of the explicit expressions in Equations (A4)-(A5), or (with some work) more elegantly directly from their integral representation. In numerical applications, when using Equations (11)-(12), it is important to recall that E(k) ∼ 1 for k → 1 − , but K(k) ∼ − ln √ 1 − k 2 ∼ − ln √ 1 − k, so that near the ring ξ = R in the ∆z = 0 plane both F 0 and F 1 diverge as however the singularity is integrable and easily treated numerically. From the theory of Green functions it is possible to prove that in cylindrical coordinates (e.g., BT08, C21) where J m are the Bessel functions of first kind and integer order m, x = (R cos ϕ, R sin ϕ, z), y = (ξ cos ϕ , ξ sin ϕ , z ), and the Hankel transform for a function f (R) of the cylindrical radius readŝ

B. GRAVITOMAGNETIC JEANS EQUATIONS IN CYLINDRICAL COORDINATES
By using the standard approach of considering velocity moments of Equation (31), and by working on cylindrical coordinates, from the velocity moment of order 0 we obtain the continuity equation and from the three moments of order 1 (respectively over v z , v R , and v ϕ ) the three momentum equations where the B field is determined by Equation (3) computed over the galaxy streaming mass density current given by the second of Equation (30). Finally, if the system is in steady state, and supported by a phase-space DF f (E, J z − RA ϕ ), from the discussion in Section 4 it follows that Equation (B1) and the last of Equations (B2) are identically verified, while the first two of Equations (B2) reduce to Equation (35).