Methods for computing giant planet formation and evolution

We present a numerical code for computing all stages of the formation and evolution of giant planets in the framework of the core instability mechanism. This code is a non-trivial adaption of the stellar binary evolution code and is based on a standard Henyey technique. To investigate the performance of this code we applied it to the computation of the formation and evolution of a Jupiter mass object from a half Earth core mass to ages in excess of the age of the Universe. We also present a new smoothed linear interpolation algorithm devised especially for the purpose of circumventing some problems found when some physical data (e.g. opacities, equation of state, etc.) are introduced into an implicit algorithm like the one employed in this work.


INTRODUCTION
For a long time, the only gas giant planets known were the four ones in our own solar system. They provided the unique observational evidence to constraint theories of giant planet formation and evolution. However, this situation changed drastically during the last decade, because of the discovery of more than 100 planets orbiting around main sequence stars (see for example the web site http://www.obspm.fr/planets), revealing that gas giant planets are very common, at least in the solar neighborhood.
It is commonly assumed that gaseous giant planets have been formed by the core instability mechanism (Mizuno 1980), which states that a solid core is formed from the accretion of planetesimals in the protoplanetary disk, followed by the capture of a massive envelope from the gaseous component of the protoplanetary nebula. The whole process may be divided in some characteristic stages (i) The accretion of solid planetesimals results in the growth of a solid core with several Earth masses. This solid core is surrounded by a tenuous gaseous envelope of very low mass.
(ii) As the solid core grows, gas is accreted at an increasing rate. At some time, gas accretion rate supersedes the accretion rate of solids.
(iii) When the mass of the envelope is comparable to the mass of the solid core, a runaway gas accretion starts. During this stage, due to the exhaustion of solids in the feeding zone of the planet, very little accretion of planetesimals occurs.
(iv) The accretion of gas is terminated, because of dissipation or tidal truncation of the nebula.
(v) The planet cools and contracts at constant mass, to its present state.
The characteristic time-scales involved in the various stages described above are very different. Steps (i) and (ii) can last for some My, but during the runaway gas accretion, a gas giant can accrete one Jupiter mass of gas during few 10 3 y. During this stage, the luminosity can rise to 10 −3 − 10 −4 L (Bodenheimer and Pollack 1986), in part due to the release of gravitational energy forced by the violent initial contraction of the envelope, during the first moments after the termination of the accretion process.
Although the core instability model is conceptually very simple, modeling it accurately is, on the contrary, diffi-cult, and usually, some particular simplifying assumptions, for each one of the different stages, have to be adopted (Bodenheimer, Hubickyj & Lissauer 2000).
At present, very few codes capable to model the formation of giant planets, incorporating the relevant physics involved in the core instability model, exist. We can mention here the ones developed by Bodenheimer and Pollack (1986) (with many further improvements), and by Wuchterl (1990).
It is the aim of the present work to describe a numerical code devised for computing the formation and evolution of giant planets in the frame of the core instability model. This code is capable to compute all the stages of the formation; subsequent detachment form the protoplanetary nebula and its final contraction and cooling. This code is a non -trivial adaption of a Henyey code tailored for computing stellar evolution in close binary systems with mass transfer (Benvenuto & De Vito 2003).
This paper is organized as follows: In Section 2 we describe the main features of our code. In Section 3 we briely describe the physical ingredients we have incorporated to our code. In Section 4 we show the performance of the code applying it to the formation and evolution of a Jupiter mass object from a half earth core mass to ages in excess of the age of the Universe. Finally, Section 5 is devoted to give our conclusions. We also present in the Appendix a new smoothed linear interpolation algorithm devised especially with the purpose of circumventing some problems found when some tabulated data (e.g. opacities, equation of state, etc.) are introduced in an implicit algorithm like the one employed in our code.

Equations of Structure and Evolution
Here, we shall briefly summarize the equations of giant planet formation and evolution to be solved by our code. As usual, we consider spherically symmetric objects, neglecting rotation and magnetic fields. In spite that it is currently accepted that giant planet formation and evolution is a phenomenon that occurs in conditions very near to hydrostatic equilibrium, we prefer to develop a full hydro code. Obviously, this is more general than an hydrostatic code. Moreover, it may happen that some planets are formed in conditions very different to those currently accepted for which hydrodynamic phenomena may be relevant.
In the conditions we are interested in, the equations of giant planet formation and evolution are (for derivation of these equations see, e.g., Clayton 1968, Kippenhahn & Weigert 1990 ii) the definition of velocity iii) the equation of mass conservation iv) the equation of energy balance v) the equation of energy transport for the radiative case and vi) the equation of energy transport for the convective case where ∇conv is the convective temperature gradient, which may be computed employing the standard mixing lenght theory (see, e.g., Kippenhahn & Weigert 1990). We have considered the differential of entropy in the form We employ the Schwarzschild criterion for the onset of convection. The total gravitational energy release due to the accretion of planetesimals L pl is given by We should remark that the planetesimal core accretion rateṀcore is not specified by these equations and remains as an input of the model, not only in its initial value but also regarding its temporal variation. In order to incorporate this energy release we introduce the rate of planetesimal energy release ε pl imposing the condition that We have found it convenient to adopt an expression for ε pl in a way that the energy release is produced near the bottom of the gaseous envelope as ε pl vanishes at menv = αMcore (where menv = mr − Mcore) and is set to zero outward. A is a constant determined by the condition given by Eq. (9): While the results are fairly insensible to the precise value of the parameter α, we shall set α = 2. δ is given by and the rest of the symbols have their standard meaning.
In dealing with the specific problem of giant planets we have to make some supplementary assumptions apart from those quoted previously. We shall handle the members of our system as spherical objects, neglecting the departure from spherical symmetry of the equipotentials (e.g. the pear -like shape of the Roche Lobe) and its evolutionary consequences. Moreover, we shall consider that the giant planet moves along a circular orbit.
As usual we shall describe the problem in Lagrangian coordinates. We shall consider the independent variable ξ, defined as In the calculations to be presented below (see Section 4), the mass of the core of the starting model is of the order of Mcore = 0.6M⊕ while the amount of gas gravitationally bound is Menv ∼ 10 −5 M⊕. Meanwhile, models of evolved planets have Mcore ∼ 20 M⊕ while Menv may be as large as Menv ∼ 3 × 10 3 M⊕. Thus, if we want to consider an interval in ξ on which we can accommodate all this evolution we need to take, say, −16 ξ 10 (see below).
As stated above, we are interested in computing the evolution of the gaseous part of the planet. In handling the core we shall simply assume it to have a constant density and neglect processes that may release energy (e.g. radioactivity). However, this may be included with minor modifications of the strategy described here. Also, the rate of core growth is an input of the model, at least in the present work (see Pollack et al. 1996 for a detailed treatment of the core growth rate).
We found it very convenient to handle radii, pressure, and temperature by means of logarithmic transformations p = ln P , θ = ln T , x = ln r, whereas lr, v are considered linearly.
For simplicity, we have written the difference equations in a centered fashion. It means that we have chosen to write a generic differential equation as a difference equation where η j+1/2 = (ηj+1 + ηj )/2, being η any quantity. The second subindex j indicates the shell of the star for which the difference equation is written. Temporal derivatives have been written in the standard backward differenced form. Employing this recipe, for example, Eq. (3) becomes Because of its definition, ξ is not a time independent coordinate. Thus, in order to write the equations taking this fact into account it is very useful to rewrite the derivative operator as Then, straightforwardly, In order to close the problem we need to define adequate boundary conditions.

