Dark Energy and Extending the Geodesic Equations of Motion: A Spectrum of Galactic Rotation Curves

A recently proposed extension of the geodesic equations of motion, where the worldline traced by a test particle now depends on the scalar curvature, is used to study the formation of galaxies and galactic rotation curves. This extension is applied to the motion of a fluid in a spherical geometry, resulting in a set of evolution equations for the fluid in the nonrelativistic and weak gravity limits. Focusing on the stationary solutions of these equations and choosing a specific class of angular momenta for the fluid in this limit, we show that dynamics under this extension can result in the formation of galaxies with rotational velocity curves (RVC) that are consistent with the Universal Rotation Curve (URC), and through previous work on the URC, the observed rotational velocity profiles of 1100 spiral galaxies. In particular, a spectrum of RVCs can form under this extension, and we find that the two extreme velocity curves predicted by it brackets the ensemble of the URCs constructed from these 1100 velocity profiles. We also find that the asymptotic behavior of the URC is consistent with that of the most probable asymptotic behavior of the RVCs predicted by the extension. A stability analysis of these stationary solutions is also done, and we find them to be stable in the galactic disk, while in the galactic hub they are stable if the period of oscillations of perturbations is longer than $0.91_{\pm0.31}$ to $1.58_{\pm 0.46}$ billion years.


