Significance of Gravitational Nonlinearities on the Dynamics of Disk Galaxies

The discrepancy between the visible mass in galaxies or galaxy clusters, and that inferred from their dynamics is well known. The prevailing solution to this problem is dark matter. Here we show that a different approach, one that conforms to both the current Standard Model of Particle Physics and General Relativity, explains the recently observed tight correlation between the galactic baryonic mass and its observed acceleration. Using direct calculations based on General Relativity's Lagrangian, and parameter-free galactic models, we show that the nonlinear effects of General Relativity make baryonic matter alone sufficient to explain this observation.


I. INTRODUCTION
An empirical tight relation between baryonic and observed accelerations in galaxies has been reported by Mc-Gaugh et al. [1, hereafter MLS2016]; larger accelerations are accounted for by the baryonic matter, i.e. there is no missing mass problem, while in lower acceleration regions, dark matter or gravitation/dynamical laws beyond Newton's are necessary. This correlation is surprising because galactic dynamics should be dictated by the total mass (believed to be predominantly dark), but instead the baryonic mass information alone is sufficient to get the total acceleration. A tight connection between dark matter and baryonic distributions would explain the observation, but such connection has not been expected. While the relation from MLS2016-including its small scatter-can be reproduced with dark matter models, [e.g. 2], the consistency of the correlation width with the observational uncertainties suggests a dynamical origin rather than an outcome of galaxy formation, since this would add an extra component to the width.
Dynamical studies of galaxies typically use Newton's gravity. There are however arguments that once galactic masses are considered, relativistic effects arising from large masses (rather than large velocities) may become important [3][4][5]. Their physical origin is that in General Relativity (GR), gravity fields self-interact. In this article, we explore whether these effects can explain the relation in MLS2016 without requiring dark matter or modifying gravitation as we currently know it. GR's selfinteraction can already account for observations suggestive of dark matter or dark energy: the universe's expansion [5], the rotation curves of disk galaxies [3,4], galaxy cluster dynamics [3]-including the Bullet cluster observation [6], and the correlation between the missing mass of elliptical galaxies and their intrinsic ellipticities [7]. It also naturally leads to the Tully-Fisher relation [8].
The article is organized as follows. We first outline the self-interaction effects in GR in Section II. Then we discuss in Section III the empirical tight dependence of baryonic mass and observed acceleration in disk (i.e. lenticular and spiral) galaxies. In Section IV, we use GR's equations to compute such correlation. These CPUintensive calculations allow us to study only a few galaxies, modeled as bulge-less disks. To cover the full range of disk galaxy morphologies, including those with significant bulge, we develop in Section V two dynamical models of disk galaxies in different, complementary ways: uniform sampling (Section V A) and random sampling (Section V B) of the galactic parameter space. In Section VI we show the results from these models, and compare them to observations. Finally, in Section VII, we summarize our findings and their importance.

