Subsections


Simulating the Capture Theory

In this chapter, the code that was developed and tested in previous chapters is used to perform preliminary investigations into two areas of the Capture Theory. The first section is concerned with the behaviour of quasi-statically collapsing bodies, as it is possible that both the protostar and the newly formed protoplanets could be in this state. The second section begins with a discussion into why quasi-statically collapsing protostars are not conducive to planetary formation. A new type of encounter, which is conducive to planetary formation, is then suggested, and preliminary simulations of this type of encounter are presented.


Quasi-Static Collapse of Spherical Bodies

An investigation into the behaviour of quasi-static collapsing bodies is particularly relevant to the Capture Theory, as both the interacting protostar and the newly formed protoplanets could potentially be in this state. If a body is to remain in equilibrium throughout its collapse, it must be capable of contracting at such a rate that it is able to generate the energy lost from radiation, and also enough additional energy to raise the temperature sufficiently to maintain the pressure balance. An upper limit on the contraction rate is set by the free-fall velocity, so that if the luminosity is greater than the rate of thermal energy generated by free-fall collapse, the body cannot possibly remain in equilibrium.


Low Mass Bodies

It is simpler to first consider low mass bodies in quasi-static collapse, as the complexities introduced by hydrogen dissociation at $2000K$, and very low opacities occurring around $1500K$, are not encountered until much later times during the collapse. In each of these simulations, an equilibrium polytropic distribution of particles is used as the starting configuration, see Section 4.2. Several equilibrium structures, with various number of particles and ratio of specific heats, were constructed, and saved for future use. These are easily scaled to any dimensions by noting that:
\begin{displaymath} T\propto\frac{M}{R} \end{displaymath} (5.1)

as derived from Equation 1.4 in Section 1.3.1.

In the limit of an infinite number of particles, it is believed that the code, as presented in previous chapters, could reliably model the quasi-static collapse under investigation, given the current assumptions concerning the state of the material5.1. However, for a finite number of particles, additional considerations must be taken into account. Suppose, incorrectly so, that as the body contracts, it retains its original polytropic profile. Under this assumption, it is straightforward to show that the surface velocity is constant throughout the collapse: Since the body is in quasi-static equilibrium, its energy is proportional to the gravitational potential energy:

\begin{displaymath} U\propto -\frac{M^{2}}{R} \end{displaymath} (5.2)

leading to
\begin{displaymath} \frac{dU}{dR}\propto\frac{M^{2}}{R^{2}} \end{displaymath} (5.3)

The rate of energy change due to radiation is proportional to:
\begin{displaymath} \frac{dU}{dt}\propto -R^{2}T^{4} \end{displaymath} (5.4)

Substituting for $T$ from Equation 5.1 and solving for $dR/dt$ leads to:
\begin{displaymath} \frac{dR}{dt}\propto -M^{2} \end{displaymath} (5.5)

For the current scenario, a numerical value for the rate of collapse can be calculated by solving the above equations more completely. However, in order to be able to perform this calculation, the effective surface temperature of the body must be known. The theoretical surface temperature of a polytrope is zero. A more appropriate temperature is that of the average temperature of the material one optical depth into the body. Indeed, the effective surface temperature of a body is usually defined as the black body temperature which corresponds to the measured luminosity and radius determined from experiment, and since the majority of the radiation leaving a body originates from one optical depth into it, the average temperature of that region is a good first approximation to the effective surface temperature.

When using SPH and this implementation of radiation transfer to model the quasi-static collapse of opaque bodies, the effective surface temperature of the body is equal to the temperature of the surface SPH particles. This is because the temperature across an isolated SPH particle is defined as being constant, see Section 3.2.7. In the limit of an infinite number of particles, the effective surface temperature would be correctly approximated by the average temperature of the particles one optical depth into the body, as this is where the majority of the power packets that leave the system would originate. However, for a finite number of particles, an incorrect value of the effective surface temperature is obtained: Consider the polytropic profiles presented in Figure 4.7 in Chapter 4. As demonstrated by the horizontal location of the arrows in this figure, the extent of the experimental radius of the polytrope depends on the number of particles used to model it. This leads to the temperature of surface particles in low resolution simulations being higher than that in high resolution simulations. The result of this is that the effective surface temperature, and hence the luminosity, is higher, leading to a larger rate of collapse.

A full solution for negating this dependence has not been investigated, however, a simple correction can be obtained by comparing the physical processes that control the luminosity of a body, and hence its rate of collapse, with the experimental processes. Experimentally, the surface leafs radiate energy freely to other parts of the system. During the radiation routine, these leafs may cool to a more realistic temperature, of the order of the theoretical surface temperature one optical depth into the body. However, the energy lost by these leafs during the radiation routine reduces the temperature of the surface particles by only a small fraction, so that at the start of the next radiation routine, the surface leafs are once more at the high temperature of the surface particles, and again lead to the overall energy of the body reducing too rapidly.

Theoretically, the luminosity at the surface of a quasi-static body, and hence the effective surface temperature, is controlled by the flux of energy through lower layers, i.e. the effective surface temperature is such that the blackbody radiation emanating from the surface is equal to the flux of energy feeding it from lower layers. This reasoning can be used to reduce the luminosity, and hence the effective surface temperature of the body, to an amount that corresponds more closely to the flux through lower layers. This is achieved by proposing that surface leafs do not actually correspond to the surface of the body, but to some layer interior to the surface. This proposal is supported by the above demonstration (Figure 4.7), that as the number of particles in the simulation is reduced, the experimental radius becomes smaller than the theoretical radius.