Introduction
With the discovery of dark energy Λ DE = (7.21 +0.83 −0.84 ) × 10 −30 g/cm 3 [1] - [3], comes a universal length scale λ DE = c/ (Λ DE G) 1/2 = 14010 800 820 Mpc, for the universe that allows for extensions of the geodesic equations of motion (GEOM). However, to be physically viable these extensions must overcome high hurdles. As outlined in [7], these hurdles include the following: Ensuring that the equivalence principle is preserved; this principle is one of the underlying principles upon which general relativity is founded, and has been experimentally verified. Requiring that the equations of motion for massless test particles are not affected; all astronomical observations-in particular, those with which the rotational velocity profiles of spiral galaxies are determined-are based on the motion of photons. Demonstrating that the extension is not prevented by attempts at showing the GEOM is the unique consequence of Einstein's field equations [4]- [6]; such proofs limit the structure of possible extensions. Finally, ensuring that effects which could have already been measured in terrestrial experiments, or observed in the motion of bodies in the solar system, are not produced; such extensions would have been automatically ruled out by experiment.
In [7] we proposed an extension of GEOM, called the extended GEOM, that satisfies these conditions. It was constructed using the dimensionless parameter c 2 R/λ DE G, where R is the Ricci scalar, and replacing the mass m of the test particle by mR[c 2 R/Λ DE G] in the Lagrangian for a test particle in general relativity. By doing so we have changed the response of the test particles to the geometry of spacetime; the worldlines of massive test particles now depend on the local scalar curvature of the spacetime. Importantly, Einstein's field equations are not changed, and thus the geometry of spacetime is still determined by them. The degree by which the worldline is changed is determined by R, which is taken to be a non-linear function of c 2 R/Λ DE G with the nonlinearity modulated by a single parameter, the power-law exponent α Λ . This exponent determines the asymptotic behavior of R for large arguments. A strict lower bound, α ΛBound , for α Λ was determined with the range of possible values for α ΛBound established in [7] by requiring that signatures of the GEOM must not have already been seen in terrestrial experiments. With reasonable choices for experimentally measurable parameters, we found that 1.28 ≤ α ΛBound ≤ 1.58 Given the scale of λ DE , it is only at galactic length scales or longer where the impact of the extended GEOM is expected to be seen, and in [8] we applied this extension to the analysis of the motion of bodies at these scales. Using a spherical model for galaxies, we calculated the density profile of a stationary galaxy given the radius r * H = 11.82 ±0.30 kpc of a typical galactic hub and the velocity v * H = 172.1 ±1.6 km/s of a typical rotational velocity curve (RVC) at this radius [8]. This r * H and v * H were determined from the observed motion of stars in 1,393 spiral galaxies [9] - [14]. The density profile for the model galaxy was determined using the extended GEOM and the following model of the RVC of the galaxy, where v H is the asymptotic velocity of the curve. The power-law exponent was set to α Λ = 1.56 ±0.10 using the Hubble length and the density of this model galaxy (the details of this analysis can be found in [8]); this value is within the bounds for α Λ found in [7]. The radius R 200 for this density profile was calculated to be 206 ±53 kpc, in reasonable agreement with observations. Importantly, σ 8 was also calculated, and was found to be 0.73 ±0. 12 , which is within experimental error of both the WMAP value of 0.71 +0.049 −0.048 [3], and the PLANCK value of 0.81 ±0.006 [15].
In [7] and [8] the focus was on using the extended GEOM to determine the properties of a stationary galaxy that has already been formed. However, if the values of R 200 and σ 8 measured are due to the extended GEOM, then the formation of galaxies must be describable, and the possible RVCs for these galaxies predictable, within this extension. Yet, given the drastic difference between galactic length scale (on the order of tens of kiloparsecs) and the length scales at which R 200 and σ 8 are relevant (on the order of a few hundred kiloparsecs and a few megaparsecs, respectively), the results of our previous paper speaks little about the formation of galaxies. Indeed, since a specific velocity curve v ideal (r) was used to begin with, it certainly cannot predict the RVCs of them. The purpose of this paper is address this lack, and to begin fulfilling these expectations. In particular, our goal here is to establish the range of possible asymptotic behaviors of the RVCs that are allowed by the extended GEOM, and to compare these predictions with observations. In [7] we showed that the energy-momentum tensor T µν for a collection of massive particles that can be treated as a fluid with density ρ, and fluid velocity u µ reduces in the nonrelativistic limit to T µν ≈ ρu µ u ν even when elements of the fluid evolve under the extended GEOM. Applying that result here, to a spherically symmetric distribution of particles rotating about a single rotational axis, we obtain a set of evolution equations Evol for the density ρ(t, r); the fluid velocity along the radial direction u r (t, r); the gravitational potential Φ(t, r); and the (specific) angular momentum L(t, r) = rv(t, r), where v(t, r) is the rotational velocity of the fluid about the rotational axis. Importantly, as our extension of the GEOM involves replacing the mass m of a test particle by mR[c 2 R/Λ DE G], and as this replacement is the same for all particles irrespective of its nature (as such the extended GEOM obeys the weak equivalence principle), the extended GEOM-and throught it the Evol -does not differentiate between baryonic and dark matter; the density ρ of the fluid is the total density of matter in the model galaxy. We have shown below that both the mass and the angular momentum of the system are conserved under this evolution. The types of galaxies that can form, and the RVCs that they can have, under the extended GEOM would then be determined by the solution of Evol for some initial distribution of mass and velocities. These evolution equations are extremely nonlinear, however, and it is doubtful that any direct attempt at solving them will yield much of use. We have taken a different approach instead.
If a choice of the initial distribution of mass and velocities results in the formation of a galaxy under the extended GEOM, then the resultant distribution of mass and velocities must result in stationary solutions-denoted by ρ ∞ (r), L ∞ , Φ ∞ , and u r ∞ -of Evol . Focusing further on galaxies where the motion of matter traces out circular orbits and Evol reduces to Eq. (11) of [8], a single second-order, nonlinear, inhomogeneous differential equation for the stationary density ρ ∞ (r) of the galaxy with the inhomogeneous term given by the angular momentum L ∞ (r) of the fluid in this limit. Importantly, solutions ρ ∞ (r) of this differential equation minimizes a stationary action S ∞ . The dependence of the structure of the galaxy on L ∞ (r)-and given that the total angular momentum is conserved, on the initial distribution of angular momentum L(0, r) of the fluid-underscores the important role that angular momentum plays in the formation of galaxies even under the extended GEOM. While it is in principle possible to choose an initial L(0, r), and then use it to determine whether a galaxy can form under the extended GEOM with this choice, and if it can, whether such a galaxy has a RVC that agrees with observations, doing so would mean evolving L(0, r) in time to L ∞ (r) using Evol . This likely is also intractable analytically.
Since a choice of angular momentum for the fluid must be made, we make this choice at the stationary limit instead of at the fluid's initial state. Using the observed properties of galaxies, we focus on a class of stationary angular momentum given by the RVC v ∞ (r) = (1 + p/q)x 2q x 2(q+p) + p/q where x = r/r * H . Here, q and p are parameters that give the asymptotic behavior of v ∞ (r) in the x 1 and x 1 limits, respectively, with the subscript denoting that we are in the stationary limit. With this choice we are able to determine whether galaxies can form under the extended GEOM, and will be able to predict their RVC. The choice itself depends only on four parameters, each of which have good physical interpretation, and each of which can either be determined (for v * H and r * H )-and thus used as inputs in the analysis-through observations, or compared (for q and p) to them. This choice is a natural generalization of v ideal (r) that is also a smooth function of r, a condition that is important for both physical and mathematical reasons. Importantly, with two free parameters and with L ∞ (r) = rv ∞ (r), v ∞ (r) can model a variety of possible angular momenta for the fluid in the stationary limit, and thus has the potential to model a variety of possible angular momentum L(0, r) at the system's initial state. It thereby defines a class of rotational velocity profiles, one for each given q and p, and importantly, the predicted values of the q and p obtained through the extended GEOM can be directly compared to observations.
To determine the values of q and p that will result in a stationary galaxy under the extended GEOM, we make use of S ∞ . Each choice of q and p results in a L ∞ (r; q, p), which in turn results in a solution ρ ∞ (r; q, p) of Evol in the stationary limit. Such a choice for L ∞ (r; q, p) need not, in general, lead to a ρ ∞ (r; q, p) that minimizes S ∞ , however. Thus, not all choices of q and p will result in the formation of a galaxy under the GEOM. To determine the values of q and p that do, we evaluate S ∞ | (ρ∞; L∞) at this ρ ∞ (r; q, p) and L ∞ (r; q, p); the resultant action then depends on the parameters q and p. Minimization of this action with respect to these parameters then gives the values of q and p that, when used in L ∞ (r; q, p), gives the ρ ∞ (r; q, p) that does minimize S ∞ . It is for these values of q and p that the extended GEOM would predict a galaxy can form. (This approach in determining q and p follows a minimization principle that is similar to the least squares and Rayleigh-Ritz variational methods for solving differential equations [16]. Like those methods, the resultant L ∞ (r) and ρ ∞ (r) obtained are an approximation of the solution of Evol in the stationary limit.) If no such q and p's can be found, then this choice of v ∞ (x) for a class of possible rotational velocity profiles of galaxies is too limited. Galaxies with a RVC given by v ∞ (x)-and likely even those approximated by it-cannot be formed under the extended GEOM.
At the end of this analysis, we find that the action S ∞ does not have one local minimum-or even a discrete number of local minima-for a distinct pair of (q, p). Rather, for each choice of q between 0.010 and 0.336 there is a p between 0.348 and 0.480 ±0.02 that minimizes S ∞ . The asymptotic behavior of the RVC in the galactic hub is thus connected with the asymptotic behavior of the RVC outside of it. This dependence between the two parameters is expected. A single galaxy is formed from a single fluid, and during its formation, fluid elements in one region will interact with the fluid elements in other regions of it. To have the structure of the galaxy inside of the galactic hub be independent of the structure of the galaxy outside of it is physically unreasonable.
That the extended GEOM predicts the formation of a variety of galaxies, each with a different density profile, is a result that is certainly consistent with observations. That the predicted RVCs for these galaxies are different is consistent with both observations and the Universal Rotation Curve (URC) proposed by Persic et. al.
In [17] Persic et. al. analyzed a homogeneous sample of 1100 RVCs of spiral galaxies, and found that only one global parameter-the luminosity of the galaxy-determines the profile of the RVC observed. To describe this dependency they proposed the universal rotation curve V U RC , which gives the velocity profile of any galaxy given its luminosity. Salucci et. al. further refined the URC model in [18], and applied it to the RVC of spiral galaxies; this refinement was then applied to dwarf spheroidal galaxies and low surface brightness galaxies in [19] and [20], respectively. One of the main results of [18] is shown in Fig. 4 of that paper, where the authors plotted an ensemble of URCs, each with a different virial mass. To show the similarity between the curves and to compare these curves to V N F W , the RVC obtained from N -body simulations of Lambda cold dark matter [21], all of the curves were rescaled and normalized to agree at the viral radius. We find that the RVCs predicted here by the extended GEOM agrees well with the curves shown in this figure.
The spectrum of RVCs predicted by the extended GEOM ranges from (q, p) = (0.010, 0.048 ±0.020 ) to (0.336, 0.387 ±0.090 ), with the median curve given by (0.172, 0.349 ±0.010 ). When v ∞ (x) is rescaled and normalized to fit the scale used in Fig. 4 of [18], we find that the ensemble of curves from V U RC is bracketed below by the (0.010, 0.048 ±0.020 ) curve and above by the (0.336, 0.387 ±0.090 ) curve; the median curve (0.172, 0.349 ±0.010 ) lies in the middle of the ensemble of URC curves, and is surprisingly close to the V N F W curve. In addition, we find that the most probable asymptotic behavior in the large x limit for the RVCs predicted by the extended GEOM has a p = 0.348, in good agreement with the profile for V N F W , which has an asymptotic power-law exponent of 0.33 + N F W with N F W < 0.1 [18].
While the minimization of S ∞ does show that stationary galaxies with L ∞ (r) can form and does predict the RVCs for these galaxies, this analysis cannot determine whether the galaxies predicted are stable under perturbations. To address this lack, we have also completed a stability analysis of the predicted galaxies by perturbing about stationary solutions of Evol . This results in a second-order, partial differential equation for first-order perturbations of the stationary radial velocity. We find that the region outside of the galactic hub is very rigid; small perturbations in the radial velocity remain small no matter the frequency of the perturbation. Within the galactic hub, on the other hand, we find that when the frequency of the perturbation is smaller than 0.267 ±0.076 to 0.47 ±0. 16 times the maximum angular velocity of the hub (corresponding to a period of 0.91 ±0.31 to 1.58 ±0.46 billion years) then perturbations in the radial velocity remains small. If, however, it is larger than this range of angular velocities then within the galactic hub small perturbations can increase exponentially with radius.
The rest of the paper is organized as follows. In Sec 2 the focus is on the evolution of fluids under the extended GEOM. The spherical model of the fluid used in this paper is presented. Difficulties in applying Evol to the formation of galaxies is pointed out, and an alternative approach using the stationary limit of Evol is proposed. Details of this approach is given in Sec 3, and the important role that the angular momentum plays is shown. A specific form for L ∞ (r) is proposed. Approximate solutions to the stationary limit of Evol are found in Sec 4 using techniques from boundary layer theory. It is then found in Sec 5 that a spectrum of RVCs is formed under the extended GEOM, and the range of this spectrum is determined. Comparisons with the URC are then made. A stability analysis of the stationary solutions is presented in Sec 6, and concluding remarks can be found in Sec 7.

