- Email: [email protected]
On the thickness and evolution of the dust layer during the formation of the solar system
Chin.Astron.Astrophys.
8 (1984)
97-104
Pergamon
Act.Astron.Sin.24 (198;?7351-361
Press.
Printed in Great Britain 0275-1062/84$10.00+,00
ON THE THICKNESS AND EVOLUTION OF THE DUST LAYER DURING THE FORMATION OF THE SOLAR SYSTEM
YUE Zeng-yuan, ZHANG Bin Department Received
of Geophysics,
1982 September
Beijing
University
6
ABSTRACT We discuss certain dynamical processes during the final stage of the sinking of the dust layer. We supposed that turbulance gave rise to a state of slow sinking (quasi-equilibrium) and evaluated the critical thickness at the onset of gravitational instability in the radial direction. We gave a precise numerical relation between 3 length-scales:c
: h* : hT = 0.02107
: 0.1592
: 1,
the first being the mean height of the dust particles at the onset of radial instability, the second being that value of the half-thickness and of the height at which the self-gravity of the dust layer is equal to the solar z-component, and the last being the longest wavelength at the onset of ring instability. We also calculated the time required for the formation of rings and found it to be far shorter than the sinking time.
1.
INTRODUCTION
Of the many questions involved in the origin and evolution of the solar system, one of great concern is that of the thickness of the dust layer, [l]. There are a number of difficulties here: (i) How should we calculate the random velocity of the dust particles? It includes not only the thermal motion, but, more importantly, we must allow for the random motion caused by turbulance in the gas, and both types of random motion depend on the size of the particles. Moreover, the velocity of the turbulant motion of the gas depends on the radial distance in a manner which is difficult to assess with any accuracy. (ii) How should we calculate the process of sticking of dust particles? It depends not only on the relative motion among the particles (hence the probability of collision), but also on the probability of sticking, which itself depends on the relative motion and the structure of the (iii) The settling of the particles. particles being a gradual, continuous process, the meaning of the term “thickness of the dust 1 aye?’ is none too clear. Let us first define a “critical thickness”: when the dust layer has settled down to this then gravitational instability in thickness, the radial direction will begin to appear.
In the usual theory of the sinking of the dust particles 11, 21, the self-gravitation of the nebular disk (a mixture of gas and dust) is neglected, because it is much smaller than the z-component of the solar gravitation. But as the layer gets thinner, the effect of this solar component gets less, while, the surface density e of the dust layer being constant, the se1 H-gravitation at the outer edge is not decreased. Therefore, when the layer has thinned to a certain value, its self-gravitation will be of the same order of magnitude as the solar z-component and become dominant. We shall call this stage the final stage of sedimentation. In a nebular disk of such huge sizes, the Reynolds number of the gas is very high. If we take the Sun-Barth distance as the characteristic length, the Kepler velocity at the Earth’s distance as the characteristic velocity and note that, here, the gas density is P 2, 1 .2 x 10T9 g/cm3, toe Hficient is PI, ant the viscosity 11% log/cm.sec, then we estimate the Reynolds number to be
(1’ At such a high Reynolds number, turbulance will inevitably be present everywhere in the nebular disk. The random turbulent motion together
98
YUE and ZHANG
with the thermal molecular motion will balance the z-component of the solar gravitation and maintain the disk at some considerable thickness. (As an example, at the Earth's neighbourhood the half-thickness will be 9.7 x lOI cm, [l], that is, only one order of magnitude less than the Sun-Earth distance). In the usual theoretical models of sinking of dust particles, we suppose particles of radii of about 10e5 cm are mixed in the gaseous disk. Because the dust particles cannot possibly have random velocities as high as the gas molecules, they cannot balance the a-component of the solar gravitation, therefore they sink to form a thin dust layer. The time-scale of sinking is about 104yr (in the region of terrestrial planets), [l]. When the dust layer has flattened so much that its self-gravitation dominates over the solar z-component, the total gravitational effect will have become far less than that at the initial stage of the sinking. Since the dust layer is now very thin, the a-component of its self-gravitationat the outer edge of the dust layer can be well approximated by the following formula (obtained on applying Grass's Theorem),
that, at every instant t during this stage, the equivalent pressure of the random motion imparted to the dust particles by turbulence is in approximate equilibrium with gravity, while the coalescence of the dust particles gradually decreases the random velocity (the equivalent pressure), so the dust layer keeps on getting thinner.
2
DYNAMICS DURING THE LATE STAGE OF SINKING. CALCULATION OF CRITICAL THICKNESS
For simplicity, we first consider the case where the half-thickness of the dust layer is already far smaller than he. In this case, the z-component of the solar gravity is negligibly small compared to the selfgravity of the dust layer. (To consider simultaneously the solar component and the self-gravity involves no difficulties of principle, and a brief discussion will be given in the Appendix). Let ap be the random velocity given to the dust particles by the turbulence. The mathematical expression for the condition of quasiequilibrium in the late stage of sinking is then
(2)
Fself,a = -2nG op,
(6)
is the surface density of the dust the solar z-component is (3) where Q
_
( )“’ *
I’
(4)
where o is the volume density of dust particlg s, the first term represents the equivalent pressure of the random motion and the second term is the z-component of the self-gravity of the dust layer, obtainable from Gauss's Theorem. Under the boundary conditions
is precisely the circular velocity of the particle under the solar gravity. The condition F,,lf a = Fez defines a height h,, 1, -
the solution of (6) is 2nCoplQ'.
(5)
Its meaning is: when the half-thickness of the dust layer is less than h*, the selfgravitation will exceed the solar a-component and becomes dominant. For the neighbourhood of the Earth, taking o % log/cm2 (according to Ref. [l], Table 8.57, we calculate that h*Q7 x 107 cm. This is 4 orders of magnitude lower than the half-thickness of the nebular disk. In other words, during the late stage of sinking, (when half-thickness
PP -
,+,scch'f. 20
where
Z.----J& aCup
(8)
If, instead of the model of continuous medium, we regard the dust particles as discrete mass points, then the same result (7) - (8) will also be given by statistical mechanics methods. We assume the gas turbulence to be isotropic, hence the random velocity ap it causes can also be regarded as approximately isotropic. For the same turbulence environment, the coagulation of dust particles will gradually decrease ap and this has two consequences: one is the continual thinning of the dust layer, the other is the
99
Evolution of bust Layer
continuous weakening of the resistance read against instability in the radial direction. 141--2iz-9’. When aP ts decreased to some critical value (18) xd 4 radial (annular) instability will set in. "Fp"> rom studies of instabilities in disk galaxies Squaring both sides of (16) and using(l7) 141, we know that, for an infinitely thin disk, the limiting condition for stability is and (18), we have
(9)
21'((1+2+~(1ZI),/In2)'ln2
Here a is the velocity dispersion, o the surface density of the stars, and K is defined by
1.
But Ml
+ 2:i,(lzl),/1at> -a+&, G
(101
hence, on introducing the non-dimensional parameter
This result is entirely applicable to dust layers, all we have to-do-is to replace a by the random velocity ap, a by up, and set
6 - 2&(lZ\)C/Ln2~
Eq. (19) is simplified into the algebraic equation
(11)
C--P
for the case of Kepler motion. To allow for the attenuating effect on radial perturbations due to a finite thickness, Eq. (9) should be modified to read
5(1+-C)- 1,
(22)
with its famous solution of golden section,
5 - +(JT-
90
(21)
1)- 0.6180.
(12)
',l.-1. XCO,
From (13), (18), (21) and (22), we find that the gravity reduction factor at the critical thickness is also precisely equal to the golden section.
where I?is the gravitational Yonstant" after allowing for the effect of finite thickness. It is related to G by, [5, 61, G - G/(1 +
l~l(lzl)/w.
d
(13)
-w-P
G Here ]kl is the radial wave number, IzI is the mean distance of the dust particles to the central plane. For the profile (7), we have (lzl)- Z&2.
5 -+(2/-i - l>- 0.61110. (24)
Noting that
(25)
(14)
Substituting (14) in (8) gives Ilp - (xGop(lZI)/1a2)~.
1
1+5
we have, from (21) and (25), (15) (lzl>,- y
Substituting (13) and (15) in (12) gives the equation to be satisfied by the "critical half-thickness" < IzI>C,
6 - zr+
2
c'& - 0.02107&.(26)
, where (27)
Since we are considering the critical thickness at the onset of gravitational instability, we take Ikl to be the wave number at this point. For an infinitely thin disk, this wave number is, [4],
24r - -,
P
*Go,
(17)
Where we have used Eq. (11). Allowing for the effect of finite thickness, we change this to
has a clear physical meaning, namely, the upper bound to the wavelength for instability in radial (annular) perturbations. (See Ref. [l]. In fact, this critical wavelength was first introduced by Toomre [7] in his study on the stability of disk galaxies). From (5) and (27), we immediately obtain the relation between the half-thicknessh* at which the self-gravity equals the solar a-component and Xt:
YUE
100
h, -
-
l
IT-
0.1592 Ir.
and
ZHANG
that happens is that the four numbers on the right side of (33) will all be multiplied by the same correction factor. Substituting (26) in (15) and using (21), we obtain the critical random veloc:ity at the onset of radial gravitational instability,
(28)
2%
Comparing (26) and (28), we have (29)
{lZl>c~ he.
This shows that, when the dust layer has flattened to the critical thickness, its self-gravity is indeed much greater than the solar z-component, and so we can neglect the latter in a first approximation. In Ref. [l], the Roche density formula was used,
PR
(34) For the surface density of the dust layer in the primitive solar nebula u we took the value from Ref. [l] (Table !:S thereof), and for the coefficient B1 in the expression for the primitive mass,
- 3.5349 - 3.534 $. (35)
M',-,%Uo and a corresponding half-thickness h, was derived from giving hR=-.
1
we took the typical value of 1.2 (Ref. [l]). We then calculated the values of ht, h,, < 1ZI >c, h, and a at various locations in the solar nebularPcdisk. The results are shown in TABLE 1 From these we derive two conclusions: (i) as soon as the random velocity given to the dust particles by the gas turbulence reaches the order of 10 cm/s! a stage of slow sinking (quasi-equilibrium)will obtain during the late stage of sinking, (ii) Before the dust layer has decreased to the thickness h r, gravitational instability in the radial direction (tending to form rings) will already have appeared in the dust layer, for < lal'c is about one order of magnitude greater than h,.
Gal
7.068 Q' This was used in Ref. [l] as an estimate for the thickness of the dust layer. Comparing (27) and (31), we have &== 0.003584 Ir. hfi = l 7.068X 4%'
(32)
The above results (26), (28), (32) can he written in the combined form h,:((Z(),:h,:lr
0.003584:0.02107:0.1592:1 (33)
-
It should be emphasized that, in deriving these precise ratios, the only assumption used is the assumption of turbulence causing a slow sinking during the late stage. The derivation does not depend on any of the factors such as the mechanism of sticking, the sticking probability, the precise magnitude of the random velocity of the dust particles etc., factors which are difficult to assess with any accuracy. Even if we consider the difference between the original mass of the Sun Mm’ and its present mass M,, (according to [I], 1’4~’ = B1M,,, 1~ Bl< 2). the above ratios will not be changed, because all
3
DYNAMICAL PROCESS OF RINGS
FORMATION OF
We now consider the process of formation of rings through instahility after the dust layer reaches the critical thickness. At this time, the layer is already very thin (near the Earth, < ]z 1>c is only 9 x lo6 cm, see TABLE l., moreover, the dust layer keeps on getting thinner), hence, for a first approximation, we shall use the model of an indefinitely thin disk in discussing the formation of rings in the horizontal
TABLE 1 Vaiues of Parameters at Various Points in the Solar Nebula (Bl= 1.2) Locatio Vercury Venus Earth -------_ Mcm) 1.66X IO’ 3.U9XlO’ 4.23X10’
6.53X10’
8.7X10’
2.64X10’
4.92X10’
6.3IxlO’
1.04X10’
3.5OXIO‘
6.51X10’
9.01X10‘
1.38~10’
1.4X10’ ---_ 1.8~10’
2..s7X10’
5.32X10’
7.40X10’
1.13X10’
3.1X10’
-------
Mcm) -viz----------(cm) P---PMcm) .,.(cm/---------WC)
Mars Jupiter Saturn Uranus Neptune
9.72
7.08
6.03
4.89
10.4
1.84X10’”
3.6X10”
5.8XlO’O
2.9X10’
5.7X10’
9.2X10’
3.88~10’
7.6~10’
1.2X10’
6.6X10’
1.3X10‘
2.08X10’
8.8
--
6.04
5.0
Evolution
of
direction. Since the random velocity of the gas (molecular thermal velocity and turbulence velocity) is far greater than that of the dust particles, if we use the usual values (Ref. [l] Tables 8.5 and 8.1) for the surface density and temperature of the gas, then its random velocity will be sufficient to suppress any gravitational instability of the gas. In other words, when gravitational instability appears in the dust layer, the gas still maintains its original state of motion. Thus, the analysis of gravitational instability of the dust layer differs from the usual case of disk galaxies in that, here, we must consider the effect of the gas on the dust layer. There are two effects: one is that,pushedalong by the gas turbulence, the dust particles will acquire a certain amount of random velocity; the other effect is that, when particle clusters move macroscopically ralative to the gas, they will experience a resistance. equations Thus, the linearized for small perturbations in the dust layer in the model of indefinitely thin disk are
au,, 6t
+
2
BUP, 1 a Qw + 7 5 (roPp1) +
I
+ Q2
-
ZQv, ___&&el-BvI__,, up0 ar
Bv,+Q~++-L[&?!%+g
av, -
Up.68
-
0,
(36)
1,
(37)
ar
-
r op. a8
at
I
4n G+,&(r).
~,(38)
(39)
Suffixes 0 and 1 refer to the base state and perturbation, respectively. ~1, ~1 are the radial and tangential perturbations in the velocity of the particle cluster, apl is the perturbation in the surface density of the dust 1 ayer , and Vl is the perturbation in the gravitational potential. Under the linear approximation, ~1 and ~1 are also the relative velocity between the particle cluster and the gas, for in the basic state, the dust layer and the gas move together around the Sun in Kepler motion, In (37) and (38)) the term containing u represents the resistance on the particle cluster by the gas, the relation between the coefficient a and the gas density thermal velocity of the gas V , radius of pg* dust particle h and density of dugt particle 6 being, according to [2], (49)
Dust
Y,
Layer
O,(r)exp(rr + i \‘kdr)~ \
-
P,(r)~p(*t+ i 1’kdr).]
v, -
y>O means a growing mode, y< 0 a decaying mode ; a purely imaginary y means neutral oscillation, a complex y means oscillations with increasing or decreasing amplitude according as Re y is positive or negative. Substituting (41) in (36) - (39), and under the fast varying phase approximation (which is an extremely good approximation here since the characteristic scale-length of the dust layer At is far less than the radial size of the nebular disk), we obtain the following setsof equations:
r+a
0
the
up, - 6, u, -
ring-shaped
perturbation
~(7 + a)’ + Q’r + 4’ a: -
~,(r)exp(rr
+
i
\‘kdr),
ap -
+
a)
0
-
0.
(43)
(44)
to estimate
the
scale of the Also we shall take ~~~~~i~~ ~~er~~~‘number’at which gravitational instability first occurs,
lkl -2kr-5.
~(7
+
a)’
time
0
Substituting +
(44) 0’7 -
and (45)
ZQ’(7 + a> -
in 0.
(43)
gives (46)
Let 7 7’-,
(46)
,8-z. Q then
1(9 + 8)’ + 9 (41)
2nGuP0J (1 -w
It is easily seen that when the resistance coefficient CY= 0, Eq. (43) reverts to the usual dispersion relation for ring-shaped density waves in gaseous disk models. When u #O, Eq. (43) is a cubic equation in y. Since, after the critical thickness is the dust particles keep on growing reached, through coagulation,ap continues to decrease This circumrstance enables further from apt. us to use the very good approximation, [1] ,
Eq.
(r)~xp(rt+ i\'kdr)s
Q/2
For a non-trivial solution, the determinant Hence we of coefficients must be zero. obtain the dispersion relation
Q
Consider
101
or 9’ + 2&a+
(8’ -
(47)
becomes 2(9 + B) 1)7j -
28 -
0.
(48)
0. (49)
YDE
102
and
Let
(54) are real and positive, and the expression gives the required root directly;
?j-x--
Eq.
ZHANG
;
(49)
ry.
then
(50) becomes
For discussing the growth rate of we are interested only in the perturbations, positive roots of (49)) hence also only in the positive roots of (51). Since (52)
~%-(~~+l)
(the case of small drag), then in (54) are complex conjegates, write it as
The boundary the positive
4(f)-(f>‘-
(53) it follows that positive root,
(51)
has
a unique
+ ($)‘]”
the the
two cases equation
is
given
~-~~‘==~(~+ll~-l)-o
by
(56)
real that
x - [- + +&y
between root of
the two terms and we can re-
is,
-t(54)
8’ _
+
(-11
+ dlzs>
(57)
~0.09017,
or There
are
two cases:
fil
B-B*-
i.e.
of
where 6 is seen that (i)
(ii) (iii)
large
drag),
defined
at
The two expressions 6= B, Small
(-
11 + Jlzs,]‘”
(58)
5 0.3003.
Combining (SO), (54) and (58)) we arrive at the final form of the relation between the growth rate of the drag coefficient:
(4)‘> I$[’ (the case
[$
drag
Large drag nQ22/@ for
limit:
then
(47).
are
n+l
asymptote: B>>l
both
It
is
terms
easily
continuous
as B+O
of
at
(61)
(62)
The n rliB curve according to (59) and (60) is shown in Fig. 1. For the gas density Pg and temperature and the dust particle radius b and TgJ . density 6, we took the data from Ref. [l] (taking 131=M,‘/Me=1.2 in calculating Since the main constituent of the gas P,) . is hydrogen, we can use the following
( <__ Fig.1
Solid curve shows the precise between n and B. Dashed curve asymptotic relation for large
relation is the B
Evolution
of
expression to calculate the thermal velocity of the gas molecules: V, -
1.24 X 10’2/?~&=.
(63)
The values of a at various locations were then found according to (40). In order to calculate the original angular velocity Q1 of the matter that eventually formed a particular planet, we must aIlow for the difference in mass between the primitive and the present sun. Specifically, for the planet with present solar distance r and r’” angular velocity a, the “mean distance and the “mean angular velocity Q”’ of its original matter are given by
103
Dust Layer
We repeat that, during the formation of the ring, the dust particles are continuing to grow, decreasing the drag coefficient c1 all the time, hence the actual growth rate is higher than shown in the Table, and the timescale To shown must also be regarded as an upper bound to the time of formation of the rings. The above calculation shows that the time required for the dust layer to break up into ring is far shorter than the time of formation (sinking) of the dust layer, which is of the order of lo4 yr. The formation of the rings provide the condition for the formation of planetesimals, and the planets were then formed following further collision and combination of the planetesimals, [I].
(64) 4
Q’r” - g,,_ From these expressions respectively Kepler motion and the conservation angular momentum, we have r* -
DISCUSSION
(65) of of
the
(66)
r/B,.
4’ - &2
(67)
Using the present values of 0 of the planets and taking Bl=l.Z, we calculated n’, hence the non-dimensional resistance coefficient 8. Then, using (60), (in fact we can use the asymptotic expression (62)), we calculated the non-dimensional growth rate n = y/G’, hence the growth rate y, and finally, the time-scale Te of the formation of ring (i.e. the time required for an e-fold These results are increase in density). shown in TABLE 2.
TABLE 2
Calculated
Values
We have evaluated the mean height < IzI >c of the dust particles at the onset of gravitational instability in the dust layers and hence proceeded to an analysis of the dynamical process of the formation of rings and an evaluation of the time-scale involved. In our calculation, we used the numerical data given in Ref. [l]. It should be pointed that the density and temperature profiles in the solar nebular disk are far from being completely solved problems. Further research on these problems should include the physical processes during the early stage of the evolution of the solar system, and the results obtained will influence the subsequent formation of planetesimals and planets. Different profiles will not only change the quantitative results. they. mav even alter the overall picture of the evolutional process.
at Various klars P--P-
Points
in the
Solar Uranus
Nebula
Jupiter
Saturn
Neptune
I’ 4.05X10-”
1.66X10-”
8.15x10-”
2.13X10-”
4.6X10-‘*
3.8X10-‘* --
3.1X10-”
YUE
104
and
ZHANG
REFERENCES [l] [Z] [B]
[4] [S] 161 171
“Taiyangxi Yanhauxue” (“Evolution of the Solar DA1 Wen-sai, System”) Vol. l., (lY7Y), Shanghai Kexue Jishu Chubanshe. Safronov, V.S., Evolution of the Protoplanetary Cloud and Formation of the Earth and the Planets, Jerusalem, 1972. Lin, C.C. and Segel, L.A., Mathematics Applied to Deterministic Prohlems in the Natural Sciences, Macmillan Publishing Co., Inc., 1974. Spiral Structure” Chinese translation C.C. Lin “Theory of Galactic by HU Wen-rui and HAN Nian-guo, (lY77), Kexue Chubanshe. Yue, Z. Y, , Geophys. Astrophys. FI!A?:~ Dyn. 20( lY82), 1 - 110. YUE Zeng-yuan, Ch7kz. Astrun. Astrwphgs. 7(1983) 186 - 194. Chinese original in Act. Astrophys. Sin. 3(1983) 103 - 112. Toomre, A., A~tt?~ph,q. ,7. 139(1964), 1217.
where
APPEND1X
E= o’l(2nC,o,0). At the critical thickness, the mean height of the dust particles < IzI >c is much less than effect exceeds h *, so that the self-gravity In the paper, the the solar z-component. latter was neglected for simplicity, and the quasi-equilihrium between self-gravity and the random velocity caused by turbulence was We now include the solar considered, and consider the quasiz-component, equilibrium hetween all three: -*i
- Q’Z -
P, dr
where I)’ = s
Y’
2nG
J:,+
= 0.
(1)
(7) condit
ions
are (8)
Let
iEd.?_ 4
then 9
’
(9)
(6)
becomes
ridy = 8+e-,,
(10)
or
(11)
(2)
.
Differentiating more, we get -#s d,! ’ dz’
The houndary
(1)
with
respect
Hence +z - 67 - e-J + c,.
to z once
_ I)*+ zncp,.
(3)
5 = &&=$ where P o is the density of dust particles in the gentral plane. Eq. (3) now becomes zd’y = 8 + e-y.
conditions
(4)
(6)
then
give
hence +
(5)
(12) (8)
(13)
C,=l,
Let y - ‘+J?p/P*),
The boundary
= d/+@y
+ 1 - ,-,)‘,a.
(14)
Using (8) once more, we obtain the solution including both self-gravity and solar a-component in the integral form 5 - +j:
(8, + pt =-‘)“”
(15)
8 (1984)
97-104
Pergamon
Act.Astron.Sin.24 (198;?7351-361
Press.
Printed in Great Britain 0275-1062/84$10.00+,00
ON THE THICKNESS AND EVOLUTION OF THE DUST LAYER DURING THE FORMATION OF THE SOLAR SYSTEM
YUE Zeng-yuan, ZHANG Bin Department Received
of Geophysics,
1982 September
Beijing
University
6
ABSTRACT We discuss certain dynamical processes during the final stage of the sinking of the dust layer. We supposed that turbulance gave rise to a state of slow sinking (quasi-equilibrium) and evaluated the critical thickness at the onset of gravitational instability in the radial direction. We gave a precise numerical relation between 3 length-scales:
: h* : hT = 0.02107
: 0.1592
: 1,
the first being the mean height of the dust particles at the onset of radial instability, the second being that value of the half-thickness and of the height at which the self-gravity of the dust layer is equal to the solar z-component, and the last being the longest wavelength at the onset of ring instability. We also calculated the time required for the formation of rings and found it to be far shorter than the sinking time.
1.
INTRODUCTION
Of the many questions involved in the origin and evolution of the solar system, one of great concern is that of the thickness of the dust layer, [l]. There are a number of difficulties here: (i) How should we calculate the random velocity of the dust particles? It includes not only the thermal motion, but, more importantly, we must allow for the random motion caused by turbulance in the gas, and both types of random motion depend on the size of the particles. Moreover, the velocity of the turbulant motion of the gas depends on the radial distance in a manner which is difficult to assess with any accuracy. (ii) How should we calculate the process of sticking of dust particles? It depends not only on the relative motion among the particles (hence the probability of collision), but also on the probability of sticking, which itself depends on the relative motion and the structure of the (iii) The settling of the particles. particles being a gradual, continuous process, the meaning of the term “thickness of the dust 1 aye?’ is none too clear. Let us first define a “critical thickness”: when the dust layer has settled down to this then gravitational instability in thickness, the radial direction will begin to appear.
In the usual theory of the sinking of the dust particles 11, 21, the self-gravitation of the nebular disk (a mixture of gas and dust) is neglected, because it is much smaller than the z-component of the solar gravitation. But as the layer gets thinner, the effect of this solar component gets less, while, the surface density e of the dust layer being constant, the se1 H-gravitation at the outer edge is not decreased. Therefore, when the layer has thinned to a certain value, its self-gravitation will be of the same order of magnitude as the solar z-component and become dominant. We shall call this stage the final stage of sedimentation. In a nebular disk of such huge sizes, the Reynolds number of the gas is very high. If we take the Sun-Barth distance as the characteristic length, the Kepler velocity at the Earth’s distance as the characteristic velocity and note that, here, the gas density is P 2, 1 .2 x 10T9 g/cm3, toe Hficient is PI, ant the viscosity 11% log/cm.sec, then we estimate the Reynolds number to be
(1’ At such a high Reynolds number, turbulance will inevitably be present everywhere in the nebular disk. The random turbulent motion together
98
YUE and ZHANG
with the thermal molecular motion will balance the z-component of the solar gravitation and maintain the disk at some considerable thickness. (As an example, at the Earth's neighbourhood the half-thickness will be 9.7 x lOI cm, [l], that is, only one order of magnitude less than the Sun-Earth distance). In the usual theoretical models of sinking of dust particles, we suppose particles of radii of about 10e5 cm are mixed in the gaseous disk. Because the dust particles cannot possibly have random velocities as high as the gas molecules, they cannot balance the a-component of the solar gravitation, therefore they sink to form a thin dust layer. The time-scale of sinking is about 104yr (in the region of terrestrial planets), [l]. When the dust layer has flattened so much that its self-gravitation dominates over the solar z-component, the total gravitational effect will have become far less than that at the initial stage of the sinking. Since the dust layer is now very thin, the a-component of its self-gravitationat the outer edge of the dust layer can be well approximated by the following formula (obtained on applying Grass's Theorem),
that, at every instant t during this stage, the equivalent pressure of the random motion imparted to the dust particles by turbulence is in approximate equilibrium with gravity, while the coalescence of the dust particles gradually decreases the random velocity (the equivalent pressure), so the dust layer keeps on getting thinner.
2
DYNAMICS DURING THE LATE STAGE OF SINKING. CALCULATION OF CRITICAL THICKNESS
For simplicity, we first consider the case where the half-thickness of the dust layer is already far smaller than he. In this case, the z-component of the solar gravity is negligibly small compared to the selfgravity of the dust layer. (To consider simultaneously the solar component and the self-gravity involves no difficulties of principle, and a brief discussion will be given in the Appendix). Let ap be the random velocity given to the dust particles by the turbulence. The mathematical expression for the condition of quasiequilibrium in the late stage of sinking is then
(2)
Fself,a = -2nG op,
(6)
is the surface density of the dust the solar z-component is (3) where Q
_
( )“’ *
I’
(4)
where o is the volume density of dust particlg s, the first term represents the equivalent pressure of the random motion and the second term is the z-component of the self-gravity of the dust layer, obtainable from Gauss's Theorem. Under the boundary conditions
is precisely the circular velocity of the particle under the solar gravity. The condition F,,lf a = Fez defines a height h,, 1, -
the solution of (6) is 2nCoplQ'.
(5)
Its meaning is: when the half-thickness of the dust layer is less than h*, the selfgravitation will exceed the solar a-component and becomes dominant. For the neighbourhood of the Earth, taking o % log/cm2 (according to Ref. [l], Table 8.57, we calculate that h*Q7 x 107 cm. This is 4 orders of magnitude lower than the half-thickness of the nebular disk. In other words, during the late stage of sinking, (when half-thickness
PP -
,+,scch'f. 20
where
Z.----J& aCup
(8)
If, instead of the model of continuous medium, we regard the dust particles as discrete mass points, then the same result (7) - (8) will also be given by statistical mechanics methods. We assume the gas turbulence to be isotropic, hence the random velocity ap it causes can also be regarded as approximately isotropic. For the same turbulence environment, the coagulation of dust particles will gradually decrease ap and this has two consequences: one is the continual thinning of the dust layer, the other is the
99
Evolution of bust Layer
continuous weakening of the resistance read against instability in the radial direction. 141--2iz-9’. When aP ts decreased to some critical value (18) xd 4 radial (annular) instability will set in. "Fp"> rom studies of instabilities in disk galaxies Squaring both sides of (16) and using(l7) 141, we know that, for an infinitely thin disk, the limiting condition for stability is and (18), we have
(9)
21'((1+2+~(1ZI),/In2)'ln2
Here a is the velocity dispersion, o the surface density of the stars, and K is defined by
1.
But Ml
+ 2:i,(lzl),/1at> -a+&, G
(101
hence, on introducing the non-dimensional parameter
This result is entirely applicable to dust layers, all we have to-do-is to replace a by the random velocity ap, a by up, and set
6 - 2&(lZ\)C/Ln2~
Eq. (19) is simplified into the algebraic equation
(11)
C--P
for the case of Kepler motion. To allow for the attenuating effect on radial perturbations due to a finite thickness, Eq. (9) should be modified to read
5(1+-C)- 1,
(22)
with its famous solution of golden section,
5 - +(JT-
90
(21)
1)- 0.6180.
(12)
',l.-1. XCO,
From (13), (18), (21) and (22), we find that the gravity reduction factor at the critical thickness is also precisely equal to the golden section.
where I?is the gravitational Yonstant" after allowing for the effect of finite thickness. It is related to G by, [5, 61, G - G/(1 +
l~l(lzl)/w.
d
(13)
-w-P
G Here ]kl is the radial wave number, IzI is the mean distance of the dust particles to the central plane. For the profile (7), we have (lzl)- Z&2.
5 -+(2/-i - l>- 0.61110. (24)
Noting that
(25)
(14)
Substituting (14) in (8) gives Ilp - (xGop(lZI)/1a2)~.
1
1+5
we have, from (21) and (25), (15) (lzl>,- y
Substituting (13) and (15) in (12) gives the equation to be satisfied by the "critical half-thickness" < IzI>C,
6 - zr+
2
c'& - 0.02107&.(26)
, where (27)
Since we are considering the critical thickness at the onset of gravitational instability, we take Ikl to be the wave number at this point. For an infinitely thin disk, this wave number is, [4],
24r - -,
P
*Go,
(17)
Where we have used Eq. (11). Allowing for the effect of finite thickness, we change this to
has a clear physical meaning, namely, the upper bound to the wavelength for instability in radial (annular) perturbations. (See Ref. [l]. In fact, this critical wavelength was first introduced by Toomre [7] in his study on the stability of disk galaxies). From (5) and (27), we immediately obtain the relation between the half-thicknessh* at which the self-gravity equals the solar a-component and Xt:
YUE
100
h, -
-
l
IT-
0.1592 Ir.
and
ZHANG
that happens is that the four numbers on the right side of (33) will all be multiplied by the same correction factor. Substituting (26) in (15) and using (21), we obtain the critical random veloc:ity at the onset of radial gravitational instability,
(28)
2%
Comparing (26) and (28), we have (29)
{lZl>c~ he.
This shows that, when the dust layer has flattened to the critical thickness, its self-gravity is indeed much greater than the solar z-component, and so we can neglect the latter in a first approximation. In Ref. [l], the Roche density formula was used,
PR
(34) For the surface density of the dust layer in the primitive solar nebula u we took the value from Ref. [l] (Table !:S thereof), and for the coefficient B1 in the expression for the primitive mass,
- 3.5349 - 3.534 $. (35)
M',-,%Uo and a corresponding half-thickness h, was derived from giving hR=-.
1
we took the typical value of 1.2 (Ref. [l]). We then calculated the values of ht, h,, < 1ZI >c, h, and a at various locations in the solar nebularPcdisk. The results are shown in TABLE 1 From these we derive two conclusions: (i) as soon as the random velocity given to the dust particles by the gas turbulence reaches the order of 10 cm/s! a stage of slow sinking (quasi-equilibrium)will obtain during the late stage of sinking, (ii) Before the dust layer has decreased to the thickness h r, gravitational instability in the radial direction (tending to form rings) will already have appeared in the dust layer, for < lal'c is about one order of magnitude greater than h,.
Gal
7.068 Q' This was used in Ref. [l] as an estimate for the thickness of the dust layer. Comparing (27) and (31), we have &== 0.003584 Ir. hfi = l 7.068X 4%'
(32)
The above results (26), (28), (32) can he written in the combined form h,:((Z(),:h,:lr
0.003584:0.02107:0.1592:1 (33)
-
It should be emphasized that, in deriving these precise ratios, the only assumption used is the assumption of turbulence causing a slow sinking during the late stage. The derivation does not depend on any of the factors such as the mechanism of sticking, the sticking probability, the precise magnitude of the random velocity of the dust particles etc., factors which are difficult to assess with any accuracy. Even if we consider the difference between the original mass of the Sun Mm’ and its present mass M,, (according to [I], 1’4~’ = B1M,,, 1~ Bl< 2). the above ratios will not be changed, because all
3
DYNAMICAL PROCESS OF RINGS
FORMATION OF
We now consider the process of formation of rings through instahility after the dust layer reaches the critical thickness. At this time, the layer is already very thin (near the Earth, < ]z 1>c is only 9 x lo6 cm, see TABLE l., moreover, the dust layer keeps on getting thinner), hence, for a first approximation, we shall use the model of an indefinitely thin disk in discussing the formation of rings in the horizontal
TABLE 1 Vaiues of Parameters at Various Points in the Solar Nebula (Bl= 1.2) Locatio Vercury Venus Earth -------_ Mcm) 1.66X IO’ 3.U9XlO’ 4.23X10’
6.53X10’
8.7X10’
2.64X10’
4.92X10’
6.3IxlO’
1.04X10’
3.5OXIO‘
6.51X10’
9.01X10‘
1.38~10’
1.4X10’ ---_ 1.8~10’
2..s7X10’
5.32X10’
7.40X10’
1.13X10’
3.1X10’
-------
Mcm) -viz----------(cm) P---PMcm) .,.(cm/---------WC)
Mars Jupiter Saturn Uranus Neptune
9.72
7.08
6.03
4.89
10.4
1.84X10’”
3.6X10”
5.8XlO’O
2.9X10’
5.7X10’
9.2X10’
3.88~10’
7.6~10’
1.2X10’
6.6X10’
1.3X10‘
2.08X10’
8.8
--
6.04
5.0
Evolution
of
direction. Since the random velocity of the gas (molecular thermal velocity and turbulence velocity) is far greater than that of the dust particles, if we use the usual values (Ref. [l] Tables 8.5 and 8.1) for the surface density and temperature of the gas, then its random velocity will be sufficient to suppress any gravitational instability of the gas. In other words, when gravitational instability appears in the dust layer, the gas still maintains its original state of motion. Thus, the analysis of gravitational instability of the dust layer differs from the usual case of disk galaxies in that, here, we must consider the effect of the gas on the dust layer. There are two effects: one is that,pushedalong by the gas turbulence, the dust particles will acquire a certain amount of random velocity; the other effect is that, when particle clusters move macroscopically ralative to the gas, they will experience a resistance. equations Thus, the linearized for small perturbations in the dust layer in the model of indefinitely thin disk are
au,, 6t
+
2
BUP, 1 a Qw + 7 5 (roPp1) +
I
+ Q2
-
ZQv, ___&&el-BvI__,, up0 ar
Bv,+Q~++-L[&?!%+g
av, -
Up.68
-
0,
(36)
1,
(37)
ar
-
r op. a8
at
I
4n G+,&(r).
~,(38)
(39)
Suffixes 0 and 1 refer to the base state and perturbation, respectively. ~1, ~1 are the radial and tangential perturbations in the velocity of the particle cluster, apl is the perturbation in the surface density of the dust 1 ayer , and Vl is the perturbation in the gravitational potential. Under the linear approximation, ~1 and ~1 are also the relative velocity between the particle cluster and the gas, for in the basic state, the dust layer and the gas move together around the Sun in Kepler motion, In (37) and (38)) the term containing u represents the resistance on the particle cluster by the gas, the relation between the coefficient a and the gas density thermal velocity of the gas V , radius of pg* dust particle h and density of dugt particle 6 being, according to [2], (49)
Dust
Y,
Layer
O,(r)exp(rr + i \‘kdr)~ \
-
P,(r)~p(*t+ i 1’kdr).]
v, -
y>O means a growing mode, y< 0 a decaying mode ; a purely imaginary y means neutral oscillation, a complex y means oscillations with increasing or decreasing amplitude according as Re y is positive or negative. Substituting (41) in (36) - (39), and under the fast varying phase approximation (which is an extremely good approximation here since the characteristic scale-length of the dust layer At is far less than the radial size of the nebular disk), we obtain the following setsof equations:
r+a
0
the
up, - 6, u, -
ring-shaped
perturbation
~(7 + a)’ + Q’r + 4’ a: -
~,(r)exp(rr
+
i
\‘kdr),
ap -
+
a)
0
-
0.
(43)
(44)
to estimate
the
scale of the Also we shall take ~~~~~i~~ ~~er~~~‘number’at which gravitational instability first occurs,
lkl -2kr-5.
~(7
+
a)’
time
0
Substituting +
(44) 0’7 -
and (45)
ZQ’(7 + a> -
in 0.
(43)
gives (46)
Let 7 7’-,
(46)
,8-z. Q then
1(9 + 8)’ + 9 (41)
2nGuP0J (1 -w
It is easily seen that when the resistance coefficient CY= 0, Eq. (43) reverts to the usual dispersion relation for ring-shaped density waves in gaseous disk models. When u #O, Eq. (43) is a cubic equation in y. Since, after the critical thickness is the dust particles keep on growing reached, through coagulation,ap continues to decrease This circumrstance enables further from apt. us to use the very good approximation, [1] ,
Eq.
(r)~xp(rt+ i\'kdr)s
Q/2
For a non-trivial solution, the determinant Hence we of coefficients must be zero. obtain the dispersion relation
Q
Consider
101
or 9’ + 2&a+
(8’ -
(47)
becomes 2(9 + B) 1)7j -
28 -
0.
(48)
0. (49)
YDE
102
and
Let
(54) are real and positive, and the expression gives the required root directly;
?j-x--
Eq.
ZHANG
;
(49)
ry.
then
(50) becomes
For discussing the growth rate of we are interested only in the perturbations, positive roots of (49)) hence also only in the positive roots of (51). Since (52)
~%-(~~+l)
(the case of small drag), then in (54) are complex conjegates, write it as
The boundary the positive
4(f)-(f>‘-
(53) it follows that positive root,
(51)
has
a unique
+ ($)‘]”
the the
two cases equation
is
given
~-~~‘==~(~+ll~-l)-o
by
(56)
real that
x - [- + +&y
between root of
the two terms and we can re-
is,
-t(54)
8’ _
+
(-11
+ dlzs>
(57)
~0.09017,
or There
are
two cases:
fil
B-B*-
i.e.
of
where 6 is seen that (i)
(ii) (iii)
large
drag),
defined
at
The two expressions 6= B, Small
(-
11 + Jlzs,]‘”
(58)
5 0.3003.
Combining (SO), (54) and (58)) we arrive at the final form of the relation between the growth rate of the drag coefficient:
(4)‘> I$[’ (the case
[$
drag
Large drag nQ22/@ for
limit:
then
(47).
are
n+l
asymptote: B>>l
both
It
is
terms
easily
continuous
as B+O
of
at
(61)
(62)
The n rliB curve according to (59) and (60) is shown in Fig. 1. For the gas density Pg and temperature and the dust particle radius b and TgJ . density 6, we took the data from Ref. [l] (taking 131=M,‘/Me=1.2 in calculating Since the main constituent of the gas P,) . is hydrogen, we can use the following
( <__ Fig.1
Solid curve shows the precise between n and B. Dashed curve asymptotic relation for large
relation is the B
Evolution
of
expression to calculate the thermal velocity of the gas molecules: V, -
1.24 X 10’2/?~&=.
(63)
The values of a at various locations were then found according to (40). In order to calculate the original angular velocity Q1 of the matter that eventually formed a particular planet, we must aIlow for the difference in mass between the primitive and the present sun. Specifically, for the planet with present solar distance r and r’” angular velocity a, the “mean distance and the “mean angular velocity Q”’ of its original matter are given by
103
Dust Layer
We repeat that, during the formation of the ring, the dust particles are continuing to grow, decreasing the drag coefficient c1 all the time, hence the actual growth rate is higher than shown in the Table, and the timescale To shown must also be regarded as an upper bound to the time of formation of the rings. The above calculation shows that the time required for the dust layer to break up into ring is far shorter than the time of formation (sinking) of the dust layer, which is of the order of lo4 yr. The formation of the rings provide the condition for the formation of planetesimals, and the planets were then formed following further collision and combination of the planetesimals, [I].
(64) 4
Q’r” - g,,_ From these expressions respectively Kepler motion and the conservation angular momentum, we have r* -
DISCUSSION
(65) of of
the
(66)
r/B,.
4’ - &2
(67)
Using the present values of 0 of the planets and taking Bl=l.Z, we calculated n’, hence the non-dimensional resistance coefficient 8. Then, using (60), (in fact we can use the asymptotic expression (62)), we calculated the non-dimensional growth rate n = y/G’, hence the growth rate y, and finally, the time-scale Te of the formation of ring (i.e. the time required for an e-fold These results are increase in density). shown in TABLE 2.
TABLE 2
Calculated
Values
We have evaluated the mean height < IzI >c of the dust particles at the onset of gravitational instability in the dust layers and hence proceeded to an analysis of the dynamical process of the formation of rings and an evaluation of the time-scale involved. In our calculation, we used the numerical data given in Ref. [l]. It should be pointed that the density and temperature profiles in the solar nebular disk are far from being completely solved problems. Further research on these problems should include the physical processes during the early stage of the evolution of the solar system, and the results obtained will influence the subsequent formation of planetesimals and planets. Different profiles will not only change the quantitative results. they. mav even alter the overall picture of the evolutional process.
at Various klars P--P-
Points
in the
Solar Uranus
Nebula
Jupiter
Saturn
Neptune
I’ 4.05X10-”
1.66X10-”
8.15x10-”
2.13X10-”
4.6X10-‘*
3.8X10-‘* --
3.1X10-”
YUE
104
and
ZHANG
REFERENCES [l] [Z] [B]
[4] [S] 161 171
“Taiyangxi Yanhauxue” (“Evolution of the Solar DA1 Wen-sai, System”) Vol. l., (lY7Y), Shanghai Kexue Jishu Chubanshe. Safronov, V.S., Evolution of the Protoplanetary Cloud and Formation of the Earth and the Planets, Jerusalem, 1972. Lin, C.C. and Segel, L.A., Mathematics Applied to Deterministic Prohlems in the Natural Sciences, Macmillan Publishing Co., Inc., 1974. Spiral Structure” Chinese translation C.C. Lin “Theory of Galactic by HU Wen-rui and HAN Nian-guo, (lY77), Kexue Chubanshe. Yue, Z. Y, , Geophys. Astrophys. FI!A?:~ Dyn. 20( lY82), 1 - 110. YUE Zeng-yuan, Ch7kz. Astrun. Astrwphgs. 7(1983) 186 - 194. Chinese original in Act. Astrophys. Sin. 3(1983) 103 - 112. Toomre, A., A~tt?~ph,q. ,7. 139(1964), 1217.
where
APPEND1X
E= o’l(2nC,o,0). At the critical thickness, the mean height of the dust particles < IzI >c is much less than effect exceeds h *, so that the self-gravity In the paper, the the solar z-component. latter was neglected for simplicity, and the quasi-equilihrium between self-gravity and the random velocity caused by turbulence was We now include the solar considered, and consider the quasiz-component, equilibrium hetween all three: -*i
- Q’Z -
P, dr
where I)’ = s
Y’
2nG
J:,+
= 0.
(1)
(7) condit
ions
are (8)
Let
iEd.?_ 4
then 9
’
(9)
(6)
becomes
ridy = 8+e-,,
(10)
or
(11)
(2)
.
Differentiating more, we get -#s d,! ’ dz’
The houndary
(1)
with
respect
Hence +z - 67 - e-J + c,.
to z once
_ I)*+ zncp,.
(3)
5 = &&=$ where P o is the density of dust particles in the gentral plane. Eq. (3) now becomes zd’y = 8 + e-y.
conditions
(4)
(6)
then
give
hence +
(5)
(12) (8)
(13)
C,=l,
Let y - ‘+J?p/P*),
The boundary
= d/+@y
+ 1 - ,-,)‘,a.
(14)
Using (8) once more, we obtain the solution including both self-gravity and solar a-component in the integral form 5 - +j:
(8, + pt =-‘)“”
(15)