Analytic solutions to the accretion of a rotating finite cloud towards a central object I. Newtonian approach

We construct a steady analytic accretion flow model for a finite rotating gas cloud that accretes material to a central gravitational object. The pressure gradients of the flow are considered to be negligible and so, the flow is ballistic. We also assume a steady flow and consider the particles at the boundary of the spherical cloud to be rotating as a rigid body, with a fixed amount of inwards radial velocity. This represents a generalisation to the traditional infinite gas cloud model described by Ulrich (1976). We show that the streamlines and density profiles obtained deviate largely from the ones calculated by Ulrich. The extra freedom in the choice of the parameters on the model can naturally account for the study of protostars formed in dense clusters by triggered mechanisms, where a wide variety of external physical mechanisms determine the boundary conditions. Also, as expected, the model predicts the formation of an equatorial accretion disc about the central object with a radius different from the one calculated by Ulrich (1976).


INTRODUCTION
The first model of a spherical symmetric accretion flow towards a central object was made by Bondi (1952). With time, it has become the key to understand different accretion processes in the universe (see e.g. Frank et al. (2002)), despite the fact that it was created for mathematical curiosity, rather than for astrophysical applications (Bondi 2005).
Bondi thought of an infinite gas cloud with a nonrelativistic gravitational central object located at the origin of coordinates. He assumed the flow to have reached a steady state. Under these assumptions Bondi was able to integrate the hydrodynamic equations for this spherical symmetric flow.
Years later, Ulrich (1976) modified Bondi's ideas by assuming all fluid particles at the border of the cloud to have a certain amount of angular momentum, following a distribution given by a rigid body rotating about the z axis. In his model, Ulrich took no account of pressure effects on the infalling gas. In other words, his analysis was ballistic, which is approximately true if the flow is supersonic and if heating by radiation and viscosity effects are negligible (cf. Mendoza et al. 2004).
The boundary conditions in Ulrich's model imply that ⋆ Now at Departamento de Astronomía, Universidad de Guanajuato, Guanajuato, México all fluid particles move in parabolas whose focus correspond to the origin of coordinates. As soon as particles arrive to the equatorial plane they collide with their symmetric counterparts and stay on it, forming an equatorial disc of radius ru made of fluid particles which rotate about the central object. Ulrich's model is usually taken as the basic model which predicts the formation of an accretion disc of particles orbiting a central object due to accretion with rotation. The solution given by Ulrich (1976) was extended by Casen & Moosman and is used (see e.g. Nagel 2007;Lin & Pringle 1990;Stahler et al. 1994) as a natural initial condition for simulations dealing with the problem of disc formation and subsequent evolution.
Both models, Bondi (1952) and Ulrich (1976) have been taken to the relativistic regime by assuming the central object to be a Schwarzschild black hole (Michel 1972;Huerta & Mendoza 2007). They have wide applications in high energy phenomena (see for example Lee & Ramirez-Ruiz 2006;Beloborodov & Illarionov 2001).
On many astrophysical situations where Ulrich's model is used, the usual assumption is that the boundary conditions for an infinite gas cloud remain valid even when it is applied to the inner core of a large accretion cloud (see e.g. Cassen & Moosman 1981;Lin & Pringle 1990;Stahler et al. 1994). This is not necessarily true, because the accretion flow at the boundary of the inner core has, in general, an inward radial velocity different from zero 1 . Also, if the angular momentum of particles follows rigid body rotation at this finite radius then, as we describe in this article, the flow may largely deviate from Ulrich's solutions, since the specific total energy (kinetic plus gravitational) is not zero and so, orbits are no longer parabolic. The relaxation of Ulrich's boundary conditions leads to a more realistic infall model with applications to a wider variety of accretion problems. With an appropriate selection of boundary conditions this could in principle modify the main results of a number of previous works in the literature dealing with various astrophysical accretion problems.
In this article we assume a central Newtonian potential produced by the central star, and a rotating gas cloud with a finite radius r0. Based on this, we construct a cylindrical symmetric flow which also predicts the existence of an accretion disc about the central object. As it will be shown, the trajectories are no longer sections of parabolas. They are general conical sections and once boundary conditions are fixed they can be of different families for each fluid particle.
In order to show the relevance of having a finite cloud radius r0 and a non-zero radial velocity vr 0 , note that this boundary conditions can be added to the standard Ulrich's ones which are defined by the mass M of the central star, an accretion rateṀ and a specific angular momentum h0 of the cloud. In other words, r0 and vr 0 characterise the modifications to Ulrich's (1976) model. The length r0 can be taken as the outer boundary of the inner region of the cloud, where its self-gravity is neglected, and the main force is given by the gravity of the star.
In real astrophysical problems this model can be applied. For example to NGC 1333 IRAS 4 region where an unusual velocity profile in an infalling envelope is proposed by Choi et al. (2004) to explain the observation. Also, it is relevant in the study of triggering star formation, because whatever the physical mechanism responsible, a whole set of combinations for r0 and vr 0 is possible. In this scenario, Ulrich's model is not applicable. There are well known mechanism for the triggering, namely stellar ejecta from winds (Foster & Boss 1996) or supernova shells (Melioli et al. 2006) as well as the ionising radiation from a star, or a collection of stars, that produces an ionisation shock front (Bertoldi 1989;Esquivel & Raga 2007). Each of these mechanisms sets particular boundary conditions for the collapse of the cloud to form a star.
The triggering is relevant in clouds with a high density of potential in sites of star formation. A first generation of stars can influence the dynamical configuration for the next generation. This mechanism is highly important when it is realized that the stars are formed in turbulent clouds (Vázquez-Semadeni et al. 2007) and that a quiescent core for star formation envisioned by Shu et al. (1987) is in some dense clusters difficult to find. A complete picture of these ideas is clearly presented by Hartmann et al. (2001).
A glimpse of the full range of boundary conditions is exemplified by Hartmann & Burkert (2007). They speculate that the rich Orion nebula cluster is formed by the collapse 1 In fact, as we will see later on, if one builds a stationary situation with zero inward radial velocity, then it would be necessary to have an infinite density at the cloud's border.
of an elliptical rotating sheet of gas with a density gradient along the major axis. A high density region is formed when one tip of the cloud collapses preferentially and so, the already formed protostars (Heitsch et al. 2008) accrete material with velocities given by the global collapse of the cloud, independent of the gravity of the collapsed object. Thus, free choosing of boundary conditions is highly important to study the regions close to the star. The resolution of these regions is beginning to be reached by computer simulations and observations.
In this context, the ballistic solution given in this paper is useful to study star and disc formation in very complex environments, where one has the need to choose boundary conditions freely and easily. However, an obvious weakness in the model is that at r0, the cloud is rigidly rotating. Thus, in the scales where this solution is valid, one either needs to invoke a magnetic field at r > r0 that is responsible to enforce rigid body rotation (Mouschovias & Paleologou 1979), or the fact that the shell feeding the star-disc system rotates rigidly during some reasonable period of time, in order to consider this process as quasi-stationary. However, there are many works in the literature (Kenyon et al. 1993;Jayawardhana et al. 2001;Whitney et al. 2003) that make such an assumption.
In Sections 2 and 3 we calculate the velocity and density fields of the flow respectively. At the end of Section 3 we recover Ulrich's (1976) model. Section 4 describes the behaviour of the flow for its different parameters and discusses its limitations. Finally, Section 5 deals with realistic astrophysical situations and the validity of our model on those scenarios.