In principle, this proposal can be implemented by surrounding the body with an additional layer of leafs, with properties equal to that of adjacent surface leafs. This virtual layer of leafs would serve to modify the energy radiated by the real surface leafs, via the attenuation process described in Section 3.2.12. The virtual layer would not radiate, and the energy absorbed in it would be discarded. In practice, a virtual layer is not required, as the energy loss associated with power packets leaving the system can simply be reduced by the attenuation factors of the surface leafs from where the power packets originate5.2. This modification reduces the collapse rate by around a factor of two for quasi-static bodies simulated with less than around $1000$ particles, although the discrepancy in collapse rate between bodies modelled with $4000$ particles and $1000$ particles is still a factor of two or so.

Another method of softening this effect, could be to reformulate the equations which govern how the particle properties are spread onto the tree structure (Section 3.2.7), so that instead of the temperature across an isolated SPH particle being constant, it falls off in a bell shaped manner similar to the density. The temperature of leafs would then automatically fall to zero at the surface of the body.

Results

In these simulations, several different mass protoplanets are considered, with masses ranging from $1 M_{J}$ to $16 M_{J}$. The minimum of one Jupiter mass is used, so that the results can be compared with those presented elsewhere in the literature. Since quasi-static collapse is being investigated, the global properties of the protoplanet evolve on a radiation timescale. Disregarding the radial motion associated with quasi-static collapse, if the SPH particles of the protoplanet are stationary, then the determined dynamic timescale controlling the time-step is of the order of the radiation timescale. However, as will be presently discussed, the protoplanets contract in dynamic quasi-static equilibrium (convective motion is present). This leads to the simulations taking a long time to complete, even when a relatively few number of particles are employed ($\approx 1000$). This is especially true of low mass protoplanets, where, in the regime under investigation in this experiment, the dynamic timescale is much shorter than the radiation timescale. In order to investigate the discrepancy in using only a small number of particles, the first experiment is a comparison of simulations employing $1067$ particles, with simulations employing $4381$ particles, for a relatively small number of time-steps. The full results of the low resolution simulations, in which many time-steps are used, are then presented.

Whilst testing the code in Chapter 4, performing comparisons between different simulations was relatively simple, as the governing equations were almost always completely scalable. However in these experiments, comparison is complicated by the inclusion of a physically realistic equation of state and material opacity, both of which depend on the temperature of the material. It was concluded that if these were kept constant during the course of a simulation, the rate of collapse of the protoplanet, once it had reached a dynamic equilibrium configuration, was approximately constant, and was proportional to the square of the mass of the protoplanet, as predicted by Equation 5.5. In deriving Equation 5.5 the assumption was that the body retained its original polytropic profile, however, this equation is equally valid for a body that retains an arbitrary profile over time. This is the case in these experiments: During the first several dynamic timescales, the protoplanet relaxes from its original polytropic profile, into a dynamic equilibrium configuration which remains approximately the same for the remainder of the simulation, as discussed below.

Comparison of Resolution

For the range of masses under consideration here, the protoplanets become opaque, and hence able to be in quasi-static equilibrium, at dimensions of around $5\,AU$. In order to make some kind of comparison, the initial average radial distance of the particles in every simulation was set to $2\,AU$. For a Jupiter mass protoplanet, this corresponds to an average temperature of particles equal to around $18\,K$, and a protoplanet that has just become opaque, $\tau\approx 10$, for the initial polytropic profile used here ($\gamma=5/3$). A protoplanet $16$ times more massive has a temperature $16$ times this, and an optical depth approximately $1000$ times as large. Figure 5.1 shows the results of this comparison of resolution experiment.

Figure 5.1: Contracting protoplanets: Comparison of resolution. Symbols are plotted every $50$ time-steps. Temperature is in units of the initial average particle temperature. Solid line, 1067 particles; Dotted line, 4381 particles; Triangles, $1 M_{J}$; Squares, $4M_{J}$; Diamonds, $16 M_{J}$.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/xmgrgraphs/... .../thesis/xmgrgraphs/compplanetsavT.eps,angle=0,width=0.9\linewidth}} \end{figure}

Since the protoplanets initially have zero velocity, they initially radiate without collapsing. This leads to the average temperature of the particles initially dropping, and subsequently to the protoplanet initially collapsing too rapidly. The result of this is that the protoplanet initially oscillates as it contracts. After several dynamic timescales, the oscillations have ceased, and the protoplanet is in a new dynamic equilibrium configuration, contracting quasi-statically at a gradually decreasing rate. This is to be expected, as the opacity increases with temperature in this regime. The graphs show that the rate of collapse in the low resolution simulations is around a factor of two larger than in the high resolution simulations. The top graph would be the mirror image of the bottom graph if the bodies were in a state of quasi-static collapse during the whole simulation, i.e. if the initial oscillations did not occur.

Dynamic Equilibrium Profiles

As previously mentioned, the protoplanets adopt dynamic equilibrium structures after several dynamic timescales. The temperature profiles of these at the end of the first experiment (average radial distance equal to $1AU$), are shown in Figure 5.2 for the high resolution ($4381$ particles) protoplanets.

Figure 5.2: Dynamic equilibrium temperature profiles for $4381$ particles at an average SPH particle radial distance of $1\,AU$. The arrows in the bottom frame indicate the velocity of the particles in the band (length of arrow) and the number of particles in the band (width of arrow).
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/xmgrgraphs/profilea2.eps,angle=0,width=0.85\linewidth}} \end{figure}

Each simulation exhibits convective motion, the more massive the protoplanet, the more violent the motion. These graphs only represent the distribution of particles at the end of the simulation. However, once dynamic quasi-static equilibrium has set in, the profiles of the protoplanets throughout the rest of the experiment, are similar to those at the end of the experiment. In the bottom frame, $M=16M_{J}$, the convective motion is so violent, that two bands of material flowing in opposite directions are clearly visible. The top band contains a large number of particles rising gradually to the surface of the protoplanet, whereas the bottom band contains a comparatively small number of particles, sinking rapidly to the centre. The length of the thin arrow corresponds to around $50$ time-steps, which, at the end of that simulation, is equal to around $4$ years. The simulations in the top two frames of Figure 5.2, also exhibit convective motion, but to a lesser degree, and the bands are more difficult to discern on a static media. By animating the profiles in time, the convective motion can easily be seen, even in the $M=1M_{J}$ simulation. For the majority of this Jupiter mass simulation, convective motion is only present in the inner regions of the protoplanet. Presumably this is due to the low temperatures in outer layers. It should be noted that for presentation on paper, these dynamic quasi-static equilibrium profiles, would probably benefit from the inclusion of small arrows on each point, corresponding to the radial velocity of the particles.

