Analytic solutions for spherical gravitational gas accretion on to a solid body

The process of gravitational accretion of initially homogeneous gas on to a solid ball is studied using the methods of ﬂuid dynamics. The ﬂuid partial differential equations for polytropic ﬂow can be solved exactly in an early stage, but this solution soon becomes discontinuous and gives rise to a shock wave. Afterwards, there is a crossover between two intermediate asymptotic self-similar regimes, where the shock wave propagates outwards according to two similarity laws, initially accelerating and then decelerating (and eventually vanishing). Lastly, we study the ﬁnal static state. Our purpose is to attain a global picture of the process.


I N T RO D U C T I O N
The problem of spherical infall of gas on to a solid object has applications in several astrophysical situations. For example, it can represent the infall of gas on to a neutron star from the surrounding cloud. Alternatively, it can be a simplified model for the formation of planetary atmospheres as a result of gas accretion on to the solid (rock) surface of a previously formed spherical core.
Aspects of this fluid dynamical problem have been treated previously in the literature. A simplified formulation with constant gravity was considered by Bisnovatyi-Kogan, Zel'dovich & Nadezhin (1972) and a self-similar solution with a shock wave was found. Self-similar solutions with (variable) gravity arising from a central point-like mass, applicable to either ordinary or neutron stars, were studied by Cheng (1977). Similarity methods favour power-law distributions. This property was amply used in the study of supernova explosion and later fallback by Chevalier (1989). An extensive study of self-similar spherical accretion with an initial power-law radial distribution of gas was provided by Kazhdan & Murzina (1994). They concluded that a realistic solution should be an interpolation between a self-similar solution with zero mass flux near the origin (of which they derived general features) and a solution with constant gravity and the correct boundary conditions at the solid surface. Chevalier (1989) and Kazhdan & Murzina (1994) assume that the incident gas flow is cold (T = 0) and, respectively, that its accretion rate decays with time with some exponent or that its density decays with radius with exponent ω < 3 (both conditions are related). A cold inflow is not a good approximation for large distances where the gas thermal energy is non-negligible with respect to its gravitational energy. On the other hand, the condition ω < 3 implies that the mass E-mail: gaite@laeff.esa.es of gas diverges at large radius, making the neglect of its self-gravity questionable.
One may try to amend these problems by endowing the incident gas with pressure. For polytropic flow, similarity demands that the initial pressure also be a power law with radius and, furthermore, determines the exponents (in terms of the polytropic exponent). This is the case studied by Cheng (1977) and it may seem somewhat unnatural. Besides, the exponent of the density distribution also leads to a diverging mass with radius and, hence, to the self-gravity problem.
In this paper, we relax the demand of having similarity at the outset and consider instead simple uniform initial conditions for the gas. In other words, the initial configuration is a solid body with a spherically symmetric mass distribution (a ball) placed in a homogeneous gas with pressure, and the problem is to analyse the subsequent evolution of this gas under the action of the gravity of the body. The evolution will consist of infall of gas with a progressive modification of the gas distribution around the body, while the gas stays homogeneous far from it. This is essentially the type of accretion treated by Bondi (1952), although the artificial imposition of stationary flow allowed him to dispense with differential equations. At any rate, a homogeneous distribution of gas is just a particular case of a power law (ω = 0) and, therefore, comparison with the results of Chevalier (1989) and Kazhdan & Murzina (1994) will be possible.
These simple initial conditions are suitable for an analytic treatment, employing the full power of the methods of non-linear fluid dynamics. These allow us to attain a clear picture of the initial nontrivial process near the surface. In particular, the exact description of the formation of a shock wave may illustrate its general features in accretion processes. Even though the evolution is not self-similar, similarity arises during the course of it, and we shall see how and why this happens. Due attention is also paid to the self-gravity question, and to the decay towards a final static state.
The paper is organized as follows. In Section 2, we introduce the relevant magnitudes in the problem and, hence, various space and time-scales, to reach a rough intuitive idea of the physics involved and to define the limits of applicability of the model. In Section 3, we introduce the fluid equations to be used. In Section 4, an exact solution is found for the initial non-trivial dynamics, which occurs near the surface of the ball. In the following section we obtain two similarity solutions valid for larger t, the first still confined within a short distance from the ball, while the second is valid for large radius. Finally, we consider the long-time asymptotic static state and discuss the results.
A note on notation: we shall frequently use the asymptotic signs ∼ and ≈; for example, f (x) ∼ g(x) or f (x) ≈ g(x) (sometimes without making the independent variable x explicit). The former means that lim f (x)/g(x) when x goes to zero or infinity, as the case may be, is finite, while the latter means, in addition, that the limit is one. On the other hand, the sign ∼ may appear in some instances with the related but looser meaning 'equal up to a numerical factor of the order of unity'.

