Relativistic corrections to the rotation curves of disk galaxies

We present a method to investigate the effect of relativistic corrections arising from large masses to the rotation curves of disk galaxies. The method employs a mean-field approximation and gravitational lensing. Applying it to a basic model of disk galaxy, we find that these corrections become important and magnified at large distances. The magnitude of the effect is sufficient to explain the galactic missing mass problem without requiring a significant amount of dark matter. A prediction of the model is that there should be a strong correlation between the inferred galactic dark mass and the galactic disk thickness. We use two independent sets of data to verify this.


INTRODUCTION
The total mass of a nearby disk galaxy is typically obtained from measuring its rotation curve and deducing from it the mass using Newton's dynamics. The rationale for this non-relativistic treatment is the small velocity of stars: v/c 1 sufficiently far from the central galactic black hole. However, the assumption that relativistic corrections are negligible may be questioned on several grounds. Inspecting the post-Newtonian [1] Lagrangian, e.g. for two masses M 1 and M 2 separated by r, shows non-Newtonian potential terms of the type G 2 M 1 M 2 (M 1 +M 2 )/2r 2 (G is the gravitational constant) that are independent of v, thus not suppressed at small v, and can be non-negligible for large enough M 1 and M 2 . These terms express the non-linear nature of General Relativity (GR), which arises from its field self-interaction: the gravitational field has an energy and hence gravitates too.
Field self-interactions are well-known in particle physics: Quantum Chromodynamics (QCD), the gauge theory of the strong force between quarks, features color-charged fields that self-interact. In fact, GR and QCD have similar Lagrangians, including self-interacting terms, as can be seen when the Einstein-Hilbert Lagrangian of GR is expanded in a polynomial form [2,3]. Field self-interaction in QCD, which causes quark confinement, exists even for static sources, as shown by the existence of numerous heavy quark bound states (in which v ≈ 0 for quarks) [4] and by classic numerical lattice calculations for v = 0 quarks [5]. This, as well as the correspondence between the respective terms of the GR and QCD Lagrangians, shows that for bodies massive enough a relativistic treatment is required regardless of their velocity. Finally, the measured speeds at the rotation curve plateaus are of several hundreds of km/s, e.g. 300 km/s (or v/c = 0.1%) for NGC 2841. They are similar to that of stars orbiting the central black hole of our galaxy and clearly display the relativistic dynamics expected in the strong regime of GR [6].
These arguments suggest that one should investigate the importance of relativistic dynamics in galaxies and how it affects the missing mass problem. From experi-ence with QCD, a non-perturbative approach is required to fully account for field self-interaction, making post-Newtonian formalism inadequate. In Refs. [2,3], a nonperturbative numerical lattice method was used. Here, we propose to approach the problem with a mean-field technique combined with gravitational lensing. There are several advantages of the approach compared to the lattice method used in [2,3]: (1) it is an entirely independent method, thereby providing a thorough check of the lattice result; (2) it is not restricted to the static limit of the lattice method and can be applied to systems with complex geometries; (3) it is significantly less CPU-intensive than a lattice calculation, and hence much faster; (4) it clarifies that the effect calculated in Refs. [2,3] is classical. The lattice approach -an inherently quantum field theory (QFT) technique -used in Refs. [2,3] may misleadingly suggest that a quantum phenomenon is involved. In fact, the classical nature of the effect is consistent with these lattice calculations being performed in the high-temperature limit in which quantum effects disappear, as discussed in Ref. [3]; (5) the lensing formalism is more familiar to astrophysicists and cosmologists, in contrast to lattice techniques with its QFT underpinning and terminology.