Evolution under the Extended GEOM
In this section we focus on the time evolution of fluids under the extended GEOM, and the use of these evolution equations in determining the formation of galaxies. We begin with a brief review of the extended GEOM as applied to individual test particles. These equations of motion are then applied to the motion of the collection of these particles that form a fluid using the energymomentum tensor for the fluid in the nonrelativistic and weak-gravity limits. A spherical model of a galaxy is then presented, and the Evol is obtained. The difficulties in using Evol to determine the structure of galaxies are pointed out, and an alternate approach using the stationary limit of Evol is proposed.

A Review of the Extended GEOM for Test Particles
As WMAP measured the pressure to energy density ratio for Dark Energy to be −0.967 +0.073 −0.072 [3]-within experimental error of the ratio expected for the cosmological constant-following [8] we identify Dark Energy with the cosmological constant, and required only that Λ DE changes so slowly that it can be considered a constant. Einstein's field equations are then where T µν is the energy-momentum tensor for matter, R µν is the Ricci tensor, Greek indices run from 0 to 3, Latin indices run from 1 to 3, and the signature of g µν is (1, −1, −1, −1). Here, we have followed [22] and taken, The extended GEOM for a test particle with mass m is obtained from the Lagrangian where for this section only x is the four-vector, x µ = (x 0 , x 1 , x 2 , x 3 ). In [7] we argued for where is chosen so that D(0) = 1. Here, α Λ is a constant, and is the only free parameter in the theory. To prevent the effects of the extension from being already seen in terrestrial experiments, we considered in [7] an experiment designed to look for anomalous accelerations through the propagation of sound waves in a gas of He 4 atoms at 4 K. Reasonable choices for experimental parameters then gives the lower bound for α Λ to be between 1.28 (for Λ DE = 10 −32 g/cm 3 ) and 1.58 (for Λ DE = 10 −29 g/cm 3 ). From Eq. (5), the canonical momentum for the particle is leading to the constraint, as expected. As then with the parametizationẋ 2 = c 2 the Euler-Lagrange equation gives It then follows from Eq. (4) that the extended GEOM for point particles is It is important to note that we have not changed Einstein's field equtions, and thus the geometry of spacetime is still given by the solution of Eq. (5). What we have done by replacing m with mR[c 2 R/Λ DE G] in the Lagrangian for a test particle in general relativity is to change the response of the motion of test particles to the geometry of spacetime. As a consequence, the worldline of the test particle is now given by the extended GEOM Eq. (14) and not the geodesic equations of motion.
Finally, in [8] we used Eq. (14) to determine the density profile of a galaxy with the velocity profile v ideal (r) given in the introduction. This was done by splitting the space around the model galaxy into three regions. While the analysis in [8] for the first two regions will change in this paper, the analysis in the third region will not. Importantly, we found that in this third region the density of a galaxy with a RVC given v ideal (r) by decreases exponentially fast at distances greater than r II = χ/(1 + 4 1+α Λ )λ DE from the center of the galaxy; the reader is referred to [8] for the details of this analysis. Since this decrease in density is not seen, the maximum distance between galaxies is 2r II ; setting this equal to the Hubble length gives α Λ = 1.56 ±0.10 .

The Evolution of Fluids under the Extended GEOM
We begin by considering a collection of particles in a region of space that can be described as a fluid. The distribution of such a fluid is given by its density ρ(x), while the four-velocity field for the fluid is given by the velocity field u µ (x). Then from Eq. (14) the four-velocity of each fluid element is given by the solution of the equation of motion, As we are interested in the formation of galaxies we work in the nonrelativistic limit. In particular, in limit ρc 2 3p we showed in [7] that by using the extended GEOM Eq. (14) the energy-momentum tensor for this fluid can be approximated as T µν ≈ ρu µ u ν , and the current density j µ ≡ T µν u ν = ρu µ is conserved: ∇ · j ≈ 0. Next, WMAP and the Supernova Legacy Survey put Ω K = −0.011 ±0.012 , and the spatial curvature is within experimental error of vanishing. The universe is essentially spatially flat. As the timescales and the lengthscales we are interested in are much shorter than cosmological scales, we approximate the scale factor in the Freeman-Lemaitre-Robertson-Walker metric to be a constant. We are therefore working in the weark gravity limit, and can take the metric to be to be g µν = η µν + h µν . Here, η µν is a flat background metric, and h µν is a small perturbation of it that has only the one nonzero component: h 00 = 2Φ/c 2 . We make this choice for g µν even though with the Λ DE term in Eq. (2) the spacetime will be significantly different from Minkowski space at scales comparable to λ DE . At 14010 800 820 Mpc, λ DE is much larger then the length scales we are interested in, however, and taking the background metric to be flat is a good approximation.
Writing the connection as Γ α µν = 0 Γ α µν + H α µν , then is the connection in the absence of matter and thus determined solely by the coordinates used, while the contribution to Γ α µν due to matter is As expected, the non-vanishing components of 0 Γ α µν are 0 Γ k ij , while for H α µν they are H 0 As u 0 ≈ c, both sides of the µ = 0 component of Eq. (15) is of order u i /c, and are negligible in the non-relativistic limit. The spatial components do survive, however, and give while mass conservation ∇ · j = 0 reduces to in the non-relativistic and weak gravity limits. Since Einstein's field equations can be expressed as As is well known, the only nonvanishing contribution to R µν in this limit is R 00 , and Eq. (20) reduces to The last term is small, however, in comparison to 4πρ, and we set it to zero from now on.