THE MODEL AND ITS VELOCITY FIELD
Let us assume that a spherical cloud with radius r0 accretes matter in a steady way towards a central object located at the origin of coordinates. All fluid particles located at r0 are taken to follow rigid body motion in the azimuthal direction, i.e. they rotate about the z axis of coordinates. At this position, particles also have a radial velocity component vr 0 .
Here and in what follows, due to the symmetry of the problem we use spherical coordinates r, θ, and φ for the radial coordinate, and the polar and azimuthal angles respectively. The particle number density n at r0 is taken to be constant with a value n0. To simplify the problem, we assume the self gravity of the gas cloud to be small as compared to the gravitational potential produced by the central object. In addition, if the gradients of pressure are small compared to the kinetic and potential energies, then all fluid particle trajectories can be approximated as ballistic. The total angular momentum of the gas cloud points in the direction of the z axis, which contains the central object. Note that when r0 → ∞ and vr 0 → 0, the model converges to the accretion problem proposed by Ulrich (1976), which is well described by Mendoza et al. (2004). In fact, Ulrich's convergence is mathematically proved at the end of Section 3.
Since the gas is being accreted from r0, it's accretion flow rate per unit massṄ at any fixed radial distance r is constant and given bẏ At r0, the distribution of specific angular momentum h is given by h = h0 sin θ0, where θ0 is the initial polar angle of the particle's trajectory.
Under all the above assumptions, the specific energy E and the specific angular momentum h are constants of motion along each particular trajectory. The specific energy E is given by where G is Newton's gravitational constant. We now introduce two dimensionless parameters µ and ν given by where ru = h 2 0 /GM is the disc's radius in Ulrich's model and E0 := GM/ru is the specific gravitational potential energy of a fluid particle evaluated at ru. The quantity µ represents the ratio of Ulrich's disc radius to the original cloud's radius. So, for example, in order to consider an infinite cloud's radius we must take the limit µ → 0. On the other hand, ν represents the initial amount of radial velocity measured in units of the Keplerian velocity v k := √ E0, at position ru. With these parameters, equation (2) can be written into dimensionless form as where the dimensionless energy ε is given by ε := 2E/E0. The trajectory of each fluid particle is contained on a plane and is given by a conic section. The origin represents one of the foci of the orbit. At a specific initial position, the particle is located at r0 and its polar and azimuthal angles are given by θ0 and φ0. Over the orbit plane, the particle trajectory is defined by an azimuthal angle ϕ which at the initial position has the value ϕ0. The trajectory of a given particle is defined by the solution to Kepler's problem (Landau & Lifshitz 1989): In the previous equation and in what follows, unless stated otherwise, distances are measured in units of ru. The eccentricity e of the orbit is given by At the border of the cloud r = r0 = 1/µ and so, substitution of this in equation (5) gives the following condition for ϕ0: Performing standard spatial rotations its easy to obtain the following formulae between the angles ϕ, ϕ0, θ, θ0, φ and φ0: Using equation (8) to rewrite equation (5) yields where We now use the fact that in spherical coordinates v φ = r sin θdφ/dt = h sin θ0/r sin θ, v θ = rdθ/dt = (dθ/dφ) (v φ / sin θ), and vr = dr/dt = (dr/dθ) (v θ /r). With this and using equation (8) the expressions for the velocity field, measured in units of v k are given by:

PARTICLE NUMBER DENSITY
In order to obtain the particle number density field we start from the continuity equation for a stationary flow, i.e.
We now integrate the previous equation over a volume consisting on a collection of streamlines. In other words, this "tube" is built in such a way that its lateral surface is bound by streamlines and its lids by spherical sections, the upper one at the cloud's border, i.e. r = r0 = 1/µ and the lower at any arbitrary radial distance r such that 0 < r < 1/µ. With this selection, equation (14) can be integrated to obtain: Taking into account the fact that the differential area element da is given by da = r 2 sin θ dθ dφ er + r sin θ dr dφ e θ + r dr dθ e φ , (16) then, equation (15) transforms to sin θ0 dθ0 dφ0 = −n vr r 2 sin θ dθ dφ, where the particle number density n is measured in units of nu :=Ṅ /4πv k r 2 u . Using the symmetry of the quantities in the azimuthal angle involved in the previous equation, we can integrate equation (17) with respect to φ from 0 to 2π, to obtain On the other hand, from equation (9) it follows that −vr sin θ sin θ0 Substitution of equations (19) and (13) into (18) leads to the required particle number density field: As a result of the collision at the equatorial plane between the streamlines coming from the northern hemisphere with those coming from the southern one, a shock front is expected to form, as well as the formation of an equatorial disc shape structure (see e.g. Ulrich 1976; Mendoza et al. 2004, and references therein).
Taking θ → π/2 in equation (9) followed by the substitution θ0 → π/2 we obtain an expression for the accretion disc radius given by: From this last equation and since naturally r d < r0, the following condition between the parameters arises: If in equations (9), (11), (12), (13) and (20) we take the parameters values µ = 0 and ν = 0 , we get the following set of equations which corresponds to the expressions calculated by Ulrich (1976).  Figure 1 shows projected streamlines for different values of the parameters µ and ν with height z and radial coordinate R := r sin θ. Each dot on a fixed curve correspond to the intersection of the real streamline and the plane φ = const as it is being swept from φ = 0 to φ = π/2. In other words, each line is the intersection of a family of real streamlines all having the same θ0, with a fixed plane φ = const. The reason of doing this is because direct projections of streamlines on a plane φ = const show false intersections of the streamlines on the obtained plots. Figure 1 shows the projected streamlines for a fixed µ and different values of ν. Figure 2 shows also projected streamlines of the flow for a disc's radius of one half of the original cloud's radius (r d = 1/2µ). The bottom panels in Figure 2 have "triangles" drawn on them. These are zones out-of-range in the model due to real intersections of the accretion streamlines. Note that this behaviour appears because the model is ballistic. We have explored numerically in which cases these intersections occur and it turns out that for r d = 1/2µ, the streamlines show intersections for ν 1. From values of ν 50 intersections appear but the out-of-range region is kept fixed. Figures 3 and 4 show density isocontours and density profiles for a fixed height measured from the equatorial plane respectively. The out-of-range region appears only in the bottom panels of Figure 3, although on the left-bottom panel it appears as a tiny region. Figure 5 shows plots of the streamlines for the case ν = 0. The streamlines on the plot have different values of the initial cloud's radius. It is interesting to note that the Note that in all these cases, there is no intersection of the flow streamlines. This appears to indicate that as long as the cloud reaches a sufficiently large radius r0, with respect to the accretion disc radius r d , the intersections disappear. Indeed, we have tested different accretion models varying r d /r0 from ∼ 0.0 to 1.0 and the previous statement turns out to be correct. As an example, a value of r d /r0 ≈ 0.08 represents the threshold of intersections of streamlines for values of ν = 10.0.