1067 Particle Experiment

The purpose of this experiment is to continue the previous experiment to much later times than possible using $4381$ particles. Figure 5.3 shows the results of these simulations.

Figure 5.3: Contracting protoplanets: 1067 particles. Symbols are plotted every $500$ time-steps. Triangles, $1 M_{J}$; Crosses, $2M_{J}$; Squares, $4M_{J}$; Circles, $8M_{J}$; Diamonds, $16 M_{J}$.
\begin{figure}\vspace{-1cm} \begin{flushright} \epsfig{figure=/home/steve/tmp/ho... ...old/thesis/xmgrgraphs/1067avT2000.eps,angle=0,width=1.0\linewidth}} \end{figure}

As well as extending the results from the previous experiment, two additional simulations, $M=2M_{J}$ and $M=8M_{J}$, are shown. The simulations were terminated when the temperature in the centre of the protoplanets reached $2000\,K$, or, as was the case for $M=1M_{J}$ and $M=2M_{J}$, the simulation was stopped as it would have taken an infeasible amount of computational run-time to complete. The main conclusion to be drawn from these graphs, is that unlike the conclusions of the investigation by Schofield (1979) and Schofield and Woolfson (1982a), bodies of less mass than around $5$ Jupiter masses, take thousands, and even tens of thousands of years to reach the point of hydrogen dissociation, rather than tens to hundreds of years.

An empirical observation, the results of which are not presented here, is that in the simulation domain of these experiments, the rate of collapse is approximately related to the opacity of the material by:

\begin{displaymath} \frac{dR}{dt}\propto -\frac{1}{\sqrt{\kappa}} \end{displaymath} (5.6)

This means that reducing the opacity of the material by an order of magnitude, increases the rate of collapse by a factor of three. This has consequences for the above conclusions if any of the assumptions concerning radiation transport are not valid5.3, or that the opacity of the material is much different to that given in Section 3.2.3. The relationship in Equation 5.6 is only valid if the body remains optically thick following a reduction in opacity.

Larger Mass Bodies

The reason for terminating the previous simulations when the temperature in the centre of the body reached $2000\,K$, is due to this first phase of quasi-static collapse coming to an end: At $2000\,K$, the released gravitational energy which usually serves to heat the gas, instead dissociates hydrogen molecules. Although the thermal energy content of two hydrogen atoms at $2000\,K$, is greater than one hydrogen molecule at $2000\,K$, the difference is small compared to the energy of dissociation, see Section 2.9. Therefore, once the temperature in the centre has reached $2000\,K$, only a small amount of further contraction is required, before the pressure in the centre of the protoplanet becomes unable to support the outer layers. At this point, a rapid collapse of the core commences, in which the density and temperature in that region increases by many orders of magnitude. Once the resulting high density core comes into thermal and gravitational equilibrium, it has dimensions two orders of magnitude smaller than the dimensions of the whole protoplanet. From the point of view of the Capture Theory, once such a body has formed it is unlikely that it will be tidally disrupted at first perihelion. Thus, from the point of view of this thesis, the details of protoplanetary evolution beyond $2000\,K$ are unimportant.

Whilst the hydrogen molecules are dissociating in the centre of the protoplanet, the opacity in that region has already dropped by three orders of magnitude5.4, due to the grains present evaporating around $1500\,K$. This means that the central regions can transfer energy efficiently to opaque outer regions. However, whilst the body is hydrodynamically supported, little energy is transfered due to the relatively small temperature gradients. As the temperature in the centre rises during core collapse, radiative energy is transported efficiently to opaque outer regions. The heating of outer layers during this core contraction phase, can lead to the surface of the protoplanet been `blown' away. The $1067$ particle $M=8M_{J}$ simulation was continued beyond $2000\,K$, and a stable core of radius $\approx 0.02\,AU$, of maximum temperature $\approx 25000\,K$ and maximum density $\approx 1\,kg\,m^{-3}$ was formed, whilst the outer layers of the protoplanet were lost.

The question of the behaviour of more massive bodies than those investigated in the previous section, is important with regards to the protostar. In early investigations of the capture mechanism, the protostar was assumed to be in a state of quasi-static collapse at the time of the interaction (Woolfson, 1964; Dormand and Woolfson, 1971), although a more freely collapsing protostar was recently envisaged (Dormand and Woolfson, 1988). The conclusion to be made from the above discussion, is that if the protostar is required to be in a state of quasi-static collapse, of the order of the dimensions originally conceived (several $AU$), its central temperature must not be greater than $2000\,K$. It is possible, however, that the formation of a high temperature and density core may support the outer layers of the protostar, and hence extend the lifetime that the protostar spends in an extended quasi-statically contracting state. Only a more physically detailed investigation could determine the structure, evolution, and lifetime of such bodies. See Bodenheimer (1993) for a review of the subject of protostar formation.

In the first Capture Theory paper (Woolfson, 1964), the protostar had a mass equal to $0.15M_{\odot}$, and a radius equal to $14.7\,AU$. If it is assumed that the distribution of material in this protostar was similar to a $\gamma=5/3$ polytrope, this results in a central temperature equal to around $1000\,K$, and a surface temperature equal to around $200\,K$. For this configuration, the combination of a relatively small optical depth across the protostar $\tau\approx 10^{4}$, and relatively high temperatures, results in the radiation timescale being comparable to the dynamic timescale. In an isolated simulation to test the evolution of this protostar, after only one time-step of $10$ years, the temperature in the centre had reduced to $850\,K$, and the temperature at the surface to $180\,K$. This rate of energy loss led to an almost free-fall collapse, with the temperature in the centre of the protostar reaching $2000\,K$ after around $50$ years, at which time the radius had reduced to $7\,AU$, and the radial velocity increased to $2\,km\,s^{-1}$. As already noted, increasing the opacity of the material increases the radiation timescale. It is not until the opacity is increased by two orders of magnitude, that a non-rotating protostar of this mass can possibly be in a state of extended quasi-static equilibrium. Rotationally supported protostars, where the primary star is more likely to capture disk material, rather than core material, have not been investigated here.

It is therefore apparent, that the original assumption of a non-rotating quasi-statically collapsing `massive' protostar, of dimensions of the order of $AU$, cannot be modelled using the current code, unless the opacity is artificially enhanced by several orders of magnitude. Interactions between rotationally supported disks and compact stars have recently been examined by Boffin et al. (1998). Their experiments showed that an interaction between a disc (with a protostar at the centre) and a star, provides a mechanism for removing energy from, or adding energy to, the orbits of the interacting stars. It was hoped that energy lost in the encounter, via the disk, would result in the two stars becoming bound. They showed that this only happened very occasionally, and thus was unlikely to be an important binary formation mechanism. However, they found that a more significant consequence of such encounters, was that they can trigger fragmentation of the disc, via tidally and compressionally induced gravitational instability, leading to the formation of additional stars in up to half of all of encounters. In their simulations, these `stars' remained bound to the protostar-disk in which they were formed.