R E L E VA N T M AG N I T U D E S A N D D I M E N S I O NA L A NA LY S I S
The initial condition (the solid ball in the homogeneous gas) is characterized by four parameters, namely, the radius R and mass M of the ball, and the pressure P 0 and density ρ 0 of the gas (we assume that the gas is perfect, inviscid and non-heat-conducting or polytropic). In addition, we have the constant of gravity G. From these five dimensional characteristic parameters, we can form two independent non-dimensional numbers, namely, the ratio of densities 4 3 πR 3 (ρ 0 /M) and the ratio RP 0 /(ρ 0 GM). The quantity P 0 /ρ 0 ∼ c 2 0 (c 0 being the speed of sound) is approximately the gas thermal energy per unit mass. On the other hand, GM/R is the gas potential energy per unit mass on the ball. Their ratio measures the relative strength of gravity in our problem. To have any significant gas infall, we must demand that the gravity of the ball dominates over the gas thermal energy, i.e. Rc 2 0 /(GM) 1. The ratio of densities is also very small. This is a necessary condition for neglecting the self-gravity of the gas, as we will do, but it is not sufficient. We shall discuss the sufficient condition after analysing the typical length and time-scales.
The basic length-scale is R, of course. We can obtain a larger length-scale by dividing by the small number Rc 2 0 /(GM): This length R marks the scale at which the gas thermal energy is similar to its potential energy. It was introduced by Bondi (1952) to define a non-dimensional radial variable. Alternatively, we may divide R by the cube root of the other small number 4 3 πR 3 (ρ 0 /M): in which we may suppress the numerical factor. 1 Obviously, this second length is the radius of a volume of gas such that its mass is similar to M. Analogously, we have a basic time-scale, namely, R 3/2 / √ G M, defined only in terms of the parameters of the ball (and interpreted as the typical time of Keplerian motion close to the surface of the ball). Dividing by [Rc 2 0 /(GM)] 3/2 we cancel the R dependence and obtain GM/c 3 0 , i.e. the time of sound propagation over the distance R (note that the time of Keplerian motion at distance R is R 3/2 / √ G M = G M/c 3 0 as well). Finally, dividing by [R 3 (ρ 0 /M)] 1/2 , we cancel both the R and M dependences and obtain 1/ √ Gρ 0 , the typical time of gravitational collapse of the gas (aside from the ball, which might act as a seed for the collapse).
We will assume that the largest typical length (the radius of a mass of gas M) is much larger than the intermediate scale R; in other words, we consider the mass of gas in the volume of radius R to be negligible in comparison with M. Then, analogously, the typical time of gas collapse is much larger than GM/c 3 0 . This is the precise condition for neglecting the self-gravity of the gas, understood in an asymptotic sense, namely, R (M/ρ 0 ) 1/3 or G M/c 3 0 1/ √ Gρ 0 . In non-dimensional form, c 6 0 /(G 3 M 2 ρ 0 ) 1, which can be restated in terms of the two previously defined non-dimensional numbers: besides that both must be very small, the ratio that is, the ratio of densities must be much smaller than the other one. Interestingly, under this condition, the parameters M and G will not appear independently in the equations of motion but only in the combination GM, so that one has only four dimensional parameters. Correspondingly, only one non-dimensional number is relevant, namely, the ratio R/R. As regards mass scales, the basic mass is M, of course. The mass of gas enclosed in the sphere of radius R, namely, M ≈ 4 3 πR 3 ρ 0 , can be interpreted as the total mass of initially bound gas. We observe that M/M ∼ R 3 ρ 0 /M 1, according to condition (3). On the other hand, this condition can also be written as M J M, introducing the gas Jeans mass M J ∼ c 3 0 /(G 3 ρ 0 ) 1/2 . This is, of course, the natural mass scale for gas collapse owing to self-gravity.
It is convenient to remark that, since (after neglecting the gas self-gravity) we still have one non-dimensional number (R/R), we can construct many length and time-scales (large or small), but they have no particular physical meaning. However, multiplying the basic length R by that number, we obtain R 2 c 2 0 /(GM) = c 2 0 /g, where g is the acceleration due to gravity on the surface of the ball. This small length-scale and its associated time-scale c 0 /g will play a role in the following.