Spherical geometry and the Evol
We now focus on spherically symmetric fluid distributions where the fluid rotates about a single rotational axis. Then using the spherical coordinates (r, θ, φ) where the zenith direction lies along the rotational axis, the velocity along the polar direction, u θ (t, r) = 0, vanishes while the density ρ(t, r), the radial velocity u r (t, r), and the rotational (azimuthal) velocity v(t, r) ≡ u φ (t, r) are functions of t and r only. For this fluid Eqs. (18), (19), and (21) reduce to where L(t, r) = rv(t, r) is the angular momentum of the fluid while is the convective derivative. Since R depends on t only implicitly through ρ, But from Eq. (24) we see that Dρ/∂t is of order u, and the last term in Eq. (22) is of order u 2 , and thus will not contribute to our analysis. Similarly, at the length scales that we are dealing with and with our interest being on the structure of galaxies, using the values of r * H , v * H , α Λ , and Λ DE given in the introduction, we find that numerically R ∼ 1 + O(10 −5 ), and we can set R = 1 in Eq. (23). Angular momentum conservation follows from Eq. (23), while Eq. (24) gives mass conservation.
Equations (22)−(25) give the set of evolution equations Evol for the fluid 1 with the solution to Evol denoted by Three out of the four equations that make up Evol are nonlinear, and as such determining a sufficient set of general boundary conditions needed to obtain a G(t, r) is nontrivial. Indeed, this nonlinearity will limit any definitive comments we can make about the existence of G(t, r), or the properties of it. For much of this paper we will be guided instead by physical principles. In particular, we expect on physical grounds that a set of initial conditions with Φ(0, r) given as the solution of Eq. (25), is needed; Evol can then be considered to be the mapping Evol : G 0 (r) → G(t, r). We also require on physical grounds that as r → ∞, the three quantities, ρ(t, r) → 0, u r (t, r) → 0, and v(t, r) → 0, must separately vanish. One purpose of this paper is to determine whether the formation of galaxies with RVCs that agree with observations is possible under the extended GEOM. As the structure of observed galaxies is essentially stationary, one approach to addressing this question would be to choose a G 0 (r), solve Evol to obtain G(t, r), and then see whether this G(t, r) evolves in the t → ∞ limit to a nontrivial, stationary solution of Evol . (It should be noted that not all choices of G 0 (r) need evolve to a stationary solution of Evol , and the limit in Eq. (30) need not exist. Note also that since G(t, r) is a solution to Evol at each t > 0, if this limit exists then G ∞ (r) is a stationary solution of Evol.) This G ∞ (t) would then be the galaxy predicted to form under the extended GEOM for this choice of G 0 (r), and its RVC could be compared to observations. However, while straightforward, there are a number of issues with this approach.
Evol gives the evolution of any initial distribution of mass and velocities in the spherical geometry. That a specific choice of G 0 (r) may not result in a stationary solution of Evol , or if it does, may not predict a galaxy whose RVC agrees with observation, does not mean that the formation of galaxies with the observed RVCs is not possible under the extended GEOM. It may simply be that the wrong G 0 (r) was chosen. On the other hand, knowing which G 0 (r) should be chosen instead is a daunting task. Indeed, given the extreme nonlinearity of Evol , determining the set of G 0 (r) for which galaxies may be formed under the extended GEOM, or proving that such a set is empty (as would be expected if galaxy formation was not possible), is a difficult task. We have instead taken a different approach, one that focuses on the stationary solutions of Evol .
If a G 0 (r) can be chosen that results in a G ∞ (r) with a RVC which is consistent with observations, then such a solution must be a stationary solution of Evol . To determine whether galaxies can form under the extended GEOM with a RVC that agrees with observation, we focus on these stationary solutions. As we show in the next section, stationary solutions of Evol are given by the solution of a nonlinear, ordinary differential equation, and do not explicitly depend on the choice of G 0 (r); evolving this G 0 (r) under the nonlinear evolution equations given by Evol can be avoided. While a stationary solution to Evol need not be stable, a stability analysis of this solution can be attained by analyzing the evolution of time-dependent perturbations about it. Such perturbations naturally linearize Evol , and their evolution is given by linear partial differential equations whose solutions are tractable. This approach of finding stationary solutions of Evol , and then analyzing the stability of these solutions is the one we have taken in this paper.

The Stationary Limit of Evol and Its Perturbation
We turn our attention to the stationary limit of Evol , and the perturbations about it. We begin by denoting the components of the stationary solution by G ∞ (r) = (0, L ∞ (r), ρ ∞ (r), Φ ∞ (r)). As we are interested in galaxies for which the trajectories of stars are nearly circular, we have taken u ∞ (r) = 0 2 . Moreover, we are in the region where 2πρ ∞ /Λ DE 1, and we may further approximate Eq. (7) as It follows that D ∞ (8πρ ∞ /Λ DE ) 1, and thus from Eq. (6), To obtain both the stationary limit of Evol and perturbations about this limit, we perturb the general solution G(t, r) of Evol about G ∞ (r) by taking r)) being the perturbation. Keeping to first order in this perturbation and separating the time-independent terms from the time dependent ones, we obtain from Evol the equations that determine both G ∞ (r) and G 1 (t, r). We begin with G ∞ (r).

Evol in the Stationary Limit
For G ∞ (r), Evol in the nonrelativistic limit reduces to where D ∞ (r) ≡ D (8πρ ∞ (r)/λ DE ). These two equations can be combined into one second-order differential equation by multiplying Eq. (33) by r 2 and taking the derivative with respect to r. Equation (34) is then used to obtain in agreement with [8].
Treating Eq. (35) as a differential equation for ρ ∞ (r) with a given source term L ∞ (r), we find that the solution to Eq. (35) minimizes the time-independent action where we have made use of Eq. (31). The integral is to r II since 2r II is the maximum separation between galaxies.

Evol for G 1 (t, r)
For G 1 (t, r), Evol reduces to in the nonrelativistic limit. As with the stationary equations, these four equations can be reduced to a single second-order differential equation.
Taking the derivative of Eq. (25) with respect to t and making use of Eq. (24), we find that whereṀ (t) is an arbitrary function of time only. Next, taking the derivative of Eq. (37) with respect to t, and making use of Eqs. (38), (39) and (41), we arrive at since ρ(t, r)u r (t, r) ≈ ρ ∞ u r 1 (t, r) to first order. The solution of this differential equation for u r 1 also minimizes an action, for a given L ∞ (r). As the current density along the radial direction j r (t, r) = ρ(t, r)u r (t, r), the mass flux through a sphere Sph(R) with radius R about the center of the galaxy is From Eq. (41), this flux is WhenṀ (t) > 0 there is a flux of mass leaving the center of the galaxy. and thus the mass of the galaxy would be increasing at the rateṀ (t) even in the limit R → 0. On the other hand, whenṀ (t) < 0 there is a flux of mass entering the center of the galaxy, and thus the mass of the galaxy would be decreasing at the rate ofṀ (t). There will thus be an essential singularity at the center of the galaxy that would either inject mass into the galaxy, or remove mass from it. In either case, the total mass of the galaxy would not be conserved iḟ M (t) = 0. As Eq. (24) asserts that mass is in fact conserved, we setṀ (t) = 0.