Method
The mean-field technique has been widely employed. In the context of gravitation, it has been used e.g. to calculate the propagation of gravitational waves near large masses [7]; to solve, by using the gravity/gauge correspondence, the non-perturbative bound-state problem of QCD [8]; or to derive the background field method introduced to quantize gravity [9]. Hence the mean-field approach is a common and diverse technique. We propose to use it together with the gravitational lensing formalism to compute the self-interaction of the gravitational fields that generates GR's characteristic non-linearities. A picture of the method would be that of a traditional gravitational field embedded in curved space-time, see Fig. 1.  1: Illustration of the background field method. The left sketch shows field lines between two masses or charges in flat spacetime. The right panel shows the field lines deformation in spacetime curved by a background field. In this article, the mass deforming spacetime is the total galactic mass and the effect on the field lines connecting two parts of the galaxy is determined by using the standard gravitational lensing formalism.
The instantaneous ("Coulomb") component of the force is given by field lines distorted by the curved space-time. We will compute here how gravity's field lines are distorted by a mass distribution in the same way as light is lensed by such distribution. Then the field line flux is computed to obtain the gravitational force including its self-interaction effects.
Electromagnetic and gravitational field lines are affected identically by curved spacetime (viz gravitational background field). In particle language, photons and gravitons follow the same null geodesics. To understand this in terms of forces, the photons and gravitions are both massless and thus have only two allowed spin projections in spite of being spin-1 (vector field) and spin-2 (tensor field), respectively. This can be understood in the following way: in Newton's gravity, a massless particle is lensed because of the energy-mass equivalence. However, the deflection is only half that predicted by GR (unless the deflected particle has no spin) because the Newtonian gravitational field is scalar and thus does not couple with the spin of a particle, leaving it unchanged. Likewise, the gravitoelectric part of the tensorial GR field interacts with a particle without affecting its spin -creating the deflection predicted in a Newtonian framework -while the gravitomagnetic part of the tensorial GR field couples with the spin, interacting by flipping it, hence doubling the deflection for particles with two allowed spin projections. Then, since the electromagnetic and gravitational fields have the same number of spin projections, the usual gravitational lensing formalism can be used directly for gravitational field lines as well.

Calculations and model of disk galaxy
The disk galaxy producing the background field is modeled as an axisymmetric homogeneous disk of surface brightness I(R, z) decreasing exponentially with pro- jected radius R and altitude z according to the characteristic radial scale h R and scale height h z , respectively. The field lines stem from the galaxy center at angle φ with respect to the z = 0 plane. We consider only field lines emerging from the galaxy center for simplicity. This approximation is justified by the large baryonic matter density near z = 0, and its fast decrease with radius; see also discussion in Section . The distortion of the field lines will be computed numerically. A small angular deviation δφ of a field line passing near a point mass M is approximately given by where h is the impact parameter and c the speed of light. For our calculation, M is the mass enclosed in a ring of thickness ∆r and height 2h. In addition, the rings at altitudes −(2j +1)h (j ∈ N) contribute to deflecting the field lines toward the z = 0 plane, while rings at +(2j + 1)h deflect the field lines away from the z = 0 plane. Nevertheless, the dominant bending comes from the rings with mid-planes at z = 0, henceforth referred to as "central rings". Fig. 2 shows a sketch of the central rings and the bending of a field line. (For simplicity, the formulae shown below are derived with a mass density decreasing exponentially with deprojected radius r. However, the final results are obtained with a more realistic mass density derived by using the Abel integral equation applied to I(R, z).) The mass M j k of a ring of radius r ≡ k∆r (k ∈ N) and of horizontal mid-height plane at altitude z = jh is: with M * tot the total mass of the galaxy (the * label signals that only baryonic mass is considered), and M k the mass of the central rings (mostly responsible for the lensing): with a ± ≡ (2k ± 1)∆r/(2h R ). This creates an angular bending for a given ring of where complete elliptical integral of the first kind. The results described in this article are for pure disk galaxies i.e. Hubble types 5 or 6. Using the method for earlier Hubble types, with prominent bulges, requires modifying the mass density distribution but does not fundamentally change the method. An example of a calculation is visualized with raytracing in Fig. 3. The figure elucidates two important facts: (1) Although the force is larger at small radii, the difference between the field lines in the right panel and those in the left panel increases with r. Hence, even if the force is strongest at small r, the consequence of self-interaction is more evident at larger r. This explains the larger "missing mass" inferred at large r where fields are weak and, seemingly paradoxically, local self-interaction is not important (δφ → 0). For example, the field lines for r 4 are straight, i.e. δφ ≈ 0, but also roughly parallel to each other in a radial plane and for small z. This implies a force with a r-dependence differing from 1/r 2 . This difference would generate a large missing/dark mass in a Newtonian analysis. In contrast to this straightforward consequence and to empirical observations, dark matter halo models naturally predict large densities of dark matter at small r, which is known as the core-cusp problem [10]. Furthermore, this explains why the gravity modifications of MOND [11] are enabled below a small characteristic acceleration.
(2) The field lines at small r display an approximately isotropic distribution, leading to the familiar 1/r 2 -dependence of gravity. At larger r and small z, indicated in Fig. 3 by a yellow band, the field lines, while still axisymmetric around z, tend to become parallel within radial planes. This leads to a ∼ 1/r behavior of the force and a logarithmic potential ln(r) within a disk of height about 0.01 in the example of Fig 3. For this thin disk, the rotation curves are flat. The field lines outside the yellow region show that gravity still acts outside the galaxy, albeit with depleted strength. This contrasts with the calculations of Refs. [2,3] in which the field needed to be set to zero outside the system. The gravitational force in the radial direction, F (r), is computed as the flux F of the gravity field through a small surface at radius r. We verified that in the case without self-interaction (δφ ≡ 0), the expected 1/r 2 dependence of the force is recovered. For an idealized galaxy with homogeneous mass distribution, the twodimensional flux in a horizontal (z=constant) plane, F h , must be unaffected by lensing because due to axisymmetry the horizontal bending from masses on one side of the field line is compensated by the bending from the masses on the other side. After verifying numerically that indeed F h ∝ 1/r, one can greatly simplify (and accelerate) the numerical calculation by generating field lines only in a vertical plane (x, y = 0, z) and recording the flux F v through a small vertical segment. Then, the force is given by The mass density is highest at small z and this is also where the effect of self-interaction is most important and consequently where the gravitational force is enhanced. Hence, the small-z parts of the disk should rotate faster than its elements at higher z. The attraction between the stars at same r but different z and gas friction, both neglected in our calculation, reduce this z-dependence of the rotation speed at fixed r. Hence we focus our calculations on a small-z disk, as it will be representative of the rotation curves. Another reason for concentrating on small z is that rotation curves are generally obtained from the gas Doppler shift, and the gas scale height is several times smaller than the stellar scale height.