F L U I D E Q UAT I O N S
Under the assumption of spherical symmetry, we are led to solve the partial differential equations of fluid dynamics in one dimension (the radial distance). These equations are non-linear and no general method of solution is available. However, given the simplicity of the initial and boundary conditions, many results can be obtained by purely analytic means, as we shall see.
We consider an adiabatic evolution of the gas or, more generally, a polytropic equation of state, P ∝ ρ γ (γ 1). Then, we have the continuity equation, the Euler equation and the thermodynamic equation: where g = GM/x 2 is the gravitational acceleration (Chevalier 1989;Kazhdan & Murzina 1994). The initial conditions are ρ(0, x) = ρ 0 , P(0, x) = P 0 and v(0, x) = 0 (x R). The boundary condition at the solid surface is v(t, R) = 0. In principle, equation (6) has the trivial solution P/ρ γ = P 0 /ρ γ 0 . Let us introduce the velocity of sound c, which satisfies Hence, the mathematical problem boils down to solving The polytropic gas flow in one dimension is amenable to powerful mathematical methods (Courant & Friedrichs 1948;Landau & Lifshitz 1987;Chorin & Marsden 1993), although the presence of gravitation complicates the problem. However, restricting our interest to a small zone near the surface, we can take the acceleration of gravity g to be constant in equation (9) (Bisnovatyi-Kogan et al. 1972). This will allow us to take advantage of those methods.

I N I T I A L S TAG E S
Initially, the gas will start falling with acceleration g, that is, with a negative velocity increasing as v ≈ −gt. However, it will be stopped at the surface of the ball, where the density and, therefore, the pressure must increase. Consequently, a wave transmitting the boundary condition v(t, R) = 0 will propagate outwards. It is therefore crucial to determine the law of propagation of waves in the present conditions. This is called, in mathematical terms, the analysis of characteristics. Once this is done, and as long as the dynamics consists of the propagation of a simple wave, one can obtain an exact solution. It is not possible here to describe in detail how this solution is obtained and our focus will be the origin of the shock wave. We refer the reader to Courant & Friedrichs (1948), Landau & Lifshitz (1987) or Chorin & Marsden (1993) for a general treatment of the theory of one-dimensional gas flow.
For the moment, we confine ourselves to a spherical shell over a surface of height much smaller than R, where we can consider g to be constant. Hence, the problem becomes that of the fall of gas on to a flat surface (the ground) and it is convenient to take this surface as the origin of x. So we must use the one-dimensional form of the vector divergence and, therefore, we must neglect the final term on the left-hand side of the continuity equation (8), writing it as The gas still unperturbed by the wave, merely falls with velocity v = −gt. In a reference system falling with it, it is at rest and, therefore, the speed of sound is c 0 . Hence, in the original reference system the wave front is located at x f = c 0 t − gt 2 /2. Then, for x > x f , the solution is trivial, and we only need to find the solution for x < x f . In the falling frame, the acceleration of gravity vanishes and the set of two equations (10) and (9) reduces to a 'gas tube problem' that can be solved by introduction of the Riemann invariants (Courant & Friedrichs 1948;Landau & Lifshitz 1987;Chorin & Marsden 1993). Let us denote by x , v and c the coordinate and variables in this frame. Initially, we have a constant state, which is preserved for The backward characteristics crossing the forward characteristic x = c 0 t transmit a constant value of the Riemann invariant J − , so the solution is a forward simple wave. Given that the velocity of sound is related to the gas velocity by c = c 0 + (γ − 1)v /2. Furthermore, the fact that the solution is a forward simple wave implies that the forward characteristics are straight lines. Hence, the propagation law of the wave v = F(x − c t) is an implicit equation for v , assuming that we can determine the function F. This is performed using the boundary condition v(t, 0) = 0. The implicit equation is algebraic and its solution is straightforward. We For a simple wave, the density is a definite function of the velocity. We can calculate it from equation (7) to be These expressions for ρ and v constitute the solution of equations (9) and (10) with the given boundary conditions. The alert reader may have noted a peculiarity of the solution found above: the coordinate of the wave front x f begins increasing but, for t > c 0 /g, changes to decreasing and, eventually, returns to the origin (free fall with initial velocity c 0 ). On the other hand, a simple analysis shows that a singularity occurs before this, namely, a shock wave. One can calculate (∂x/∂v) t , and its null locus defines a line in the xt-plane, namely, x = [2c 0 + (γ − 1)gt] 2 /(8γ g). Initially, the corresponding x is larger than x f , so it is unphysical, but, for t = 2c 0 /[(γ + 1)g], it meets the wave front; that is, the falling gas overtakes the gas just below it and a shock wave develops henceforth. From that moment onwards, the solution (11) and (12) becomes multivalued and ceases to be valid (see Fig. 1). Although it is possible, in principle, to study the situation in which a shock wave propagates, the analysis cannot be performed in terms of a simple wave and it is far more complicated.
Of course, the foregoing analysis is valid as long as everything happens in a thin shell, which implies that the typical length to R R (the already mentioned condition of gravity-dominated evolution).