The Capture Mechanism

The initial state of the approaching protostar is unknown. Originally it was postulated that the protostar was in a state of slow quasi-static collapse, so that even though the cross-section of the encounter may be small, the available time for a suitable interaction is large. As a result of simulating such encounters, a new type of encounter is suggested, in which the protostar is large in extent, but rapidly collapsing. It is discussed why the original premise is not conducive to planetary formation, and why the latter premise is.

Initial Orbital Components

In most encounters, the simulation commences with the protostar outside of the Roche limit of the Sun. This is so that the tidal forces exerted on the protostar are small, so that the protostar can initially be approximated as being spherical. Also under these conditions, a simple two body approximation can be used to determine the initial position and velocity of the protostar relative to the Sun. To simplify analysis of results, these simulations are performed in the reference frame of the Sun. Many procedures are possible for defining the orbit. In these simulations, it is convenient to specify the eccentricity, $e$, the perihelion distance, $q$, and the initial Sun-protostar distance, $d$. Figure 5.4 shows a typical initial configuration for $e>1$.

Figure 5.4: Orbital components of a typical encounter.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/startinga2.eps,angle=0,width=0.25\linewidth}} \end{figure}

In this diagram, the mean anomaly, $\theta$, and the semi-latus rectum, $p$, are also shown, where:
\begin{displaymath} p=q(1+e) \end{displaymath} (5.7)

The polar equation of the conic section is:
\begin{displaymath} d=\frac{p}{1+e\,cos(\theta)} \end{displaymath} (5.8)

Since the protostar approaches in an arbitrary direction, the analysis can be simplified by orientating the axes of the encounter, such that the orbital motion of the protostar lies in the $x-y$ plane, with orbital symmetry about the x-axis, as is the case in Figure 5.4. The initial coordinates of the protostar are then:
\begin{displaymath} {\setlength\arraycolsep{2pt} \begin{array}{cl} x_{0}=&d\cos(\theta)\\ y_{0}=&-d\sin(\theta) \end{array}} \end{displaymath} (5.9)

The initial velocity components of the protostar are obtained by differentiating the position components in Equation 5.9, with respect to time:
\begin{displaymath} {\setlength\arraycolsep{2pt} \begin{array}{cl} \vspace{8pt} ... ...{d^{2}\dot{\theta}}{p}\Big\{e+\cos(\theta)\Big\}} \end{array}} \end{displaymath} (5.10)

where $\dot{d}$ in the above derivation has been replaced by:
\begin{displaymath} \dot{d}=e\frac{d^{2}\dot{\theta}}{p}\sin(\theta) \end{displaymath} (5.11)

which was obtained by differentiating Equation 5.8 with respect to time. The magnitude of the specific angular momentum, $h$, of the protostar about the Sun, is given by:
\begin{displaymath} h=d^{2}\dot{\theta} \end{displaymath} (5.12)

and also by:
\begin{displaymath} h=\sqrt{\mu\,p} \end{displaymath} (5.13)

where $\mu=G(M_{\odot}+M_{P})$, due to the calculation being performed in the reference frame of the Sun, where $M_{P}$ is the mass of the protostar. Solving for $\dot{\theta}$, and substituting this into Equations 5.10, the initial velocity components of the protostar become:
\begin{displaymath} {\setlength\arraycolsep{2pt} \begin{array}{cl} \vspace{8pt} ... ...p}\Bigg)^{\frac{1}{2}}\Big\{e+\cos(\theta)\Big\}} \end{array}} \end{displaymath} (5.14)

The approach distance, $d$, is usually defined as some fraction, $f$, of the Roche limit:
\begin{displaymath} d=1.5f\Bigg(\frac{M_{\odot}}{\rho_{P}}\Bigg)^{\frac{1}{3}} \end{displaymath} (5.15)

where $\rho_{P}$ is the maximum density in the protostar. For quasi-statically collapsing protostars, see Section 5.2.2, a value of $f$ between $1$ and $2$ is usually used. For a more freely collapsing protostar, see Section 5.2.3, values of $f$ between $0.5$ and $1.5$ are commonly used.