Inner Boundary Conditions
We shall assume that the high density core has a constant density ρcore which we have fixed at a value of ρcore = 3 g cm −3 . Furthermore, we have neglected any energy release coming from the core apart from the one due to its growth by accretion of planetesimals. Thus, in this case, the inner boundary conditions to be applied at the bottom of the envelope (where Mr = Mcore) are given by and Finally, the inner boundary condition for velocity is given by core growth given by

The Formation Stage
In the formation stage, as usual, we shall consider the external radius of the planet R as the minimum of the accretion radius Racc and the Hill radius R Hill defined as (where c is the local velocity of sound) and R Hill = a M 3Mstar Then, the planetary radius is respectively. Physically, the accretion radius is the place at which the molecular velocity equals the escape one, while the Hill radius corresponds to the equivalent radius of a sphere with a volume equal to that of the Roche lobe of the planet. Usually, it has been considered that at R the temperature and density of the planet correspond to those of the protoplanetary nebula denoted as T neb and ρ neb respectively. From a numerical point of view we have found it very convenient to impose the boundary conditions in a different way. Let us consider, as discussed above, the possibility of extending the grid far beyond the planetary edge. If we introduce of some kind of softening of the gravitational potential due to the presence of a limiting physical agent, this softening makes the gradient of the gravitational potential to drop to zero near the planetary radius and thus, we shall have, apart from hydrodynamical effects, a constant pressure. But, as the gradient of temperature is proportional to the pressure gradient, we shall also have a constant temperature region. Consequently, the density in such a region will have also a flat profile. In these conditions we can impose the physical boundary conditions corresponding to the planetary surface at the outermost point of the grid chosen to be located far outside the planet. In other words, handling the outer boundary conditions in the way described above we set T = T neb and P = P neb for a value of the independent coordinate ξ = ξ edge far larger than that corresponding to the planet. For an adequate softening of the gravitational potential we shall have T = T neb and P = P neb for values of ξ from ξ edge up to the corresponding to the actual planetary surface ξ surf , where ξ surf = ln (M/Mcore − 1). When the planet begins to undergo a noticeable growth, R will increase. In accounting this effect we have only to change the radius in the function adopted for the quoted softening, but the planetary boundary conditions will be automatically fulfilled.
Specifically, we have introduced a restricted -threebody spherically -averaged gravitational potential which is the standard one, multiplied by the factor We apply this factor up to a predetermined fraction of the planetary radius r ζR. For larger values of r we have found it convenient to use a Fermi -like function asking for continuity of the function and its first derivative at the point r = ζR. This gives A = 2 1 − ζ 3 and β = 1 − ζ 3 / 6ζ 2 . We have tested values of ζ in the range ζ = 0.90 − 0.99 finding that the global properties of the evolution the planet are fairly insensitive to the value of ζ. Also, from the numerical point of view it produces very good convergence if we take care in defining a grid dense enough near r = R.
If we consider strictly flat profiles for T and P at the interval ξ surf ∼ < ξ ξ edge , this treatment will allow the planetary envelope to grow at a rates that can be arbitrarily high (of course, the precise rate of growth of the planet will be solution of the equations). However, if the protoplanetary nebula has a low density gap near planetary surface, this will naturally impose a upper limit for the gas accretion rate of the envelope. Preliminary calculations, beyond the scope of the present paper, indicate that the scheme we present here is fairly adequate for considering migration, the occurrence of a gap in the gas distribution in the protoplanetary disk and also non constant planetesimal core accretion rates.