S I M I L A R I T Y S O L U T I O N S
A one-dimensional gas flow problem may be formulated with similarity variables [in which the fluid partial differential equations become ordinary differential equations (ODEs)] only when the initial and boundary conditions are sufficiently simple (Sedov 1982). Certainly, this is not our case, for we can construct independent length and time-scales (see Section 2). Nevertheless, it is commonly observed that non-linear equations that do not have similarity solutions at the outset develop them in an intermediate asymptotic regime, that is, a regime between two very different scales where both scales can be neglected (Barenblatt 1996). We will actually find two independent intermediate asymptotic regimes.
First, we shall study the self-similar solution that develops from the exact solution found in Section 4, i.e. with constant gravity. The corresponding ODEs have been formulated by Bisnovatyi-Kogan et al. (1972). However, we shall derive a slightly different (though equivalent) system of equations more amenable to analytic study.
Afterwards, we proceed to consider the variation of gravity with radius, connecting with the similarity solutions of Chevalier (1989) and Kazhdan & Murzina (1994). These solutions are for an initial density distribution given by ρ = Bx −ω , but a constant density is just a particular case. The most interesting solutions are those with zero mass flux near the origin, named 'of type 3' (Kazhdan & Murzina 1994), since they are expected to match the self-similar solution with constant gravity.