The Role of Angular Momentum
In the passage from G(t, r) to G ∞ (r), Evol reduces from four equations determining four functions to two equations determining three functions; Evol thus becomes a system of underdetermined differential equations in the stationary limit. This can be seen explicitly in Eq. (35) where ρ ∞ (r) is only determined once L ∞ (r) is known. This underdeterminacy is expected, and is consistent with observations. Suppose instead that Evol reduces in the stationary limit to a system of three differential equations for the three non-zero components of G ∞ (r). This system of differential equations would be complete, and could then be solved without reference to the initial conditions G 0 (r); they would only have to satisfy the same the boundary conditions at r = 0 and r → ∞ that are required of G(t, r). The solutions of these differential equations will include three arbitrary constants that would then be determined by these boundary conditions. It would then follow that the stationary galaxies predicted by the extended GEOM would all be the same irrespective of the choice of initial condition G 0 (r). This certainly is not what is observed.
That the set of differential equations given by Evol is underdetermined means that G ∞ (r) depends indirectly on the choice of initial conditions G 0 (r). A choice of G 0 (r) will, after evolving with Evol , give a L ∞ (r) that will, through the solution of Eq. (35), give ρ ∞ (r) as well. Indeed, the dependence of Eq. (35) on L ∞ (r), and its role as the driving term for determining ρ ∞ (r) underscores the important role that angular momentum plays in the formation of galaxies. However, since the evolution of G 0 (r) to get a L ∞ (r) would require the solution of Evol , and since such a solution would also give ρ ∞ (r) in the first place, one can question the usefulness of focusing on the stationary solutions of Evol and Eq. (35). In the end, it comes down to the choice of angular momentum for the fluid, and when this choice is made.
We can certainly make this choice at t = 0 by choosing a specific G 0 (r), and this choice may then result in a L ∞ (r) in the stationary limit after evolution under Evol . We can also make this choice at the stationary limit by choosing a L ∞ (r) directly, with the expectation that, since the angular momentum of the system is conserved, there exists some choice of G 0 (r) that will give this L ∞ (r) after evolving with Evol . By making the choice at the stationary limit we circumvent the difficulty of solving the nonlinear partial differential equations in Evol . Moreover, with an appropriate choice of L ∞ (r) we will also be able to determine whether it is possible for the extended GEOM to form galaxies with RVCs that agree with observation. This choice of L ∞ can be determined using S ∞ and the fact for given a stationary L ∞ the stationary density ρ ∞ must minimize S ∞ .
The v ∞ (r) considered here in Eq. (1) is a natural generalization of v ideal (r). It has a well-defined asymptotic behavior in both the x → 0 and x → ∞ limits, and the two behaviors are smoothly joined together. It depends on four parameters, r * H , v * H , q, p, and as q and p give the power law behavior of v ∞ (r) in the asymptotic limits x → 0 and x → ∞ respectively, each parameter has a definite physical interpretation. Two of the parameters, r * H and v * H , are set by observations, and are considered fixed. The remaining two parameters, q and p, are considered to be free, and forms a two-dimensional parameter space P. As each chosen q and p will give a different RVC, this v ∞ (r) describes a class of velocity profiles, each similar in form, and, since by construction v ∞ (x = 1) = v * H , all of whom can be compared with one another. We will, in addition, require that v ∞ (x) → 0 as x → 0 and x → ∞; this in turn requires that q > 0 and p > 0. Newtonian gravity, on the other hand, would set q = 1 and p = 1/2; we follow suit and limit q ≤ 1 and p ≤ 1/2. The parameter space is thus restricted to the strip P = {(q, p) : 0 < q ≤ 1, 0 < p ≤ 1/2}. The main focus of this paper is determining the region of P for which galaxies with velocity profiles v ∞ (x) can be formed under the extended GEOM.
We choose L ∞ (r) = rv ∞ (r). Then each q and p will give a specific L ∞ (r; q, p), and through the solution of Eq. (35), a density profile ρ ∞ (r; q, p) for a galaxy in the stationary limit. But while each choice of q and p may ultimately result in a ρ ∞ , such a ρ ∞ need not minimize S ∞ . By evaluating S ∞ at ρ ∞ (r; q, p) and L ∞ (r; q, p), and minimizing the action with respect to the parameters, we are able to determine the values of q and p that do produce a L ∞ (r; q, p) and a ρ ∞ (r; q, p) which minimizes the action. It is for these values of q and p that a stationary galaxy with a RVC given by v ∞ (x) can be formed under the extended GEOM. Importantly, if no such q and p can be found, then such galaxies could not form under the extended GEOM.

Stationary Solutions
We now turn our attention to finding solutions to Eq. (35). This is possible due to the two drastically different length scales in the theory, r * H and λ DE . A straightforward use of them results in a small parameter that allows for a perturbative solution of Eq. (35). This parameter multiplies the highestorder derivative of the differential equation, however, and thus an application of perturbation theory results in a reduction of the order of this differential equation. The perturbation theory is therefore singular, and thus techniques from boundary layer theory (see Ch. 9 of [23]) will have to be applied. However, there are such significant differences between Eq. (47) and its boundary conditions, and the differential equations and boundary conditions analyzed in [23] that the analysis outlined and the terminology used in [23] cannot be directly applied here. Rather, they will instead serve as guidance for our analysis. Indeed, like the boundary layer analysis [23], perturbative solutions for our differential equation must be found in two or more regions of space, and a consistent solution only exists if these regions overlap. In our case, we will find both the leading and the first-order perturbation solutions of Eq. (47) within the galactic hub and within the galactic disk, and then we will determine for which v ∞ (x) these two regions overlap. Unlike the differential equations considered in [23], however, one of the boundary conditions for Eq. (47) is in the limit x → ∞, and thus use of the outer-and inner-terminology used in [23] would be confusing at best, and we do not use it here. More importantly, with r * H and λ DE we have two different length scores, and this will allow for a different approach to finding the inner-limit solution than that described in [23]. We begin at the r * H length scale.

The Solution in the Galactic Hub Region
In this region we make use of the length scale r * H , and take x = r/r * H . Then Eq. (35) becomes where The values for r * H , v * H and λ DE given in the introduction are used to obtain = 6.131 × 10 −3 after evaluating χ at α Λ = 1.56 ±0. 10 . Using as an expansion parameter, we first take the outer limit → 0 [23], and expand Υ(x) to first order in 2 : Υ(x) = Υ 0 (x) + 2 Υ 1 (x). Equation (47) gives for the leading term and like the nonlinear Carrier equation [23] the reduction of order in Eq. (47) results in an algebraic equation that is easily solved. Indeed, the resultant equation in this limit is algebraic to all orders in the perturbation expansion. In particular, it gives for the first-order perturbation The resulting density is Equation (52) is valid for values of x for which 2 |Υ 1 (x)| < |Υ 0 (x)|, or equivalently, when E hub (x) ≡ 2 |Υ 1 (x)/Υ 0 (x)| < 1. This condition establishes the region R hub of space for which Eq. (52) is valid. We will see below that R hub encompasses the region around r = 0, and thus corresponds to the galactic hub. Indeed, this can be seen directly by setting = 0 in Eq. (52); ρ ∞ (r) then reduces to the Newtonian result.

The Solution in the Galactic Disk Region
In this region we make use of the length scale λ DE , and takex ≡ r/(χ 1/2 λ DE ). Then definingȳ(x) ≡ 8πρ ∞ (x)/Λ DE , and Eq. (35) becomesF Forx and expand Eq. (35) to first order inx 2p The leading termȳ a is given by The solution to Eq. (57) was found in [8]; the details of how this was done can be found there. Here, we will only need the followinḡ where For the first-order perturbationȳ 1 , on the other hand, Eq. (35) gives The extent of the region R disk where this perturbative expansion is valid is determined by the condition As we will see below, this region excludes the point r = 0, and thus corresponds to the galactic disk.
Equation (60) is straightforwardly solved to givē where ν = |Σ α Λ +1 | − 1/4, j c (x) = 1 x 5/2 cos(ν log x), and j s (x) = 1 Herex 0 is a point in the intersection of the regions R hub and R disk . In the next section we will determine the region of P where this intersection is nonempty. For now, we will assume that we are working in this region. We then make use of the following expansion to evaluate the integral where the constants are determined by requiring the density to be smooth at the transition point x 0 =x 0 /x H between the regions R hub and R disk , As for the particular solution ρ disk part (x), when x c < x 0 , while when x c x Z n pq +2 (70)