The Evolution Stage
In handling the outer boundary conditions at the evolutionary stage we have to take into account the irradiation from the central star. This irradiation has a non negligible effect on the evolution of the giant planet and should be especially important for the case of small planets at advances evolutionary stages (like Saturn) or in the case of extrasolar giant planets orbiting very close to the central star (see Section 1).
Notice that irradiation has an obvious non radial nature. However, in order to work in spherical symmetry we shall consider that this energy is distributed uniformly in the whole surface of the planet. Otherwise we would be forced to change our whole treatment from the very beginning. Let us remark that in the case of the giant planets of our Solar System, rotation is very fast making irradiation of the planetary surface to be approximately uniform. In more general conditions we should expect some kind of fluid circulation from the irradiated hemisphere to the opposite one driven by the presence of a temperature gradient. While this effect should be of fundamental relevance in computing the spectra of irradiated planets, it seems that assuming an uniform irradiation is acceptable for our purposes.
In the frame of these approximations, we define Lirr as the fraction of the energy irradiated by the parent star absorbed by the illuminated hemisphere of the planet (see, e.g., Guillot 2001) Here L * is the luminosity of the parent star and A is the Bond albedo of the planet. Regarding the strategy for handling outer layers integrations, we shall introduce a method that represents a generalization from the one presented by Kippenhahn, Weigert & Hofmeister (1967) devised for the case of stellar evolution. They proposed to divide the log L − log T ef f (where L is the luminosity and T ef f is the effective temperature) plane in rectangle triangles, employing outer layers integrations at the vertexes of a given triangle in which falls the point corresponding to the actual log L − log T ef f values of a given model. They compute outer layer integrations for the conditions at a each vertex and get the outer boundary conditions to be applied to the model by means of a two dimensional linear interpolation.
In stellar evolution, luminosity, effective temperature and radius are related by the well-known relation in which all these quantities are positive. However, in the case of an irradiated planet we have (see, e.g., Guillot 2001) where Lirr is given by Eq. (27). If we assume a constant Lirr, the planet would cool down asymptotically to an effective temperature T ef f 0 given by However, for most of the recently discovered extrasolar giant planets, the parent star is undergoing core hydrogen burning (the so called main sequence stage; see, e.g., Kippenhahn & Weigert 1990). At main sequence, stars suffer from a increase in its luminosity 1 . Consequently, even without migration, radiation increases as a function of time. Thus, it is possible to have a situation in which the planet reaches an equilibrium situation at a finite age. From that moment on, as irradiation increases monotonically, the planet will begin to absorb energy from its parent star, i.e., L becomes negative. This situation does not happen in stellar evolution, and force us to modify Kippenhahn et al. (1967) strategy.
In order to maintain this kind of strategy we are forced to divide a plane defined by logarithmic axes (this is very convenient because the enormous variations in the outer conditions of the giant planets during their whole evolution). Thus, obviously, these axes must correspond to positively defined quantities. In the case of the evolution of giant planets including irradiation from its parent star we shall divide the plane log R − log T ef f in rectangle triangles. Assuming a value for the radius we immediately compute a value of Lirr (see Eq. (27)), and then, by means of Eq. (29) we get the corresponding value of L. Then, we are in conditions to perform a standard outer layers integration which can be accomplished by standard methods (e.g., Runge -Kutta).
In computing the structure of the outermost layers of the planet we have employed Eqs. (1), (3), (5), and (6) neglecting terms containing temporal derivatives (which means that we have neglected the inertia and heat content of these layers). Consequently, these layers have a constant luminosity value. As in deeper layers, we have considered the fully non-ideal EOS (see Section 3), non-adiabatic convection and imposed the Schwarszchild criterium for the onset of convection.
As it is usual in stellar evolution calculations, these layers are integrated taking the total pressure as the independent variable. Regarding the fraction of mass we include in these outer layers integrations, it is advisable to keep the amount of mass in the outer layers integration MOL as small as possible in order to avoid neglecting a significant gravitational energy release due to contraction. In the case of our code we have taken MOL/M planet ∼ < 10 −3