Constant gravity
In this case, we have one less parameter, because the parameters R and GM only appear in the combination g = GM/R 2 . However, we can still form the two independent non-dimensional variables gt/c 0 and gx/c 2 0 , so we will have to remove another parameter to have a similarity solution. If we assume that the initial gas temperature is very low, then c 0 → 0 and the only possible non-dimensional variable is ξ = x/(gt 2 ) (or a function of it). (We measure the distance x from the ground.) This corresponds to the intermediate asymptotics To take into account the presence of shock waves and the consequent dissipation, we consider the full set of equations (5), (6) and (10). We can express them in non-dimensional form by introducing Some straightforward algebra then yields Upon making the change of variables (corresponding to changing to the falling reference frame), the ξ −1 term in the second equation vanishes and the system of ordinary differential equations simplifies tõ Bisnovatyi-Kogan et al. (1972) do not make the change to the falling reference frame, and restrict themselves to solving the equations by a series expansion near the ground. Instead, after the change of reference, we can follow the general (and more powerful) methods exposed by Sedov (1982). In particular, we note that in terms of the variable τ = lnξ , the previous equations constitute a system of autonomous non-linear ODEs, namely, One can study this system of equations with the methods of nonlinear ODEs, that is, study its singular points, phase portrait, etc. However, given the homogeneity properties of these equations with respect top and r, it proves convenient to introduce the variablẽ θ =p/r . It is (in the falling frame) the non-dimensional form of the temperature for a perfect gas with pressure P and density ρ, i.e. T = µ (x 2 /t 2 ) θ, where µ is the molar mass and T is measured in energy units. Hence, The solution of this equation provides the relation between the 'velocity'ũ and the 'temperature'θ,θ(ũ). Substituting it back into equation (21), we have an ODE forũ(τ ), which is immediately solved by a quadrature. Analogously, one solves for r (τ ).
Let us see whether we have gathered enough information to determine a unique solution of equation (24). Given thatθ → 0 as ξ → ∞, this must be a fixed point. There are three singular points of equation (24) withθ = 0, namely, withũ = 0, 1 or 2, respectively. If we want to have a finite v as x → ∞, we must take the point (0, 0). So we know that the solution of equation (24) that we seek must end at the origin. On the other hand, the singular point (2, 0) corresponds to the boundary conditions at x = 0. However, the solution that departs from this point does not end at the origin (see Fig. 2). This does not mean that there is an inconsistency, as there can be discontinuities in the solution. Indeed, in most self-similar solutions the boundary conditions can only be fulfilled if there is a point of discontinuity (Sedov 1982). Therefore, recalling the exact solution of Section 4, we observe that a shock wave arises and propagates upwards. Then we must consider discontinuities associated with the presence of a shock wave. The conditions for conservation of mass, momentum and energy across a shock surface read (Courant & Friedrichs 1948;Landau & Lifshitz 1987;Chorin & Marsden 1993) where subindices refer to the values on either side of the shock and w is the enthalpy (or the appropriate thermodynamic function for general γ ); for a perfect gas, w = γ P/[ρ (γ − 1)]. These equations hold in a coordinate system in which the shock surface is at rest. On the other hand, the location of the shock must be at fixed ξ , say ξ s . Consequently, the shock wave velocity is i.e. u s = 2. Then, r 1 (u 1 − 2) = r 2 (u 2 − 2), valid in both the rest and the falling reference frames. Using the falling frame and eliminating r 2 /r 1 between the first and second equations, We have arrived at an involutive mapping of the plane (ũ,θ ) as the relation between the variables at either side of the shock discontinuity. Since we have the condition that as x → ∞ the solution goes to (0, 0) and (in the falling frame) these values hold all the way down to the shock surface, we just need the image of the origin under the mapping. This yields the point Therefore, we need the solution of equation (24) that goes through this point. We can see in an example (Fig. 2) that the strand of this solution that we need ends at the point (2, 0), hence satisfying the boundary condition at x = 0. 2 Now, substituting the function just computed back into equation (21), we obtain dτ = γθ(ũ) − (ũ − 2) 2 −θ(ũ)(2 + γũ) +ũ(2 − 3ũ +ũ 2 ) dũ, from which we can deduce the interval of τ between the points (4/(γ + 1), 8(γ − 1)/(γ + 1) 2 ) and (2, 0) by integration. The result for γ = 7 5 (the adiabatic index of a perfect diatomic gas) is 0.088 0104. Hence, we obtain the quotient between the corresponding values ofξ . This quotient gives the location of the shock relative to the ground: recovering the original notation that distinguishes both reference frames,ξ s /(1/2) = exp(0.0880104) = 1.09200 ⇒ ξ s = 0.0919995/2 = 0.0459997, so that the coordinate of the shock is x s = 0.045 9997gt 2 .
It is easy to compute the post-shock u in the rest frame: u = [(1.092/2)(5/3) − 1)]/0.046 = −1.957. Hence, the post-shock velocity is v = (x s /t) u = ξ s ugt = −0.09gt. Comparing with the preshock value v = −gt, we see a great reduction. The corresponding kinetic energy is dissipated as heat. The similarity properties of the velocity and its jump at the shock are shown in Fig. 3.