Quasi-Static Collapsing Protostar Interactions

Encounters with quasi-statically collapsing protostars can be broadly classed into two categories: Those in which the protostar contracts slowly, and those in which the protostar contracts rapidly. For the purpose of this classification, a slowly contracting protostar, is one in which the timescale of the interaction is much shorter than the contraction timescale, e.g. such a protostar placed well outside the Roche limit of the Sun, $3$ to $5$ times say, will contract insignificantly, by less than $\approx10\%$ of its radius, before it has reached perihelion within the Roche limit. A rapidly contracting protostar, is one in which the timescale of the interaction is much longer than the contraction timescale, e.g. a protostar of this type placed on the Roche limit of the Sun, collapses at such a rate that the Roche limit recedes faster than the protostar approaches perihelion. In such encounters, the protostar is stable to disruption at perihelion. Clearly the limit of a slowly contracting protostar, is one in which the protostar does not contract at all. The limit of a rapidly contracting protostar, is one in which the protostar ceases to behave quasi-statically, and approaches a contraction rate close to free-fall collapse.

Protostars which collapse at a moderate rate, such that when placed on the Roche limit of the Sun, they are just disrupted at perihelion, tend to produce results similar to those obtained using slowly contracting protostars, rather than rapidly contracting ones. Little attention was paid to either moderately contracting protostars or rapidly contracting protostars. This was mainly because their lifetime in a quasi-statically contracting state is relatively short, such that the question of their prior evolution before the encounter becomes relevant. Also, the definition of a rapidly contracting protostar, requires that its initial approach distance be within the Roche limit of the Sun for disruption to occur. For very large extended protostars in the first stages of formation, this is a valid proposition, as there seems no reason why the formation sites of such bodies, may not reside within the Roche limit of the star they are to interact with. However, for protostars that are assumed to preexist before the encounter, which is necessarily the case for quasi-statically contracting protostars, this is not a valid proposition. This analysis led to an investigation of the evolution of protostars in the early stages of formation, in the presence of a tidal field from a more massive body, as described in Section 5.2.3.

In order to investigate the possibility of planetary formation, as a result of an interaction with a quasi-static, slowly contracting, extended protostar, at least one of the following assumptions have to be made.

  1. Radiation transfer is negligible, and can therefore be omitted.
  2. The opacity is two or three orders of magnitude larger than thought.
  3. The protostar is hundreds of times less massive than the Sun.
These assumptions are based on the conclusions of the investigation into quasi-statically collapsing bodies in Section 5.1. If an assumption is only weakly adhered to, the rate of collapse of the protostar could be such that it is stable to disruption at perihelion. Even if the protostar collapses slowly enough that it would have been disrupted at perihelion, if its central temperature reaches $2000\,K$ before tidal disruption has occurred, the ensuing collapse of the core would result in it almost certainly becoming stable to disruption at perihelion.

The simulation domain of point 1, has been thoroughly investigated during the course of this thesis, at a time when radiation transfer had yet to be included in the code, and indeed, the null result from those experiments were the primary reason for the subsequent investigation and development of a radiation transport algorithm. In the limit of a large opacity, point 2 reduces to point 1. The option of using low mass protostars, point 3, is perhaps the most physically realistic of the possible proposals, but, as shall be presently discussed, it is even less conducive to planetary formation than the first two.

In all of the above proposals, radiation transfer is largely unimportant before tidal disruption has occurred, as the radiation timescale is much longer than the dynamic timescale, and hence also the interaction timescale, for slowly contracting protostars. However, the properties of the material subsequent to tidal disruption is unknown, and therefore, whilst performing simulations under the assumptions of point 2 and point 3, radiation transport was included. It was hoped that with the inclusion of radiation transport, one of the following two possible sequences of events would occur. The first of these proposed that since the material of the protostar is in a state of collapse before disruption, a small amount drawn from it would continue to collapse in isolation. Failing this, the second sequence of events proposed that following tidal disruption, if the material in the filament expanded, it would cool hydrodynamically. This expansion, and consequent cooling, would reduce the optical depth of the material significantly, such that further cooling by radiation was possible. If either of these sequence of events occurred, dynamical fragmentation within the filament, leading to protoplanetary condensations forming, might occur.

Quantitatively, the above three assumptions all result in simulations which demonstrate the same behaviour, from which a conclusion could be made: If a filament of material is captured from a slowly contracting protostar, the filament enters orbit as a disk, or perhaps as rings of material, in which protoplanetary fragmentation does not occur. The reason for this null result was subsequently investigated. This investigation began with questioning why it seemed that neither of the proposed hypotheses had operated - it was expected that if the first hypothesis did not operate, then the second would. With regards to the first hypothesis, the assumption that material drawn from a collapsing body continues to collapse in isolation, is not valid for bodies collapsing in quasi-static equilibrium. For such bodies, Equation 5.1 applies at all times. Consequently, material drawn from the protostar will initially expand, due to its temperature and dimensions being of the order of those of the protostar, but its mass being much smaller.

After it had become apparent that the first sequence of events does not occur for quasi-statically collapsing bodies, it followed that the second sequence of events probably does. Indeed, the filament of material drawn from the protostar, does initially expand and cool, and hence lose further energy by radiation. However, the expansion is driven mainly by tidal forces from the Sun, rather than the internal hydrodynamical forces of the filament. This tidal expansion increases the Roche limit of the material in the filament, so that for this type of experiment, even though the filament may initially be receding from the Sun, it never moves outside its Roche limit, i.e. even in the limit of the material cooling to absolute zero, the filament is unable to fragment to form protoplanets.