II. SELF-INTERACTION EFFECTS IN GENERAL RELATIVITY
Field self-interaction induces the non-linearities of GR. The phenomenon is neglected when Newton's law of gravity is used, as typically done in dynamical studies of galaxies or galaxy clusters. However, such a phenomenon becomes significant once the masses involved are large enough. Furthermore, it is not suppressed by low velocity-unlike some of the more familiar relativistic effects-as revealed by the inspection of the post-Newtonian equations [9]. In fact, the same phenomenon exists for the strong nuclear interaction and is especially prominent for slow-moving quark systems (heavy hadrons). It produces the well-known quark confining linear potential.
The connection between self-interaction and nonlinearities is seen e.g. by expanding the Einstein-Hilbert field Lagrangian of GR into a polynomial form: (1) where g µν is the metric, R µν the Ricci tensor, M the system mass and G is the gravitational constant. In the natural units ( = c = 1) used throughout this article, [G] = energy −2 . The polynomial is obtained by expanding g µν around a constant metric η µν of choice, with ϕ µν ≡ g µν − η µν the gravitational field [see e.g. 10]. The brackets are shorthands for sums over Lorentz-invariant terms [4]. Field self-interaction originates from the n > 0 terms in Eq. (1), distinguishing GR from Newton's theory, for which Lagrangian is given by the n = 0 term. One consequence of the n > 0 terms is that they effectively increase gravity's strength. It is thus reasonable to investigate whether they may help solve the missing mass problem. In fact, it was shown that they allow us to quantitatively reproduce the rotation curves of galaxies and dynamics of clusters without need for dark matter, also providing a natural explanation for the flatness of the rotation curves [3].
The framework used in these studies is analogous to that of Quantum Chromodynamics (QCD, the gauge theory of the strong interaction). GR and QCD Lagrangians are similar and both contain field self-interaction terms. In QCD, their effect is well-known as they are magnified by the large QCD coupling, typically α s 0.1 at the transition between perturbative and strong regimes [11]. In GR, self-interaction becomes important for non-negligible values of the coupling GM/L (L is the system's characteristic scale, used here to form the dimensionless coupling needed to heuristically assess self-interaction's importance), typically for GM/L 10 −3 [4]. In QCD, a major effect of self-interaction is a stronger binding of quarks, resulting in their confinement. In GR, self-interaction increases gravity's strength, which can explain the missing mass problem [3].
A key feature for this article is the suppression of self-interaction effects in isotropic and homogeneous systems [3]: • In a two-point system, large GM/L or α s values lead to a constant force between the two points (and none elsewhere), i.e. the string-like flux-tube, well-known in QCD.
• For a homogeneous disk, because of the symmetry, the flux collapses only outside the disk plane. The resulting force between the disk center and a point in the disk at distance r decreases as 1/r. This automatically produces flat rotation curves if the disk density decreases fast enough, as it does for disk galaxies.
• For a homogeneous sphere, the force recovers its usual 1/r 2 behavior since the flux has no particular direction or plane of collapse.
This symmetry dependence has led to predict a correlation between the missing mass of elliptical galaxies and their ellipticity [7]. It also offers a simple explanation for the relation reported in MLS2016: since disk galaxies contain a central high-density bulge that is usually nearly spherical [12], self-interaction effects are suppressed there by the near-spherical symmetry. Hence, in the region of higher accelerations, the 1/r 2 behavior of gravity assumed in galaxy dynamic analyses is actually valid. In the disk itself, the absence of spherical symmetry enables field self-interaction and thus, the dynamics in lower acceleration region is non-linear. If the self-interaction is not accounted for and Newton's law is instead used to analyze the data, a mass discrepancy arises.

III. BARYONIC MASS-ACCELERATION DEPENDENCE
The correlation between the radial acceleration traced by rotation curves (g obs ) and that predicted by the observed distribution of baryons (g bar ) reported in MLS2016 was established after analyzing 2693 points in 153 disk galaxies with varying morphologies, masses, sizes, and gas fractions. A good functional form fitting the correlation is: where g † is an acceleration scale, the only parameter of the fit. In the remainder of the article, we show that the observed correlation may entirely be due to the nonlinear GR effects which are neglected in the traditional Newtonian analysis. In the next section, we use the results of direct calculations of rotation curves of actual galaxies [3]-modeled as bulge-less disks-to show that when the galactic bulge of the actual galaxy is negligible, the calculation yields a relation that agrees with the empirical correlation from MLS2016. In the two subsequent sections, we develop models to include the effect of bulges and to account for the variation of morphology of disk galaxies.

