ARTICLE IN PRESS

Planetary and Space Science 55 (2007) 1299–1309
www.elsevier.com/locate/pss

Hydrostatic, isothermal planetary models from the pointof view of dynamical systems: ‘‘Determination of explicit
formulas for the point of critical mass’’
Johannes Schönke
Arbeitsgruppe für nichtlineare Dynamik, Institut für Theoretische Physik der Universität Bremen, Germany
Received 11 September 2006; received in revised form 2 April 2007; accepted 13 April 2007
Available online 19 April 2007

Abstract
We investigate a planetary model in spherical symmetry, which consists of a solid core and an envelope of ideal and isothermal gas,embedded in a gaseous nebula. The model equations describe equilibrium states of the envelope. So far, no analytical expressions fortheir solutions exist, but of course, numerical results have been computed. The point of critical mass, above which no more static
solutions for the envelope exist, could not be determined analytically until now. We derive explicit formulas for the core mass and the gas
density at the core surface, for the point of critical mass. The critical core mass is also an indicator for the ability of a core to keep its
envelope when the surrounding nebula is removed, because at this point, the core’s inﬂuence extends up to the outer boundary at the Hill
radius.
r 2007 Elsevier Ltd. All rights reserved.
Keywords: Planet formation; Solar system; Autonomous dynamical systems; Point transformations

0. Introduction
In recent years, planet formation has raised a lot of newquestions due to observational advances like the discovery
of extra-solar gas giants. In the nucleated instability
hypothesis, gas envelopes are formed around planetarycores, which in turn are built up by accretion of km-sized
solid bodies. The envelope structure has been investigatedwith hydrostatic models, e.g. by Ikoma et al. (2001),
Mizuno (1980), Papaloizou and Terquem (1999), Stevenson (1982), Wuchterl (1993). There was found an uppermass limit for static envelopes—the critical mass—above
which no static solutions exist for given nebula conditions.
The dependence of this critical mass on nebula conditions
or material properties (opacities) is a quite complex
problem. In the isothermal idealization used in this article,
we can make an analytical approach, giving us deeper
insights into the basic features of core–envelope structuresand especially the critical mass.
Tel.: +49 6221 3282970.

E-mail address: schoenke@ita.uni-heidelberg.de.
0032-0633/$ - see front matter r 2007 Elsevier Ltd. All rights reserved.
doi:10.1016/j.pss.2007.04.002

As there are several deﬁnitions of the critical mass andthe respective critical core mass (CCM) (cf. Wuchterl,1991), we have to specify what we mean with ‘‘criticalmass’’. We follow the deﬁnition introduced by Mizuno(1980) where the CCM is deﬁned as the ﬁrst maximum of
M c ðM tot Þ. For a given total mass M tot an eigenvalue
problem is solved, giving the respective core mass M c . Outof a series of increasing total masses a relation M c ðM tot Þarises. If this relation has a maximum, a contradiction
occurs. Seen as a quasi-static time evolution, the core mass
is a constantly growing structure but the maximum showsthat there is no static solution for a core mass above themaximum, so this core mass is called the critical one after
Mizuno (1980), hereafter called ‘‘classical’’ CCM. Thisdeﬁnition makes no statements about the further evolutionof a critical conﬁguration, only the static possibilities are
ruled out. It should be noted that the described procedure
to get the CCM is done for ﬁxed nebula parameters, i.e. forgiven outer density and temperature. This fact leads toanother description of the phenomenon, also used in thisarticle (cf. Fig. 1), presenting a manifold of possible staticcore–envelope structures in a parameter set of core mass

ARTICLE IN PRESS
1300

J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

Fig. 1. Contour plot of envelope mass at the Hill radius, function of core mass M c and gas density rcs at core surface. Parameters are given in the text.

M c and gas density (or pressure) at the core surface rcs(Broeg, 2006; Pečnik and Wuchterl, 2005). In this description we ﬁnd the point of CCM if we ﬁx the outer
temperature and follow a line of constant outer density inthe direction of increasing total mass. We will reach at least
one maximum in the core mass, corresponding exactly tothe classical CCM. A nice visualization how the M c ðM tot Þplot can be understood as a projection out of the manifold
is shown in Broeg (2006). It shows that the CCM is closely
related to the highest core mass value of the so-called‘‘region IV’’ in the described parameter set, our point O in
Fig. 1. The direct connection to our work gives theisothermal investigation by Pečnik and Wuchterl (2005).
Besides the local CCM, deﬁned there as the outer density
depending classical CCM, they introduce a global CCM
which is independent of the outer density and corresponds
again to our point O. In fact this global CCM is theclassical CCM for the highest possible outer density
permitting hydrostatic equilibrium. The ﬁnding that theclassical CCM is quite independent of the outer density(Mizuno, 1980) can be explained with the strong decreaseof the outer density for increasing core mass beyond theglobal CCM, the point O, because of the strong gravitational inﬂuence the core gets. So the global and classical
CCMs do not differ signiﬁcantly since all outer densities
are realized in a short core mass range. But this behaviourchanges for small orbital distances (Ikoma et al., 2001) or
fully convective envelopes (Wuchterl, 1993). For small
orbital distances this happens because the above-mentioned fast drop in outer density beyond the point O isslowed down due to increasing gas temperature. Naturally
the global CCM has the smallest value under all possible
CCMs with their respective nebula densities, since it isrelated to the highest outer density.