Variable gravity
If we consider heights of the order of the planet radius or larger, we have to return to the full continuity equation (4) and reset g = GM/x 2 in the Euler equation (5) (x is again the distance from the centre). We can find similarity solutions in the intermediate asymptotics R x R, in terms of the non-dimensional variable ξ = x 3 /(GMt 2 ). Hence, we can express the continuity and Euler equations as The energy equation becomes These three equations look similar to those corresponding to constant gravity, and they can also be reduced to two differential equations for u(ξ ) and θ(ξ ), respectively, (Cheng 1977;Kazhdan & Murzina 1994). Unfortunately, they are more difficult to analyse: the change of variables that transformed the constant-gravity equations into an autonomous system is not available. In fact, the free-fall (pressureless) problem is now non-trivial, although it can be solved.
If p = 0, equation (30) decouples, becoming an equation for just u(ξ ), namely, To implement the PDE initial condition, at ξ → ∞, it is convenient to make the change of variableξ = ξ −1 , transforming the equation The point (0, 0) is a singular point of the ODE: a solution is unspecified, unless we introduce an additional condition. The initial condition v = 0 implies that v = √ G M/xξ −1/2 u goes to zero aŝ ξ → 0. To impose this, it is best to linearize and solve the ODE around (0, 0), obtaining so that the singular point is a node. Since limˆξ →0ξ −1/2 u = C, the particular solution of the non-linear equation (33) in which we are interested is the only one that satisfies u (0) = −1 (C = 0).
The full free-fall solution can be obtained in parametric form using the Lagrangian formalism (see Appendix A). It reads where α ∈ (0, π/2). The corresponding graphs for u(ξ ) and r (ξ ) are shown in Fig. 4. This free-fall solution is to be distinguished from the simpler freefall solution considered by Cheng (1977) or Kazhdan & Murzina (1994), namely, u = − 2 / ξ , which corresponds to vanishing initial energy or, in other words, to the initial configuration consisting of gas particles at rest at infinite distance. Note that the dominant term of the solution (34) is proportional to √ 1/ξ and coincides with their free-fall solution for C = − √ 2, being then an exact solution of the non-linear equation (33). Furthermore, for other non-null values of C, we have different self-similar asymptotic solutions, with a given initial velocity proportional to √ G M/x. On the other hand, we can derive u = − 2 / ξ from the parametric solution in the limit ξ → 0 (corresponding to α → π/2 and meaning small distance or long time), since the initial energy of the gas particles becomes negligible in this limit.
The complete similarity solution consists of the free-fall solution and an inner solution with pressure matching across a spherical shock wave. The boundary conditions of the PDE at x = 0 determine the form of the inner solution of the two differential equations for u(ξ ) and θ(ξ ) near ξ = 0. It must be a solution of 'type 3', with one free parameter (Kazhdan & Murzina 1994): where u 0 and θ 0 are related by [3(γ − 1)β + 3γ − 4]u 0 = 2βθ 0 and β = − 1 3 + 2 3 √ (γ − 1)/γ . The match with the previous free-fall solution across the shock determines the free parameter. There is no straightforward way to find the shock coordinate ξ s and a sort of 'trial and error' procedure is necessary (Chevalier 1989;Kazhdan & Murzina 1994). An approximate value can be obtained by matching the expressions (38), namely, ξ s = 0.013.
Both the pre-shock and the post-shock gas velocities are given by v(x s , t) = (GMξ s ) 1/3 u(ξ s )t −1/3 , for the respective u values. Hence, the similarity mass flow across the shock isṁ = 4πx 2 ρv = 4πG Mρ 0 ξ s u(ξ s )r (ξ s )t, which increases with time. This may seem paradoxical, since the velocity decreases as t −1/3 , but the increase in shell mass arising from the spherical geometry compensates for it.
We note that the shock velocity experiences at x s ∼ R a crossover from v s = 2x s /t = 2 ξ s gt to v s = 2/3 (x s /t) = 2/3 (GMξ s ) 1/3 t −1/3 (with different ξ s ), turning from accelerating to decelerating. The previous solution is valid for x s − R R, whereas the present solution is valid for x s R. Once the second solution takes hold, given that the velocity of sound for the far-away gas is c 0 , the shock disappears after its velocity becomes subsonic. This occurs for (GMξ s ) 1/3 t −1/3 ∼ c 0 ⇒ t ∼ GM/c 3 0 and x s ∼ R.