Consistent Solutions
The region where ρ hub ∞ is valid is given by R hub = {x : E hub (x) < 1}; this region is simply connected, and Eq. (35) exists when the two regions overlap, R hub ∩ R disk = ∅ [23]. It is only in this case that a x 0 can be chosen, and the arbitrary constantsC cos andC sin in the homogenous solution to Eq. (35) be determined. The focus of this subsection is on determining both R hub and R disk , along with the region P ∩ ⊂ P for which their intersection is not empty: Considering first the limit x → 0, we find that Then x hub L = 0 for q ≤ α Λ /(α Λ +1) and q = 1, while when α Λ /(α Λ +1) < q < 1, The largest that x hub L can be is 2.69 × 10 −4 or roughly 3.18 pc. In the x → ∞ limit, on the other hand, ρ disk are found numerically using the following process. For x hub U , a value of x trial is chosen, and E hub (x trial ) is calculated. If E hub (x trial ) > 1, the value of x trial is decreased while if E hub (x trial ) < 1, it is increased. With this new value for x trial , E hub (x trial ) is again calculated, and the process is repeated until sufficient accuracy is achieved. This final x trial is then set equal to x hub U . A similar process is used to determine x disk L . Fig. 1 A 3D plot of (∆x)∩ with respect to q and p. Note the region in red where (∆x)∩ < 0.
In this region R out ∩ R in = ∅.
The results of this calculation is shown in Fig. 1. Both x hub U and x disk L were calculated to an accuracy of 0.0001 starting at q = 0.01 and continuing in 0.1 increments from q = 0.1 to q = 1.0. Similarly, p starts at 0.01, and increases in increments of 0.01 until p = 0.50 is reached. We find that (∆x) ∩ > 0 everywhere except for the red triangular-shaped region shown in the figure. This region is bounded by the lines q = 1.0, p = 0.50, and a curve that starts at (1.0, 0.33) and ends at (0.3, 0.50). Outside of this triangular region the regions R hub and R disk overlap, and there is a consistent boundary-layer solution to Eq. (35).
We emphasize that while (∆x) ∩ < 0 in the triangular-shaped region, this does not mean that there are no solutions to Eq. (35) in this region of P. All that can be concluded is that the singular perturbation analysis that divides space into only two regions cannot be applied. A consistent solution may be possible when a third, intermediate region is introduced to interpolate between the two regions, for example. This region would be given by the solution to Eq. (35) obtained by setting the inhomogeneous term equal to the term proportional to 2 . However, such solutions would depend more on the detail behavior of v ∞ (r) in the transition region between the galactic hub and the disk-and thus on how this transition is modeled-than on the asymptotic properties of the RVC. Moreover, we have found values of q and p do that minimize the stationary action, and they lie far outside of the triangular region. As our focus is on the asymptotic behavior of RVCs, the two-region, boundary-layer solution is sufficient for our purposes.
The values of (∆x) ∩ in P ∩ range from 0.0014 to 0.6339, and, given that 6(∆x) ∩ < (x hub U + x disk L )/2, are quite small when compared to either x hub , the size of (∆x) ∩ is such that this choice of x 0 does not have much of an impact on our analysis. Nevertheless, since we will find in the next section that the action is dominated by the behavior of v ∞ (x) in the galactic disk, we choose x 0 = x hub U to maximize the contribution of the galactic hub to the action.

A Spectrum of Rotational Velocity Curves
In the region P ∩ the boundary-layer method gives as the solution of Eq. (35). Such a solution must also minimize the action S ∞ , however. In this section we will determine the values of q and p that do, and in doing so, determine the RVCs that can form under the extended GEOM. We will find that a continuous range of RVCs is possible, and this spectrum of RVCs is in agreement with the URC. When evaluated at L ∞ and ρ ∞ , the action breaks up into two pieces, corresponding to the solutions in the regions R hub and R disk . We begin with the region R hub .

The Action in the Region R disk
While in Sec 4.2 we usedȳ andx, in this section we find it more convenient to use y = (Λ DE /8πρ * H )ȳ and x. Then Expanding Eq. (36) aboutȳ a , and keeping terms linear in y 1 , This action naturally breaks up into two additional pieces, S disk , with the first piece consisting of the first two terms in Eq. (79). When S disk-asym ∞ is evaluated at the solution Eq. (58), they can be integrated to give after Eq. (59) is used. For the second piece consisting of the third and fourth terms in Eq. (79), after an integration by parts and making use of Eq. (60), it reduces to Making use of Eq. (63) again, this last integral becomes for x c < x 0 , while for x c > x 0 ,