The above conclusion applies broadly to all classes of slowly contracting quasi-static protostars. It is possible, however, that given a sufficiently fast encounter, material drawn from the protostar will recede rapidly from the Sun and out of the Roche limit, and thus be able to fragment. Although protoplanets formed from such an encounter, would not be captured by the Sun. The failure of this type of simulation to produce protoplanets, can be explained in terms of the Jeans criteria and the Roche limit. When a slowly contracting protostar in quasi-static equilibrium enters the Roche limit of the Sun, its thermal and gravitational energies are in balance, and hence the whole of the protostar only just satisfies Jeans' criteria. Material drawn from the protostar is less than Jeans critical, and therefore initially expands, resulting in the Roche limit of the material, centred on the Sun, also expanding. Whilst the captured material is moving to aphelion, it is continually expanding due to the differential tidal forces acting on it from the Sun, and as a consequence, it never moves outside of the Roche limit. For small mass protostars, this effect is even more pronounced, due to the large ratio of the initial size of the Roche limit, to that of the initial radius of the protostar.

As an example of the behaviour exhibited by this type of encounter, consider Figure 5.5.

Figure 5.5: An example of a quasi-statically collapsing protostar encounter, with parameters equal to that of Woolfson's 1964 paper: $M=0.15M_{\odot}$, $R=14.7\,AU$, $q=3R$, $e=1$. All $11191$ SPH particles, represented as points, are shown at intervals of 50 years. The circle represents the Roche limit of the protostar, calculated from the central density of the protostar at $t=0$. At $t=250$ years the Roche limit of the material in the filament is around $200\,AU$.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/results/ses/classic/196411191fig.eps,angle=0,width=1.0\linewidth}} \end{figure}

The parameters of this simulation are similar to those which gave Figure 1.9 in Chapter 1, which was taken from the first Capture Theory paper (Woolfson, 1964). Even though the existence of such a slowly collapsing protostar is questionable, this experiment serves as a comparison between methods. More physically realistic quasi-static encounter simulations, in which the mass of the protostar is reduced, or the opacity of the material is increased, give very similar results. Option $1$, omitting radiation transport, was implemented here. The $11191$ particle protostar initially has an equilibrium profile of a $\gamma=5/3$ polytope, a mass of $M=0.15M_{\odot}$, a radius $R=14.7\,AU$, a perihelion distance of $q=3R$, and an eccentricity of $e=1$.

Although these types of experiments do not lead to the formation of planets from the direct collapse of material in the filament, they do result in a considerable fraction of the mass of the protostar (up to around $20\%$) being captured by the Sun. In this particular simulation, one ring of material, of mass equal to around $15$ times that of Jupiter, with eccentricity equal to around $0.5$, and semi-latus rectum of around $20\,AU$ was formed. In simulations in which the viscosity is reduced by an order of magnitude, or the resolution increased substantially, little velocity mixing occurs, and the captured material retains the orbital parameters of the filament at the time it left the protostar. This leads to several rings of material forming about the Sun, within which more conventional theories of planetary formation could possibly operate.


Near Free-Fall Protostar Interactions: Tidal Induced Fragmentation Followed by Capture

This type of interaction was conceived as a result of the investigation in the previous section. To summerise those conclusions, the captured filament does not fragment because the material drawn from the protostar is less than Jeans critical. This leads to the material in the filament initially expanding, and as a result, it never moves outside of the Roche limit of the Sun. An additional factor to be considered when tidal forces are present, hitherto undiscussed, is that material in an elliptical orbit becomes pinched together as it approaches perihelion, and stretched apart as it approaches aphelion. Since in the previous type of encounter, the material leaves the protostar just before aphelion, this leads to an additional expansion of the material. Before the investigation into this new type of encounter commenced, the last line of investigation for the previous type of encounter, focused on examining the nature of this pinching effect, as it was hypothesised that it may result in the density of the material at perihelion, becoming sufficiently enhanced, that fragmentation may be possible. The conclusion from this small investigation, was that although the density increased significantly at perihelion, the corresponding reduction in the material-Sun separation at that distance, meant that the material was still well within the Roche limit.

This new type of encounter, in similarity to the previous type, assumes a preformed Sun, with a luminosity and mass equal to its current day values. However, the protostar is in an earlier stage of formation than previously envisaged. For model purposes, the protostar is initially spherical, with a uniform temperature and density, and is stationary in its frame of reference. Hence the mass of the protostar, along with its initial radius, temperature and orbital components, are simulation parameters. Although none of the former assumptions can strictly be true, they serve to contract the problem somewhat. Indeed, the protostar would probably have a finite rate of collapse, it would certainly not be spherical, and it would probably have a temperature gradient across it as a result of heating from local sources, e.g. from the approaching star.

Figure 5.6: A schematic representation of tidal induced fragmentation followed by capture. (a) The protostar approaches on a hyperbolic orbit, (b) As it collapses, it deforms into an egg shape, (c) The whole protostar is stretched into an arc shaped filament of material at perihelion, (d) As the filament leaves perihelion it straightens up, (e) The filament fragments to produce several protoplanetary condensations, (f) In parabolic encounters, upto half of the protoplanets are captured into high eccentricity orbits.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/picmoda1.eps,angle=0,width=0.97\linewidth}} \end{figure}

Figure 5.6 shows the typical characteristics of this type of encounter. As the protostar approaches the Sun, Figure 5.6(a), it begins to collapse under self gravity, Figure 5.6(b). Under a given range of, as yet unknown, starting conditions, the rate of collapse of the protostar is sufficiently slow, that it does not form a core stable to disruption at perihelion, and sufficiently fast, that, aided by the pinch effect, it forms an extended arc shaped filament of Jeans critical, but within the Roche limit, material at perihelion, Figure 5.6(c). As the arc of material leaves perihelion, it straightens up somewhat, Figure 5.6(d), and due to the speed of the encounter, and the fact that the material was Jeans critical at perihelion, it fragments into several protoplanets, Figure 5.6(e), some of which are captured into orbit about the Sun, Figure 5.6(f). Thus in this type of encounter, the whole protostar becomes the filament, with upto half of its mass being captured under suitable conditions. Because of the characteristics of this type of encounter, it is referred to as ``tidal induced fragmentation followed by capture''.