Being a main assumption of this article, the isothermal
approximation has to be discussed. As all attempts tojustify this assumption involve a lot of vague estimations
about opacity, luminosity and other temperature gradientrelated parameter values (e.g. Pečnik and Wuchterl, 2005),
we choose a rather pragmatic way to show the applicabilityof the isothermal assumption. Clearly, for high densities
especially in the compact near core parts of the envelope,signiﬁcant temperature gradients appear. But in theparameter range of interest for this work, i.e. fordensities up to r  1  103 kg m3 (cf. the ‘‘height’’ of thepoint O in Fig. 1), we justify the approximation bycomparing our results with respective non-isothermal
investigations. Critical density (cf. Section 4.1) as well as
the global CCM (cf. Section 4.6) agrees up to a factor of 2
with the respective non-isothermal results. So the isothermal investigation provides more than just qualitative
answers.
For a review of the whole topic cf. Wuchterl et al. (2000),for a detailed discussion of the critical mass cf. Wuchterl(1991).
1. The model and its equations
Our model contains a solid, rigid core, with mass M c andmean density Rc , so that we get the core radiussﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
3 3M crc ¼.
4pRc
The envelope starts at rc with a given density rðrc Þ :¼ rcs(gas density at the core surface). To describe the mass
distribution MðrÞ and density proﬁle rðrÞ in the envelope,
we need the following equations. The force density balanceARTICLE IN PRESS
J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

1301

Fig. 2. Examples for radial density proﬁles for the four different regions in Fig. 1.

for the hydrostatic equilibrium is

2. The numerical results so far and the questions they raiserP ¼ rrf,

The system (3), (4) has previously been examined withnumerical methods, cf. Pečnik and Wuchterl (2005),
Schönke (2005). The integration
was alwaysﬃ stopped at
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
the so-called Hill radius rH ¼ a 3 MðrH Þ=3M S , where a isthe planet’s orbital distance from the central star, MðrH Þ
the mass at the Hill radius (i.e. the total mass M tot ) and M S
the mass of the star. rH approximates the region ofgravitational inﬂuence of the planet (cf. Section 4.4). Fig. 1
shows the envelope mass at the Hill radius M env ¼
M tot  M c , for a variety of values of M c and rcs(a ¼ 5:2 AU, T ¼ 123 K, M S ¼ M , Rc ¼ 5500 kg m3 ,mg ¼ 3:32  1027 kg (molecular hydrogen)).
We can distinguish four different regions. The crosses inFig. 1 mark the M c and rcs values, for which we see a radial
density proﬁle in Fig. 2, one for each region. The densityrðrH Þ at distance r ¼ rH from the planet connects to thenebula density, so if we embed the planet in a certain
nebula, we have to choose the right proﬁle. Now, we give avery short summary of the nature of the different regions.
Region I (‘‘earth-like planets’’) has massive cores comparedto the mass of the envelope, the gravitational potential isdominated by the core and the density decreases exponentially with radius (like on earth with our atmosphere as theenvelope). Region II (‘‘gas giant planets’’) has high rcs
values, so there is so much mass in the envelope near thecore, that self-gravity of the gas sets in and forms acompact envelope; for radii further out, the density drops
very fast to vacuum. Proﬁle II in Fig. 2 is set constant tolog r ¼ 50 (perfect vacuum) where integration gives
smaller values. Region III (‘‘asteroids’’) has small cores,
which hardly inﬂuence the surrounding gas so that theenvelope is quasi-homogenous, i.e. rðrc Þ  rðrH Þ. Region
IV (‘‘protoplanets’’) has massive and extended envelopes,

where P is the gas pressure and f the gravitational
potential. In radial symmetry and with
rf ¼ f g ¼

GM r
,
r2 r

with f g the speciﬁc gravitational force, we getdPðrÞGMðrÞrðrÞ
¼
.
dr
r2

(1)

In addition, we need the equation of state for the ideal gas
in the isothermal case1PðrÞ ¼

kB T
rðrÞ;
mg

T ¼ const,(2)

with mg the atomic mass of one gas particle. Introducing(2) in (1) we obtain
Gmg MðrÞrðrÞ
drðrÞ
¼
.
dr
r2
kB T

(3)

Mass density rðrÞ and mass distribution MðrÞ are related bydMðrÞ
¼ 4pr2 rðrÞ.
dr(4)

Together, Eqs. (3) and (4) represent a nonlinear system ofordinary differential equations, whose solutions are uniquely determined given the two initial conditions Mðrc Þ ¼
M c and rðrc Þ ¼ rcs .
1

Up to r  1  103 kg m3 , the assumption of ideal and isothermal gas isa valid approximation. For an extensive discussion of the applicability ofideal gas and the isothermal assumption cf. Pečnik and Wuchterl (2005).

ARTICLE IN PRESS
J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

1302