The Stage in Between
A non trivial problem is the way in which the planet detaches from the protoplanetary nebula. In our code we have made the simplest hypothesis: when the object reaches a prefixed amount of gravitationally bounded matter we abruptly change the boundary conditions from those of the protoplanetary nebula to those described in the previous subsubsection corresponding to the evolutionary stage. Needless to say, this represents a gross oversimplification to the actual physical situation. This should be regarded as completely unsatisfactory in the case that we are particularly worried about the behavior of the planet at moments very close to detachment. However, for the present work our main concern is to reach an accurate description of the whole behavior of the system. In this sense we consider that out treatment suffices for our purposes in an initial approach to the problem.
As it will be described below, the planet reaches a very short lived high luminosity stage which we interpret to be essentially unobservable. Also, the maximum in the effective temperature of the planet is reached very soon after detachment from the protoplanetary nebula. A very important point is to investigate the dependence of the shape of the evolutionary track of the planet upon the way it de-taches form the protoplanetary nebula employing more realistic models than the one assumed in this paper.

the overall problem
Taking into account all the details described above, we handle the equations in a fully implicit way by means of the Henyey technique as presented in Kippenhahn et al. (1967). We found very fast convergence of the code in most of the considered conditions with the remarkable exception of the moment of detachment from the protoplanetary nebula.

THE PHYSICAL INGREDIENTS OF THE CODE
In the present version of the code we have incorporated the equation of state (EOS) presented by Saumon, Chabrier & van Horn (1995). This EOS has been especially devised for computing low mass objects like brown dwarfs and giant planets. It represents a very detailed treatment for hydrogen plasma an another less detailed for helium. While the treatment performed in computing the Saumon et al. (1995) EOS is detailed enough for the purposes of this work, we found some serious numerical problems in handling it as a part of an iterative scheme like the code presented in this paper. As it can be expected, the density of the plasma is a rather smooth function throughout the whole interval of temperatures and pressures, however, quite contrarily quantities that represent second derivatives of the free energies (e.g. CP , ∇ ad and δ) are functions with steep variations For formation stages and low temperatures we have considered grain opacities given by Semenov et al.(2003) for the case of low densities while for higher densities have been taken from Pollack, McKay, & Christofferson (1985). For evolutionary stages we considered the opacity data given by Guillot (1999). For temperatures above 10 3 K we considered Alexander & Ferguson (1994) molecular opacities, which are available up to T 10 4 K, for higher temperatures we considered opacities given by Rogers & Iglesias (1992). Notice, that at such conditions we expect the interior of giant planets to be in convective equilibrium. Because of this reason, radiative opacities are of little importance in such a regime.
The above described difficulties found in handling CP , ∇ ad and δ are even more serious in the case of opacities. This is so especially at very low temperatures at which dust gives large contributions.
While the EOS is reasonably well established, it seems this not to be the case for opacities. In a recent work, Podolak (2003) reexamined the opacity due to grains finding it to be significantly lower than earlier estimates. Remarkably, the values of opacities are among the most important ingredients in determining the critical mass of the core for reaching the core instability. Lower opacities than the ones considered here would produce lower critical core mass values that, in turn, would be formed faster alleviating the well-known problem of the timescales of formation of giant planets in the frame of the core instability mechanism (see also Section 1).
In the near future we plan to incorporate the few very low energy nuclear reactions (particularly deuterium burning) suffered by very massive giant planets (with masses M ∼ > 14 MJup; see, e.g., Burrows, et al. 1995).