IV. DIRECT CALCULATIONS
The rotation curves of several disk galaxies were computed in [3] based on Eq. (1) and using numerical lattice calculations in the static limit [4]. These calculations neglect the galactic bulge and approximate a spiral galaxy with a disk featuring an exponentially-falling density profile. They were carried out for nearly bulge-less Hubble types 5 and 6 galaxies (NGC 2403, 3198 and 6503), and for Hubble types 3 and 4 galaxies (NGC 2841, 2903 and 7331), which have moderate bulges. Using these results, we can compute the total acceleration g SI stemming from baryonic matter and including GR's field selfinteraction-analog of g obs from MLS2016. Plotting it versus the Newtonian acceleration g N obtained from the same distribution of baryonic matter, but ignoring GR's self-interaction-analog of g bar from MLS2016-one obtains the results shown in the top panel of Fig. 1. The curves for types 5 and 6 galaxies agree well with the observed correlation, thereby providing an explanation for it in bulge-less galaxies. However, the curves for types 3 and 4 galaxies, while qualitatively following the correlation, overestimate g SI and lie on the edge of the observed distribution. That the empirical correlation is reproduced only for bulge-less galaxies supports that 1) the correlation from MLS2016 is explainable by GR's self-interaction without requiring dark matter or modification of the known laws of Nature, and 2) at large acceleration, i.e. typically for small galactic radii, the bulge reduces the value of g SI since self-interaction effects cancel for isotropically distributed matter.
Although based directly on the GR's L, the lattice approach is limited since it is computationally costlyrestricting the study to a few galaxies-and it applies only to simple geometry-limiting the studied galaxies to late Hubble types. To study the correlation from MLS2016 over the wide range of disk galaxy morphologies, we developed two models based on: 1) the 1/r gravitational force resulting from solving Eq. (1) for a disk of axisymmetrically distributed matter; and 2) the expectation that GR field self-interaction effects cancel for spherically symmetric distributions, such as that of a bulge, restoring the familiar 1/r 2 force. These dynamical models are described in the next section.