which get self-gravitating beyond some radial distance (thedensity decrease changes from exponential to 1=r2 , cf. thebend in the proﬁle IV of Fig. 2), but are not compact yet.
At a speciﬁc rcs value, the near core parts of the envelope
become compact in a local sense, but do not inﬂuence thewhole structure.
The question is how to make a clear statement about
where the transition occurs between, e.g. region III and I.
Precisely, if we assume the border between III and I in Fig.
1 to be vertical: which minimal mass does a core need, so
we can call it an ‘‘earth-like planet’’ rather than an‘‘asteroid’’? Another question concerns the problem where
we can ﬁnd the so-called point of critical mass (O in Fig. 1),
deﬁned to be the point where the four regions in Fig. 1
meet? This point should lie on the upper end of the borderline between III and I, of course, and its rcs should be equal
to the speciﬁc value, for which objects in IV get locallycompact near the core. But which value is it? The
procedure so far was to plot something like Fig. 1 and totell the points from this graphic. But this is unsatisfactory,
because we do not know what happens, if we change thetemperature or other important variables, like the orbital
distance (which will also change the temperature), thedensity of the core or the atomic weight of the gas. What
happens at these border lines and at O in particular? Whysuch a sudden change in the qualitative behaviour? To
answer these questions we need a new representation of thesystem (3), (4) and an analytical expression for O,
depending on quantities like the above-mentioned ones.
In the following, we delineate a way to achieve these goals.
3. Mathematical preparations
This section is a rather technical part. It discusses thetransformation of the original system into an autonomous
one, i.e. a system which does not contain the independent
variable explicitly. In this form, we can represent all
solutions in terms of a phase portrait which gives us adeeper insight into the system’s properties. The application
to the physical problem starts with Section 3.3. The reader
not interested in the mathematical reasoning may skip theother sections.

(7)changes (5) and (6) into
dy
C yz
¼
,
dv
3 v4=3

(8)

dz
¼ y.
(9)
dv
If we differentiate (9), z00 ¼ y0 , we can introduce this into
(8) and change the system of two ﬁrst-order differential
equations into one second-order equation
d2 z
C
dz
(10)
¼  v4=3 z .
2
dv
3
dv
This is a generalized Emden–Fowler equation, i.e. a special
case of s€ ¼ Atn sm s_l with A; n; m; l 2 R. In our casen ¼ 4=3, m ¼ 1, l ¼ 1 (and A ¼ C=3),2 no analyticalsolution is known so far. (For other combinations of n; m; l,there exists information about solutions, mainly in parametric form. In addition, a general procedure for arbitraryn; m; l is known, which reduces this kind of equations to an
Abel equation of the second kind, but this again is notanalytically solvable cf. Polyanin and Zaitsev (2003).)
However, a particular solution can always be given which
has no free constants of integration and describes theasymptotic behaviour of the general solution. This solution
will play an important role in the following. More details
about generalized Emden–Fowler equations are given in
Polyanin and Zaitsev (2003).
3.2. Reduction to an autonomous dynamical system via pointtransformations
To extract the asymptotic behaviour of z at large valuesof v, we make the ansatz3
z~ ¼ va z;

v~ ¼ ln v,



d2 z
d~z d2 z~
a2
þ
¼
v
aða

1Þ~
z
þ
ð2a

1Þ
,
dv2
d~v d~v2

First, we introduce dimensionless variablesMðrÞ ¼ M  zðxÞ;v ¼ x3 ) dv ¼ 3x2 dx ¼ 3v2=3 dx,(11)where a is a free parameter. Using (11) together with


dz
d~z
¼ va1 a~z þ
,
(12)
dv
d~v

3.1. Transformation to a generalized Emden–Fowler
equationr ¼ r0 x;

2=3 1=3

C ¼ 62=3 p 1 =3GM  Rc mg =3kB T  4:259  102 ð ðRc =½kgm3 Þ1=3 ðmg =½uÞ=ðT=½KÞÞ, (u ¼ 1:66  1027 kg (atomic
mass unit)). The transformation of the independent
variable x into a kind of ‘‘volume-like’’ variable v in theformrðrÞ ¼ Rc yðxÞ,where M  ¼ 5:976  1024 kg is thepearth
mass, ﬃ Rc is theﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
mean density of the core and r0 ¼ 3 3M  =ð4pRc Þ. So now,
(3) and (4) becomedyyz¼ C 2 ,
dx
x

(5)

dz
¼ 3x2 y,
dx

(6)

(13)

in (10), we obtainva2 ½   ¼ C 2a7=3
d~z
v
a~z þ
z~.
3
d~v

(14)

2
The factor A is unimportant, it can always be changed to unity with thetransformation s ! A1=ð1mlÞ s~ (if m þ la1).
3
That (11) extracts the asymptotic behaviour is not obvious, it will turnout later that it is just the case.

ARTICLE IN PRESS
J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

As a is a free parameter, we can choose it in such a way thatthe exponents of v in (14) are equal on both sides of theequation, i.e. demandinga  2 ¼ 2a  73 ) a ¼ 13.
With (14), we have our requested autonomous equationd2 z~ 1 d~z C d~z
2C¼
z~ þ z~  z~2 .
2
3 d~v 3 d~v
9
9
d~v

(15)

From Eq. (15), it is straightforward to construct theautonomous dynamical system. But some care must be~ which istaken in the choice of our second variable y,needed to form a ﬁrst-order system again. The simplest
choice, y~ ¼ d~z=d~v, would lead us to negative y~ for some yrange, as can be seen with (9) and (12). A much better
alternative, showing analogies to the z transformation (11)(as we will see later), isd~z
1¼ y~  z~,
d~v
3

(16)

we obtain
dj¼ c  j,
ds

)

(17)

which keeps y~ always positive (use (9) and (12) to provethat). Introducing (16) and (17) into (15), we getdy~ 1~  C z~Þ.
¼ yð2
d~v 3

(18)

3.4. Connections to homology invariants
A comparison of our variables with the well-known
homology invariants U and V, given by (e.g. Kippenhahn
and Weigert, 1990)

