On the formation of our solar system and many other protoplanetary systems observed by ALMA and SPHERE

TL;DR

This paper proposes an analytical model based on the Lane-Emden equation with rotation to explain the formation of planets in protoplanetary disks, accounting for observed gaps and early core formation in our solar system and others.

Contribution

It introduces a boundary-condition-independent analytical solution that explains the formation of planetary cores through oscillations in protoplanetary disks, aligning with recent observations.

Findings

The model predicts density maxima where planets can form.

Applied to our solar system, it explains the early formation of at least eleven planetary cores.

The approach is being extended to observed disks by ALMA and SPHERE.

Abstract

In view of the many recent observations conducted by ALMA and SPHERE, it is becoming clear that protoplanetary disks form planets in narrow annular gaps at various distances from the central protostars before these protostars are actually fully formed and the gaseous disks have concluded their accretion/dispersal processes. This is in marked contrast to the many multi-planet exoplanetary systems that do not conform to this pristine picture. This major discrepancy calls for an explanation. We provide such an explanation in this work, based on analytical solutions of the cylindrical isothermal Lane-Emden equation with rotation which do not depend on boundary conditions. These ``intrinsic'' solutions of the differential equation attract the solutions of the Cauchy problem and force them to oscillate permanently. The oscillations create density maxima in which dust and planetesimals are trapped and they can form protoplanetary cores during the very early isothermal evolution of such protoplanetary nebulae. We apply this model to our solar nebula that formed in-situ a minimum of eleven protoplanetary cores that have grown to planets which have survived undisturbed to the present day. We are also in the process of applying the same model to the ALMA/DSHARP disks.

Equations

$$\Omega (R)={\Omega }_0\cdot f(x),$$

$$P=c_0^2\cdot \rho ,$$

$$\frac{1}{\rho }\frac{dP}{dR}+\frac{d\Phi }{dR}={\Omega }^{2}R,$$

$$\frac{1}{R}\frac{d}{dR}R\frac{d\Phi }{dR}=4\pi G\rho ,$$

$$\frac{1}{x}\frac{d}{dx}x\frac{d}{dx}\ln \tau +\tau =\frac{{\beta }_0^2}{2x}\frac{d}{dx}\left(x^2{f}^2\right),$$

$$R_0^2\equiv \frac{c_0^2}{4\pi G{\rho }_0}=\frac{c_0^2}{2{\Omega }_J^2}.$$

$$\tau (x)=\frac{{\beta }_0^2}{2x}\frac{d}{dx}\left(x^2{f}^2\right),$$

$$\frac{d}{dx}x\frac{d}{dx}\ln \tau \equiv 0,$$

$$\ln \tau =\ln \frac{{\beta }_0^2}{2}-\ln x+\ln \frac{d}{dx}\left(x^2{f}^2\right),$$

$$\frac{d}{dx}x\frac{d}{dx}\ln \frac{d}{dx}\left(x^2{f}^2\right)=0.$$

$$\frac{d}{dx}\left(x^2{f}^2\right)=A{x}^k,$$

$$\tau (x)=\frac{{\beta }_0^2}{2}\cdot A{x}^{k-1},$$

$$f(x)=\frac{\sqrt{A\cdot g(x)+B}}{x},$$