Results

Only a very small subset of the possible simulation domain has been investigated to date. In this section, an example encounter is presented, in which the initial conditions are suitable for this mechanism to operate. There are several factors which limit the choice of initial conditions. Currently, encounters in which material of the protostar passes through the Sun cannot be simulated. This is because the Sun interacts with SPH particles via a point mass force, due to the correct procedure for softening the force, without affecting the orbital elements of the particles, having yet to be decided. Spatial resolution forms another computational limitation: The number of SPH particles in the simulation must be such, that the resolution in the filament is fully resolved at all times. This has consequences for the mass and initial size of the protostar, as these both to some extent control the dimensions of the filament. If the filament is extremely long and thin, as is the case of filaments formed from protostars with large initial radii, then the possibility for fragmentation is severely reduced when the filament is not sufficiently resolved, due to the large smearing of the gravitational force, and less than $\approx3N_{N}$ particles contained in potential condensations5.5. Protostars with initially large radii, form thinner and longer filaments than those with small radii, as they have further to contract before the material in the filament becomes opaque, and hence able to support itself hydrodynamically. Since the mass is a contributing factor to the opacity of the material, it also influences the dimensions of the filament. A full understanding of the relationship between all parameters in the system is probably not possible analytically. Further experimentation will reveal more into the nature of this type of encounter.

An Example Encounter: Without Solar Radiation

When investigating this type of encounter, it is useful to perform each simulation twice, once with solar radiation, and once without, in order that its influence maybe assessed. In this simulation, the protostar has a mass equal to $M=M_{\odot}/20$, an initial radius equal to $R=100\,AU$, and an initial temperature equal to $T=25\,K$ such that the material is initially in virial equilibrium. The encounter is parabolic, $e=1$, with perihelion distance equal to $q=60\,AU$, and approach distance equal to $f=0.68$ of the Roche limit, which corresponds to $d=450\,AU$ and a mean anomaly of $\theta=137^{\circ}$. In an attempt to increase the resolution, so as not to restrict the fragmentation process, the number of neighbours is set to $N_{N}=25$ when simulating encounters of this type. Figure 5.7 shows the results from this simulation.

Figure 5.7: An example encounter without solar radiation. Initial simulation parameters: $M=0.05M_{\odot}$, $R=100\,AU$, $q=0.6R$, $e=1$, $d=450\,AU$, $T=25\,K$. All $11113$ SPH particles, represented as points, are shown at intervals of 400 years. The circle represents the Roche limit of the material, calculated from the maximum density in the filament at $t=800$ years.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/results/ses/modern/exp1/allon1fig.eps,angle=0,width=1.0\linewidth}} \end{figure}

After $200$ years, the protostar has cooled to around $10\,K$, and the central regions have began to collapse freely. By $1600$ years, the filament has fragmented into six condensations, of which three, and possibly four, are captured by the Sun. The mass of these range from two to ten times that of Jupiter. Figure 5.8 shows enlarged areas of interest of this simulation.

Figure 5.8: Enlarged areas of interest of the encounter in Figure 5.7. (a) A thin slice in the x-y plane, with particles hued according to density. Radius of circle associated with each particle, $0.4h$, (b) Filament at perihelion. The circle centred on the Sun represents the Roche limit of the material in the filament with greatest density, (c) Condensations forming in the filament, (d) Captured protoplanets.
\begin{figure}\begin{center} \begin{tabular}{c@{\quad}c} \epsfig{figure=/home/st... ...$5\times 10^{-8} $\ & 150 & 910 \\ \hline \end{tabular}\end{center}\end{figure}

Figure 5.9: Fragmentation of the filament into protoplanets. Shown approximately every $50$ years.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/results/ses/modern/exp1/filamentfig.eps,angle=0,width=1.0\linewidth}} \end{figure}

In the table in Figure 5.8, the third column, $\rho_{max}$, denotes the maximum density of the material present, the fourth column, $\bar{T}_{max}$, denotes the mean maximum temperature. This is calculated by discarding particles with uncharacteristically high temperatures, and taking an average of the remaining particles in the high temperature region. This procedure is adopted, due to particles occasionally developing temperatures around a factor of two different from neighbouring particles, when radiation transport is included. The fifth column, $N_{S}$, denotes the number of time-steps since the beginning of the simulation.

Figure 5.8(a) shows the configuration of the protostar just before perihelion. Its collapsing motion, together with tidal pinching, is responsible for this shape. For future reference, when solar radiation is not included, the temperature distribution of the protostar is similar to the density distribution at this time. The fragmentation process of Figure 5.8(c) is of particular interest, and this is demonstrated in greater detail in Figure 5.9. The timescale for fragmentation is seen to be of the order of $200$ years for this encounter. The mass of the captured protoplanets range from three to nine times that of Jupiter, the subsequent evolution of which are considered in Section 5.2.7.

An Example Encounter: With Solar Radiation

The parameters of this simulation are exactly the same as those in the previous simulation, the only difference being that solar radiation is included. Figure 5.10 shows the results from this simulation.

Figure 5.10: An example encounter with solar radiation. Initial simulation parameters: $M=0.05M_{\odot}$, $R=100\,AU$, $q=0.6R$, $e=1$, $d=450\,AU$, $T=25\,K$. All $11113$ SPH particles, represented as points, are shown at intervals of 400 years. The circle represents the Roche limit of the material, calculated from the maximum density in the filament at $t=800$ years.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/results/ses/modern/exp1withsolar/allon1fig.eps,angle=0,width=1.0\linewidth}} \end{figure}