3.3. A physical approach to the problem
With the information accumulated in the last sections,
we have the opportunity to review the problem in a clear
way. We deﬁne the dimensionless quantityjðrÞ :¼

Gmg MðrÞ
,kB T r(19)

which is obviously a global gravitational potential4 divided
by the thermal energy, andcðrÞ :¼

Gmg
4pr2 rðrÞ,
kB T

(20)

which is a bit more complicated to describe. With (4), we
see that c / dM=dr, which is the differential mass increase.
But c is also something else, namely a local gravitational
potential (again divided by the thermal energy). It is ameasure for the gravitational interaction of neighbouringmass shells and therefore an indicator for self-gravity. We
will see that the whole structure of the envelope can beunderstood by investigating the interplay of the twopotentials j and c. Introducing (19) and (20) into (3)and (4) and then substituting the independent variable
s :¼ lnðr=rc Þ,
4

It is convenient to have j as well as c positive by deﬁnition.

d ln M 4pr3 r¼,
d ln r
M

(21)

(24)

d ln P r GM¼,
(25)
d ln r P rshows some similarity. Since we consider an ideal gas, i.e.
r=P ¼ mg =ðkB TÞ ¼ const (cf. (2)), a comparison of (24)and (25) with (19) and (20) leads toV :¼ 

j ¼ V;

Eqs. (16) and (18) represent our dynamical system.

(22)

dc
¼ cð2  jÞ.
(23)
ds
The system (22), (23) has the same structure as (16), (18),i.e. it is an autonomous dynamical system. The advantage
now is that we have a clear picture of the quantities j andc, in contrast to z~ and y~ in the last sections (but as a matterof fact, z~ and j or rather y~ and c are identical up to afactor).

U :¼
d2 z~ dy~ 1 d~z dy~ 11

¼
 y~ þ z~,
¼
2
~
d~
v
3
d~
v
d
v
3
9
d~v

1303

c ¼ UV .

(26)

So j and c are homology invariants as well, because theproduct UV is still an invariant quantity. In fact all
homology invariants lead (together with (21)) to autonomous systems of differential equations. Maintaining
homology, i.e. to ensure a kind of similarity between
different solutions of a system, requires the spatial
independence of the evolution of the invariants.
With the variables U and V we obtain a much more
complicated system than (22) and (23). Our equations
contain only one non-linearity (the term c j in (23)),
whereas in the U  V variables we would get quadratic
terms. No combination of the invariants was found leadingto ‘‘simpler’’ expressions than (22) and (23).
4. Properties of the system
There are two important critical lines in the system (22),(23) (the third critical line, c ¼ 0 is trivial: r ¼ 0), namely
c¼j½with dj=ds ¼ 0,

(27)

j¼2

½with dc=ds ¼ 0,(28)

describing extrema of j and c, respectively. These lines
together with the phase portrait of the system are shown in
Fig. 3. All information about the radial proﬁles MðrÞ aswell as rðrÞ is contained in the ﬁgure. Given initialconditions M c , rcs and rc (via Rc ), we obtain from (19)and (20) an initial point (jc ; cc ) in Fig. 3.

ARTICLE IN PRESS
J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

1304

Fig. 3. Phase portrait for the system (22), (23), with the critical lines
c ¼ j and j ¼ 2.

As we are interested in positive values of j and c, themost noticeable object in Fig. 3 is the global attractor at
the point ðj; cÞ ¼ ð2; 2Þ, which follows from dj=ds ¼
dc=ds ¼ 0 at this point, i.e. the two critical lines c ¼ j
and j ¼ 2 intersect there, and the fact that every initialpoint ðj; cÞ 2 Rþ  Rþ is drawn to it for s ! 1 (ignoringthat we stop at the Hill radius sH / ln rH ). With (19) and(20) we see immediately that ðj; cÞ ¼ ð2; 2Þ ¼ const is
related to the proﬁles r / 1=r2 and M / r, which weresaid to be typical for the self-gravitating, noncompact partsof proﬁles in region IV (cf. Section 2). That means that all
proﬁles in all regions would have such a part, if their Hill
radius were not reached or they had not become compact
before. In general, every gas sphere, with or without a core,
would evolve into this asymptotic structure, for big radii.5
In these envelope parts the gas interacts with its local
vicinity only. The solution at ð2; 2Þ is the particularsolution, mentioned in Section 3.1. Our transformations(19), (20) extracted the asymptotic behaviour r / 1=r2 andM / r from the original variables r and M.

the mean density, local self-gravity has the same inﬂuenceon the gas as the total mass inside the actual sphere. For
c4j, the local self-gravity of the gas dominates the global
gravitational potential. Here, we see a fundamental
conclusion for all isothermal core–envelope conﬁgurations
in radial symmetry. Applied to the core surface, withhrðrc Þi ¼ Rc , the critical density rcrit is nothing else but the
speciﬁc rcs value in region IV, mentioned in Section 2,
because with the local self-gravity beginning to overrule theglobal potential at the core surface, we just have theproperty which distinguishes between a core-gravity resp. a
self-gravity dominated core surface, which was said to be
the basic difference between planets below resp. above thepoint X in Fig. 1, Section 2. Therefore rcrit is one
coordinate rO of the point X!
R(30)
rO ¼ c .
3
The beauty of this important result arises from its
incredible simplicity! There is a speciﬁc density, for which
the gas at core surface starts to show self-gravity effects.
This critical density is reached, when it is exactly one thirdof the mean density of the rigid core (which can havearbitrary density distribution). For this transition point, we
ﬁnd yet another change in the stability behaviour, as a new
instability sets in and ampliﬁes the already existing
instability (cf. Schönke, 2005). The comparison of (30)with isothermal numerical results shows exact agreement(e.g. Pečnik and Wuchterl, 2005). Non-isothermal investigations by Broeg (2006), including detailed EOS andenergy transfer, contain characteristic values for gas
pressure and temperature at the core surface for thecorresponding point O. Using the respective EOS,
we get the critical density at the core surface for thatpoint. Interestingly, this critical density has the samevalue for different states, i.e. different temperaturesand pressures at the core surface (due to, e.g. different
nebula parameters). This is a clear agreement with our
result that the critical density exclusively depends on thecore density. The non-isothermal value is approximately
31000 kg m3 , whereas we get 5500
 1800 kg m3 .