Minimization of S ∞ (ρ∞,L∞)
When evaluated at ρ ∞ the action becomes The first two terms depend on q and p indirectly, through x 0 . Since the integrand in the first term is integrable, and because we chose x 0 = x hub U ∼ 1 − 5, they do not contribute appreciably to the action. It is the third term, with its dependence on x II = 2.264 × 10 5 , q, and p, that dominates the behavior of S ∞ (ρ∞,L∞) , and will determine the values of q and p that minimizes it. To emphasize this, and to isolate the dependence of S ∞ (ρ∞,L∞) on these parameters, we define Note that the dependence of S(q, p) on P ∩ is dominated by the x 2 II y 1 (x II ) term in Eq. (81). This x 2 II y 1 (x II ) in turn is dominated by two terms, one coming from the particular solution to Eq. (60), which is proportional to (x c /x II ) Z 0 pq , and the other coming from the homogeneous solution to Eq. (60), which is proportional to (x c /x II ) 1/2 . For S(q, p) to be small, the homogeneous solution must dominate, and thus it is in the region of P ∩ where Z 0 pq ≥ 1/2-which in turn requires p ≥ 1/4-that the minima of S(q, p) will be found. The precise values for q and p that minimizes S(q, p) are determined numerically using the following process. We first determine x hub U to an accuracy of 10 −6 for q starting at 0.01 and continuing in 0.1 increments from 0.1 to 1.0, and for p starting at 0.01 and continuing in 0.01 increments to 0.50. We then set x 0 = x hub U , and calculate S(q, p) for these values of q and p. The result is a two-dimensional surface above P ∩ . We find that there is a slight concavity in the surface in the rectangular region of P ∩ bounded by the lines q = 0.01, q = 0.40, p = 0.34, p = 0.50. For each choice of q in this region there is a p(q) such that ∂S ∂p (q,p(q)) = 0.
The solution to Eq. (86) gives p(q) as a function of q, and thus defines a curve on the base space P ∩ . The lifting of this curve to the surface S(q, p) gives the curve S(q) = (q, p(q), S(q, p(q))) on which the action is local minimum for each choice of q. Fig. 3 Graph of the dependence of p(q) on q along S(q). Notice that for much of the graph p(q) is a weak function of q. Notice also that the uncertainty in p(q) grows rapidly as q → 0.336.
To determine this function p(q) and the curve S(q), we determine x hub U to an accuracy of 10 −6 for q starting at 0.200 and continuing in 0.002 increments until 0.210 is reached, and for p starting at 0.345 and continuing in 0.001 increments until 0.355 is reached. These x hub U are then used to calculate S(q, p), and for each q the p that minimizes S(q, p) is determined along with the value of S(q, p) at this point. Up to N = 15000 terms in the series in Eqs. (68) − (70) and (82) − (83) is used in calculating S(q, p). As expected, the collection of these points form a curve on the action surface S(q, p). We then follow this curve along values of q that are less than 0.200 and along values of q that are greater than 0.210 in increments of 0.002 until we reach a q ∈ P ∩ for which a minimum of the action cannot be found. The result of this calculation is shown in Figs. 2 and 3. Figure 2a is a graph of the minimum action curve S(q) above a region of the base parameter space P ∩ ; not shown is the surface on which this curve lies. The projection of this curve onto the S(q, p) − q plane is shown in Fig. 2b, while the projection of the curve onto the S(q, p) − p plane is given in Fig. 2c. Notice the asymptote for the curve shown in Fig. 2b and 2c. Figure 3 shows the graph of p(q) versus q, and is the projection of S(q) onto the q − p plane. Each point (q, p(q)) on the curve gives a L ∞ (x) that results in a density ρ ∞ that minimizes S ∞ . Thus, each point (q, p(q)) on this curve gives the density and, through v ∞ (x), the RVC of a galaxy that can form under the extended GEOM. Notice that p(q) is nearly flat for most of the values of q shown, and thus many of these galaxies will have RVCs that have similar asymptotic behavior in the galactic disk, while at the same time have very different behavior in the galactic hub. This can be seen explicitly in Fig. 4.
By focusing only on the large x asymptotic behavior of the RVC, we use the increments 0.001 by which p was increased when determining S(q) as a bin size ∆p, and determine the probability of finding a galaxy with an asymptotic power-law exponent between p and p+∆p. This is done by simply counting the number of RVCs with a value for p between p and p+∆p without regard to their corresponding values of q; the resultant histogram is shown in Fig. 4. Notice that the most probable value of p-with 19 out of the 164 possible RVCs, or 11.6%-is 0.348. The median of this distribution of curves is at p = 0.349 ±0.01 ; the corresponding RVC has a (q, p) = (0.172, 0.349 ±0.01 ). In addition, 50% of the possible RVCs have a p ≤ 0.355 while 95% of the curves have a p ≤ 0.404.
Shown also in Figs. 2 and 3 are the estimated uncertainties in determining S(q, p(q)) and p(q). Notice in particular the large increase in uncertainty in p(q) as q → 0.336; this is precisely the location of the asymptote for S(q). The largest contribution to the uncertainties is due to the power-law exponent α Λ . While the uncertainties in r * H , v * H , and λ DE also contribute to the uncertainties in S(q, p(q)) and p(q), including these contributions to the uncertainties in any reliable manner would require determining x hub U to an accuracy much higher than 10 −6 ; we did not do so here. Instead, we focus on the uncertainty due to α Λ by increasing the value of α Λ to α Λ + ∆α Λ with ∆α Λ = 0.01. A new curve S(q) was then calculated, and the graph p(q) determined. The uncertainty in S(q, p(q)) was calculated from the difference in S(q, p(q)) due to this change in α Λ ; the uncertainty in p(q) was calculated in a similar way. While large, this is the smallest ∆α Λ that could be used without increasing the accuracy in x hub U to beyond 10 −6 . For these reasons we caution that the uncertainty shown in Figs. 2 and 3 is an estimate.
After analyzing a homogeneous sample of 1100 RVCs, Persic et. al. [17] found that the profile of the RVC for a galaxy is determined by a single parameter, the luminosity of the galaxy. They further showed that these profiles could be described by a single function of this luminosity. This work was further refined by Salucci et. al. [18] where they expressed the square of the URC as the sum of two terms, This V U RCD gives the stellar contribution to the URC, while V U RCH gives the dark matter component. To demonstrate the self-similarity of the URC, and to compare the URC to V N F W , an ensemble of URCs, each with a different virial mass M vir , was plotted in Fig. 4 of [18]. This was accomplished by rescaling x opt = r/R opt → x vir = r/R vir , where R opt and R vir are the optical and virial radii, respectively, and normalizing both the ensemble of URCs and the V N F W so that all the curves agree at x vir = 1. The similarity between these curves becomes readily apparent.
By rescaling x → x scale = x/8 so that the maximum of v ∞ (x) now occurs near the location of the maxima of the URCs in Fig. 4 of [18], and rescaling v ∞ so that v ∞ (x scale = 1) = 1, we have added the spectrum of RVCs predicted by the extended GEOM to this graph. The result is shown in Fig. 5. The two ends of the p(q) graph in Fig. 3 correspond to (0.010, 0.480 ±0.020 ) and (0.336, 0.387 ±0.090 ); all the predicted RVCs are bracketed below by the (0.010, 0.480 ±0.020 ) curve and above by the (0.336, 0.387 ±0.090 ) curve. These two curves, shown in red in Fig. 5, are superimposed on Fig. 4 of [18] along with the median RVC curve given by (0.172, 0.349 ±0.010 ); this median curve is shown in blue in Fig. 5. The two extreme RVCs, (0.010, 0.480 ±0.020 ) and (0.336, 0.387 ±0.090 ), also bracket the ensemble of URCs from [18]. While the (0.336, 0.387 ±0.090 ) curve lies significantly above the highest URC curve shown, the uncertainty p for this curve is both very large and is the highest of the extended GEOM RVCs. The median RVC curve (0.172, 0.349 ±0.010 ) of the extended GEOM lies also in the middle of the ensemble of URC graphed in Fig. 4 of [18]. While the two extreme curves from the extended GEOM does not approach the V N F W RVC as closely as the URC curves, in the region 2 ≤ x ≤ 18 Salucci et. al. have shown that the Burkert and NFW profiles can be approximated as Fig. 5 Graphs of the extended GEOM RVC superimposed on Fig. 4 of [18]. The two graphs in red is given below by the (0.010, 0.048 ±0.020 ) curve, and above by the (0.336.0.387 ±0.090 ) curve. They bracket the ensemble of UHCs (grey-scale lines) and V N F W (solid black line) from [18]. Also plotted is the median (0.172, 0.349 ±0.010 ) curve in blue.
where N F W < 0.1 3 , and for x near 18, V U RCH (x) ∼ x −0.33− N F W . The histogram in Fig. 4 shows that the most probable asymptotic behavior of a RVC predicted by the extended GEOM has a p = 0.348. Such a RVC would have the asymptotic behavior v ∞ ∼ x −0.348 , in good agreement with the V N F W , and the ensemble of URC curves. While the extreme curve (0.010, 0.480 ±0.020 ) does not have a p that is within N F W of 0.33, 95% of the predicted curves have a p ≤ 0.404, and is within N F W of 0.33.

Stability Analysis
We now turn our attention to the stability analysis of the stationary solutions found in the previous section. The equation determining the perturbation u r 1 (t, r) of the radial velocity in the stationary limit was found in Sec 3.2. Instead of working with u r 1 (t, r), however, we work with the radial current density j r (t, r) = ρ(t, r)u r (t, r) ≈ ρ ∞ (r)u r 1 (t, r) since u ∞ (r) = 0, and through it, the flux of mass through a sphere Sph(r) of radius r, Equation (42) then becomes Using the WKB approximation to order , we find that where If ε 2 h < 0, then M 2 0 (x) > 0, and µ(t, x) will be an exponential function of t, and an oscillatory function of x. The situation is more complicated if ε 2 h > 0. While µ(t, x) will always be an oscillatory function of t, it will be an exponential function of x when M 2 0 (x) < 0, and an oscillatory function of x when M 2 0 (x) > 0. We are interested in the case when µ(t, x) is an oscillatory function of both t and x, and therefore bounded. Given that in R hub we have x < x hub U , we can always choose a ε h for which this is true. The function v 2 ∞ (x)/3x 2 is monotonically decreasing when q ≤ 1. Moreover, when q < 1, v 2 ∞ (x)/3x 2 → ∞ as x → 0. As x < x 0 for all x ∈ R hub , Thus, to ensure that This in turn imposes an upper limit to the frequency, of oscillation for the mode. As ξ(x) is then assured to be an oscillatory function, we define as the local wavenumber for the oscillation with λ(x) = 2π/k(x) being its corresponding wavelength. When x → 0, When, however, x → x − 0 and ε h → ε max h , we may approximate Here, This P (z) ∼ v * H 2 /c 2 , and is small compared to the terms proportional to z 2 and ν 2 terms in Eq. (105). We thus treat the P (z)ξ term as a perturbation, and solve Eq. (105) perturbatively by taking ξ = ξ 0 + ξ 1 . Then for z = 2(1 + 3α Λ )ε d z, These are Bessel's equations of imaginary orderν = (1 + α Λ )ν. Using the same terminology and notation in [24], we find Here,z 0 ≡ 2(1 + 3α Λ )ε d /ȳ a (x 0 ) 1/2 , and G(s,z) is the Green's function, In the limit x → ∞,z → ∞, and F iν (z) ∼ cos(z − π/4)/ √z while G iν (z) ∼ cos(z − π/4)/ √z ; ξ 0 thus dies off as 1/ √z . For ξ 1 (z), we first note that E hub ∼ y 1 /ȳ a ∼z 2ȳ 1 . Then since P (z) ∼z 2 E disk . From Sec 4.2, y 1 (z) consists of a particular solution y P and a homogenous solutionȳ h . The particular solution behaves as y p (z) ∼ 1/z 2(1+α Λ )(1+p) for largez, and as suchz 3ȳ p ∼ 1/z 2(1+αΛ)(1+p)−3 . For the integral Eq. (110) to converge at largez, 2(1 + α Λ )(1 + p) − 3 > 0 or p > (1 − 2α Λ )/2(1 + α Λ ). This is always true on P. Next, for the homogenous solution, y h (z) ∼ 1/z 5(1+α Λ )/2 , so thatz 3ȳ h (z) ∼ 1/z (5α Λ −1)/3 . The contribution bȳ y h (z) to the integral also converges. Thus, ξ(z) is well-behaved everywhere.