Galaxy Model and systematic studies
The physical parameters of the galaxy model and their chosen nominal values are: (1) the galaxy total mass M * tot = 3 × 10 11 M . This is a typical baryonic mass of a large spiral galaxy; (2) the scale height h z . The stellar scale height of disk galaxies typically obeys h * z ≈ 0.1h R , while the gas scale height h gas z is smaller: h gas z ≈ (0.02-0.04)h R . For the parameter h z , we use a typical value of the gas scale height, h z ≈ 0.03h R , for two reasons: first, the high-density gas produces the stars. Hence, the initial star distribution has a smaller h * z than the one observed in mature galaxies, as demonstrated by the smaller h * z of young blue stars compared to the h * z of older stars. Although h * z increases with time due to dispersion from two-body interactions, angular momentum conservation and the smallness of the vertical speed v z acquired from two-body interaction, v z v θ , impose that the rotation speed of the stars remains nearly constant. The rotation speed is thus determined by the force and density well within a small-z disk before star dispersion occurs. The second reason is the drag of the lighter higher-z layers mentioned in Section ; (3) the radial scale h R . The present radial scale for spiral galaxies is typically h R ≈ (0.5-5) kpc. However, h R increases with time: in the Milky Way, for instance, h R for old stars is about twice smaller than that of young stars [12]. For the same reason as for h R , the initial (smaller) radial scale must be chosen. We thus take h R = 1.5 kpc as a typical initial value.
We varied these parameters within the ranges 10 11 M ≤ M * tot ≤ 10 12 M , 0.01 kpc ≤ h z ≤ 0.5 kpc and 0.5 kpc ≤ h R ≤ 5 kpc. The effect of self-interaction mostly depends on M * tot , h z and only to a lesser extent on h R . Furthermore, in actual galaxies, the smooth decrease of the density with r and h is an average dependence. Measurements of HI gas density distribution show that the actual gas density varies by as much as several times the average. It can be assumed that the stellar density fluctuates similarly. Since averaging is a linear operation, average quantities may not be adequate inputs for non-linear systems. Thus, to investigate the effect of density fluctuations, we randomly varied each M j k by using a Gaussian distribution centered on the value that M j k would have assuming an exponential decrease, and of various widths, with tails truncated symmetrically so that M j k > 0. We determined that the effect is negligible, including the dependence on the width.