THE CASE OF A JUPITER -MASS OBJECT
As a first application of the code presented above, we shall present the results we have found employing the above -described code corresponding to the formation of a Jupiter mass object of solar composition at a fixed orbit of 5.2 AU. We should remind the reader that it is not our aim here to present state -of -the -art model of the formation of Jupiter but to show the results the numerical scheme presented above is capable to produce. At this point we should remark that many of the characteristics of the solution we shall describe below are also present in previous calculations (see, e.g. Bodenheimer, et al. 2000).
With the aim of simplicity, we shall assume a constant rate of accretion of planetesimals 10 −6 M⊕ y −1 , a constant density for the core ρcore = 3 g cm −3 and start with a model with 0.6 M⊕ which has attached a gravitationally bound gaseous envelope of ≈ 10 −6 M⊕. Regarding the physical conditions at the protoplanetary nebula, we shall assume that they remain constant with a density ρ neb = 10 −10 g cm −3 and temperature T neb = 100 K. Also, we shall employ the adiabatic temperature gradient in convective zones.
The computation comprised about ten thousand models resolved in about two thousand mesh points. The main characteristics of the numerical solution are presented in Figs. 1 -8 and in Table 1. We shall include the solar irradiation assuming a constant luminosity of 1L .
In Fig. 1 we show the mass of the core and the total planetary mass as a unction of time. Notice that the final growth of the gaseous envelope occurs in a fairly short timescale as expected for the core instability mechanism.
In Fig. 2 we show the luminosity of the planet as a function of time. Notice that there exists a very sharp peak in the luminosity due to gravitational energy release driven by the violent contraction of the outermost layers of the gaseous envelope. Previous to that moment, luminosity was rather proportional to t 2/3 because we assumed a constant core growth rate. After the peak, luminosity decays nearly exponentially up to the moment in which solar irradiation becomes important. At these advanced stages cooling noticeably slows down.
In Fig. 3 we depict the effective temperature of the planet as a function of time. In this plot we have only included the evolutionary stages in which a physically plausible effective temperature can be defined (see Eq. (29)). Remarkably, the maximum in the effective temperature is reached almost immediately after detachment from the nebula.
In Fig. 4 we show the evolutionary track of the model. After detachment the object begins to evolve to higher effective temperatures and lower luminosities. After maximum effective temperature it begins to cool down approximately on a constant radius track. Qualitatively, this part of the track resembles the ones corresponding to very low mass white dwarf stars. This is as it should be expected for a an object with a semi-degenerate interior.  Since then on we considered constant mass evolution. We should remark that the mass of the core as a function of time is an input whereas the mass of the gaseous envelope is solution of the equations of structure and evolution (see Section 2). The starting model has a total mass of 0.6 M ⊕ and a very tiny amount of bounded gas. As core grows, about 20 My later the mass of the planet has undergone an appreciable growth and now half the total mass is in the core and the other in the gaseous envelope. From this moment on there occurs the runaway instability and very soon the planet reaches its assumed final mass value.
gaseous envelope of the planet. During fromation temperature and density increase monotonically up to the moment of the beginning of the runaway growth of the envelope. At these conditions, the bottom of the envelope reaches a maximum density and the decompress approximately at constant temperature. As consequence of the end of the formation stage, the bottom of the envelope begin to compress, reaching a maximum temperature of 6.8 × 10 3 K at an age of 32.4 My. Notice that the maximum temperature reached by the planet is far lower than the necessary for the occurrence of deuterium burning. In Figs. 6-7 we show the profiles of density and temperature for the formation and evolution stages of the planet. Notice that for most of the stages included in these plots, the profiles are very similar in each figure. This is because at formation stages the densities attained at the planetary interior are so low that non ideal effects are of minor relevance. So, the EOS is well accounted for by a simple ideal, non- Here, formation (dotted lines) and evolution stages (solid line) are clearly differentiated (see Fig. 1). Formation corresponds to times before the runaway gas accretion instability, which corresponds to moments just before the flash-like luminosity peak. Detachment from the protoplanetary nebula is assumed to occur when the object has a total gravitationally bound mass of 1 M J up . The evolutionary stages at constant mass correspond to moments after the luminosity flash. Filled square indicates the position of Jupiter, while in the inset we compare the computed luminosity of the Jupiter model at an age equal to that of the Solar System with the observed one with its corresponding error bar (notice that in the inset the horizontal scale is linear).  Here we have only considered evolutionary stages. Notice that the maximum of the effective temperature occurs just after detachment from the nebula (see Table 1). Filled square indicates the position of Jupiter, while in the inset we compare the computed effective temperature of the Jupiter model at an age equal to that of the Solar System with the observed one with its corresponding error bar (notice that in the inset the horizontal scale is linear).    . Temperature and density increase monotonically up to the moment of the beginning of the runaway growth of the envelope. As consequence of the end of the formation stage, the bottom of the envelope begins to compress, reaching its maximum temperature (see Table 1). The here computed conditions for the present Jupiter are represented with a solid square.
degenerate EOS. Notice, however, an important differece: while the density profiles evolves to higher values in most of the evolutionary stages, temperature profiles do not. This is so at the final evolutionary stages of the object at which it undergoes a global cooling. We also show, in Fig. 8 the evolution of the gaseous envelope in the thermodynamic plane. Finally, for completeness, we show in Fig. 9 the evolution of the radii of spheres containing a fixed amount of mass We should remark that, in spite that in this computaton we did not set any upper limit to the growth of the envelope massṀenv (see § 2.3.1 for discussion), we found it has beeṅ Menv 10 −3 M⊕/y even during the runaway growth.
In order to compare with observational data corresponding to the present characteristics of Jupiter, we should quote that at the present age of the Solar System of 4.55 × 10 9 y, the luminosity of Jupiter is log LJup/L = −8.660 ± 0.004 while its effective temperature log (T ef f )Jup = 2.0948 ± 0.001 (Pearl & Conrath 1991). In these data the refledted light is not included. In Figs. 2 -4 we have included insets in which we included the corresponding computed evolutionary curve together with the observations together with its error bars. We notice that at the age of the Solar System the Jupiter model is slightly smaller than the real planet. While the effective temperature shows a good agreement with observatinos, the model reveals itself as somewhat underluminous. We consider these results as satisfactory taking into account the simplifying assumptions we made in this computation.