V. DYNAMICAL MODELS
To circumvent the limitations of the direct lattice calculation, we constructed two elementary models for disk galaxies. They both compute the acceleration including self-interaction, g SI , and the Newtonian acceleration due to the baryonic matter, g N . If self-interaction is responsible for the correlation between g obs and g bar shown in MLS2016, then the same correlation should exist between g SI and g N . Both g SI and g N are computed at a set of radii r, from the galactic center to its outermost parts. This is carried out for all of the galaxies-with their parameters sampling observed correlations reported in literatureselected by the two models.
The modeled galaxies have two components: a spheroidal bulge and a larger disk. Both contain only baryonic matter following the light distribution, i.e. there is no dark matter and gas is either neglected or follow the stellar distribution.
The bulge is modeled with the projected surface brightness Sérsic profile [13] used in [12]: I b (R) = I e 10 −bn[(R/Re) 1/n −1] , where R is the projected radius, I e is the surface brightness at the half-light radius R e , n is the Sérsic parameter and b n ≈ 0.868n − 0.142 [14]. The internal mass density ρ b (r), where r is the deprojected radius, is computed from the surface brightness by numerically solving the Abel integral. Since GR's self-interaction effects cancel for isotropic homogeneous distributions, the potential in the bulge has the usual Newtonian form, Φ b (r) = GM enc b (r)/r, where M enc b (r) is the mass of the bulge enclosed within a sphere of the deprojected radius r.
The disk is modeled with the usual radial profile I d (R) = I 0 e −R/h , where I 0 is the central surface brightness, h is the disk scale length, and possible effects from the disk thickness are neglected. Again, the corresponding mass density ρ d (r) is computed from the Abel integral. Self-interaction in a homogeneous disk leads to a potential Φ d (r) = G M enc d (r) ln(r), with M enc d (r) the disk mass enclosed within a radius r, and G the effective coupling of gravity in two dimensions, which depends on the physical characteristics of the disk [4]. Rather than being computed ab-initio as in Section IV, which would demand a separate lattice calculation for each galaxy, G is obtained by requiring the continuity of the acceleration at the transition r t between the two components.
The parameters characterizing a galaxy-the bulge and disk masses, M b , M d (from which ρ b,0 and ρ d,0 are obtained, respectively), R e , n and h-span their observed ranges for S0 to Sd galaxies [12,15,16]. There are known relations between these parameters [12,16]: We randomly sample R e and M b from the ranges of observed values to obtain the remaining galactic parameters-h, n, M d -through Eqs. (3)-(5). Thus, there are no adjustable parameters in our models. We stress that the accuracy of the empirical relation Eqs. (3)- (5) is not critical to this work, their purpose being only to provide reasonable values of the galactic parameter space we select. While the simplicity of our models would make it of limited interest for investigating the intricate peculiarities of galaxies, such simplicity is beneficial for the present study: no numerous parameters nor phenomena (e.g. baryonic feedback) are needed for adjustment to reproduce the correlation from MLS2016. That the correlation emerges directly from basic models underlines its basic nature.
The two dynamical models introduced in the remainder of this section share the above description. From here, they differ in two aspects. The first is in how the observed correlations in Eqs.  [12,17]. The correlations are applied strictly, i.e. without accounting for the scatter seen in actual data, since systematically spanning the observed typical ranges for the parameters contributes to the width of the correlation reported in MLS2016, and accounting for such scatter would partly double-count and thus overestimate the width.
Inside the spherical bulge-dominated region (r < r t , where r t is the bulge-to-disk transition radius, defined later in the text; denoted by subscript r < r t ), the selfinteraction cancels, and the two gravitational accelerations are the same: where M enc d (r) = 2π r 0r ρ d (r)dr.
In the disk-dominated region (r > r t ; denoted by subscript r > r t ), numerical lattice calculations indicate that self-interaction leads to a collapse in the gravitational field lines [3,4]. The bulge density there is less significant than that of the disk, but is still present. The total acceleration is: The Newtonian acceleration g N,r>rt in this region has the same form as in the bulge, given in Eq. (6). G is determined by requiring the accelerations to match at r t : g SI,r<rt (r t ) = g SI,r>rt (r t ). Thus, G = G/r t , by construction.
The accelerations in the bulge and disk regions are smoothly connected using a Fermi-Dirac function centered at r t = 2R e and of width r t /2: F (r) = 1/ 1 + e 2(r−rt)/rt . Therefore, the acceleration with selfinteraction is: g SI (r) = F (r)g SI,r<rt (r) + (1 − F (r)) g SI,r>rt (r), (10) while the Newtonian acceleration is: The dependence of the model result on the choice of width value is minimal: abruptly transitioning between bulge and disk, i.e. using a step-function rather than F (r), yields quantitatively similar results. The small dependence on the functional form for the transition is supported by the agreement between Models 1 and 2 which use different methods for the transition, as we discuss next.

B. Model 2: Random Sampling of the Galactic Parameter Space
For Model 2, in order to systematically sample the galaxy parameter space, we simulated many galaxies whose randomly sampled parameters simultaneously satisfy Eqs. (3)-(5). A candidate galaxy is generated by first randomly sampling a distribution of R e , and then using it to randomly sample the observed correlations in Eqs. (3)-(4) to obtain h and n. These are then combined with Eq. (5) to find ρ b,0 and ρ d,0 , thereby completing the parameter set for a single candidate galaxy. The candidate galaxy is then accepted if its parameters satisfy all of the correlations to within one standard deviation. An example of a set of sampled galaxies is shown in Fig. 2.  [16]. The best χ 2 fit to the observed data is denoted with a solid line, and one dex in dashed lines. Red circles denote generated galaxies.
The acceleration in the bulge and disk regions are connected with a step-function H(x) = 1 for x < 0, and 0 otherwise, such that at r t , the acceleration is kept continuous by the proper choice of G . Since actual transitions between bulges and disks are blurred, and to mitigate a possible bias in choosing a specific value of r t , r t is randomly sampled from a normal distribution centered at 2R e , with standard deviation of R e and symmetrically truncated tails. Models 1 and 2 use different methods for the bulgedisk transition. The agreement between the two models suggests that they are indifferent to a particular choice of transition method. The acceleration including selfinteraction is g SI (r) = H(r − r t )g SI,r<rt (r) + (1 − H(r − r t )) g SI,r>rt (r), (12) where g SI,r<rt (r) and g SI,r>rt (r) are as in Model 1, given in Eqs. (6)-(9), and g N (r) is as in Eq. (11).
In Model 2, we also modeled the effect of the bulge being spheroidal rather than spherical by introducing a polar dependence: ρ b (r, φ) = ρ b (r)(1 − cos 2 φ), with ρ b (r) the spherical bulge density used in Models 1 and 2. This refinement did not noticeably change the results, thereby further proving their robustness.