3 kg mThe reason for the generally smaller non-isothermal value
seems to be a non-trivial question which cannot beanswered so far.

4.1. The critical density
4.2. The minimum core mass
The ﬁrst critical line (27) describes the points where an
important event takes place: the strength of c becomesequal to the strength of j. Introducing (19) and (20) into
(27), we get a critical density
rcrit ðrÞ ¼

MðrÞ hrðrÞi
,
¼
4pr3
3

(29)

where hrðrÞi denotes the mean density inside the sphere ofradius r. Whenever the local density reaches one third of5
Because of the outer boundary (Hill radius), only proﬁles in region IV
have this structure.

If we deﬁne a speciﬁc, asymptotic energy
 :¼

2kB T
,
mg

(31)

the second critical line (28) can be identiﬁed with (using (19))
GMðrÞ
¼ .
(32)
r
The term ‘‘asymptotic energy’’ should become clear after
looking at the asymptotic solutions of the trajectories of oursystem: As we get closer to the attractor, GM=r tends towards

ARTICLE IN PRESS
J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

1305

 (i.e. j ! 2). If GM=r4 (j42), a gas particle feels astronger force Fg ¼ GMmg r=r3 towards the centre than in
the asymptotic case, whereas the force is lower if GM=ro(jo2).
Again applied to the core surface (M ¼ M c , r ¼ rc ), thecondition (32) leads us to a minimum core mass



 

3 1=2   3=2
3 1=2 2kB T 3=2
min
Mc ¼
¼
. (33)
4pRc
G
4pRc
Gmg
At a core mass of M min
the envelope at the core surfacec‘‘feels’’ like it is far away from the core and evolves in theasymptotic way. Core masses smaller than M min
havecpractically no inﬂuence on the envelope, the gas reacts onthe ‘‘lacking potential’’ even with a small density excess at
the core surface. In general, with a value of M minc , a core
starts to influence the envelope, it is the ‘‘core mass’’
analogue to rO . The counterpart of this border line,
representing M minc , is not so obvious in Fig. 1 because theplot shows the envelope mass at the Hill radius, and theinner envelope structure is hidden in this global property.
Nevertheless M min
can be found in form of the linecconnecting the minima of the contour lines in region III atlogðM minc =M  Þ ¼ 2:85. These minima correspond to
maxima in the envelope mass M env . For M c oM minc ,M env increases with M c due to the growing Hill radiuscaused by the increasing (but not yet important) core mass.
For M c 4M minc , the core forces the density to decrease
faster near the core to maintain hydrostatic equilibrium,
and M env decreases.
4.3. Oscillations in the asymptotics
The outer envelope parts of the objects in region IV are agood example for the interplay between the potentials jand c. An interesting feature of these envelopes, related tothe asymptotic solutions and not mentioned in Section 2, isa damped oscillation ﬁrst observed by Pečnik (2001). The
density proﬁles rðrÞ, e.g. are oscillating around 1=r2 andreach it for r ! 1 (this can hardly be suspected from Fig.
2, the picture is too small). An example for such a density
proﬁle in region IV is shown in Fig. 4. To discuss thisfeature, we investigate the vicinity of the attractor, i.e.
perform a linearization of the phase portrait at the criticalpoint. The two eigenvalues of the linearized system are
pﬃﬃﬃ
l1;2 ¼ 12ð1 i 7Þ,(34)which show us, that j andpcﬃﬃﬃ behave like a damped
oscillator, i.e. / exp½ð1=2 i 7=2Þs. It means that there
is an alternating up and down of global and local
gravitational energy, as a proﬁle evolves outwards. Thisis obviously a critical situation for the envelope: it
periodically starts some kind of self-interaction (i.e. selfgravity), but it is always a local phenomenon, instead of acompact state (which would be global). This behaviour is ahint that these kinds of envelopes may not be stableobjects, which has indeed been found (cf. Schönke, 2005).

Fig. 4. Example for the damped oscillation of density proﬁles in region
IV. The core mass is 1022 kg.