RESULTS
The result of the calculation for the nominal galactic parameters is shown in Fig. 4 along with best fits to the fluxes using the form F v = 1/r + α(r/r 0 ) n and F ≡ F v F h = F v /r, with r 0 ≡ 1 kpc setting the scale for the fit parameter α, and n = 1, 0 or -1. For field lines perfectly parallel in a vertical plane, n = 0. That n > 0 signifies that the field lines tend to converge at large r.
Since the force depends significantly on z, we averaged it over h z . Incidentally, this average is equivalent to a calculation done for z center = 0.2 h z , or to reducing M * tot by 19%.

Rotation curves
Once field self-interaction is effectively included in the calculation of F (r), the Newtonian kinematical formalism can be used to obtain the rotation curves. The gravitational force can be parameterized as F (r) = Gm 1 m 2 (1/r 2 + nαr n−1 /r n 0 ) with α and n determined from the fit to the flux, see Fig. 4. This and the radial equilibrium condition for a body in circular orbit at distance r yield the rotation speed v θ (r) = The dotted blue line shows the fit Fv = 1/r + 0.0042r/r0, with r0 ≡ 1 kpc. Bottom: total flux F(r) (red line). Its best fit is F = 1/r 2 + 0.0042/r0 (dotted blue line). The dashed black line is the Newtonian 1/r 2 expectation where field selfinteractions are not accounted for.
GM (r)[1/r + α(r/r 0 ) n ]. Here, M (r) is the mass enclosed within deprojected radius r and determined by using the Abel equation applied to I(R) ∝ e −R/h R , with R the projected radius. The resulting v θ (r) is shown in the top panel of Fig. 5 (red line) and displays the plateau that is typically observed at large r.
We note that to derive v θ (r), the force is approximated as if M (r) were concentrated at r = 0. This is exactly true only for spherically symmetric matter distributions and F ∝ 1/r 2 (Newton's shell theorem). However, the sharp density peaking near the galaxy center makes this approximation acceptable. We checked that for a disk galaxy with an exponentially decreasing density and for F ∝ 1/r 2 , applying Newton's shell theorem yields comparable results to an exact calculation for r h R . In any case, this is irrelevant in practice since most galaxies have an approximately spherical bulge and a high density nucleus. Moreover, the observations of interest, the rotation curve plateaus, occur at radii larger than several h R . For the same reasons, one may approximate that all field lines originate from r = 0 rather than along r = 0.

Effective dark mass profile
The effective dark mass profile can be obtained straightforwardly: in a dark matter framework F =  Fig. 4. The red line is the rotation curve accounting for General Relativity's field self-interaction. The black line is without self-interaction (Newtonian case). The dotted line is the quadratic difference of the two. This represents the missing/dark mass that would need to be introduced in an analysis employing Newton's law -rather than General Relativityin order to recover the plateau seen in the observed rotation curves. Bottom: mass contained within r. Each curve corresponds to one of the rotation curves of the top panel analyzed using Newton's law. The red and dotted lines include field self-interaction and are thus effective masses in which the effect of self-interaction is effectively expressed as an equivalent mass. The black line provides the actual galactic mass.
GM m/r 2 and v θ (r) = [G(M * (r)+M dark (r))/r] 1 /2 , where M * (r) and M dark (r) are the baryonic and dark masses enclosed within r, respectively. The equivalent dark mass profile is then M dark (r) = rv 2 θ (r)/G − M * (r). This is shown by the dotted line in the bottom panel of Fig. 5. Its rise with r illustrates the discussion in Section that the consequence of self-interaction is more evident at larger r despite negligible local self-interactions there.

PREDICTIONS AND VERIFICATIONS
The model shows that the effect of field self-interaction depends mainly on the total baryonic mass M * tot and the scale height h z . This suggests that the dark mass M dark inferred from a galaxy rotation curve correlates with the galaxy's h z . We check this prediction by using two different sets of data [14,15] that provide M d and disk characteristic scales. The first set, from Sofue [14], provides only the radial scale length h R . However, the approximate proportional relation h z = h R allows us to perform the analysis (a correlation analysis is independent of the value of ). The second set, from Martinson et al. [15], The squares represent the data from Ref. [14] and triangles from [15]. M dark /M * disk correlates with hz. The dashed and dotted lines give the best fit to Refs. [14] and [15], respectively, using a ah b z form. The plain line is the lensing calculation described in this article.
provides both h z and M d . Fig. 6 shows M dark /M * disk vs. h z for both sets (set [15] was recombined in 200 pc bins of h z for clarity). We normalized M dark by the disk baryonic mass M * disk to cancel the expected dependence of GR's non-linearity with the baryonic mass. The predicted correlation is clearly visible. Fitting it with M dark /M * disk = ah b z yields the best fit parameters b = −1.48 ± 0.11 for set [14] and b = −1.25 ± 0.14 for set [15]. Although the galaxies constituting the two sets are distinct, their two values of b agree well and are clearly non-zero, indicating a correlation. We can also quantify the degree of correlation by calculating the Pearson linear correlation coefficient between ln(M dark /M * disk ) and ln h z . It is ρ = −0.70 and ρ = −0.69 for sets [14] and [15] respectively, again indicating clear correlations between h z and M dark /M * disk . Normalizing M dark to the total baryonic mass M * disk+bulge yields similar result but with lower values of b: −1.10 ± 0.13 and −1.5 ± 0.30 for Refs. [14] and [15], respectively. This is because no causal connection is expected between h z and M * bulge . Also shown in Fig. 6 are the results of our lensing approach calculated for M * tot = 0.9 × 10 11 M (the approximate average of the baryonic masses in set [14]), h R = 1.5 kpc, and a range of values for the thin disk/gas vertical scale 12 kpc ≤ h gas z ≤ 200 kpc. In Refs. [14] and [15], the inferred dark mass is estimated at R 200 which ranges from tens to hundreds of kpc. Hence, we computed our effective dark mass for a radius of 100 kpc. We rescaled our h gas z so that its range matches that of Refs. [14] and [15]. The result, shown by the plain line in Fig. 6, agrees well with the experimental data, with a χ 2 similar to that of the best ah b z fit.

SUMMARY
The values of the galactic masses and characteristic distances suggest that field self-interaction, a feature of general relativity (GR), needs to be included in studies of galaxy dynamics [2,16,17]. Field self-interaction increases gravity's strength compared to the Newtonian expectation. The effect will become noticeable in systems with large enough masses. In a Newtonian analysis of such systems, the effect would be misinterpreted as a missing mass (dark matter). Furthermore, the similarities between the GR's Lagrangian [3] and observations linked to dark matter and dark energy [18] on the one hand and on the other the QCD's Lagrangian and hadron structure phenomenology, offer another compelling reason to investigate GR self-interaction as an explanation for the universe's dark content. Finally, the exclusion by direct searches of most of the natural phase space for WIMP and axion candidates, and the emergence of dark energy from GR self-interactions [18], make this explanation of the missing mass problem more plausible.
In this article, we presented a new approach to compute the effects of GR's self-interaction based on a meanfield technique and the formalism of gravitational lensing. The findings that self-interaction effects are important for galaxy dynamics and lead to flat rotation curves agree with that of Refs. [2,3] in which a different method (numerical lattice calculation of path integrals) was used to account for field self-interaction. The present method is faster, applicable to any mass distribution, and its formalism is more familiar to the field of astronomy than the path integral formalism customary in particle physics. However, the method is less directly based on GR's equations than the path integral approach. The resulting calculations yield flat rotation curves for disk galaxies and a straightforward explanation of why the missing mass discrepancy worsens at large galactic radii. In contrast, dark matter halo models must be tuned to the specific density profile of a given galaxy to produce a flat rotation curve (disk-halo conspiracy), and these models naturally predict larger dark matter density at small radii (core-cusp problem). The present approach predicts a correlation between the missing mass inferred in galaxies and their vertical scale length. We verified this correlation with two separate data sets.