The Stability of Stationary Solutions
From Sec 6.2 we see that in R disk small perturbations in the current flux µ will remain small, and the stationary solutions to Evol found in Sec 5.3 are thus stable. The situation is more complicated in R hub , however.
We see from the analysis in Sec 6.1 that µ be an oscillatory-and thus bounded-function of both t and r when ε 2 h ≥ 0 and 0 ≤ ε h < ε max h . This only occurs when the frequency of oscillations of the perturbation ω < ω * H v ∞ (x 0 )/x 0 . For the stationary solutions found in Sec 5.3, 0.267 ±0.076 ≤ v ∞ (x 0 )/x 0 ≤ 0.47 ±0. 16 . These stationary solutions are therefore stable in R hub as long as the period of oscillations for the perturbations is longer than T max with 0.91 ±0.31 ≤ T max ≤ 1.58 ±0.46 billion years. Such perturbations have a maximum wavelength of ∼ 1.5r * H to ∼ 7r * H .

Concluding Remarks
The choice of v ∞ (x); the direct connection between the parameters used in its construction and observations; and the ability of v ∞ (x) to model a wide range of rotational velocity profiles, have allowed for those profiles that are consistent with the extended GEOM to be determined. Indeed, while each point in P ∩ corresponds to a different RVC, and while each RVC may give a L ∞ (r) that results in a solution ρ ∞ (r) of the stationary Evol , it is only along the curve (q, p(q)) in P ∩ for which a stationary solution that minimizes the action can be found. As each ρ ∞ (r) obtained from the stationary Evol would correspond a galaxy with a RVC given by v ∞ (x), it is therefore only galaxies with RVCs given along this curve that will be formed under the extended GEOM. This spectrum of allowed RVCs is consistent with the URC, and given that the URC is constructed through the observations of the velocity profiles of 1100 spiral galaxies, it is consistent with observations as well. In fact, the two extreme RVCs predicted by the extended GEOM bracket the ensemble of V U RC shown in [18], while the median curve has a form similar to both the V U RC and V N F W . Moreover, the asymptotic behavior of URCs and V N F W is in good agreement with that of the RVC with the most probable p predicted by the extended GEOM. Importantly, we have also shown that these stationary solutions in the galactic disk are stable under perturbations, while in the galactic hub they are stable as long as the period of oscillations of the perturbation is longer than 0.91 ±0.31 to 1.58 ±0.46 billion years; these perturbations have wavelengths shorter than ∼ 1.5r * H to ∼ 7r * H . When comparing the graphs of the RVC obtained using the extended GEOM with the URC in Fig. 5, it becomes readily clear that v ∞ (x) may be too simplistic in the transition region between the two asymptotic limits x → 0 and x → ∞. This is borne out by the shallowness of the minima in S(q, p); we would expect that the more accurate the choice of v ∞ (x) is, the deeper the minima will be. There are, in fact, many ways of smoothly joining together the asymptotic behavior of the rotational velocity profile in the two limits x → 0 and x → ∞. Given this, it is likely that future progress in determining the RVCs that can form under the extended GEOM using the stationary-solution approach presented here will come from making a different choice in v ∞ (x) as much as, or even more than, from more accurate numerical calculations.
We have used a spherical model for our galaxy with the fluid rotating about a single axis. Moreover, the values for the parameters r * H and v * H used here were obtained through observations of the motion of stars in spiral galaxies. As such, the results we have obtained can be most directly applied to the formation of spiral galaxies. It applicability to the formation of other types of galaxies, such as those analysed in [19] and [20], is still an open question, and is a topic of future research.
The extension of the GEOM we have considered here replaces the mass of a test particle m by mRc 2 R/Λ DE G in the Lagrangian for a test particle in general relativity. By doing so we have changed the response of the motion of test particles to the geometry of spacetime; the worldline of the test particles is now determined by the extended GEOM, and not the GEOM. Einstein's field equations are not changed, and the geometry of spacetime is still determined by the solution of them. Importantly, this approach does not differentiate between baryons and dark matter, and does not change the worldlines of massless particles. It is for these reasons that we were able to show in [7] that the extended GEOM is not excluded by the deflection of electromagnetic waves by the Sun, or through the advancement of the perihelion of Mercury.
Another approach to modifying gravity, called modified gravity theories in general, takes a different approach. The focus of these approaches is to change general relativity itself. Examples of such theories include the Jordan-Brans-Dicke theory where the gravitational constant is replaced by a scalar field, a scalar-tensor theory where the cosmological constant is replaced by a scalar field, and f (r) theories where the Ricci scalar in the Hilbert action is replaced by a function f (R) of it (see [25] for a review). Such modifications of general relativity inherently changes Einstein's field equations, and as such the resulting geometry of spacetime. Importantly, the response of test particles to this geometry is not changed, and the worldline of the particle is still determined by the GEOM. In particular, both the worldlines of massive and massless particles are effected, and as such Solar system tests of general relativity place stringent limitations on such theories and the introduction of screening mechanisms are needed (see [26] for a review and [27] for an application of the screening).
While the choice of L ∞ (r), and the construction of v ∞ (x) was driven by observations and the requirement that all the parameters used in v ∞ (x) have a definite physical interpretation, they are nevertheless choices. They were made with the expectation that there are choices of initial conditions G 0 (r) which, when evolved to the stationary limit by Evol , will indeed result in both the L ∞ (r) chosen and the stationary solution ρ ∞ (r) resulting from this choice. Whether such expectation is born out is a question that requires, if not solving the Evol , then at the least a perturbative analysis of it near the stationary limit. This analysis is also a focus of future research.

Declarations
Funding and Conflicts of Interest: The author did not receive support from any organization for the submitted work, nor does the author have any financial or proprietary interests in any material discussed in this article.