Out of curiosity, we calculate the wavelength of theoscillation. Frompﬃﬃﬃ (34), we tell the (dimensionless) wavenumber k ¼ 7p=2
ﬃﬃﬃ and therefore the wavelengthLs ¼ 2p=k ¼ 4p= 7 pﬃﬃ4:75, so we get the ‘‘real’’ wavelength Lr ¼ eLs ¼ e4p= 7  115, which is not a length6 but
tells us, if we started at an arbitrary radius r0 , that we went
through one period after reaching 115 r0 .
4.4. The Hill radius
Most of the structures in Fig. 1, like the shapes of thedifferent regions and their borders, discussed in Section 2,
are strongly related to the choice of the outer boundary of
the envelope. Especially the border lines between region III
and IV as well as I and III, and therefore the CCM M O ,depend on this choice, in contrast to rcrit or M min
c , whichare not inﬂuenced by the outer boundary, but represent
points of fundamental changes concerning the inner
structure of a core–envelope conﬁguration. We use theHill radiuspﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
rH ðMðrÞÞ ¼ a 3 MðrÞ=3M S ,(35)
as the outer boundary of the envelope. rH describes theregion of gravitational inﬂuence of the planet. In fact, it
presents the oscillating sphere of the surface where thegradient of the effective potential7 vanishes. The key
problem of the incorporation of rH is its dependence onMðrÞ (cf. (35)) which in turn depends on rH .
4.5. The border line between III and IV
The way to understand the border line structure of
region IV is to ask the following question: Suppose we are6

Due to the autonomous equations we ‘‘lost’’ the length scale, thesystem is invariant under scaling of r.
7feff ¼ fstar þ fplanet þ fcentrifugal in the corotating system.

ARTICLE IN PRESS
1306

J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

exactly at the Hill radius. If we make another radial step droutwards, do we get so much more mass dM, that the newHill radius rH ðM þ dMÞ equals our new positionrH ðMÞ þ dr? This borderline case deﬁnes the boundary
between region III and IV. The answer to the above
question is either ‘‘No, we get even more mass, so the newrH is bigger than our new position’’, if we are in region IV,or ‘‘No, we do not get enough new mass, so the new rH issmaller than our new position’’, if we are in region III. It isclear that this explains the sudden increase in M env at theboundary between the two regions: In region IV we can go
further out, because rH þ dr is still inside the new Hill
sphere, whereas in region III we have to stop, because
otherwise we would leave the region of gravitational
inﬂuence of the planet.
So we have to compare the differential mass increases atthe Hill radius. With (35), the critical increase isdM 9M S 2
¼ 3 rH .
drH
a

(36)

The general mass increase at rH is (using (4))
dM
¼ 4pr2H rðrH Þ.
dr rH

(37)

(38)

which can be integrated easilyGmg M c 1 1rðrÞ ¼ rcs exp

,
r rc
kB T



Gmg M c 11
.
rH :¼ rðrH Þ ¼ rcs exp
rH rc
kB T
After eliminating rc and rH we get8
"
!

 #
Gmg ð3M S Þ1=3
4pRc 1=3
rH ¼ rcs exp
M 2=3
.
c
kB T
a
3

(40)
For M c ! 0 we get the critical central density of a gas
sphere, above which it can build up an extensive envelope
9M S
,
(41)
4pa3
which is called r̄, because it is the mean density hrðrH Þi ofan object, cf. (29) and (35) (all models have the same r̄ forﬁxed a and M S ). It is also the highest possible outer densityof any envelope. Actually, the outer density equals r̄ forevery point on the border line of region IV (towards region
III and region II) in the approximation (38). There are no
static solutions for outer densities higher than r̄, thereforer̄ is a critical nebula density, a density analogue to theCCM.

ðM c ! 0Þ ¼
r̄ :¼ rIII;IVcs

4.6. The critical core mass

To ﬁnd an expression for the outer gas density rH :¼ rðrH Þ,
we can take advantage of the fact that there exists
an analytic solution for the regions I and III, which is aperfect approximation, as long as we stay out of region II
or IV. As a matter of fact, we use an approximation to ﬁnd
out where the approximation itself breaks down! As inregion I and III the core mass is always at least ahundred times larger than the envelope’s mass, we cansay that the mass is constant and equal to the core mass.
This procedure, known as the Roche approximation,
simpliﬁes (3) to
Gmg M c rðrÞ
drðrÞ
¼
,
dr
k B T r2

and the border line between region III and IV to be
!
"
#

Gmg
9M S
4pRc 1=3 ð3M S Þ1=3
III;IV
2=3
rcs ðM c Þ ¼
exp

.
Mc
4pa3
kB T
3
a

(39)

If we introduce (39) into (37) and then equate (36) and (37),
we ﬁnd"!
 #Gmg ð3M S Þ1=3
9M S
4pRc 1=3
2=3
¼ 4prcs exp
Mc
,
a3
kB T
a
3
8
If a40:1 AU (with M S ¼ M  , Rc ¼ R ), then rH brc and the term
ð3M S Þ1=3 =a is negligible in all the forthcoming formulas.