$$g(x)\ \equiv\ \left\{\begin{array}[]{cc}x^{k+1}/(k+1)\,,&\ {\rm if}\ \ \ k\neq-1\\ \ln x\,,&\ {\rm if}\ \ \ k=-1\end{array}\right.\ ,$$

$$\frac{1}{x}\frac{d}{dx}x\left[\frac{d\psi }{dx}\right]=\frac{1}{x}\frac{d}{dx}x\left[\frac{1}{2}{\beta }_0^2\cdot x{f}^2\right].$$

$$\frac{c_0^2}{{\rho }_0}=4\pi G{R}_0^2=9\times 10^{16}{\text{leavevmode}\mathrm{cm}}^5{\text{leavevmode}\mathrm{g}}^{-1}{\text{leavevmode}\mathrm{s}}^{-2},$$

$${\rho }_0=9\times 10^{-10}\left(\frac{T}{\overline{\mu }}\right)\text{leavevmode}\mathrm{g}{\text{leavevmode}\mathrm{cm}}^{-3},$$

$${\rho }_0\ge 4\times 10^{-9}\text{leavevmode}\mathrm{g}{\text{leavevmode}\mathrm{cm}}^{-3}.$$

$${\Omega }_J=4\times 10^{-8}\text{leavevmode}\mathrm{rad}{\text{leavevmode}\mathrm{s}}^{-1}.$$

$${\Omega }_0=1.64\times 10^{-8}\text{leavevmode}\mathrm{rad}{\text{leavevmode}\mathrm{s}}^{-1}.$$

Full text

1 Introduction

1.1 Planet Arrangements in Multi-Exoplanetary Systems

Discoveries of exoplanets started with Mayor & Queloz (1995) and Marcy & Butler (1996) and in the subsequent 24 years we ended up finding thousands of planets, many of them in multi-exoplanetary systems (www.exoplanets.org). None of the systems containing several planets resembles our own solar system. With this statement, which has also been made by others in the field (e.g., Izidoro & Raymond, 2018; Raymond et al., 2018), we take exception to the many layman articles claiming that exoplanetary systems have been found that are similar to our own (Science Daily, 7/4/2003; Universe Today, 8/24/2010 and 11/27/2013; Space.com, 7/25/2012; Cosmos, 8/3/2015; The Verge, 2/22/2017). The reason for such “fake news” is that professional astronomers have not voiced their opinions loudly, and the reason for their reluctance is probably that the “truth” is utterly disturbing: once again, we are facing a severely unpalatable conclusion, that our solar system is special, and we all know that this cannot be right.

Raymond et al. (2018) recently tried to compile and study some exoplanetary systems that, like our own, have a Jupiter-mass planet at a distance of about 5 AU. In the end, they concluded that our solar system is “weird” and it does not resemble any of the systems currently known. They and Izidoro & Raymond (2018) then argued that our solar system has been sculpted by multiple rare processes that could not successfully combine in any other known system. This premise is unacceptable because it reinforces the point made above about our solar system being special. Furthermore, as we shall see below, it has been invalidated by the recent ALMA/DSHARP and SPHERE observations (ALMA Partnership, 2015; Macías et al., 2018; Avenhaus et al., 2018; Clarke et al., 2018; Keppler et al., 2018; Guzmán et al., 2018; Isella et al., 2018; Zhang et al., 2018; Kudo et al., 2018; Long et al., 2018; Pineda et al., 2018, to name a few).

The right approach to the problem is not to find ways to justify what we all know to be false; instead, we should be trying to understand (a) why protoplanets did not migrate in the solar nebula (unlike in most exosystems), and (b) why we have not yet observed some exosystems that resemble our own. We believe that we already have the answers to these two questions, answers that are not known to many or they have been overlooked by the few.

1.2 The Recent ALMA and SPHERE Observations

The answers to the above questions cannot probably come from planet searches because the sample of nearby stars accessible to observations is severely distance and volume limited. In addition, these systems are fully formed, the gaseous disks are gone, and there has been plenty of time for instabilities to cause migration and destruction of any pre-existing planets in well-organized arrangements, such as the planets in our solar system. On the other hand, ALMA and SPHERE observations are subject to the same limitations to a much lesser degree, they are observing much younger systems, and they can map out young protoplanetary disks with resolutions down to 1 AU. The beautiful images shown in ALMA Partnership (2015), Avenhaus et al. (2018), and Clarke et al. (2018) show relatively thick disks orbiting young or not-fully-formed protostars with planets having already formed and having carved out many annular gaps in the gas.

We point out that the radial brightness maps for CI Tau shown in Clarke et al. (2018) and those shown by ALMA/DSHARP (Andrews et al., 2016; Ruane, 2017; Dullemond et al., 2018; Favre et al., 2018; Guzmán et al., 2018; Huang et al., 2018; Isella et al., 2018; Pérez et al., 2018; Zhang et al., 2018) resemble the radial density profiles that we have also obtained for our theoretical model nebula, which we describe below.

1.3 Isothermal Lane-Emden Equation with Rotation

We have used a new technique to solve analytically the Lane-Emden equation (Lane, 1869; Emden, 1907) for an isothermal, self-gravitating gas in cylindrical symmetry. The analytical solutions are singular and intrinsic to the equation itself, as they do not depend on or reproduced by any particular boundary conditions. In simple terms, these solutions are favored by the differential equation which has no regard for any boundary conditions that may be imposed externally by the Cauchy problem. When the Cauchy problem is solved numerically, its physical solutions are unable to match the preferred intrinsic solution (because of the imposed boundary conditions), but they are attracted to it, and they are forced to oscillate permanently about this fundamental solution. Such mandatory oscillations develop both in the radial density profile and in the gravitational potential. The model incorporates in its input a large family of differential rotation profiles, and differential rotation dictates the spacing and radial extent of the multiple local potential wells. These are approximately equidistant only for the case of uniform rotation. The model lacks scale, but we can remedy this problem in solar-system modeling by fitting the locations of the planets to successive radial potential minima and by setting the distance of the third minimum from the center equal to 1 AU (planet Earth).

Applied to the very early isothermal solar nebula, this model predicts that planets ought to form safely inside consecutive potential minima from which dust and pebbles could not escape, unlike dissipative gas that continues to inflow straight toward the center to form the protosun. The central mass is a lot lower than 1$M_{\odot}$, so we predict that protoplanetary cores ought to form long before the protosun. This prediction, also made by Greaves & Rice (2010), has recently received strong support from the observations of the still-forming star TMC1A (Harsono et al., 2018) that is only about 0.1 Myr old.

One objection to this model has been that it is cylindrical. We do not believe this is a problem because the ALMA-SPHERE disks are not geometrically thin Keplerian disks (see Lee et al., 2017, 2018; Zhang et al., 2018) since the protostars have not fully formed and the gaseous self-gravitating disks are relatively thick and massive. In addition, we also subscribe to an old notion put forth by Sir James Jeans (1914) last century: “All the essential physical features of the natural three-dimensional problem appear to be reproduced in the simpler cylindrical problem, so that it seems legitimate to hope that an argument by analogy may not lead to entirely erroneous result.” 111Jeans (1914) was solving the cylindrical Lane-Emden equation with uniform rotation when he made these remarks. The solution turned out to be a zeroth-order Bessel function oscillating about the intrinsic solution $\tau=\beta_{0}^{2}$ (see below). It took us more than 100 years to do better than that and solve analytically differentially rotating Lane-Emden systems.

Another objection to this oscillatory model stems from a misunderstanding of its application regime. The isothermal nebular model described below does not apply to the present state of the solar system or the ALMA disks that have undergone evolution for millions of years and the gas has been heated by various processes (Tohline, 2002, and references therein). It only applies to the early isothermal evolution of their nebulae and long before any heating processes and ionization set in. For example, we derive in § 2.6.1 maximum typical densities of molecular hydrogen and neutral helium of $\sim 4\times 10^{-9}{\rm\leavevmode\nobreak\ g}{\rm\leavevmode\nobreak\ cm}^{-3}$ in the core of our model solar nebula, values that are roughly consistent with the extended isothermal phase of the solar nebula (Tohline, 2002). Our view then is that extrasolar protoplanetary cores were formed during the same very early isothermal phase of evolution, and the solar-system planets and the observed dark gaps in present-day ALMA disks have survived ever since, as there were no resonant interactions or self-gravity induced instabilities in these disks capable of destroying their pristine, well-organized structures in the subsequent 0.1 Myr - 4.5 Gyr.

1.4 Outline

In § 2, we describe the intrinsic solutions of the Lane-Emden equation and the resulting model of the solar nebula. In § 3, we apply the model to the dark disk gaps observed by ALMA and believed to be the sites of already formed planets. In particular, we search for unstable dark gaps in mean-motion resonances and/or in a Titius-Bode arrangement. These empirical phenomena weigh heavily on to the long-term stability of such still-forming exoplanetary systems. Finally, in § 4, we summarize our results. More detailed modeling of ALMA disks is also currently under way.

2 Intrinsic and Oscillatory Solutions of the Isothermal Lane-Emden Equation with Rotation

2.1 Isothermal Equilibrium Models

We consider the axisymmetric equilibria that are available to a rotating self-gravitating gas in the absence of viscosity and magnetic fields. We adopt cylindrical coordinates ($R,\phi,z$) and the assumption of cylindrical symmetry ($\partial/\partial z=0$) which lets us ignore $z$-dependent gradients and reduces the problem to one dimension, the distance $R$ from the rotation axis. This technique has become common practice in studies of rotating, self-gravitating, fluid disks (e.g., Goodman & Narayan, 1988; Christodoulou & Narayan, 1992; Christodoulou et al., 1996) because it simplifies the stability analyses of effectively two-dimensional modes of disturbance. In what follows, we are interested in equilibrium structures that describe the physical conditions close to the midplane of a geometrically thick protoplanetary disk, so the assumptions $\partial/\partial\phi=0=\partial/\partial z$ allow us to tackle the problem by solving ordinary differential equations (ODEs).

We further adopt a rotation law of the form

[TABLE]

where $x\equiv R/R_{0}$ is a dimensionless radius and the scale length $R_{0}$ will be specified in eq. (6) below. Furthermore, $\Omega(R)$ is the angular velocity, $\Omega_{0}$ is the value of $\Omega$ at some fixed radius, and the dimensionless function $f(x)$ for differential rotation is generally an arbitrary function of $x$. For centrally condensed models, it is convenient to choose $\Omega_{0}=\Omega(0)$ and the regularity condition $f(0)=1$.

Finally, we assume an isothermal equation of state of the form

[TABLE]

where $P$ is the thermal pressure, $\rho$ is the gas density, and $c_{0}$ is the constant isothermal sound speed.

2.2 The Lane-Emden Equation With Rotation

Axisymmetric and cylindrically symmetric, nonmagnetic equilibria for a perfect fluid are described by the equation of hydrostatic equilibrium

[TABLE]

where the gravitational potential $\Phi(R)$ satisfies Poisson’s equation

[TABLE]

where $G$ is the Newtonian gravitational constant. Combining eqs. (1)-(4) and using the definition $x\equiv R/R_{0}$, we find a second-order nonlinear innonhomogeneous ODE that can be cast in the form

[TABLE]

where  $\tau\equiv\rho/\rho_{0}$,  $\rho_{0}$ is the maximum density or a fixed cutoff density for singular/annular models,  $\beta_{0}\equiv\Omega_{0}/\Omega_{J}$,  $\Omega_{J}^{2}\equiv 2\pi G\rho_{0}$,  and

[TABLE]

The term $\Omega_{J}$ represents the gravitational (Jeans) frequency and the dimensionless rotation parameter $\beta_{0}$ measures centrifugal support against self-gravity; in general, $0\leq\beta_{0}<1$, since the gas is also partially supported by pressure gradients in the radial direction.

Eq. (5) reduces to the classical isothermal Lane-Emden equation in the absence of rotation ($\beta_{0}=0$).222And for a flat rotation curve of the form $f(x)=1/x$. In the following subsection, we derive analytically a class of particular solutions of the general problem (eq. (5) with an enormous family $f(x)$ of differential rotation profiles).

2.3 Differentially Rotating Solutions

When the right-hand side (hereafter RHS) of eq. (5) is nonzero (i.e., when $\beta_{0}\neq 0$ and $f(x)\neq 1/x$), the property of scale invariance is lost from all cases of interest (uniform rotation, power-law rotation, etc.), irrespective of the prescription chosen for the differential rotation function $f(x)$.333Eq. (5) with a nonzero RHS is scale invariant only for $f(x)=\sqrt{A\ln x+B}/x$, where $A$ and $B$ are arbitrary constants. This case can be solved by transforming the scale-invariant ODE to its autonomous form, but there is no need to do so; the same solution is obtained easier by the method described in this subsection. Then, eq. (5) has no special symmetry associated with it, and this is probably why some interesting features that we describe below have gone unnoticed for so long.

The RHS of eq. (5) is not merely a rotation-dependent correction term to the classical isothermal Lane-Emden equation. The introduction of rotation changes the properties of the ODE to such a large extent that the known nonrotating solutions (Stodólkiewicz, 1963; Ostriker, 1964) cannot guide the effort to find rotating equilibrium solutions. (This will be evident in Figs. 2 and 3 below.) In fact, it is the functional form of the RHS that determines now the preferred (intrinsic) solutions of the ODE: By equating the last two terms of eq. (5), we can write down an entire class of particular equilibrium solutions, viz.

[TABLE]

provided that

[TABLE]

also holds true. Using eq. (7), we write

[TABLE]

and substituting this form into eq. (8) we find an ODE for all the differential-rotation laws $f(x)$ that satisfy eq. (8) identically and make eq. (7) a family of exact solutions of the Lane-Emden equation with rotation:

[TABLE]

This third-order linear ODE can be readily integrated to yield the following results:

[TABLE]

implying that

[TABLE]

and that

[TABLE]

where $A$, $B$, and $k$ are arbitrary integration constants and

[TABLE]

implying that  $dg/dx=x^{k}$  for all values of $k$. With so many free parameters ($A$, $B$, and $k$) in the differential-rotation profile, these solutions can easily become a theorist’s playground. Here we highlight just a few interesting cases:

1. (a)

Parameter Constraints.—Eq. (12) shows that $\tau(x)>0$ only for $A>0$. This constraint also limits the physical values of $k$ when $B\leq 0$ in eq. (13); for example, $k\geq-1$ when $B=0$. This limitation can be easily circumvented by implementing composite rotation profiles with $B>0$ (see item (d) and § 2.5 below). 2. (b)

Monotonically Decreasing Profiles.—Eq. (12) shows that $\tau(x)$ is a decreasing function of $x$ for $k<1$. The same condition is sufficient to also make $f(x)$ a decreasing function of $x$ provided that $B\geq 0$ in eq. (13). 3. (c)

Uniform Rotation.—For $A=2$, $B=0$, and $k=1$, eq. (13) reduces to $f(x)=1$ and the equilibrium density (eq. (12)) then is  $\tau(x)=\beta_{0}^{2}={\rm constant}$. Note that this constant cannot be adjusted freely, and this is the reason for the oscillatory density profile in a uniformly rotating polytropic model found by Jeans (1914) and later by Robe (1968). (This solution is effectively a zeroth-order Bessel function and it oscillates in a regular manner about the $\tau=\beta_{0}^{2}$ intrinsic solution provided by the uniform rotation.) 4. (d)

Composite Profiles with $B>0$.—Steep density profiles with $k<-1$ can be obtained by selecting $B>0$ and by incorporating a central core region in uniform rotation. Even more complex equilibrium profiles can be constructed by combining two or more density power laws (see § 2.5). 5. (e)

Asymptotic Regime.—For $k<-1$ and $B>0$, eq. (14) shows that $g(x)\to 0$ as $x\to\infty$ and eq. (13) then exhibits the asymptotic behavior $f(x)\to\sqrt{B}/x$. Therefore, all steep density profiles with $k<-1$ and $\tau(x)\propto x^{k-1}$ approach a flat rotation curve ($\Omega R\to$ constant) from above at large radii, independently of the value of $k$. Such a flat rotation curve has been recently observed in the Class 0 young protostellar system HH 211-mms in Perseus (Lee et al., 2018).

2.4 Physical Interpretation

From the perspective of the physics that dictates the above profiles, the solutions (7) of the Lane-Emden equation (5) describe a class of differentially rotating self-gravitating equilibria in which $z$-gradients are neglected and the radial gradient of the gravitational acceleration is balanced exactly by the radial gradient of the centrifugal acceleration at every radius $x$. This occurs because, in the isothermal Lane-Emden equation, we have gone to second order by taking an extra derivative on the components of the equation of hydrostatic equilibrium. The balance of gradients can be seen, most easily, by substituting eq. (7) into the one-dimensional Poisson’s equation ${\bf\nabla}^{2}\psi=\tau$, where $\psi\equiv\Phi/c_{0}^{2}$ is the normalized potential; the result is

[TABLE]

In this equation, the bracketed terms are the gravitational and centrifugal accelerations, respectively. This type of balance is different than the hydrostatic balance commonly discussed between the magnitudes of these two accelerations in rotating gravitating systems; and the power-law density solutions are borne out of this conformance of the two gradients. In the isothermal gaseous case of interest here, a power-law density profile satisfies naturally the condition that the radial variation of the enthalpy gradient  $\rho^{-1}(dP/d\ln R)$  be zero (see eq. (8)) and so the pure power-law profile is not at all influenced by the radial variation of the pressure gradient—it becomes an exact intrinsic solution of eq. (5).

2.5 Composite Equilibrium Profiles

Uniform rotation produces nearly equidistant peaks in the density profile. Such equidistant locations are also found in the inner three planets of our solar system and in the outer planets beyond Jupiter, and they are responsible for invalidating the Titius-Bode rule (Christodoulou & Kazanas, 2007, 2017). A uniformly rotating inner core is also necessary to avoid the singularity of the power-law density profile at $x=0$. We construct then composite equilibria that incorporate an inner and an outer flat region in the density profile. Fig. 1 shows three such intrinsic solutions444Intrinsic singular solutions that do not depend on boundary conditions have been found in many other ODEs and they have often been called “trivial.” As we shall see here, they are extremely important in that they fully determine the behavior of the solutions of the Cauchy problem when boundary conditions are externally imposed. of eq. (5) for slopes of $k=-1.5,-2$, and $-3$.

In Fig. 2, we show a numerical solution of the Cauchy problem with the usual boundary conditions for a centrally condensed model: $\tau(0)=1$ and $[d\tau/dx](0)=1$. The integrations were performed with Matlab (Shampine & Reichelt, 1997; Shampine et al., 1999). The rotating intrinsic solution and the nonrotating solution are also plotted for comparison purposes. The two rotating solutions share the same differential rotation profile (for more details, see Christodoulou & Kazanas, 2007). The following features stand out in the figure:

• (a)

At small radii, the numerical (Cauchy) solution starts out very close to the nonrotating solution.

• (b)

As soon as it crosses below the intrinsic solution, it gets attracted, turns around, and tries to match it. But the Cauchy solution cannot match the intrinsic solution because of the imposed boundary conditions. Thus, it is forced to oscillate permanently about the intrinsic solution, and this is how the many density maxima are created.

• (c)

Density maxima are approximately equidistant in the inner and outer regions, where the intrinsic density profile is flat, although the equilibrium disk is still differentially rotating.

2.6 A Model of the Early Isothermal Solar Nebula

In Fig. 3, we show the same model optimized to fit the locations of the planets in our solar system. Planets are expected to be found at density maxima, where the gravitational potential has minima in which dust and pebbles are expected to be trapped during the entire evolution of the solar nebula. The free parameters of the model and their best-fit values are: $k=-1.5$, $\beta_{0}=0.41$ (or, equivalently, the inner core radius $x_{1}=0.82$ AU), and the outer radius $x_{2}=11$ AU. These are determined by setting the location of the third density maximum (Earth) to 1 AU, which also fixes the radial scale length of the model to $R_{0}=0.022$ AU ($\simeq 4.7$ solar radii). This best-fit model of the solar nebula is slowly rotating and thus stable to nonaxisymmetric perturbations. Its rotation parameter $\beta_{0}=0.41$ is below the critical value555This critical value is obtained from the $\alpha$-parameter criterion for stability of rotating, self-gravitating, gaseous disks $\alpha\leq\alpha_{*}=0.35$ (Christodoulou et al., 1995) by combining $\beta_{0}\equiv\Omega_{0}/\Omega_{J}$ with $\alpha\simeq\Omega_{0}/(\Omega_{J}\sqrt{2})$ to get $\beta_{*}\simeq\alpha_{*}\sqrt{2}\simeq 0.50$. of $\beta_{*}\simeq 0.50$ for the onset of nonaxisymmetric instabilities (see also § 2.6.2 below).

2.6.1 Equation of State and Threshold Density for Planet Formation in the Solar Nebula

Using the scale length of the disk ($R_{0}=0.022$ AU) in eq. (6), we can write the equation of state for the gas as

[TABLE]

where $c_{0}$ and $\rho_{0}$ are the local sound speed and the local density in the inner disk, respectively. For an isothermal gas at temperature $T$,  $c_{0}^{2}={\cal R}T/\overline{\mu}$,  where $\overline{\mu}$ is the mean molecular weight and ${\cal R}$ is the universal gas constant. Hence, eq. (16) can be rewritten as

[TABLE]

where $T$ and $\overline{\mu}$ are measured in degrees Kelvin and ${\rm\leavevmode\nobreak\ g}{\rm\leavevmode\nobreak\ mol}^{-1}$, respectively.

For the coldest gas with $T\geq 10$ K and $\overline{\mu}=2.34{\rm\leavevmode\nobreak\ g}{\rm\leavevmode\nobreak\ mol}^{-1}$ (molecular hydrogen and neutral helium with fractional abundances $X=0.70$ and $Y=0.28$ by mass, respectively), we find that

[TABLE]

This value is comfortably larger than the well-known threshold for planet formation in the solar nebula ($\rho_{*}\simeq 10^{-9}{\rm\leavevmode\nobreak\ g}{\rm\leavevmode\nobreak\ cm}^{-3}$; see, e.g., Lissauer, 1993) and it implies that the conditions for planet formation were already in place, at least in the inner disk, during the early isothermal phase (Tohline, 2002) of the solar nebula.

2.6.2 Rotational State of the Solar Nebula

Using the characteristic density $\rho_{0}$ of the inner disk (eq. (18)) in the definition of  $\Omega_{J}\equiv\sqrt{2\pi G\rho_{0}}$,  we can determine the Jeans frequency of the disk:

[TABLE]

Then, using the model’s value $\beta_{0}=0.41$ in the definition of  $\beta_{0}\equiv\Omega_{0}/\Omega_{J}$, we can determine the angular velocity of the uniformly-rotating core ($x_{1}\leq 0.82$ AU), viz.

[TABLE]

For reference, this value of $\Omega_{0}$ corresponds to an orbital period of 12 yr. In the present solar system, that would correspond to a Keplerian orbit with semimajor axis  $a=5.2$ AU. Thus, the core of the solar nebula was rotating about as slowly as Jupiter is presently revolving around the Sun. This slow rotation of the core provides another indication that the solar nebula was not prone to nonaxisymmetric, self-gravity induced instabilities.

3 Application to the Dark Gaps Discovered by ALMA

3.1 Preliminaries

The above numerical solutions of the Cauchy problem have enjoyed more than ten years of obscurity as viable models of the solar nebula because of two main reasons, both linked to planet migration: (a) Virtually all of the exoplanetary systems discovered so far show gaseous giants that have migrated toward their central stars, having destroyed in the process any pre-existing, well-organized, inner planets. (b) In the solar system, an event called the Late Heavy Bombardment (LHB) has been hypothesized to have occurred $\sim$3.8 Gyr ago in an effort to explain cratering on the moon’s surface. Migration of the gaseous giants and resonant interactions could explain such an event.

Both of these objections appear to be invalid nowadays:

• (a)

Recent ALMA, SPHERE, and Keck/NIRC2 observations show planets forming in circular orbits in young protoplanetary disks. We point out that no images show dark crossovers between the annular gaps, which would be an indication that migration is taking place (ALMA Partnership, 2015; Andrews et al., 2016; Ruane, 2017; Macías et al., 2018; Avenhaus et al., 2018; Clarke et al., 2018; Keppler et al., 2018; Guzmán et al., 2018; Isella et al., 2018; Zhang et al., 2018; Dullemond et al., 2018; Favre et al., 2018; Huang et al., 2018; Pérez et al., 2018; Kudo et al., 2018; Long et al., 2018; Pineda et al., 2018; van der Marel et al., 2019).

• (b)

The LHB has been recently disputed both on observational grounds and by simulations (Reyes-Ruiz et al., 2015; Kaib & Chambers, 2015; Zellner, 2017; Nesvorný et al., 2017; Nesvorný, 2018; Clement et al., 2018). The unfortunate fact that this hypothetical event was successfully modeled by many past simulations invalidates all of these previous modeling efforts altogether. Researchers in the field should really take notice of this fact, as unpleasant as its realization may be.

So it seems that the tide has turned and the answer to question (a) in § 1.1 appears to be that the protoplanetary cores in our solar system formed in situ, safely inside gravitational potential minima in which dust and pebbles were trapped. None of these minima were in mean-motion resonance (MMR) (because of the particular radial density/rotation profiles) after the formation of the solar system concluded, and there were no significant perturbations to stir up the system. Thus, extensive planet migration does not appear to have occurred in our solar system, except for the limited radial excursions that the cores may have undertaken inside their local gravitational potential minima (of order of 1-3 AU for the outer gaseous giants).

Finally, we note that our model’s intrinsic density profile argues against the Nice model (Tsiganis et al., 2005; Morbidelli et al., 2005; Gomez et al., 2005) for the another reason as well: the giant protoplanets needed extended gravitational potential wells (of a few AU) in order to accumulate their large solid cores. Had they all formed in a compact configuration in an outer low-density region of the disk, there would have been insufficient amounts of planetesimals available for the process. In the model of Fig. 3, the density drops by one order of magnitude in the region currently occupied by the gaseous giants, and beyond 10 AU, a slender annulus with a width of a few AU contains barely enough dust/planetesimals to form one giant core.

Based on the recent results produced by ALMA and SPHERE, we conclude that protoplanets in configurations that resemble the planet distribution in our solar system may not be found orbiting around nearby stars (a distance/volume limited sample anyway); instead our chances are much better with circumstellar disks around very young protostars. In such cases, we should check whether any dark annular gaps (a signature of planet formation) are located in MMRs; and whether the Titius-Bode rule plays any role in such young disks. Planets forming near MMRs are a source of future instability and the Titius-Bode rule (Christodoulou & Kazanas, 2017) is important in that observers use it in order to predict the probable locations of any yet undetected planets in multi-planet exosystems.

We note however that our investigation of MMRs continues to be valid in systems in which additional undetected planets have formed, although our study below will be incomplete. On the other hand, the analysis of the Titius-Bode rule becomes invalid in the presence of additional undetected planets. With these issues in mind, we proceed to analyze five well-organized ALMA disks with multiple circular, dark, well-organized gaps.

3.2 No Mean-Motion Resonances in the ALMA Disks

1. HL Tau.—We investigated the dark rings of HL Tau (ALMA Partnership, 2015) for posible MMRs. We calculated all ratios $(a_{i}/a_{j})^{3/2}$ ($i\neq j$) between the semimajor axes of the gaps and no pair of orbits appear to be in MMR. One ratio deviates from exact low-order commensurability by 7.4% and the remaining ratios deviate by more than 10%. Thus, this young solar system appears to be a survivor for ages to come.

2. CI Tau.—Clarke et al. (2018) have reported that the dark gaps of this system are not in a resonant configuration. Based on the relatively high disk mass (0.92 $M_{\odot}$), they speculated that migration may still take place in this environment. Of course, this is not going to occur, if the planets are residing in gravitational potential minima of the gas distribution.

3. AS 209.—The latest ALMA observations (Guzmán et al., 2018; Huang et al., 2018; Zhang et al., 2018) show seven dark annular gaps in this system. We calculated all ratios $(a_{i}/a_{j})^{3/2}$ ($i\neq j$) between the semimajor axes of the rings and no pair of orbits appear to be in MMR. One ratio deviates from exact low-order commensurability by 5.2% and the remaining ratios deviate by more than 12%. Thus, this young solar system also appears to be in a nonresonant configuration.

4. HD 163296.—Isella et al. (2018) show four dark rings in this system. Elliptical fitting of the isophotes has produced precise measurements of their semimajor axes (9.96, 44.77, 86.61, and 140.62 AU). No pair of orbits appear to be in MMR. Two ratios deviate from exact low-order commensurability by 5.0% and 6.9% and the remaining ratios deviate by more than 31%. This young solar system also appears to be in a nonresonant configuration.

5. TW Hya.—ALMA and Keck/NIRC2 observations (Andrews et al., 2016; Ruane, 2017) indicate the presence of dark gaps at 1, 24, 41, 47, and 88 AU. Once again, no pair of orbits appear to be in MMR. A possible exception is the period ratio between the second and fifth dark gaps that may be thought to be in a 7:1 MMR (a deviation of 2.1%, as opposed to all other deviations that are larger than 14%). In this case as well, we cannot conclude that the system is in a resonant configuration.

We conclude that the above protoplanetary systems appear to be stable and they are not in danger of suffering resonant interactions/instabilities after the gas disks are gone. In this sense, these systems appear to be similar to our own solar system, not only because of their well-organized structures, but also because of their long-term stability properties. Thus, we argue that a population of protoplanetary systems roughly similar to the solar nebula has now been found by the latest high-resolution observations of ALMA/DSHARP disks; and this identification provides an answer to question (b) in § 1.1 above.

3.3 No Titius-Bode Arrangements in the ALMA Disks

We also examined the dark rings of the above systems to find out if any successive pairs obey approximately the Titius-Bode rule. In such a case, that would imply that the differential rotation and radial density profiles would be similar to those of the solar nebula. The results show that none of these disks has similar physical characteristics to the solar nebula.

We calculated a geometric progression of the form $a_{i}=\sqrt{a_{i-1}\cdot a_{i+1}}$ and also an arithmetic progression of the form $a_{i}=(a_{i-1}+a_{i+1})/2$ for a uniformly rotating profile, where $a_{i}$ is the semimajor axis of the $i^{\rm th}$ orbit. The smallest geometric deviations between neighboring gaps range from 6% to 39% in the five systems listed in § 3.2. Furthermore, a uniform rotation profile in the inner regions is also ruled out for the two systems each of which shows seven dark gaps (the smallest deviations are 8.3% for AS 209 and 15% for HL Tau, respectively).

But of course there could be more planets forming in these systems that are yet undetected. On the other hand, there is no known universal mechanism (neither one presently hypothesized) that would impose more or less the differential rotation profile of the solar nebula to these disks. Therefore, the above results arguing against the Titius-Bode rule also testify to the expected diversity of radial rotation/density/temperature profiles in young planet-forming protoplanetary nebulae such as the observed ALMA/DSHARP disks (§ 1).

4 Summary

The latest high-resolution observations from ALMA/DSHARP and SPHERE have produced images of protoplanetary disks in which 3-7 planets have already formed (§ 3.2) in well-organized and non-interacting dark gaps, a signature of planet formation very much unlike the picture which has been produced from exoplanet searches (§ 1.1). There is no evidence of violent evolution in these systems, no dramatic planet migrations, no dominant nonaxisymmetric features, no self-gravity induced instabilities, and no stirring from streaming/radial motions in the gas. Unlike in extrasolar systems, nebular evolution has not concluded yet, the gas disks are still present and, for the most part, they are not in Keplerian rotation (see, e.g., Lee et al., 2018), although some form of differential rotation is expected to be pervasive in these systems anyway.

We are then facing an enormous discrepancy between the end-product of planetary formation and the very early stages of ALMA-observed protoplanetary nebulae that occur in relatively massive gaseous disks still in the process of forming their central stars. This discrepancy deserves an explanation.

As far as we know, there currently exists only one physical model of well-organized, ordered planet formation in gaseous isothermal disks. This model was described in § 2. It is based on a formal mathematical property of nonlinear differential equations, such as the isothermal Lane-Emden equation with rotation that is applicable to our context. Such equations admit “intrinsic” (often called trivial) solutions that are independent of boundary conditions. Trivial or not, these solutions are preferred by the equations themselves and they determine all other physically relevant solutions (the Cauchy problem). This is now understood as follows: When boundary conditions are externally imposed, the resulting Cauchy problem does not conform to the preferences of the differential equation and produces new (physical) solutions that are either repelled from or attracted to the intrinsic solutions (Christodoulou et al., 2016a, b; Katatbeh & Christodoulou, 2018). Cauchy solutions attracted to an underlying intrinsic solution have no choice but to oscillate permanently about it, since the chosen boundary conditions forbid an exact match.

The Lane-Emden equation with rotation that is relevant to the solar nebula is one of these equations whose Cauchy solutions are attracted to the intrinsic solution and they are forced to oscillate permanently about it. The corresponding equilibrium density profile develops density maxima (i.e., gravitational potential minima) which can trap dust and planetesimals very early on in the isothermal phase of evolution of real protostellar disks. Protoplanets then form in relative safety inside these potential minima and long before the central stars are formed (see also Greaves & Rice, 2010).

Virtually all protoplanetary disks observed by ALMA and SPHERE and during the DSHAR Project exhibit these characteristics, well-organized annular dark gaps with planets orbiting in circular orbits and in the complete absence of gravitational interactions, radial migrations, or self-gravitating instabilities (§ 3). This robust and novel mathematical model argues strongly against violent planet formation scenarios such as those speculated in the past (planet destruction via migration of gas giants, the Nice model, and instabilities caused by planet crossing of MMRs and/or disk self-gravity). We believe that future observations will only add support to this simple but fundamental physical picture of early planet formation in protoplanetary disks during their early phase of isothermal evolution. We are in the process of modeling the ALMA disks to learn about the physical conditions of these planet-forming young systems.