CONCLUSIONS
The goal of this paper has been to present a numerical code intended for computing all stages of formation and evolution of giant planets in the frame of the core instability mechanism. This code is based on the standard Henyey technique usually employed in stellar evolution calculations.
Perhaps the key point of this method is the kind of change of variables given in Eq. (13), particularly considering the total mass in units of the core mass. If we naively choose a more straightforward independent variable like, e.g., total mass we would face serious troubles at computing the temporal derivative of the entropy during formation stages. With this kind of natural coordinate we would have an outwards migration of the envelope. Thus, in order to compute the difference of entropy between consecutive models at the same envelope mass element we would be forced to introduce interpolations. But any interpolation introduce numerical noise. This noise may be irrelevant at some stages but it is catastrophic during the runaway growth, preventing its calculation.
In testing the code we have computed the formation and evolution of a Jupiter mass object from a 0.6 M⊕ core mass to ages in excess of the age of the Universe. While this should not be considered a state-of-the-art calculation (because of the many simplifying assumptions, see Section 4), it shows that the general structure of the code works fairly well in simulating the whole process of formation and evolution of giant planets.
We also present a new smoothed linear interpolation algorithm devised especially with the purpose of circumventing some problems we found when some physical data is introduced in an implicit algorithm like the one employed in this work. This has been very important in allowing us to incorporate the detailed physics of the problem. These ingredients have some dramatic discontinuities due to very important physical reasons (in opacities: the grain evaporation; in the equation of state: the plasma phase transition). While the motivation in developing the method has been very specific, in our opinion it may be a valuable method of interpolation in more general conditions. We want to thank the anonymous referee for his (her) constructive criticism that enabled us to largely improve the original version of the present work. OGB was partially supported by project FONDAP 15010003. AB acknowledges the financial support of AGNPCyT, PICT 03-11044. Figure 6. The evolution of the logarithm of the density as a function of the transformed mass scale given by Eq. (13). We represent the profiles of formation (evolution) models in dashed (solid) lines. Here, the included models correspond to those whose main characteristics are given in Table 1. Initially, formation models have a very tiny amount of matter attached to the solid core and, consequently, the fall of density from the conditions at the bottom of the envelope to the ones at the protoplanetary nebula occur in a very narrow mass interval. As core grows, the density profile becomes less steep. This is especially so at runaway conditions. When the object detaches from the protoplanetary nebula it undergoes a global monotonic contraction up to the final computed model. In the profiles corresponding to evolution stage, we have not included the part correponding to the outermost layers integrations (  We have included the same models as in Fig. 6 and lines have the same meaning as there. Notice the strong resemblance of the temperature profiles an the density ones (see Fig. 6). This is due to the fact that during the whole formation stage, the equation of state of the gas is essentially that corresponding to a ideal, non degenerate gas. An important difference occurs at the last computed stages. Due to the global cooling down of the planet, the final plotted profile is cooler to the previous one. Here we also excluded outer layers integrations from the figure while the heavy line corresponds to the computed conditions for the present Jupiter. formation and evolution. In spite of the very specific motivation in constructing this interpolation method, this is very general and may be interesting in other applications. Because of this reason, hereon we shall describe it in detail.
In many numerical simulations of astrophysical interest, it is very common the neccesity of interpolating tabulated data, representing some physical quantity. Usually, these quantities are computed by means of complex numerical procedures. Even if the codes that produce the tabulated data were available, using them as subroutines would result impractical and wasteful.
Good examples of such a case are found in Section 3 where physical quantities such as the equation of state of the plasma, radiative opacities, etc., in planetary interiors must be obtained as a function of pressure (or density) and temperature from a given table. Although these tables use to be dense enough so as to provide good interpolated values, the numerical derivatives of many interpolation algorithms exhibit discontinuities which often make it difficult the numerical solution of the problem (i.e. in the case of implicit schemes like the one presented in this paper, the convergence of iterative algorithms usually turns out to be cumbersome). Even if convergence is attained, notice that smooth derivatives are neccesary in order to prevent spurious oscillations around steep gradients (Dorfi 1997). Simple approximations such as cubic splines (Press et al. 1986) are unreliable, because they may exhibit strong oscillations between the tabulated values. Such spurious oscillations may be even stronger in the derivative of the function inhibiting convergence of the main iterative loop.
In principle, there exists a large number of ways to accomplish such a task. One possibility is to employ rational spline interpolation or the Akima (1970) interpolants. In this Appendix we shall present an alternative method to represent tabular data, which is not an interpolation but a kind of "smooth representation" which offers continuous, analytic derivatives. This method is very simple, being thus an alternative to be used in numerical simulations where the smoothness of the involved functions is of prime importance.
This Appendix will be organized as follows: In Section A1 we present the method and its main properties. Some specific weight functions are presented and studied in Secion A2. Section A3 is devoted to the study of the effect of the smoothing here proposed on the spectrum of frequencies of the interpolation. In Section A4, some examples and comparisons with the most standard procedures are displayed. The last Section of this Appendix ( § A5) will be devoted to general comments on the properties of the method.

A1 The Algorithm
Let us suppose that we have a table of values yn, n = 1, . . . , N of the function y(x), for given xn, not necessarily evenly spaced. A simple way to interpolate a value of y(x) for a non tabulated x (xn < x < xn+1) is by means of a linear interpolation where y n = (yn+1 − yn)/(xn+1 − xn). The linear interpolation is continuous throughout the table but its numerical derivative is a step piecewise function, thus exhibing discontinuities at the tabulated xn. Because of the reasons detailed in Section A, such discontinuities in the derivative are higly undesirable. In what follows we shall present an algorithm specifically designed to avoid this problem in the frame of linear interpolations which can also be straightforwardly extended to interpolations of higher degrees. The central idea of our algorithm is to represent the function y(x) by means of an appropriate weighting of the linear piecewise interpolated values yint with a given weight function. This weight function must be chosen in order to assign high weights to values of yint near the interpolation region, and considering with less relevance those values away this point. Being θ(x, xn, xn+1) the weight function, we define the representation of y(x) as θ(x, xn, xn+1) must fulfill some conditions. It must be positive defined and smooth and if x is far from the interval xn, xn+1, then θ(x, xn, xn+1) → 0. If we require continuous first derivatives of the smooth representation of the tabular data, it is straightforward to show that the first derivative of the weight function θ(x, xn, xn+1) must also be continuous. A natural way of fulfilling all these requirements is to define θ(x, xn, xn+1) as with the exception of the first and the last intervals, for which where w(x, x ) is a function that must verify the normalization condition Notice that w(x, x ) may be chosen in several ways. From a practical point of view, however, Eq. (A3) should be solvable analitically. The values of y θ (x) obtained with our smoothing algorithm can be interpreted as the weighted sum of the contributions of each piece of the tabulated data.
It is important to notice that performing the smoothing in the way we propose, the resulting function does not pass exactly by the tabulated data. In the case that we choose is the Dirac function, integrating it in the corresponding interval we get the step function, which in turn makes us to recover the original linear interpolation.
If the error in the linear interpolation is then, the error in the smooth representation e θ (x) is en(x) θ(x, xn, xn+1).
In the case that θ(x, xn, xn+1) is choosen to be a narrow -peaked function around x, the error in the smooth representation should be comparable to the corresponding to the standard linear interpolation. However, in the case that θ(x, xn, xn+1) were not narrow -peaked, it seems not easy to give a useful bound to the error. Nevertheless, this case is of little interest, simply because the representation of the function would be far from the original values.

A2 Some particular cases of weight functions
A possible choice for w(x, x ) is the rational function, whose expression is with This function is bell-shaped, with its maximum at x = x and with a characteristic width given by the free parameter Γ. As required, Ln(x − x ) has analytical primitive. We will denote the integral In spite that the analytic integral exists for arbitrary value of n, values of n > 2 will result in complicate expressions, expensive from a numerical point of view. Moreover, studying the example proposed in Section A4 we found no obvious advantage in chosing n = 3. Thus, here on, we shall concentrate on cases n = 1, 2.
The expression for yL n (x) is For the case of n = 1 we have whereas, for the case of n = 2 we have In the above expressions

A3 The Spectrum of Frequencies of the Interpolation
The smoothing operation should filter out only those high frequencies present in the tabulated data which could make difficult the convergence of the employed algorithms. Otherwise, there is the risk to obtain stable solutions but very different (or even meaningless) from the true one. Thus, let us study the problem of the modification of the spectrum of frequencies of the linear interpolation when we apply the process of smoothing we propose in this Appendix. We shall adopt the definitions for the Fourier Transform and Antitransform as stated in, e.g., Duff & Naylor (1966) Applying these definitions to the algorithm presented in Section A2, the Fourier transform of the smooth representation can be written as where Θ(s, xj, xj+1) is the Fourier transform of the weight function θ(x, xj, xj+1) in the corresponding interval. Now, let us study the particular case of the weight functions presented in Section A2. The Fourier transform of the weight function in this case is The last integral has a value independent of the interval. For the particular cases of n = 1, 2 their values are π/Γ exp (−sΓ) and π/(2Γ where Θstep(s, xj , xj+1) is the Fourier transform of the step function in this interval. It is clear that the smoothing of the linear interpolation acts as a lowpass band filter. Thus, for the cases of n = 1, 2 frequencies for which ωΓ ∼ > 1 are strongly damped.

A4 A Numerical Example
In this section we shall present an example of how our proposed algorithm works in the case of being employed to interpolate a definite function. For such purpose, we have chosen the values of radiative opacities relevant for the physical conditions attained in the outer envelope of giant planets (Guillot, et al. 1994). In particular, the selected curve correspond to element abundances corresponding to Jupiter at a density of ρ = 10 −14 g cm −3 . The curve is defined by 17 points non evenly spaced. The maximum and minimum intervals are 0.180 and 0.039 respectively. In order to test the algoirithm we presented above, we shall employ the two weight functions considered in Section A2 for cases n = 1, 2 assuming values of Γ similar to those corresponding to the table spacing.
We show in Fig. A1 the results of applying the algorithm we propose for the case of n = 1 for the data set cited above. In doing so we have assumed values for the free parameter Γ of Γ = 0.03, 0.06, 0.12, 0.24. The results presented in Fig. A1 clearly indicate that the method discussed in the present work is poor and not recommendable for the case of n = 1. The overall behaviour of the method could be improved restricting the interval of x but this would imply the introduction of highly undesirable boundary effects. Because of this reason, we shall not consider the case n = 1 any further.
Let us now consider the case n = 2. The results of the interpolation of the function are depicted in Fig. A2 . From inspection of this figure it is quite clear that in this case the overall behaviour of the interpolation largely supersede that corresponding to the case of n = 1. Only for the largest value of Γ there are some significative differences near the absolute minimum of the function.
For comparison, we have also included in such figures results corresponding to interpolation of the data at the same Figure A1. The interpolation of the opacity data values (represented by filled circles by means of the weigth function defined by Eq. (A12) with values of Γ = 0.03, 0.06, 0.12, 0.24. The larger the Γ value, the worse the interpolated value. Also the results predicted by the Akima (1970) spline algorithm are included with short dashed lines. It is clear that the weight function Eq. (A12) provides a poor smoothed interpolation. This is so mainly because of its asymptotic behaviour (see text for more details).
points performed with the Akima (1970) cubic spline algorithm which constructs a piecewise cubic polynomial curve close to a manually drawn curve. Notice that for small values of Γ, comparable to the minimum spacing of the present data tabulation, our algorithm with n = 2 produces results almost indistinguishable from those corresponding to the Akima spline.
In Fig. A3 we show the results corresponding to the derivative of the function for case n = 2 and the set of values of Γ employed in Figs. A1-A2 together with those given by the Akima spline. In spite that the Akima spline gives a fairly good interpolation as shown in Fig. A2, it is clear that it produces a derivative of the function appreciably more oscillating than the one obtained with the smoothing here proposed, even for the smallest considered Γ values. This is, as discused above, a highly desirable property at handling interpolation algorithms in iterative schemes.

A5 Some General Remarks on the Algorithm
In this Appendix we have presented a simple smoothed linear interpolation algorithm intendend to represent tabulated functions. We presented a general treatment, and studied two particular cases in which we considered specific weight factors resulting from the integration of adequate rational Figure A2. The interpolation of the opacity data values (represented by filled circles) by means of the weigth function defined by Eq. (A13) for the same values of Γ as in Fig. A1. The larger the Γ value, the worse the interpolated value. The weight function Eq. (A13) provides a very good smoothed interpolation that, for the case of low Γ values is very similar to that given by the Akima spline algorithm (short dashed lines).
functions in each interval of the tabulated data (see Section A2).
In order to test the performance of our algorithm, we have chosen a particular data set and computed interpolated values for different values of the free parameter Γ. We have found that, as expected, the results of the interpolation are closer to the original linear interpolation as Γ → 0, while the function is smoother the larger Γ is. For the case n = 1 results are poor, while for n = 2 the algorithm works nicely when we assume Γ values not larger than data spacing. Moreover, our technique provides interpolated values very similar to those given by the Akima spline interopolation ( Fig. A2) but smoother derivatives. (Fig. A3).
As our main motivation in developing this alghorithm was to have a reliable inetrpolation technique adequate to employ as a routine in implicit relaxation calculation, we see an interesting property of this technique as compared to that presented by Akima. Let us suppose that we are trying to solve a given system of equations by an iterative technique that fails to converge. Then, in the frame of this algorithm, we can, in principle adopt a very large Γ value that will provide a very inexact but smooth interpolation. After we get convergence we can compute a sequence of solutions of the iteration with lowering the Γ value in several steps. In such a way we can embbed the initial artificial solution to produce a realistic one. In the case of the Akima spline this Figure A3. The derivative of the interpolation of the opacity data values assuming the weigth function defined by Eq. (A13) for the same values of Γ together with the derivative given by the Akima spline interpolant. Notice that Akima algorithm gives strongly oscillating derivatives is not possible, simply because there is no free parameter available to be handled in an equivalent way.
Also, it should be noticed that the algorithm presented in this Appendix can be straightforwardly generalized for representing multidimensional functions. This is so if we choose weight functions in separate variables. This paper has been typeset from a T E X/ L A T E X file prepared by the author.