With the knowledge of the critical density at the coresurface rO (cf. Section 4.1), it is now straightforward to
calculate the CCM M O with the help of (40). We just needto ask for the M c value where the borderline between III
and IV crosses the critical density rO , i.e. we have to solve
the equation rIII;IVðM O Þ ¼ rO ¼ Rc =3, givingcs0 13=2
kB T
4pa3 Rc
ln
B Gm
C
27M S
g
B
C
M O ¼ B(42)
C .
@ 4pRc 1=3 ð3M S Þ1=3 A3aWith M O we found the second and ﬁnal coordinate of thepoint of critical mass O! We can think of M O as the end of
a process, which began with the minimum core mass. At
M min
the core started to inﬂuence the envelope near theccore, whereas at M O the core interacts with the wholeenvelope up to the Hill radius. For the parameters like in
Fig. 1, we get a CCM of logðM O =M  Þ  1:17, in
excellent agreement to numerical results from otherisothermal investigations (Pečnik and Wuchterl, 2005).
The obvious discrepancy to non-isothermal publications(Mizuno, 1980; Ikoma et al., 2001), which give values up to
the order of 10M  , has three main reasons. As mentioned
before, M O is the global CCM, i.e. the CCM for the highest
possible outer density r̄ and therefore the minimum CCM
value for all possible nebula densities (cf. Introduction).
But if we look into Fig. 1 in Ikoma et al. (2001), we canactually ﬁnd the respective M O value. There, the classical
CCM is shown at 0.1 AU for different nebula densities up
3
to 1000 rH(rH
0  0:67 kg m
0 is the MMSN density,
Hayashi et al., 1985) which is close to but already above
the (isothermal) maximum outer density of r̄ ¼ 0:42 kg m3
at that orbital distance (cf. (41)). For 1000 rH
0 , i.e. the veryright end of the graph, the CCM is given as

ARTICLE IN PRESS
J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

1307

M critc  0:55M  , where (42) with the respective parameters
gives M O ¼ 0:38M  , which is not such a big deviation
anymore. The second reason for the still remaining
difference is of course the isothermal assumption giving
naturally smaller values for the CCM because the constanttemperature is certainly a lower limit for the planet
especially in the high density parts close to the core.
Upcoming analytic investigations with polytropic relations(to be published soon) will show signiﬁcantly increased
values for M O . Last but not least there are deviations due
to non-ideal effects of the gas not considered in this article.
4.7. The border line between I and III
As we mentioned before, M O also represents the borderline between the regions I and III, i.e.
M I;III¼ M O,
c(43)

so an object with M c 4M O can be called an ‘‘earth-likeplanet’’ and its core holds the envelope in a ‘‘ﬁrm grip’’. If
we deﬁne a density contrastrd :¼ cs ,(44)
rH
we see with (39) that in region I and III, d is constant for aconstant core mass. At M c ¼ M O , we ﬁnd
dO ¼

4pa3 Rc rO
.
¼
r̄
27M S

(45)

Therefore, the ability to keep an atmosphere when the
nebula is removed can also be described with a minimum
density contrast of dO . If dodO , we have just an ‘‘asteroid’’.
4.8. Properties of region IV
Proﬁles in region IV end up at ðj; cÞ ¼ ð2; 2Þ sooner orlater, as we mentioned before. The total mass and the outerdensity do not vary much in the whole region, as can beseen in Fig. 1. If we evaluate (19) and (20) for ðj; cÞ ¼ ð2; 2Þ
at r ¼ rH , we obtain the mean values for M env (the heightof the ‘‘island’’ in Fig. 1) and rH in region IV
2kB Ta 3=2
IV
hM env i ¼
ð3M S Þ1=2 ,
(46)
Gmg
hrIV
Hi¼

3M S r̄¼ .
4pa3 3Fig. 5. Summary of the characteristic quantities, discussed in the previous
sections, shown in the M c  rcs parameter space (cf. Fig. 1).

4.9. A brief summary of the formulas
For an overall review, Fig. 5 collects the important
quantities calculated so far. A comparison of its lines andpoints to the numerically computed structures in Fig. 1
shows excellent agreement. Using the abbreviations
C :¼ ð4pRc =3Þ1=3 and D :¼ ð3M S Þ1=3 =a ¼ ð4pr̄=3Þ1=3 andE :¼ kB T=ðGmg Þ we can summarize the derived formulas9
R
rO ¼ c ,
3
9M S
,
4pa3
 3=2
2E
M min
¼
,
c
C
r̄ ¼



 3=2
E
r
ln O
,
MO ¼
CD
r̄


C  D 2=3
M
rIII;IV
ðM
Þ
¼
r̄
exp
.
c
cs
c
E

4.10. Dependence of the critical core mass on the orbital
distance(47)

These values are exact at the point ðM c ; rcs Þ ¼ ðM minc ; rO Þ,of course. The results are still correct for non-isothermal
planets if only the outermost envelope parts are isothermal,
as comparisons with Broeg (2006) show. The temperature
to be inserted in (46) is then of course the one from theisothermal outer part, i.e. the nebula temperature. Aninteresting fact is that hrIV
H i corresponds to the ‘‘mean
density inside the orbit’’ (M S distributed over a sphere withradius a), yet another property of the asymptotic solutions.

It is interesting to investigate the dependence of M O on
the different parameters, especially the orbital distance a,because it usually changes the temperature, too. If we use acommon temperature model for a solar nebula, e.g. by
Hayashi et al. (1985)
1=2

 
T
L 1=4
a
¼ 280
,
(48)
½K
½L 
½AU
9

If abð9M S =4pRc Þ1=3 , then CbD and D can be neglected.

ARTICLE IN PRESS
1308

J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

Fig. 6. The critical core mass M O (thick lines) and the minimal core mass
(thin lines) as functions of orbital distance for two different
M mincmolecular weights. The temperature dependence TðaÞ is assumed to be(48). The molecular weight is just shifting the functions (cf. (42) and (33)).