T H E S TAT I C S TAT E
Since the shock wave eventually vanishes, leaving the gas with more entropy than in the initial state, we expect that it must reach a stationary state in the long run. The condition that v(x) = 0 (from the continuity equation and the boundary condition) implies that the stationary solution is actually static. The hydrostatic equation, with P ∝ ρ γ and g constant, has the solution: where the subscript 'b' denotes values at the surface of the ball. Actually, this formula describes well the lower region of the present atmosphere of the Earth. Nevertheless, it is unphysical when taken over the entire range of x (for γ > 1, that is, for a non-isothermal distribution): ρ γ −1 (which is proportional to the temperature) decreases linearly with height, becoming negative for some height. It is natural that there cannot be a stationary state with constant gravity, given that nothing can prevent the free fall of the gas. In contrast, if we take variable g = GM/x 2 , In this case, it is more convenient to take the reference at infinity, where the density and pressure keep their initial values, so that As expected, R gives the sphere out of which the gas can be considered gravitationally unbound (since its thermal plus gravitational energy is positive) and, consistently, the scale of crossover to the asymptotic (x → ∞) density ρ 0 .
For small radius x R, we can neglect the unity in the righthand side of equation (41), so that we have power laws. In particular, the density can be written as by introducing n according to γ = 1 + 1/n. It coincides with Cheng's hydrostatic solution, after substituting ρ 0 R n = (G M) n γ −n ρ n+1 0 P −n 0 . However, Cheng's 0 subscript indicates arbitrary reference values, given that ρ n+1 P −n is constant. Therefore, we see that it is convenient to take the asymptotic values as a reference, as long as they are non-vanishing, R being the scale where ρ and P take those values. Moreover, we can let those asymptotic values vanish, by taking R → ∞ while the product ρ 0 R n has a finite limit. On the other hand, since the ratio between the surface values of temperature or density and the initial values is given by 1, we may prefer to choose the former values as a reference. Then we can analogously take the limit R → 0 and ρ b → ∞, while keeping ρ b R n constant. Since we have pure power laws when R/R → ∞, the reference is arbitrary.
We must remark that the law (41) gives a density contrast ρ − ρ 0 that decreases too slowly with the radius, yielding a divergent accumulated mass of gas as x → ∞. In fact, the exact form of the hydrostatic equation involves the gravity of the gas and, hence, is equivalent to the equation that describes stellar structure, namely, the Lane-Emden equation (Chandrasekhar 1967). If the bound mass of gas is much smaller than the mass of solid material, that is, if R (M/ρ 0 ) 1/3 , the solution of this equation is well approximated by equation (41), while for x > (M/ρ 0 ) 1/3 the decay is faster than the decay ∼1/x given by it. One can then approximate the mass of accreted gas by the total mass of gas inside the sphere with radius R, which is larger than the initially bound gas mass, M, because a portion of the gas that is initially outside the sphere falls inside it. The detailed calculation in Appendix B shows that the ratio of this portion to M is of the order of unity for γ 4 3 (in particular, it takes the value 2.65 for γ = 7 5 ). We have seen in the previous section that the spherical shock vanishes when it reaches a radius ∼R. This is consistent with the fact that the gas outside a sphere of radius R is immune to the gravitational influence of the planet. On the other hand, since the shock vanishes for t ∼ GM/c 3 0 , this implies that the similarity solution becomes useless for t > GM/c 3 0 , marking the crossover to the static state. Nevertheless, it is necessary to remark that the static state is an intermediate asymptotic regime not only in space, as was commented above, but also in time, it being always required that t 1/ √ Gρ 0 .