Only two condensations are formed, one at each end of the filament, and thus it is clear that the inclusion of solar radiation, in this encounter at least, is an important factor. The complete mechanism leading to this behaviour, is not yet understood, but it seems that additional heating from the Sun, leads to additional hydrodynamic support in the filament, and consequently to a decrease in the ratio of the length of the filament to its width. The filament thus fragments into fewer condensations, as predicted by Jeans, and as discussed in Section 1.3.1. Figure 5.11 shows enlarged areas of interest of this simulation.

Figure 5.11: Enlarged areas of interest of the encounter in Figure 5.10. Particles in the thin slice plots (a&b) are hued according to temperature, with the radius of the circle associated with each particle equal to $0.4h$. (a) Protostar becoming elongated as it approaches perihelion. (b) Filament at perihelion, (c) Condensations forming in the filament, (d) Captured protoplanet.
\begin{figure}\begin{center} \begin{tabular}{c@{\quad}c} \epsfig{figure=/home/st... ...$4\times 10^{-8} $\ & 170 & 970 \\ \hline \end{tabular}\end{center}\end{figure}

Figure 5.11(a) shows the configuration of the protostar just before perihelion. It is clear that the Sun is responsible for a temperature gradient developing across the protostar. At perihelion, Figure 5.11(b), the temperature gradient is still present, but to a lesser extent. The solar heating is responsible for the temperature and density in the filament at perihelion, being around one and a half times as hot and diffuse as that at this point in the previous simulation. Figure 5.11(c) shows the filament beginning to condense at both ends, to form one bound protoplanet, and one unbound protoplanet. The mass of the unbound protoplanet is around $5$ times that of Jupiter. The mass of the captured protoplanet is around $8$ times that of Jupiter, the subsequent evolution of which is considered in the next section.


Evolution of Protoplanets

In these simulations, the scale of the interaction is such that the captured protoplanets have initial orbital periods of around $3000$ years, and mass between $3$ and $8$ times that of Jupiter. As concluded in Section 5.1, the more massive of these should reach the point of hydrogen dissociation in around $500$ years, with the least massive taking around $10$ times longer. Although the protostar is initially stationary in its frame of reference, the angular momentum imparted to it during the course of the interaction, results in the protoplanets containing a finite angular momentum in their spin. A procedure for assessing the magnitude of this angular momentum, has yet to be included into the visualisation and analysis programs, but should be trivial to implement.

Quite often in this type of encounter, two protoplanets are formed bound to each other, such that they orbit about their centre of mass, and the centre of mass orbits about the Sun. The protoplanets tend not to remain bound however. In simulations where the bound protoplanets are of low mass, or the radiation is turned off so that their orbits can be followed for longer, the protoplanets become unbound subsequent to their first perihelion passage with the Sun. In one such simulation, one protoplanet entered a considerably larger orbit, at the expense of the size of the other.

Initially the protoplanets are roughly spherical in shape, due to their revolution rate being small at these dimensions. However, due to conservation of angular momentum, as the protoplanets contract in size, their revolution rate increases. Eventually the protoplanets are spinning at such a rate, that material starts to be left behind in a disk, and the near spherical core reduces in mass. Although this seems a logical explanation why experimentally the core is seen to lose mass to a disk as it contracts, this mechanism has not been investigated here, and could possibly be due to numerical artifacts, perhaps related to artificial viscosity5.6.

The subsequent evolution of the protoplanets captured in these simulations, are shown in Figure 5.12.

Figure 5.12: Protoplanetary evolution. The scales are spatially aligned. All $\approx4500$ SPH particles in the top frame (without solar radiation), and $\approx2000$ SPH particles in the bottom frame (with solar radiation), are shown.
\begin{figure}\centering {\framebox{\epsfig{figure=/home/steve/tmp/hold/thesis/r... ...4.7\times10^{-5}$&0.5\\ \hline \end{tabular}\end{center}\normalsize\end{figure}

The top frame of this figure is from the simulation without solar radiation, and the bottom frame is from the simulation with solar radiation. The scales are the same size in each frame, and are located at the same position in space relative to the Sun. In the table in this figure, $M$ is the approximate mass of the near spherical core in units of a Jupiter mass, $T$ is the maximum temperature, $\rho$ is the maximum density, and $R$ is the approximate radius of the core. All SPH particles are plotted every $150$ years. By $3050$ years without solar radiation, and by $4050$ years with solar radiation, the most massive protoplanet has reached the point of hydrogen dissociation, and the simulation cannot be continued further.

During these simulations, the more massive of the protoplanetary cores reduce in size by an order of magnitude, and consequently, their density increases by approximately three orders of magnitude. At the end of the simulation, each of the cores contain around two thirds of their original mass. The remainder of the mass forms a disk of material in orbit about the core. Each core should remain in thermal and gravitational equilibrium during its contraction. Analysis of the table in Figure 5.12 reveals that to within the accuracy of the measured quantities, this is approximately so:

\begin{displaymath} \frac{TR}{M}=100\pm20\,K\,AU\,M^{-1}_{J} \end{displaymath} (5.16)

where the above quantity was discussed in Section 1.3.1 and Section 5.1.1, and is true for all cores at all times. It is possible that the least massive of protoplanets, Figure 5.12(a), will not be sufficiently compact at perihelion to withstand tidal disruption from the Sun. An estimate of its Roche limit can be calculated from its density at the end of the simulation:
\begin{displaymath} R_{L}\approx 1.5\Bigg(\frac{2\times10^{30}}{7.7\times10^{-6}}\Bigg)^{\frac{1}{3}}\approx4\,AU \end{displaymath} (5.17)

The perihelion distance of the orbit of this protoplanet is unknown. However, since the original perihelion distance of the protostar was $60\,AU$, it is unlikely that the perihelion distance of the protoplanet will be an order of magnitude less than this, and hence its core will more than likely be stable to disruption.

Front Page
Stephen Oxley 2002-01-19