The Case Against Dark Matter and Modified Gravity: Flat Rotation Curves Are a Rigorous Requirement in Rotating Self-Gravitating Newtonian Gaseous Disks

By solving analytically the various types of Lane-Emden equations with rotation, we have discovered two new coupled fundamental properties of rotating, self-gravitating, gaseous disks in equilibrium: Isothermal disks must, on average, exhibit strict power-law density profiles in radius $x$ on their equatorial planes of the form $A x^{k-1}$, where $A$ and $k-1$ are the integration constants, and"flat"rotation curves precisely such as those observed in spiral galaxy disks. Polytropic disks must, on average, exhibit strict density profiles of the form $\left[\ln(A x^k)\right]^n$, where $n$ is the polytropic index, and"flat"rotation curves described by square roots of upper incomplete gamma functions. By"on average,"we mean that, irrespective of the chosen boundary conditions, the actual profiles must oscillate around and remain close to the strict mean profiles of the analytic singular equilibrium solutions. We call such singular solutions the"intrinsic"solutions of the differential equations because they are demanded by the second-order equations themselves with no regard to the Cauchy problem. The results are directly applicable to gaseous galaxy disks that have long been known to be isothermal and to protoplanetary disks during the extended isothermal and adiabatic phases of their evolution. In galactic gas dynamics, they have the potential to resolve the dark matter--modified gravity controversy in a sweeping manner, as they render both of these hypotheses unnecessary. In protoplanetary disk research, they provide observers with powerful new probing tool, as they predict a clear and simple connection between the radial density profiles and the rotation curves of self-gravitating disks in their very early (pre-Class 0 and perhaps the youngest Class Young Stellar Objects) phases of evolution.