VI. RESULTS
Direct lattice calculation and the two dynamical models allow us to compute the accelerations for a set of galaxies whose parameters follow the typical observed ranges for disk galaxies. The acceleration including nonlinear self-interaction (g SI ) is plotted in Fig. 1 versus the acceleration computed with the same baryonic mass distribution but assuming Newtonian gravity (g N ). This is compared to the observed correlation between g obs and g bar reported in MLS2016. The top panel shows the results for the direct calculation, the middle panel the results of the Model 1, and bottom panel for Model 2. Our computed correlations agree well with the empirical observation, without invoking dark matter or new laws of gravity/dynamics. We formed the residual distribution between the result of Model 2 and the best fit to the observational data. It is shown in the embedded figure in the bottom panel of Fig. 1. (Since Model 1 samples the galactic phase space uniformly rather than using normal distributions, its residual has little meaning and we show only that of Model 2.) The residual distribution needs not to be exactly gaussian since it is formed using the best fit to a different distribution than that of Model 2, the best fit form being adapted to this different distribution. However, the residual distribution is approximately gaussian, with a width of about 0.09 dex, comparable but sightly smaller than the 0.11 dex observed in MLS2016. This is expected since our models employ idealized mass profiles that do not incorporate substructures such as arms in the disk. These density variations create fluctuations in a given rotation curve that increase the width of the observed residual.

VII. DISCUSSION AND CONCLUSION
Our findings support the possibility that GR's selfinteraction effects increase the gravitational force in large, non-isotropic mass distributions. When applied to disk galaxies, the increased force on the observed matter transposes to the missing mass needed in the traditional Newtonian analyses. We have thus proposed a plausible explanation for the correlation between the luminous mass in galaxies and their observed gravitational acceleration shown in MLS2016. That this correlation is encapsulated in our basic, parameter-free, models indicates its fundamental origin.
The explanation proposed here is natural in the sense that it is a consequence of the fundamental equations of GR and of the characteristic magnitudes of the galactic gravitational fields, and in the sense that no fine tuning is necessary. This contrasts with the dark matter approach that necessitates both yet unknown particles and a fine tuning in galaxy evolution and baryon-dark matter feedbacks [see e.g. 2]. We used several approaches that are quite different, thus leading to a robust conclusion.
The work presented here adds to a set of studies that provide straightforward and natural explanations for the dynamical observations suggestive of dark matter and dark energy, but without requiring them nor modifying the known laws of Nature. This includes flat rotation curves of galaxies [3], the Bullet cluster [3,6], galaxy cluster dynamics [3], and the evolution of the universe [5]. The Tully-Fisher relation [8] also finds an immediate explanation [3]. There are compelling parallels between those observations and QCD phenomenology, e.g. the equivalence between galaxies' Tully-Fisher relation, and hadrons' Regge trajectories [3,4], plausibly due to the similarity between GR's and QCD's underlying fundamental equations. The fact that these phenomena are well-known for other aspects of Nature that possess a similar basic formalism; the current absence of natural and compelling theory for the origin of dark matter (supersymmetry being now essentially ruled out); and the yet unsuccessful direct detection of a dark matter candidate or its production in accelerators despite coverage of the phase-space expected for its characteristics; all support the approach we present here as a credible solution to the missing mass problem.