STREAMLINES AND DENSITY PROFILES
From Figures 6 and 7 it is clear that the particle number density does not accumulate that much in the origin as ν increases.

ASTROPHYSICAL APPLICATIONS
From the point of view of Spectral Energy Distribution (SED) modelling of protostellar cores, one of the main ingredients required for a radiative transfer calculation is a specification of the density profile in the core. Many authors use Ulrich's model to achieve this (see for e.g. Kenyon et al. 1993;Butner et al. 1991;Adams & Shu 1986;  Whitney et al. 2003). A radiation transfer simulation requires a lot of information: luminosity of the star, composition, abundance and optical properties of the material in the core, in terms of its location with respect to the star. Thus, a full analysis of the SED produced by the use of the density given in this article is beyond its scope, but must be consider in the future.
In order to give a qualitative idea of the difference in the results, using either Ulrich's density or the density calculated in Section 3 of this article, we note that less density means more luminosity arriving from the object. The main argument is that the optical depth, which is proportional to the surface density calculated including all the material along a line of sight, is an element that give us an idea of the SED produced by the core.
For a disc viewed pole-on, in Figure 9 we present the axisymmetric surface density Σ as a function of the radial coordinate, for Ulrich's model (Σu) and from this paper (Σa for vr 0 = 0 and Σ b for vr 0 = 5 × 10 4 cm/s, according to the results of Hennebelle et al. (2004)) . This comparison is done for 3 representative cases of r0 and Ω = h0/r0, such that ru = 20 AU is the same for all. We take three different sets of boundary conditions: (a) r0 = 300 AU, Ω = 10 −11 s −1 , (b) r0 = 3000 AU, Ω = 10 −13 s −1 , and (c) r0 = 30000 AU, Ω = 10 −15 s −1 . Note that in all three cases, h0 is kept fixed since for all models the radius ru has been chosen to have a fixed value. The other common parameters on these models are given by the mass of the central star M = 1 M⊙ and the  Figure 4. Both panels show density profiles for a fixed azimuthal angle, i.e. φ = const and a fixed height z = 0, 0.1, 0.2, 0.3, 0.5 measured from the equatorial plane. The figures were drawn using values of ν = 0.64, 1.62 respectively for a fixed value of the radius of the disc r d = 1/2µ. Note that for a height z = 0 the particle number density diverges at r = 0, ±0.5. Distances are measured in units of the cloud's radius r = 1/µ. mass accretion rateṀ = 10 −6 M⊙ yr −1 . The most typical astronomical case is (b) in accordance with the observations of Benson & Myers (1989); Jijina et al. (1999).
In Figure 9, one can see that Ulrich's surface density Σu is smaller than the one calculated using our model, in the inner regions of the core (few times r d ). The difference between both curves increases as r0 decreases. Thus, the flux of radiation in the typical frequencies produced in this region are larger in the former than in the latter. Of course, a full analysis of this problem awaits an appropriate simulation, but since case (b) is quite common in accordance to observations, we can see that the modification to Ulrich's model is quite important. From the figures it follows that even if vr 0 has a null value, our model shows significant deviations as compared to the ones predicted by an Ulrich approach.
The last part of this section is intended to show the relevance of having flexible boundary conditions for the accretion flow studied in this article. The col- lapse of non-rotating clouds was studied by Larson (1969); Penston (1969);Hunter (1977). As an example, Whitworth & Ward-Thompson (2001) built a model, in order to fit lifetimes and accretion rates for class 0 protostars (this stage ends until half the final mass of the central object is assembled). In this model, the density differs from the singular isothermal sphere (SIS), cf. Shu (1977), because in the former, the density has an inner flat profile. Moreover, observations of protostellar cores (Ward-Thompson et al. 1999) suggest that initial conditions for protostellar collapse depart from SIS.
At this point, it is clear that the boundary conditions in a collapsing core depend on the environment. Dynamical collapses are present in regions like ρ Ophiuchi and Perseus clusters, where the collapse and subsequent star formation is triggered by external agents ( Figure 6 has. All plots have density as a function of the equatorial radial coordinate R for different fixed heights z = 0, 0.15, 0.3, 0.5. In all figures, the plot at z = 0 has density divergences at R = 0, ±1. shock wave. In regions like Taurus it is possible that cores collapse by themselves without an external perturbation and so, the mass accretion rate is closer to the value calculated for the SIS (Shu 1977).
Triggered star formation by an external mechanism is throughfully studied by Hennebelle et al. (2004). This work presents numerical simulations of a rotating core, with high spatial resolution, for the study of the formation of a disc. The relevance of the perturbation is measured with their pa- rameter φ, which is the ratio between the time-scale in which the internal pressure doubles, to the sound-crossing time of the core. A high value of φ, represents a slow compression (spontaneous collapse); the opposite is found with a small value, i.e. fast compression (strongly externally triggered).
The difference between both cases is clearly seen in the size-scale of the disc. In the first one, the disc evolves to a two-armed spiral pattern without fragmentation. Ring formation is the outcome usually found in the second case. Most of the material reaching the centre has too much angular momentum and so, it moves outwards and establishes a dense ring (Nagel 2007). This configuration promotes fragmentation. Despite the fact that in this case Hennebelle et al. (2004) found no central protostar, there must be situations where certainly it is present. In conclusion, either spontaneous or externally triggered collapse can be studied with the ballistic solution presented in this paper.
Both, Hennebelle et al. (2004) and our work start from a core in solid body rotation. The former evolves to a core of differential rotation, where a rapid compression implies higher velocities and densities in the outer parts of the disc. This fast compression intuitively means that a larger vr 0 is present. Indeed, if in our model we set a large vr 0 , then higher velocities and densities are obtained in the same region. For, substituting θ = π/2 and θ0 in terms of r from equation (9) in equation (13) it is found that for a given r, the velocity vr and particle number density n increase with an increasing ν. This conclusion, expected intuitively, reinforces the idea that the solution presented in this paper can be taken for studies of disc formation with an analytical approach.
The robustness of this kind of solutions is that, depending on the particular astrophysical situation, we can choose a ballistic solution with appropriate boundary conditions. For example, Throop & Bally (2008) studied numerically the accretion of material into a star-disc system moving through the gas in a cluster. An appropriate solution, should help to study small scales, where the star's gravity represents the main contribution. For the case of star formation in molecular clouds, Bate & Bonnell (2005) describe the accumulation of material as "competitive accretion". The matter comes from different places with different boundary conditions. In this case, ballistic solutions are able to give glimpses of the physics behind the formation of star-disc systems.

CONCLUSION
We have constructed an analytic accretion flow that expands Ulrich's (1976) model in the sense that the radius of the rigid body rotating cloud is finite and an allowance for a radial accretion velocity at the cloud's border is assumed. When the radius of the cloud tends to infinity and the radial velocity of the input flow at infinity goes to zero, the model converges to that of Ulrich.
We have also shown that for real astronomical systems, where typical sizes of the clouds are ∼ 30000 AU and for which the initial inward velocity is null, the differences with Ulrich's model is quite noticeable. The problem with this hypothesis is that to achieve a fixed accretion rate, then the density at the cloud's border needs to be infinite. The difference between both models becomes even greater if one assumes an initial inward velocity different from zero. Since this assumption makes the density not to diverge at the clouds border, it is the more reasonable model to take for a real astronomical system.
We are developing a relativistic model with all these properties in order to generalise the model constructed by Huerta & Mendoza (2007). These results will be the subject of a subsequent article.