Introduction
A large number of observations, mostly in the 21 cm emission line of neutral hydrogen, have firmly established that the rotation curves of spiral galaxy disks do not exhibit a Keplerian falloff, in fact, most of them remain flat or slightly increasing as far away from the centers as they can be observed (Freeman 1970;Rogstad & Shostak 1972;Roberts & Rots 1973;Bosma 1978;Rubin 1980;Rubin et al. 1980;Bosma 1981a,b;Rubin et al. 1982;Bahcall & Casertano 1985;van Albada & Sancisi 1986;Kent 1987;Begeman 1987;Persic & Salucci 1988;Begeman 1989;Persic & Salucci 1990;Carignan et al. 1990;Casertano & van Gorkom 1991;Broeils 1992;Persic & Salucci 1995;Persic et al. 1996;Salucci & Persic 1997). The radial scales over which the neutral hydrogen disks can be observed reach out to ∼100 kpc in the largest spiral galaxies and since the rotation curves remain flat, it was postulated by many researchers that some unseen extended mass distribution ought to exist all the way out to hundreds of kiloparsecs from the galaxy centers. Thus the Dark Matter Hypothesis was born, and soon the news leaked out to the rest of the physics community and intensive and extensive searches for dark matter particles and fields boomed into existence. On the other hand, some researchers who certainly felt uncomfortable with this new "aetherial" hypothesis proposed that Newtonian gravity should instead be modified at galactic scales and beyond, in order to solve the problem of the fast rotation of HI galaxy disks (Milgrom 1983a,b,c;Tohline 1983;Felten 1984;Sanders 1984Sanders , 1986Mannheim & Kazanas 1989;Christodoulou 1991;Mannheim & O'Brien 2011, 2012. The gas in spiral galaxies is distributed in centrally concentrated, vertically thin disks. For this reason, it was expected that the rotation curves had to turn over at some intermediate radius and begin a decline that would be indicative of the absence of substantial amounts of matter at large radii. This view about the luminous matter is nowadays considered so settled and clear that it has made its way into introductory Astronomy textbooks that compare and contrast the kinematics of spiral galaxies to the kinematics observed in our solar system (the Keplerian falloff mentioned above). In this work, we show that this elementary perception is quite naive and totally wrong because it ignores the influence (in fact, the dominance) of pressure and enthalpy gradients in self-gravitating gaseous disks. Personally, we believe that it amounts to a blunder because, before the results described below and even back in the 1980s when all this was unfolding, we had important clues that pointed out the importance of gas pressure in determining the equilibrium structures of gaseous disks. For instance, the sound-crossing time at 10 kpc in a 10 4 K cold galaxy disk is 1 Gyr, a value that lies to within 1.5-3 of the rotation timescales at 10 kpc in all galaxies with rotation speeds of 100-200 km s −1 . Gallagher et al. (1984) noted that the star formation histories of a sample of irregular and spiral galaxies did not indicate the presence of dark matter in the low-mass galaxies of the sample. But the most obvious clue was the existence of the Mestel (1963) disk, a centrally concentrated potential-density pair with a flat rotation curve (see also Jalali & Abolghasemi 2002). The Mestel disk has always been considered just a toy model (Binney & Tremaine 1987) and the situation did not improve when Schulz (2012) showed that the finite Mestel disk requires significantly less mass to produce the flat rotation curve. The main argument against the universal adoption of the Mestel disk has been the absence of a physical law or reason that would make galaxy disks assume this specific surface density profile and rotation law.
The above argument is not based on solid reasoning. When nature shows us that she has widely adopted a specific property (the flat rotation curves in galaxy disks), Aristotelian Logic dictates that we should search for a new law or reason, in order to understand the universality of this property and establish its physical meaning; not to create ghosts (particles and fields), aethers, and new forces that effectively facilitate our aversion to confronting the facts.
We do not claim that the Mestel (1963) disk is the answer to establishing the universality of flat rotation curves in galaxy disks; only that it has always been a telling clue that gravity does not pull the strings and is not in control in gaseous self-gravitating disks. Furthermore, we have solved the full Newtonian problem and we now know precisely how such universal rotation curves emerge in spiral galaxy disks. The resolution of this ubiquitous problem is the subject of this paper. Before we can delve into the physics of the problem, we need to correct some common misconceptions that appear in the theory of second-order differential equations and which also have made their way into the textbooks. We do so in § 2. Then, in § 3, we revisit the theory of rotating Newtonian isothermal gaseous-disk equilibrium models and we calculate analytically the mean shapes of their density profiles and their rotation curves.
The results match precisely the shapes of the rotation curves of spiral galaxy disks with no additional assumptions of any kind. So these results make a strong case against both dark matter and modified gravity and their implications have far-reaching consequences all the way to cosmology. For completeness, we describe in § 4 polytropic models that also demand monotonically increasing rotation curves because they are subject to the same physical principles. These models are also applicable to very young protoplanetary disks (certainly to pre-Class 0 disks and possibly to the youngest Class 0 non-Keplerian disks). Finally, we conclude with a discussion of all the pertinent issues and our results in § 5.

Second-Order Differential Equations and the Cauchy Problem
In mathematical physics, the trivial solutions of the various second-order differential equations are commonly ignored as being uninteresting; and too much attention is paid to the Cauchy problem in determining arbitrary constants as opposed to the internal properties of the equations themselves that have no regard for externally imposed conditions of any type. Both of these practices are damaging as they work to hinder our efforts toward solving the physical problems described by the differential equations in the first place. Such practices are relevant to all linear and nonlinear second-order equations of physics, so we can discuss and clarify the various issues involved by using any well-known equation. We choose to make use of the Bessel differential equation in this section.
The Bessel equation (Watson 1922) has regular solutions that are called Bessel functions of order m and a trivial solution y = 0. The trivial solution is not singular as it can be obtained from the regular solutions by an appropriate choice of the arbitrary constants. Nevertheless, its name indicates that y = 0 is of no interest at all. It is also well-known that the Bessel functions all oscillate about the x-axis (Watson 1922), but this statement is grossly inaccurate and obscures the truth: the regular solutions oscillate about the trivial solution which just happens to coincide with the x-axis in this case. We demonstrate this important point by solving numerically an inhomogeneous m = 0 Bessel equation of the form along with nonsingular boundary conditions that attempt to initially push the regular solutions away from the new trivial solutions y = K: y(1) = 1, y ′ (1) = −10 for K = 5; and y(2) = 1, y ′ (2) = 10 for K = −5. The results are shown in Fig. 1. Both regular solutions have nothing to do with the x-axis; instead, they turn around and settle into oscillations that clearly occur about the new trivial solutions y = K. This behavior can be demonstrated for all linear second-order equations of mathematical physics with oscillatory regular solutions (Christodoulou et al. 2015) and for (non)linear equations of the Lane-Emden type (Robe 1968;Schmitz & Ebert 1986;Schneider & Schmitz 1995;Christodoulou & Kazanas 2007). The lesson to be learned is that the so-called trivial solutions of second-order equations are not at all trivial. They are in fact favored by the differential equations themselves which have no regard for the externally imposed boundary conditions. Thus, we will heretofore call these solutions the intrinsic solutions that are preferred and demanded by the equations themselves, irrespective of the externally imposed Cauchy problem.
When the Cauchy problem is solved, as in Fig. 1, the externally imposed boundary conditions are usually at odds with the underlying equation and the regular solutions cannot match the favored intrinsic solution. As a result, the regular solutions are forced by the equation itself to oscillate about the intrinsic solution as soon as they intersect this favored solution the first time. Thus, the intrinsic solutions act as attractors of the regular solutions which, in turn, are forced to always stay near and around the more dominant intrinsic solutions. We view this behavior as a triumph of the differential equation (and its intrinsically favored solution) over the Cauchy problem (and the particular solution it strives to produce).
This striking behavior remains intact in at least some nonlinear second-order equations.
Very clear examples in which rotation is involved can be found in Schneider & Schmitz (1995) and Christodoulou & Kazanas (2007). Here we provide two additional nonlinear examples of the dominance of intrinsic solutions over regular solutions drawn from Lane-Emden equations (Lane 1869-70;Emden 1907) in the absence of rotation. The isothermal Lane-Emden equations take the form where D is the dimensionality of space (D = 3 in spherical coordinates and D = 2 in cylindrical coordinates). Singular and regular solutions have been obtained in many applications (Stodólkiewicz 1963;Ostriker 1964;Schmitz & Ebert 1986;Binney & Tremaine 1987;Christodoulou & Kazanas 2007) and they are all nonoscillatory. 1 The reason for this is quite obvious: eq. (3) does not have an intrinsic solution because e y = 0. Why the latter condition precludes an intrinsic solution will become clear in the next section, where we describe a procedure for obtaining intrinsic solutions.
In stark contrast, the polytropic Lane-Emden equations take the form where n > 0 is the polytropic index, and it possesses the intrinsic solution y = 0. Although few analytic solutions are known (Chandrasekhar 1939), numerical integrations show that, depending on n, this equation has both oscillatory and nonoscillatory solutions. For eq. (4), we have derived a precise criterion for the existence of oscillatory solutions and we will describe it in a future publication. This criterion predicts that for D = 2 (cylindrical form), all solutions with odd integer n-values are oscillatory; while for D = 3 (spherical form), only the n = 1 and n = 3 integer-n solutions are oscillatory. Numerical integrations (using the physical boundary conditions y(0) = 1, y ′ (0) = 0) easily confirm these results (see Figs. 2 and 3). The reason for the existence of nonoscillatory solutions is that for the corresponding choices of n, the differential equation is not a harmonic oscillator (Christodoulou et al. 2015). This is also true for the modified Bessel equation (Watson 1922) that is known to possess only nonoscillatory solutions. Its real solutions cannot be oscillatory 2 and they are then prohibited from intersecting the intrinsic solution more than once (Christodoulou et al. 2015).

Isothermal Self-Gravitating Newtonian Gaseous Disks
In what follows, we use the arbitrary scaling constants R o and ρ o to normalize the disk radius R and density ρ(R), respectively. We thus define the dimensionless radius x ≡ R/R o and density τ (x) ≡ ρ(R)/ρ o . Velocities V (R) are also normalized consistently by the constant where G is the Newtonian gravitational constant, in which case we also define the dimensionless rotation velocity v(x) ≡ V (R)/V o . The same scaling also applies to the sound speed C o of the gas which in this section is a constant, i.e., the dimensionless sound speed is The cylindrical isothermal Lane-Emden equation (Lane 1869-70;Emden 1907) with rotation can then be written in dimensionless form as This equation describes the radial (x) equilibrium of a rotating, self-gravitating, gaseous disk or cylinder in which the gas obeys the isothermal equation of state is the dimensionless pressure of the gas. Eq. (5) is valid exactly for infinite cylinders and to a high degree of approximation in the equatorial (symmetry) planes of disks (see the Appendix). This latter point has been demonstrated convincingly by the calculations of Hayashi et al. (1982), Schmitz & Ebert (1986), and Narita et al. (1990). In particular, the last two investigations of finite disks uncovered equatorial density profiles that were strictly oscillatory under proper boundary conditions, just as was predicted by the analysis of § 2 above. Hayashi et al. (1982) and Schmitz (1986) studied also the stability of such equilibria and found that, except for the very flattened disks and the nearly spherical configurations, the intermediate models are stable. The very flattened disks of Hayashi et al. (1982), in particular, were unstable to ring formation that causes their equatorial power-law density profiles to become oscillatory, in agreement with the numerical solutions of eq. (5) obtained by Christodoulou & Kazanas (2007).
Despite the exact analytic results of the researchers quoted above, an objection has been raised over the years concerning the validity of using cylindrical coordinates to study axisymmetric, vertically thin disks rather than just infinite cylinders; and this is also the major sticking point in the present work, so it needs to be addressed in detail. We defer the analysis of the Lane-Emden equations for vertically thin disks to an Appendix because it is much easier to follow the derivations after the intrinsic analytic solution has been obtained for cylinders as follows.
Christodoulou & Kazanas (2007) described a procedure for obtaining the intrinsic solution of eq. (5): If we equate the last two terms: then this is an intrinsic solution provided that the rest of the equation (the radial variation of the logarithmic gradient of the enthalpy) vanishes: Eqs. (6) and (7) form a system in which v(x) is totally dependent on τ (x). First we solve eq. (7) to obtain the radial density profile: and then we solve eq. (6) to determine the rotation curve of the intrinsic solution: where The solution contains 3 free parameters, the integration constants A, B, and k. Parameter B sets the vertical scale of the rotation curve v(x), so we can choose B = 1 in what follows. The density profile τ (x) is a power law with index k − 1, while k can be thought as the power-law index of the surface density profile σ(x) ∝ x k (Christodoulou & Kazanas 2007).
Figs. 4-6 show the shapes of the rotation curves obtained from eqs. (9) and (10) for various choices of the constants A and k. The results are scale-invariant (a property of power-law density profiles), so the radial scale is arbitrary (while, on the other hand, B = 1 scales v). It is not surprising (see § 3.1 below) that, in these Newtonian models, one sees most of the shapes of the "flat" rotation curves observed in spiral galaxies.
The only shapes missing from the figures are those of the falling rotation curves (Casertano & van Gorko 1991). Apparently, in some compact galaxies, the above equilibrium profiles did not endure (the sound-crossing time at 30 kpc for C o = 10 km s −1 is 3 Gyr, so it seems that there has been enough time to achieve equilibrium); perhaps because of interactions with nearby galaxies; or because the "external" gravity of the massive bulges eliminates the intrinsic solution (see § 3.1 and footnote 4 below). But even in these objects, the falloff is slower than Keplerian which implies that gas pressure is still fighting to establish its own preferred profiles. This must be the case since the intrinsic solutions are favored by the equilibrium differential equation itself (see § 2) in the absence of other external forces. In the two galaxies observed by Casertano & van Gorkom (1991), certain segments of the rotation curves are flat and that indicates to us that the radial self-gravitational equilibrium has fallen apart at some locations but not (yet) everywhere in these disks.
The above results can be summarized as follows: The derived density profiles are simple power laws in radius and the rotation curves are flat or slightly increasing at large radii (eqs. [9] and [10]) irrespective of the value of the index k. Spiral galaxies have always been fitted with exponential density profiles, thus we do not know which values of k occur in nature. Galaxy profiles will have to be fitted again, but the payoff this time will be substantial: when k is determined from observations, the large-scale rotation curve (away from the center) will also be obtained independently from the velocity measurements. Thus, the observational results will be tested for consistency within the same data set in each case. This also holds true for protoplanetary disks in their early isothermal or adiabatic (see § 4 below) phases; but first we need to find such purely self-gravitating disks in a pre-Class 0 YSO (Young Stellar Object) stage, and the 21 cm HI line in emission or absorption may give us a chance (Kamp et al. 2008).

Physical Interpretation
But how can such power-law density profiles produce and support flat or slowly increasing rotation curves? The problem has always been that the centrifugal force remains too high in the outer regions of the disk where Newtonian gravity weakens substantially. How do these equilibrium models get around this discrepancy? The answer to these questions was given in Christodoulou & Kazanas (2007) and we repeat it here: The Lane-Emden eq. (5) is a second-order differential equation. As such, it respects but does not rely solely on force balance. (As will be readily seen in eq. [12] below, force balance is guaranteed by the way that the specific enthalpy of the gas is operating.) This second-order equation describes locally the detailed radial (d/dx) variation of the logarithmic gradients (d/d ln x) of the potentials involved in the struggle to reach equilibrium. So, initially, it is the log-gradient of the enthalpy (with help from rotation) that sets out to oppose the log-gradient of the gravitational potential. This competition can be seen by rewriting eq. (5) in the form where h(x) ≡ dp/τ and ψ(x) are the dimensionless enthalpy per unit mass and the Newtonian self-gravitational potential, respectively.
And in that case, how and why does the average intrinsic solution (eqs.
[8] and [9]) come into existence? The answer to these questions is even more illuminating: The intrinsic solution was derived above by imposing two separate conditions, that the centrifugal force should match the gravitational force (see eq. [6] and the last two terms in eq. [11]): while, at the same time, the log-gradient of the enthalpy should retire from the competition by assuming a constant profile (see eq. [7] and the leading term in eq. [11]): This occurs only for a power-law density profile (eq. [8]) and for a specific rotation profile (eqs. [9] and [10]), and these profiles become internal properties characteristic of the equilibrium disk. Stated more simply, up until now people believed that rotation was not related to the structure of the equilibrium disk and that they could adopt any arbitrary rotation profile for gaseous disks. We see now that this is not true, the radial density and rotation profiles are strongly coupled and uniquely determined through the intrinsic solution discussed above.
The only surprise in this narrative is the unique way that the differential equation (or the disk) finds to promote and establish the above intrinsic solution: rather than trying to simultaneously balance the variations of the potentials involved at every single radius (not possible because the corresponding timescales vary widely at different radii), the disk assumes gradually a logarithmic specific-enthalpy profile determined from the solution of eq. (13): So it is the thermodynamic potential h(x) of the gas that becomes logarithmic, and not the gravitational potential that people have been trying to make it so for nearly 50 years! Naturally, the disk establishes such a logarithmic profile because this law guarantees precise force balance (eq. [12]) at all of its equatorial radii. The action of h(x) unfolds in the physical disk from inside-out over timescales of the order of the local sound-crossing time R/C o . So the global equilibrium becomes complete after a time ∼ R max /C o , where R max is the outer radius of the disk.
We understand physically the preceding results in the following manner: Self-gravity is a long-range force (by Gauss's law, the gravitational potential at radius x depends on the entire mass interior to x) and it cannot adjust its potential in the disk to effect local changes to the density distribution. In other words, the density is a source term in the Poisson equation that determines the gravitational potential, but not the other way around. In stark contrast, enthalpy is a local potential whose action depends only on the local behavior of the pressure and the density of the disk. When h(x) assumes its logarithmic radial profile (eq. [14]), it dictates that the local density adjust according to the local log-pressure gradient (τ ∝ dp/d ln x). By doing that, the enthalpy uses the density in order to modify the sourcing of both the gravitational potential (via the Poisson equation ∇ 2 ψ = τ ) and the centrifugal potential (via eq. [6]). The result of this tactic is a rotation law that is entirely dependent on the distribution of enthalpy (see eq. [6]) that does not feed back (v does not enter in eq. [7]). The rotation so produced is capable of balancing gravity at all radii all by itself (eq. [12]), and the enthalpy retires from the struggle for equilibrium (eq. [7]) having implicitly won the competition at every radius.
It is important to keep in mind that the enthalpy wins the struggle because the two equations of the intrinsic equilibrium solution (eqs.
[12], [13]) and the Poisson equation (∇ 2 ψ = τ ) do not allow for feedback loops and counter-sourcing by self-gravity or rotation. Thus, the state of rotation, the density profile, and the local self-gravity have all been manipulated unilaterally by the action of the enthalpy.

Concluding Remarks
The above equations give us a new probe into gaseous astrophysical disk systems (eqs. [6] and [7]). For spiral galaxy disks, this probe ought to routinely confirm the above-described density-rotation coupling, as the observations already exist. At the same time, we are full of anticipation about what we can learn from the very early phases of purely self-gravitating protoplanetary disks (before the central protostars form and dominate the dynamics) for which the rotation curves have not been measured with accuracy yet, but the radial density profiles in Class 0 YSOs have been obtained (Williams & Cieza 2011). Our prediction, of course, is that, if found, such early disks (preferably earlier than Class 0) will be observed to have "flat" rotation curves. 3

Polytropic Self-Gravitating Newtonian Gaseous Disks
The cylindrical polytropic Lane-Emden equation (Lane 1869-70;Emden 1907) with rotation can be written in dimensionless form as where n > 0 is the polytropic index and the dimensionless constant sound speed c o was defined for ρ = ρ o . (In general, the square of the sound speed c 2 (x) ≡ dp/dτ varies as τ 1/n across the medium.) This equation describes the radial (x) equilibrium of a rotating, selfgravitating, gaseous disk or cylinder in which the gas obeys a polytropic equation of state p ∝ τ 1+1/n . As in § 3, eq. (15) is valid exactly for infinite cylinders and to a high degree of approximation in the equatorial (symmetry) planes of disks (see the Appendix). This latter point is supported by the calculations of Schmitz & Ebert (1987), Schmitz (1988), and Schneider & Schmitz (1995) who studied also the stability of thin-disk and cylindrical equilibria and found large regions of the parameter space with stable models for all values of n > 1; and a sizeable region in which flattened disks with power-law density profiles were unstable to ring formation that causes their profiles to become oscillatory, just as was predicted by the analysis of § 2 above.
We repeat the procedure outlined in § 3 in order to obtain the intrinsic solution of 3 In the case of the Class 0 young system L1527 (Tobin et al. 2012(Tobin et al. , 2013, the rotation curve is not flat, but we did not catch this 0.3 Myr old system early enough (the sound-crossing time at 100 AU for a 10 K cold gas with C o = 0.3 km s −1 is 1500 yr). On the other hand, Yen et al. (2015a,b) report several other Class 0 YSOs whose rotation is not Keplerian outside of the inner few AU. eq. (15): If we equate again the last two terms: then this is an intrinsic solution provided that the rest of the equation (the radial variation of the logarithmic gradient of the enthalpy) vanishes: Eqs. (16) and (17) form a system in which v(x) is totally dependent on τ (x). First we solve eq. (17) to obtain the radial density profile: and then we solve eq. (16) to determine the rotation curve of the intrinsic solution: v where A > 0, n > 0, k < 0, Ax k ≥ 1 (i.e., x ≤ A −1/k ), and the upper incomplete gamma function is defined as The solution contains 4 free parameters, the integration constants A, B, k, and the polytropic index n. Parameter B sets the vertical scale of the rotation curve v(x), so we can choose B = 1 in what follows.
Figs. 7-10 show the shapes of the rotation curves obtained from eq. (19) for two polytropes with n = 1.5 and n = 3 and for various choices of the constants A and k. As in the isothermal case of § 3, the rotation profiles are slowly increasing or flat with radius x. In this case however, we need to obey the condition x ≤ A −1/k (z ≥ 0 in eq. [20]) in the calculation of the gamma function and so the rotation curves terminate when x reaches its maximum value. Two basic trends are noted in the figure captions as well: (a) for fixed n, the curves rise more steeply for steeper values of the index k < 0; and (b) for fixed k, the curves become flatter for higher values of the polytropic index n > 0.

Physical Interpretation
The polytropic Lane-Emden equation with rotation and its intrinsic solution assume the exact same forms as in the isothermal case ( § 3.1) when the polytropic equation of state p ∝ τ 1+1/n is used to introduce the specific enthalpy h(x) ≡ dp/τ . Therefore, the fundamental equations discussed in § 3.1 also apply to polytropic models with finite radial extent and the enthalpy plays the exact same role in manipulating the source terms of the gravitational potential and the rotational potential.
The fact that the thermodynamic potential h(x) operates locally in the exact same manner (i.e., h ∝ ln x) in polytropic and isothermal equilibria helps us correct another common misconception: Since the time that the results of Hayashi et al. (1982) came to light (they studied isothermal self-gravitating gaseous disks that oddly exhibited power-law density profiles and "flat" rotation curves in equilibrium), it has been often stated that isothermal disk models tend to exhibit flat rotation curves because they happen to have a "special" mass distribution (i.e., their specific angular momentum is proportional to the mass interior to radius x, or their mass grows linearly with x). This is not true. The above results show without a doubt that there are no special equilibrium models; and that the isothermal and the polytropic equilibria are both subject to the same fundamental physics at the local level, where the thermodynamic potential h(x) operates and dominates when the only gravity it faces is the self-gravity of the disk. As for the validity of using the cylindrical coordinate system with its "special" ∇ 2 operator, this issue is addressed in detail in the Appendix.
Furthermore, as we have seen above, a flat rotation curve does not imply and does not need a linearly increasing mass distribution (see also the Appendix). Even if the underlying density profile is steeply decreasing in radius (as the light is in spiral galaxy disks), a "flat" rotation curve is a requirement in the equilibrium solutions derived in this section and in § 3. Therefore, one can assume a constant mass-to-light ratio and build a Newtonian selfgravitating galaxy disk model with a steeply declining light distribution that will still be required to have a rising rotation curve.

Concluding Remarks
The polytropic intrinsic solution (eqs. [16] and [17]) gives us a new probe into nonisothermal astrophysical disk systems. In particular, protoplanetary disks undergo various early phases of adiabatic evolution (Fig. 2 in Tohline 2002). We predict that, if found, such disks will be observed to have the same fundamental characteristics (eqs. [18] and [19]) irrespective of the polytropic index n appropriate for each adiabatic phase. In fact, observations of the density profiles, fitted with eq. (18), may be able to determine, not only the profile constants k and A, but also the value of n, thereby deriving the equation of state of the gas independently from theoretical models. The only problem is that we need to find such systems very early in their development (see footnote 3 above), and this is a very difficult task (Kamp et al. 2008;Williams & Cieza 2011;Tsitali et al. 2013).

Discussion
In this paper, we have investigated the Lane-Emden equations with rotation that describe the equilibrium structures of rotating, self-gravitating, gaseous disks and cylinders ( § 3 and 4; see also the Appendix). We have obtained new analytic singular solutions that we call intrinsic solutions because they are dictated and favored by the differential equations themselves with no regard to physical boundary conditions that are externally imposed and that shape up the regular solutions of the equations (the Cauchy problem). In § § 2-4, we have effectively shown that second-order differential equations (both linear and nonlinear) are at odds with the Cauchy problem because the equations show a strong preference for their own intrinsic solutions that, for inhomogeneous equations of the Lane-Emden type, cannot be obtained by solving the boundary-value problem (hence they are singular). The regular Cauchy-type solutions are then attracted to and forced to oscillate about the intrinsic solutions, which means that they do their best to match those dominant solutions. This results in "regular" oscillatory density and rotation profiles whose averages are precisely the underlying intrinsic solutions (Christodoulou & Kazanas 2007;Christodoulou et al. 2015). In this sense, the differential equations succeed in imposing their preferences to the Cauchy problem. This, by itself, is an important conclusion that has ramifications beyond astrophysics for the theory of second-order differential equations of mathematical physics.
The intrinsic solutions are very much related to the so-called trivial solutions of differential equations (Christodoulou & Kazanas 2007). We now understand that there are no trivial solutions, in fact such solutions of second-order equations are quite dominant (see § 2): In many cases of interest, knowing the trivial solution of an equation implies that we know the average behavior of all the regular solutions that depend on various types of boundary conditions but, nevertheless, end up oscillating about the intrinsic solution, provided that the differential equation is a harmonic oscillator (as the ordinary Bessel equations in § 2; the nonisothermal Lane-Emden equations without rotation in § 2; and the Lane-Emden equations with rotation in § § 3-4).
The mean density and rotation profiles that we have derived analytically in § § 3 and 4 are dominated by natural logarithms and power laws. This is the implicit reason that Marr (2015) has recently suceeded in matching the shapes of the rotation curves of a sample of 37 spiral galaxies by using the log-normal probability distribution (and no dark matter or modified gravity) to describe the density profiles in the equatorial planes of the disks. This surface density distribution (eq. [2] in Marr 2015) is equivalent to a variable power law of the form x k(x) in normalized radius x in which the index k(x) varies across the disk as where s 2 is the variance of the distribution, a free parameter to be fitted for each spiral galaxy model. In principle, a slowly varying k(x) is permitted by our analytic solution because the specific enthalpy h(x) is a strictly local potential function (eq. [13] in § 3.1); and if the power-law index k has to vary radially in order for the equilibrium disk to obey some other fundamental law (e.g., as Marr states, the total entropy of the overall configuration should be maximized), then such adjustment may occur on timescales determined by the local sound speed.
Returning to the astrophysical context, the intrinsic solutions of the various types of the Lane-Emden equation with rotation have, for the first time, succeeded in explaining the flat rotation curves of spiral galaxy disks without the need of invoking dark matter or modified gravity. Flat rotation curves are a rigorous requirement of Newtonian gravity in gaseous self-gravitating astrophysical disks. The only fair way to describe this result (see § § 3, 4) is that Sir Isaac Newton (1687) is vindicated once again, and our searches for dark matter in the universe and our attempts to modify Newtonian gravity (references are listed in § 1) have sadly amounted to just a "wild-goose chase." In fact, since the flat rotation curves of spiral galaxies can now be explained at such a fundamental level, the massive observational results collected over the years must be considered as yet another test that Newtonian gravity has successfully passed on scales of ∼ 10 − 100 kpc and at nonrelativistic velocities. This is an impressive achievement when compared to previous tests conducted on and limited to scales no larger than that of our solar system.
Our results render the Dark Matter Hypothesis unnecessary on galaxy scales. This removes the largest pillar of this hypothesis (flat rotation curves have remained to this day the "strongest piece of evidence" in favor of dark matter, but not any longer). But this "aetherial" hypothesis is not about to roll over and die without a fight. The next areas of confrontation will be larger than galaxy scales and cosmology. We are very much encouraged from a recent report of the absence of dark matter on larger than galaxy-disk scales: Magain & Chantry (2013) analyzed 25 gravitational lenses and found that their mass determinations indicate the absence of extended dark matter haloes all the way out to distances comparable to the Einstein ring (the separation between lensed quasar images).
On the other hand, we do not anticipate any serious problems materializing in cosmology because we do not believe that there currently is any credible observational evidence in favor of dark matter or modified gravity on those largest scales. In fact, some results that argue against the necessity for dark matter on various scales have timidly begun to appear (López-Corona 2015;Lane et al. 2015;Moni Bidin et al. 2012. For these reasons, here is how we approach the issue of dark matter now: Christiaan Huygens presented his theory of elastic longitudinal light waves propagating in "aether" to the Paris Academy of Sciences in 1678 and published his views a few years later (Huygens 1690). That year, the physics world entered a Dark Age that lasted nearly 200 years, until the genious of Michelson & Morley (1887) finally showed that the universe is not filled with aether. It seems that we also live in another Dark Age, the "Dark Matter Dark Age," that commenced in the 1970s when K. C. Freeman (1970) and others reported that the rotation curves of spiral galaxies were not falling with radius. By all accounts, the current Dark Age has lasted for nearly 50 years; and the sooner we get out of it, the better for our understanding of the large scales of the universe around us.
The analytic solutions that we have derived in this work should also find applications in the field of protoplanetary disk research ( § § 3.2 and 4.2), especially if very young, purely self-gravitating disks could be found in the future (Kamp et al. 2008). At present, only observations of Class 0 YSOs are widely available (Williams & Cieza 2011), but these systems are not young enough and do not have flat rotation curves; their protostars have formed and they are changing the dynamics and kinematics of the disks. 4 We are encouraged however by the report of Tsitali et al. (2013) who discovered an increasing rotation curve in the inner 2000 − 8000 AU of the "first hydrostatic core" candidate Cha-MMS1.
We thank Joel Tohline for feedback and guidance over many years; and John Marr, Hsi-Wei Yen, and Earl Schulz for many fruitful discussions, especially those concerning the observations of astrophysical self-gravitating disks. DMC was supported in part by NASA grant NNX14-AF77G.

A. Vertically Thin Disks Versus Infinite Cylinders
When applied to the equatorial planes (Z = 0) of thin disks, eqs. (5) and (15) do not include the specific enthalpy term where z ≡ Z/R o . This term is small for cold gases since it scales as the sound speed squared c 2 o . Nevertheless, strong objections have been raised about the validity of this approximation.
Here we address such objections as follows.
The analysis presented in § § 3 and 4 shows that the specific enthalpy in cylinders assumes a logarithmic radial profile (eq. [14]). Now, pressure is an isotropic force and its nature is to push spherically out in self-gravitating gases. Therefore, eq. (14) suggests that, in axisymmetry, the specific enthalpy of the gas ought to assume a spherical form where r is the spherical radius normalized by R o . This behavior is artificially suppressed in cylinders by dropping the dependence on z from the equations. But, in principle, it cannot be suppressed for thin disks, so this term should at least be quantified in the equatorial planes of disks: (a) Applying the radial (x) component of the cylindrical Laplacian in eq. (A2) and then setting z = 0, we confirm eq. (14); that is, this term of the Lane-Emden equations vanishes on the equatorial plane, as was found in the intrinsic solutions of § § 3 and 4.
(b) Combining eqs. (A1) and (A2), we find that This is the term that was ignored in the force balance described by eq. (12). But it makes only a minor contribution to the force balance over all radii where c o << v(x). Specifically, it modifies eqs. (6) and (16) to where ℓ < 0 is a proportionality constant of order c 2 o . For our models, ℓ = (k − 1)c 2 o for isothermal disks and ℓ = nkc 2 o for polytropic disks. Since ℓ < 0 (because k < 0), the last term in eq. (A4) opposes self-gravity. By integration of this equation, we find that its contribution to the rotation profile v 2 (x) is logarithmic: where B is the integration constant. We see now that B and the dimensionless cylindrical "mass" function m(x) = τ xdx make the largest contributions to the rotation curve. This dependence (v 2 ∝ B + m) differs conclusively from the conventional thinking that is responsible for the acceptance of dark matter; that v 2 /r = Gm/r 2 , thus a constant v requires m ∝ r.
We should elaborate on the B + m term in the above analysis. Hidden in this term are the two different Newtonian models that we have discovered and that require rising rotation curves in disks and cylinders: (1) In all isothermal models with k ≥ −1 (Fig. 4) and in all polytropic models (Figs. 7-10), the indefinite integral in eq. (A5) is positive definite. Then v 2 ∝ B + m(x), and the rotation curves rise at all radii as more mass is added by the integration.
(2) On the other hand, in isothermal models with k < −1 (Figs. 5 and 6), the indefinite intergal in eq. (A5) is negative definite (although the definite integral for the mass is still positive). Then, in effect, v 2 ∝ B − |m(x)|, and the rotation curves approach asymptotically the constant value v = √ B as m(x) decreases with x. Because τ (x) is a very steeply decreasing function of x in such cases, then m(x) → 0 from below quite fast, and that makes the rotation curves appear quite flat over most radii in Figs. 5 and 6.
Returning to the pressure term (ℓ ln x) in eq. (A5), if it is included in v 2 (x) in future models, this term has the potential to drive the rotation curves down at large radii because it is negative for x > 1. This behavior can occur only in disks since this term has been suppressed in infinite cylinders. So, unlike disks, cylindrical filaments can never exhibit falling rotation curves in equilibrium (unless of course their age is smaller than the soundcrossing time). Our analysis in § § 3 and 4 was carried out without this term in order to bring out the physics of such Newtonian systems. Nevertheless, the pressure term in eq. (A5) or some similar approximation may be of interest to researchers planning to remodel the rotation curves of spiral galaxies. But it does not appear to be necessary to modeling the rotation of pre-Class 0 protoplanetary disks because such starless disks are subject to extended infall of matter that ought to create cylindrical quasi-equilibria to a good approximation.
This preprint was prepared with the AAS L A T E X macros v5.2.  (2) subject to the following boundary conditions: in the K = 5 case, we use y(1) = 1 and y ′ (1) = −10; in the K = −5 case, we use y(2) = 1 and y ′ (2) = 10. In both cases, the regular solutions are forced to oscillate and stay near the dominant intrinsic solutions y = K (dashed lines).  (4) for integer values of n from 1 to 7 and subject to the usual boundary conditions y(0) = 1 and y ′ (0) = 0. Only the solutions for n = 1 and n = 3 oscillate about the intrinsic solution y = 0 (dashed line).