D I S C U S S I O N
We have performed a thorough analysis of a simple model for gas accretion on to a solid body. We have seen that there are two relevant length scales, namely, R and (M/ρ 0 ) 1/3 , the relative values of which determines whether the dynamics consists of accretion or collapse: to neglect self-gravity, it is required that R (M/ρ 0 ) 1/3 . We have not restricted ourselves to similarity solutions. We have divided the space dimension (the radius x) into a near zone (near the surface of the ball, where x − R R and g is constant) and a far zone (x R radius and an exact static state for very long times. These solutions provide a good intuitive understanding of the entire dynamics. The crucial feature is the formation and propagation of a shock wave until its dissipation. The shock wave arises over a short time ∼c 0 /g and grows rapidly. As long as its propagation is strongly supersonic, we can consider the gas to be cold (c 0 = 0) or, in other words, the dynamics to be driven by gravity. The shock wave experiences a crossover between one stage of strong dissipation and acceleration to another of slowing down, owing to decreasing gravitational energy input, until it eventually vanishes, leaving behind dense and hot accreted gas in a quasi-static distribution.
The process of accretion begins at t = 0 and continues until the end, when v → 0, but it can be considered finished when t ∼ GM/c 3 0 . Since the velocity experiences a strong reduction at the shock, with a consequent increase in the density, it is sensible to define the accretion rate as the value of the mass flowṁ there. We have seen that this value is proportional to t during the second selfsimilar regime, so the amount of mass being accreted only begins to decrease at the end of it. The total accreted mass in the static state is not very large if γ 4 3 (γ = 4 3 for the blackbody opacity limit), in spite of the fact that the corresponding ratio of surface density to ρ 0 is very large {[R/(n R)] n }.
Although an adequate polytropic index accounts for some form of heat conduction or radiation, the polytropic static state for a heat-conducting or radiating gas is unstable against further cooling. Therefore, we must remark that the true static state, from a strictly theoretical thermodynamic point of view, is only achieved by cooling until reaching an isothermal state. The corresponding isothermal density distribution is much more concentrated near the surface of the ball (it is exponential), leading to a far larger accreted mass. However, the time to reach this isothermal distribution would be, arguably, larger than the gas collapse time.
We may wonder how our picture would be modified by different initial conditions, namely, by a non-homogeneous density distribution (e.g. a power law). Gravity only has an effect inside the sphere of radius R, in which the gas rapidly becomes supersonic and approaches free fall, until the shock wave arrives. The free-fall solution (35)-(37) is independent of the value of ρ 0 , which can be a function of x (see Appendix A), but the form of the shock wave and the post-shock gas distribution will change.
We may also consider an initial non-homogeneous temperature distribution. Near the body surface, we can make a linear approximation to it. If its slope is sufficiently smaller (in absolute value) than the static value provided by equation (39), we can still consider it isothermal and predict that the fall of gas will produce a shock wave. Nevertheless, we can envisage distributions close to the static one such that the gas falls gently on to the surface without ever giving rise to a shock wave (such as, for example, the distribution corresponding to the final stage of our process, after the vanishing of the shock wave). Presumably, these distributions lead to little accretion. Now, I briefly estimate the numerical scales suitable for the application to neutron star accretion or planet accretion. For a neutron star of 1 M with R 10 4 m, g 10 12 m s −2 . The speed of sound (proportional to the square root of the temperature) can be large, say c 0 10 000 m s −1 . Then, c 0 /g 10 −8 s and c 2 0 /g 10 −4 m, R 3/2 / √ G M = √ R/g 10 −4 s, but R 10 12 m and GM/c 3 0 10 8 s. It is not clear how to determine ρ 0 . A 'typical' galactic value would give a collapse time similar to star formation times, i.e., of the order of a million years. Considering a planet placed in the protoplanetary nebula, we may take for g the Earth value of 10 m s −2 and, for example, c 0 300 m s −1 (its current value in the low at-mosphere of the Earth); the initial typical scales of the problem, namely, the time c 0 /g and the distance c 2 0 /g, take values of 30 s and 9 km, respectively. The basic time-scale is √ R/g 800 s (for an Earth-like planet with R 6000 km). We further estimate R G M/c 2 0 5 × 10 6 km and GM/c 3 0 5 × 10 7 s 1 yr. A plausible value for the density of the protoplanetary nebula is ρ 0 ∼ 10 −11 g cm −3 . We see that its collapse time 1/ √ Gρ 0 10 9 s is sufficiently larger than GM/c 3 0 to neglect self-gravity.

AC K N OW L E D G M E N T S
I thank A. Mancho, J. M. Martín-García, J. Pérez-Mercader and M. P. Zorzano for conversations, A. Domínguez for comments on the manuscript, and R. Gómez-Blanco for reading a preliminary version of the manuscript. This work is supported by a 'Ramóny Cajal' contract and by grant BFM2002-01014 of the Ministerio de Ciencia y Tecnología.
Inverting (A7) we obtain x/x 0 as a function of ξ , which upon substitution in (A4) yields its self-similar form. This is the solution of the ODE that satisfies the appropriate initial condition (see Section 5.2). The method exposed above works in general for free self-similar motion in an external gravitational field. In the present case, one can calculate the integral in equation (A6) by the change of variables σ = cos 2 ϕ. This allows one to obtain the solution of the differential equation (A1) in parametric form: let and u = − cos −3 (α) α + sin(2α) 2 sin α, where α ∈ (0, π/2). It is easy to verify that this parametric solution fulfils the initial condition lim ξ →∞ u = −1.

A P P E N D I X B : AC C R E T E D M A S S I N T H E S TAT I C S O L U T I O N
Let us calculate the gas mass enclosed in the sphere of radius R in the final static state: 1.5 2 2.5 3 5 10 15 20 Figure B1. Accreted mass ratio as a function of n (γ = 1 + 1/n).
This integral is convergent for R → 0 if γ > 4 3 . We can easily express it in non-dimensional form: Now, it is convenient to make the change γ = 1 + 1/n. The resulting integral can be computed in terms of a hypergeometric function: The value of grows almost linearly, with a small slope, up to near the pole n = 3 (γ = 4 3 ), as Fig. B1 shows. In particular, it takes the value 2.650 99 for γ = 7 5 (n = 2.5). The pole shows the divergence of the integral as s → 0. Then, it is not possible to make R/R = 0 near the pole. The n = 3 integral yields Its value is smaller than 5 even for R/R as large as 10 9 . This paper has been typeset from a T E X/L A T E X file prepared by the author.