with L as the luminosity of the central star, we caneliminate T and get a relation M O ðaÞ / a3=4 ½lnðk aÞ3=2
with a constant k depending on M S and Rc (cf. (42)). M O ðaÞ
and M minc ðaÞ are shown in Fig. 6 for two molecular weights(and L ¼ L , other parameters like in Fig. 1). Forlogða0 =½AUÞ  2:07, M O equals M min
c . This happens
because due to the decreasing orbital distance, r̄ increases
and gets close to rO , the needed density contrast at theCCM becomes very small (dO  1) and so does M O . In
contrast the minimum core mass increases continuouslywith decreasing orbital distance (or increasing temperature). For aoa0 it is not possible to dominate an envelopewith the core potential alone. Only compact objects (e.g. inregion II with rcs 4rO ) keep their envelopes when the
nebula vanishes.
5. Conclusions
We investigated an isothermal core–envelope structure inradial symmetry as a model for a planet embedded in a
nebula around a star. The underlying system of differential
equations was transformed into an autonomous one, so wewere able to represent a great variety of different solutions,i.e. radial proﬁles for mass and density, in a twodimensional phase portrait with the trajectories as thesolutions. This phase portrait is dissipative and contains apoint attractor, so all trajectories tend towards the same
limit. Thus all radial proﬁles exhibit an asymptotic‘‘singular proﬁle’’ for r ! 1, which has the form r /1=r2 or rather M / r.
The transformed variables represent two kinds ofgravitational potentials, one is the global potential j andthe other one is the local potential c describing self-gravity.
The properties of the interplay between j and c contain alot of useful information.

An important result is the critical gas density above
which self-gravity effects overrule the global gravitational
potential. This happens when the local density is one thirdof the mean density inside the actual sphere. At the coresurface this corresponds to one third of the mean coredensity and therefore gives us the density value rO of thepoint of critical mass. If the gas density at the core is higher
than rO , self-gravity starts to inﬂuence the envelope in a
global way. For the core mass we found the analogue to rO
in the minimal core mass M min
c . Cores more massive than
M min
have
a
signiﬁcant
inﬂuence
on the envelope.
c
Given a boundary for the region of gravitational
inﬂuence in form of the Hill radius, we found the CCM
M O , marking the point where the core potential inﬂuences
the whole Hill sphere. In the evolution of the protoplanetary disk, the value of M O decides whether or not a planet
is able to retain its gaseous envelope as the disk gas
vanishes. Planets with cores above M O will keep theirenvelope, while others will not (if not compact throughself-gravity). The gas density analogue to M O is the borderline between region II and IV. As this line is not vertical, cf.
Fig. 1, the critical gas density at the core surface, which
ensures the transition to a compact state (i.e. entering
region II), depends on the core mass. The bigger the core
the smaller is the needed gas density (the core ‘‘supports’’
the compactiﬁcation process).
The formulas derived for rO and M O describe the point ofcritical mass and therefore the important point where there
are no more static states of the envelope in general and no
more core growth driven, quasi-static evolution is possible.
The investigation of the behaviour of mass proﬁles MðrÞ atthe Hill radius led us to an explicit expression rcs ðM c Þ of theborder line between ‘‘asteroids’’ in region III and ‘‘protoplanets’’ with extensive envelopes in region IV. For vanishingcore mass, this relation gives the critical central density of agas sphere, above which extensive envelopes are built up.
With the help of the asymptotic solutions we were able todetermine the mean values of envelope mass and outerdensity for region IV. These quantities are quite constant
there, because for r ! rH all proﬁles reach the unique
asymptotic solution. Minor variations in the quantities at
rH are due to damped oscillations of the proﬁles near theasymptotic solution. These oscillations were investigatedthrough linearization of the vicinity of the attractor in thephase portrait, including the calculation of the eigenvalues
and the determination of the typical wavelength.
Acknowledgements
I want to thank Günther Wuchterl for the funny andinteresting discussions and for showing me, what is really
important when doing science; Bojan Pečnik for merging
his thoughts with mine in times when something about our‘‘island’’ seemed pretty strange; Dr. Herlt from the ITP of
the physics department at the University of Jena, for a
fruitful discussion about point transformations; Peter H.
Richter, who encouraged me to write this treatise, although

ARTICLE IN PRESS
J. Schönke / Planetary and Space Science 55 (2007) 1299–1309

I might do something different now and who gave me some
useful hints for the discussion of the dynamical system;
Werner M. Tscharnuter for some remarks, concerning
connections to the homology invariants; the two unknownreferees for helping to improve the clarity of the paper.
References
Broeg, C., 2006. Dissertation, Universität Jena.
Hayashi, C., Nakazawa, K., Nakagawa, Y., 1985. Protostars and Planets
II, p. 1100.
Ikoma, M., Emori, H., Nakazawa, K., 2001. ApJ 553, 999.

1309

Kippenhahn, R., Weigert, A., 1990. Stellar Structure and Evolution.
Springer, Berlin.
Mizuno, H., 1980. Prog. Theor. Phys. 64, 544.
Papaloizou, J.C.B., Terquem, C., 1999. ApJ 521, 823.
Pečnik, B., 2001. Master Thesis, University of Zagreb.
Pečnik, B., Wuchterl, G., 2005. A&A 440, 1183.
Polyanin, A.D., Zaitsev, V.F., 2003. Exact solutions for ODE’s, second ed.
Chapman & Hall/CRC, London/Boca Raton.
Schönke, J., 2005. Diploma Thesis, Universität Jena.
Stevenson, D.J., 1982. PSS 30, 755.
Wuchterl, G., 1991. Icarus 91, 39.
Wuchterl, G., 1993. Icarus 106, 323.
Wuchterl, G., Guillot, T., Lissauer, J.J., 2000. Protostars and Planets IV,
1081.