{'fraction_non_alphanumeric': 0.0410291887289466, 'fraction_numerical': 0.019982438868637417, 'mean_word_length': 4.522186306200411, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Approach to Collapse of Molecular Clouds

The dense molecular cloud cores that form stars, like other self-gravitating objects, undergo bulk oscillations. Just at the point of gravitational instability, their fundamental oscillation mode has zero frequency. We study, using perturbation theory, the evolution of a spherical cloud that possesses such a frozen mode. We find that the cloud undergoes a prolonged epoch of subsonic, accelerating contraction. This slow contraction occurs whether the cloud is initially inflated or compressed by the oscillation. The subsonic motion described here could underlie the spectral infall signature observed in many starless dense cores.


Introduction
The central process in star formation is the gravitational collapse of a dense, molecular cloud core. Such objects are well studied empirically, through a variety of techniques. Less than half of observed dense cores already contain embedded stars, with the fraction increasing at higher mean column density and volume density (Beichman et al. 1986;Jessop & Ward-Thompson 2000). The remaining, so-called starless, cores exhibit a range of density contrasts and temperature gradients (e.g., Tafalla et al. 2004;Crapsi et al. 2007). Those of higher density do not differ in their gross properties from cores with stars . Observations, then, have not yet provided a clearcut answer to the salient question: How does the collapse process actually begin?
One important clue is the low level of internal motion observed within dense cores. This fact was established early, from the relatively narrow linewidths of optically thin molecular tracers such as NH 3 (Myers & Benson 1983). The linewidths reveal gas motion that is largely subsonic, in contrast to the turbulent or wave-like character of the cores' external environments (Goodman et al. 1998). Other studies bolster the view from spectroscopy. Alves et al. (2001) used near-infrared extinction mapping to reconstruct the density structure of the starless core B68. After azimuthal averaging, the density profile is close to a theoretical one for an isothermal sphere in balance between self-gravity and gas pressure (Ebert 1955;Bonnor 1956). Thus, the object appears to be in dynamical equilibrium, as expected when internal motion is subsonic. Similar results have since been found for a number of other starless cores (Tafalla et al. 2004;Kandori et al 2005). The inferred density profiles are consistent with earlier ones reconstructed through submillimeter emission maps (Ward-Thompson et al. 1999).
Real dense cores are not perfect spheres, in part because they derive additional support from the anisotropic force associated with the interstellar magnetic field (Crutcher 1999). In any case, the fact that clouds begin collapse from a near-equilibrium state is significant theoretically. What is the nature of this equilibrium? Clearly, the initial state cannot be dynamically stable with respect to small perturbations, as no collapse would ensue. 1 On the other hand, it is also unlikely to be dynamically unstable; it is difficult to imagine the prehistory of an object that could have veered away sharply from equilibrium at any time. The most natural initial state just prior to collapse for a star-forming dense core is one of marginal stability. While the object is in force balance, its fundamental mode of oscillation has zero frequency. Working again within the idealized framework of isothermal spheres, our paper addresses two specific questions. How does a marginally stable, isothermal cloud evolve with time? Can we account in a simple and compelling way for the tendency of such an object to collapse?
Collapse calculations have a long history, and we are by no means the first to recognize the particular importance of the marginally stable isothermal sphere as an initial state. Hunter (1977), Foster & Chevalier (1993), Ogino et al. (1999), and Aikawa et al (2005) all followed the collapse of such a configuration. However, these authors were concerned with developments at a relatively advanced stage. The subtle question of how the cloud first evolves away from equilibrium was effectively bypassed through assuming that the cloud was initially slightly overdense with respect to equilibrium, thereby guaranteeing eventual collapse. Our calculation focuses entirely on the transition issue. We show, using the tools of perturbation theory, that marginally stable clouds enter a protracted phase of slow, but accelerating inward contraction. We do not follow the evolution deep into collapse, a task that could be pursued using more standard methods.
Our study of this early contraction phase holds more than purely theoretical interest. A large fraction of starless cores show convincing signs of inward motion (Williams et al. 1999;Gregersen & Evans 2000;Lee et al 2001;Schnee et al. 2007). The spectroscopic signature of infall is an asymmetric emission line, often self-absorbed, that is skewed toward the blue. Inferred speeds are well below thermal, but greater than those associated with gravitational settling me-diated by ambipolar diffusion of the magnetic field (as calculated, e.g., by Ciolek & Mouschovias 1994). Our accelerating contraction is an attractive candidate to explain this intriguing, and perhaps pivotal, set of observations.
In Section 2 below, we present our method of solution. After nondimensionalization of the dynamical equations, we introduce a perturbation expansion that allows us to separate out the first-order, oscillatory motion from the second-order displacement that progressively evolves with time. Sections 3 and 4 are devoted to both analytical derivations and numerical results. In the first, we display the cloud's normal modes of oscillation, including the critical one of zero frequency. In the second, we trace the bulk contraction of a cloud subject to the zero-frequency mode. Finally, Section 5 discusses the possible connection to the infall signature of starless cores, and indicates fruitful extensions of this work.

Physical Assumptions
Our task of following early cloud evolution is facilitated by adopting a simplified physical picture. These simplifications are by now traditional, but it is important to revisit them as new observations arise and theoretical understanding grows. Thus, our employment of an isothermal equation of state, with a gas temperature that is constant in both space and time, is technically inconsistent with recent observational studies of starless cores. The inferred temperature of L1544, a starless core in Taurus that exhibits infall, decreases from 12 K at the outer edge to 6 K near the center (Crapsi et al. 2007). The central temperature is depressed by the partial shielding of dust grains from external starlight; these grains are thermally coupled to the gas through collisions. However, the region that is effectively shielded comprises little of the total mass, and collapse calculations that account for the finite temperature gradient show it to have a relatively minor effect (Keto & Field 2005).
Another traditional assumption, and one ostensibly even more radical, is that our model cloud is spherical. It has long been known that the projected molecular-line emission contours of starless cores, like dense cores in general, are more accurately elliptical, with mean axis ratios of about 2:1 ). The deprojected, or intrinsic, distribution of shapes must be inferred statistically from observations of a large sample. Under the condition that the underlying structure be axisymmetric, a prolate configuration best matches the data (Myers et al. 1991;Ryden 1996). Relaxing this restriction, Goodwin et al. (2002) found that a triaxial configuration is preferred (see also Jones et al. 2001). This result is puzzling. As we stated earlier, a departure from spherical symmetry could partially be explained by additional magnetic support. However, Galli (2005) has demonstrated that there are no triaxial magnetostatic structures, at least for scale-free configurations. It may be, as Galli suggests, that the observed shapes reflect distortion created by internal oscillations, a phenomenon we ourselves shall invoke shortly. In any event, our adoption of a spher-ical geometry, along with our neglect of magnetic support, can only be viewed as a convenient, first approximation. The essential finding of accelerating contraction from marginal stability should continue to hold within a more complete description of the equilibrium, magnetostatic state.
Within the last decade, a number of researchers have questioned the very existence of equilibrium structures. They have made the assertion, based on numerical simulations, that dense cores are transient, fully dynamic entities. Such studies, recently reviewed by Ballesteros-Paredes et al. (2007), treat an interior portion of the larger, parent cloud as a computational box filled with gas, with or without an embedded magnetic field. Turbulence is simulated by stirring the gas. This energy input is eventually offset by the unavoidable numerical dissipation present in the calculation. Once such a steady state is attained, external stirring is turned off and self-gravity turned on. Sufficiently overdense structures collapse on themselves, in apparent imitation of the star formation process. 2 It is intriguing that these structures resemble, in their masses, sizes, and geometric aspect ratios, observed dense cores. Even triaxial structures have been obtained (e.g., Offner et al. 2008). However, the key observational distinction between dense cores and their surroundings is their subsonic (or, more properly, sub-Alfvénic) internal velocity. This fact implies, as we have stressed, that the entities, unlike their numerical surrogates, are in force balance. Moreover, the large observed fraction of starless cores indicates a correspondingly protracted epoch in a typical dense core's lifetime prior to collapse Jessop & Ward-Thompson 2000;Kirk et al. 2005). If these structures indeed condense from a turbulent medium, their initial growth may be retarded by a relatively strong ambient magnetic field (Gálvan-Madrid et al. 2007).

Dynamical Equations and their Nondimensionalization
We consider a spherical cloud if fixed mass M and isothermal sound speed a. The equation of state relating the internal pressure P and density ρ is The outer edge of the cloud is set by the condition that P falls to some external value P e , also constant. In spherical symmetry, any matter outside the boundary exerts no gravitational force on the cloud. Hence we do not need to consider this medium in detail, beyond the fact that it exerts a fixed pressure. In particularly, we need not assume (unrealistically) that the external gas has a relatively high temperature and low density, as is frequently done in collapse calculations (e.g., Foster & Chevalier 1993). We shall see that the cloud boundary, as just defined, shrinks during the evolution. This finding suggests that the distinction between "cloud" and "external medium" is somewhat artificial, and that inward motion is spatially extended. Indeed, the observed infall signature of starless cores is striking for its broad extent (Myers et al. 2000).
In Eulerian coordinates (r, t), the equation of mass continuity is where u(r, t) is the fluid velocity. This velocity satisfies momentum conservation: Here, m is the mass enclosed within any radius: As the total cloud mass is fixed, we shall find it more convenient to switch to Lagrangian coordinates (m, t). The mass continuity equation (2) becomes Here, we have replaced the velocity u by the time-derivative of the (now dependent) variable r: and have again suppressed all subscripts in the partial derivatives. The momentum equation (3) transforms to ∂ 2 r ∂t 2 = −4 π r 2 a 2 ∂ρ ∂m − G m r 2 .
Finally, equation (4) for m is replaced by Equations (5), (7), and (8) constitute our basic dynamical equations. We may further cast them into nondimensional form using the three constants G, a, and P e . We define nondimensional versions of m, r, ρ, and t:m We further replaceρ by a new nondimensional variableφ, defined through After dropping the tilde, our dynamical equations become We will also be using equation (6) relating u and r. If we setũ ≡ u/a, then this equation remains the same nondimensionally.

Perturbation Expansion
Equations (14)- (16) are to be solved subject to the inner boundary condition r(0, t) = 0 and the outer one φ(M, t) = 0. The latter is just the requirement, expressed in nondimensional language, that the gas pressure at the edge equal the fixed value P e . We must also specify the cloud's initial configuration. This state is itself slightly perturbed from true equilibrium, under the influence of the oscillation mode. We are thus motivated to introduce a perturbation expansion of the dependent variables about the equilibrium state. We write: Here, ǫ is a small, free parameter. The physical significance of ǫ is that it will set the amplitude of the oscillation mode present at the outset (see §4 below). Quantities with the 0 subscript refer to the time-independent, equilibrium state.
We are now in a position to separate out the hierarchy of cloud motions. We substitute our series expansions into equations (14) -(16) and equate like powers of ǫ. At the lowest (zeroth) order in ǫ, both sides in the equation of mass continuity vanish. However, equations (15) and (16) yield, to the same order, the structural equations of the underlying, equilibrium state: When we equate terms proportional to ǫ, all three equations give a non-trivial result: These zeroth-and first-order equations are sufficient to describe small oscillations about a general equilibrium state. However, we are interested in that unique equilibrium for which the fundamental oscillation frequency vanishes. The perturbed cloud will not be perfectly static, but will evolve slowly, through second-order effects. Thus, we also need to equate terms in the dynamical equations proportional to ǫ 2 . We find The two acceleration terms A 1 and A 2 in equation (25) are

Equilibria
We consider first equations (19) and (20), that describe the equilibrium state. These equations are subject to the boundary conditions r 0 (0) = 0 and φ 0 (M ) = 0, to be applied at the center and outer edge, respectively. More accurately, the second condition defines the boundary of the equilibrium configuration. As is well known, there exists a one-parameter family of such structures, each characterized by its own center-to-edge density contrast. We denote this contrast as ρ c , since it is identical to the nondimensional central density in our formulation.
Because of the inner boundary condition, the righthand sides of both equations are singular at the origin. To start our numerical integration, we therefore expand both r 0 (m) and φ 0 (m) in the appropriate power series: The first two coefficients are and higher ones are determined recursively.
Based on our numerical integration, the dashed curves in Figure 1 are density profiles for states with ρ c = 5 and ρ c = 30. Note that we have plotted ρ 0 ≡ exp(−φ • ) as a function of r 0 , in the conventional manner. The middle, solid curve corresponds to ρ c = 14.04, for which M = 4.19. This is the famous Bonnor-Ebert state. As we shall verify shortly, it is the starting configuration of interest for the present problem. The other two curves represent outer limits to the empirical fitting of density profiles for most starless dense cores (Kandori et al 2005).

Normal Modes
Equations (21)-(23) describe, to linear order, internal motion of the equilibrium cloud. All such motion can be decomposed into a series of normal modes. The fundamental, or breathing, mode is generally expected to have the largest amplitude, but higher harmonics may also be present (as well as non-radial oscillations; see Keto et al. 2006).
To obtain the full set of spherical normal modes, we first solve equation (23) for the density perturbation: We substitute this expression into the momentum equation (22) and obtain a partial differential equation for r 1 (m, t): The coefficients in this equation are all functions of the chosen equilibrium state.
Before proceeding with the solution of equation (34), we note that the additional first-order equation (21) has not been used in the derivation. In fact, this relation is redundant, and already contained implicitly in equation (33). To see this, take the time derivative of the latter: = r 2 0 e −φ• ∂ 2 r 1 ∂t ∂m + 2 r 0 e −φ• dr 0 dm which is just equation (21). Here we have used the time-independence of all equilibrium quantities, and have also employed equation (20).
Returning to equation (34), we seek oscillatory solutions for r 1 (m, t). Thus, we set Since equation (34) is linear in r 1 , the coefficient A is an arbitrary constant. After substitution, we find an ordinary differential equation for ξ(m): The boundary conditions are now ξ 1 (0) = 0 and φ 1 (M ) = 0. Using equation (33), the latter may be transformed into a condition on ξ 1 (m): We note that equation (37) has a regular singular point at the origin. Once more, we start the integration through a power series development: The constant c 0 may be chosen arbitrarily. The next coefficient is and successive ones may be similarly found.
For any central density ρ c of the equilibrium cloud, there exists a sequence of ω 2 -values such that ξ 1 (m) meets the two boundary conditions. When ρ c = 14.04, the lowest value of ω 2 is zero. The corresponding ξ 1 (m) is the oscillation mode of interest. If we denote the fundamental, first, and second harmonics as ξ 10 , ξ 11 , and ξ 12 , then the corresponding ω 2 -values are 0, 8.37, and 24.3. Figure 2 plots these three normal modes, again using r 0 as the independent variable. In all cases, we have set the coefficient c 0 in equation (39) to unity.
Note finally that we may use the equilibrium relations, equations (19) and (20), to rewrite equation (37) as This equation is of the Sturm-Liouville type. As expected physically, all of its eigenvalues ω 2 are real. Distinct eigenmodes are orthogonal: Since only the zero-frequency oscillation will be of interest in the following discussion, we will soon revert to the simpler notation ξ 1 (m) for the fundamental mode, and assume the standard normalization c 0 = 1 for this function. After also setting the coefficient A in equation (36) to unity, we will thus be making the identification r 1 (m, t) = ξ 1 (m).

Fundamental Equation
Once we have selected the oscillatory mode of the equilibrium state, the second-order equations (24)-(28) describe additional motion. This motion is, of course, critical when the underlying oscillation has zero frequency, but we shall first keep the discussion more general. As in our derivation of the normal modes, it is convenient to solve for the density perturbation (now φ 2 ) from the mass-radius relation, equation (26): where we have used equation (33) for φ 1 . We substitute both this latter relation and equation (43) for φ 2 into the righthand side of the momentum equation (25). We thus derive our fundamental partial differential equation for the displacement r 2 (m, t): where We again remark that it has been unnecessary to invoke the continuity equation (24). The explanation, as before, is that this equation yields no new information. To see this, take the time derivative of equation (43) and multiply through by e φ• : After using equation (33) for φ 1 and equation (35c) for ∂φ 1 /∂t, we also have = 2 r 0 ∂r 1 ∂m ∂r 1 ∂t + r 4 0 e −φ• ∂r 1 ∂m ∂ 2 r 1 ∂t ∂m + 4 r 1 ∂r 0 ∂m ∂r 1 ∂t + 2 r 0 r 1 ∂ 2 r 1 ∂t ∂m , where we have again utilized equation (20). Combining equations (46) and (47b) then gives which is equivalent to equation (24).
Except for the inhomogeneous term F , equation (44) for r 2 (m, t) is identical to equation (34) for r 1 (m, t). It is precisely because of this extra term that the second-order motion cannot be described as a normal mode of oscillation. Once we single out the zero-frequency, fundamental mode for the first-order motion, F can be written purely as a function of m: The boundary conditions on r 2 (m, t) are the usual ones: r 2 (0, t) = 0 and φ 2 (M, t) = 0. From equation (43), the outer boundary condition is more usefully recast as

Method of Solution
Solving the partial differential equation (44) requires that we first specify the full initial state of the cloud. Assuming the fundamental mode has zero frequency, i.e., that r 1 = ξ 1 , we still need to set the functional form of r 2 (m, 0). As our fiducial initial state, we demand that the cloud's density profile be just that resulting from the normal mode acting on the equilibrium state. That is, we set φ 2 (m, 0) = 0. With this condition, equation (17) then gives the physical interpretation of the parameter ǫ as the nondimensional amplitude of the initial φ-perturbation.
The vanishing of φ 2 (m, 0) does not mean that r 2 (m, 0) = 0. Instead, the initial r 2 -profile follows by setting φ 2 = 0 in equation (43) and specializing to r 1 = ξ 1 : We solve this ordinary differential equation for r 2 (m, 0). After again noting the regular singular point at the origin, we begin the numerical integration through a power-law expansion: The expansion coefficients are readily found: The top curve of Figure 3 displays the calculated r 2 (m, 0), again as a function of r 0 .
We also need to specify the initial velocity ∂r 2 /∂t (m, 0). An interesting, and physically relevant situation is when the cloud is perfectly static, i.e., when More realistically, the cloud has slowly evolved from some earlier state, perhaps through the effect of ambipolar diffusion. In the absence of a more complete evolutionary picture, we shall adopt the simplest, zero-velocity, initial condition.
Equation (44) may be solved numerically through the method of characteristics. For this purpose, it is convenient to adopt r 0 as the independent, spatial variable. Using the connection between r 0 and m in equation (20), the fundamental equation (44) becomes where F is now written as The reason for changing Lagrangian variables from m to r 0 is that the characteristics of equation (55) are simply dr 0 dt = ± 1 .
In light of our nondimensionalization, we see that small disturbances in the cloud travel at the isothermal sound speed, an intuitively appealing result. We solve equation (55) by propagating the partial derivatives ∂r 2 /∂r 0 and ∂r 2 /∂t along the (r 0 , t) grid, as detailed in the Appendix. The starting value of ∂r 2 /∂t is everywhere zero, according to equation (54). As we advance in time, we also incorporate both the inner boundary condition r 0 = 0 and the outer one of equation (50b). The latter is now more conveniently written as

Numerical Results
Integration of the partial differential equation (55) for r 2 (r 0 , t) shows that this function globally decreases for t > 0. That is, all mass shells contract, although the cloud was neither expanding nor contracting in its initial state. If this state were slightly compressed from equilibrium, the parameter ǫ would be negative. If slightly expanded, ǫ would be positive. In either case, the second-order effects governing subsequent evolution are proportional to ǫ 2 . Thus, the cloud inevitably contracts. This contraction begins slowly, but accelerates with time. To illustrate the effect, Figure 3 displays r 2 (r 0 , t) for t = 0, 1, and 2. We see how any initial, slight inflation (assuming ǫ > 0) is everywhere reversed by t = 1. By t = 2, mass shells roughly midway from the center are clearly contracting the most rapidly.
The time span shown in the figure exceeds that associated with free-fall collapse. The latter is conventionally taken to be where ρ c is the cloud's central density. (Here, we are temporarily reverting to dimensional variables.) Our nondimensional timet is given by equation (12), so that Here, ρ e ≡ P e /a 2 . For the Bonnor-Ebert initial state (ρ c /ρ e = 14.04), the numerical coefficient in equation (60) is 1.95. Thus, the last profile shown in Figure 3 corresponds to 2 × 1.95 = 3.90 free-fall times.
We are primarily interested, not in the displacement of mass shells, but in the contraction speed. From equation (18), the nondimensional velocity (i.e., its value relative to the sound speed a) is given by ǫ 2 ∂r 2 /∂t. The actual velocity of each mass shell thus depends on the value of ǫ, which we do not know a priori. Keto et al. (2006) have recently matched molecular-line profiles from the starless core B68 by assuming the object is undergoing a nonradial oscillation mode, with a dimensionless amplitude of 0.25. Using this figure as rough guide, we provisionally choose ǫ = +0.2 and display the resulting velocity profiles in Figure 4. Here, the times are identical to those in Figure 3. However, the independent spatial variable is now the full radius r, as given by equation (18).
We see that the velocity of all mass shells is negative, with a magnitude that increases with time. Thus, contraction is indeed accelerating. By the last time shown, the largest contraction speed within the cloud is about 0.2 times the sound speed. Interestingly, the cloud edge also has negative velocity, suggesting that contraction is spreading into the exterior region.
Closer to the cloud center, it appears from Figure 4 that the velocity at any time increases linearly with radius. This feature of the evolution was found by Foster & Chevalier (1993) in their numerical simulation. We may derive the result analytically by examining the fundamental equation (55) in this limit. Consider first the inhomogeneous term F . Now the series expansion for ξ 1 (m) in equation (39) tells us that ξ 1 = e 0 r 0 to lowest order, where e 0 is a constant. From equation (56), F has three terms that diverge as r 0 approaches zero. Their sum is which vanishes. However, the rest of equation (55) has coefficients that still do diverge. Close to the origin, the equation reduces to Here we have used the fact that, since r 2 = 0 at the origin for all times, we may neglect the second time derivative on the left side of equation (55). Since equation (62) itself holds for all times, we may take its time derivative and obtain the analogous relation for the velocity v ≡ ǫ 2 ∂r 2 /∂t: The non-divergent solution to this equation is that v is indeed proportional to r 0 . Figure 5 shows the evolution of the cloud's density profile ρ(r, t) over the same time interval as in Figures 3 and 4. Since ǫ was assumed to be positive, the cloud begins in a slightly inflated state, with a central density of only 5.15. However, subsequent contraction drives up the density, which reaches a central value of 15.5 by t = 2. Since all internal velocities are still subsonic at this time, the density profile is consistent with a cloud that, to first order, is in hydrostatic balance. Force balance will, of course, be strongly violated in the future, as the cloud enters a state of true collapse.
The perturbative nature of our calculation limits us to describing the initial transition phase. Our calculation is only valid to the point where second-order terms in the expansions become comparable to their first-order counterparts. For ǫ = +0.2, the maximum absolute value of ǫ 2 r 2 (m, t) in equation (18) becomes equal to the maximum value of ǫ r 1 (m, t) at t = 1.9. Thus, the calculation is still marginally valid for the t = 2 profiles shown in Figures 3, 4, and 5, but not beyond those.

Alternative Initial Condition
It is important to verify that the evolutionary results shown thus far are not sensitive to the detailed initial state chosen. Recall that our fiducial state represented a pure first-order perturbation of the density, in that we set φ 2 = 0 for all mass shells. This choice was convenient, but arbitrary. We could have chosen any initial density profile consistent with the boundary condition φ 2 (M, t) = 0.
With this φ 2 -profile, the outer boundary condition is satisfied, but the central φ 2 is now unity. Since ρ = exp (−φ), the initial cloud is less dense than before. If we insert equation (61) into equation (43), we may again integrate an ordinary differential equation for r 2 (m, 0). As shown in the left panel of Figure 6, the r 2 -profile is thus inflated relative to the previous initial state, with a maximum fractional difference of 43 percent. Nevertheless, subsequent contraction brings the cloud to a very similar configuration. The right panel of Figure 6 displays the new r 2 -profile at t = 2. It differs from the one obtained using the original starting state by at most 8 percent. The profile of velocity, ǫ 2 ∂r 2 /∂t, is also very close to the previous one, with the maximum speeds differing by less than 1 percent at t = 2.

Discussion
In the picture introduced here, star-forming dense cores undergo a prolonged phase of contraction before their ultimate collapse. The contraction is slow because the cloud's self-gravity is still nearly counteracted by the outward pressure gradient. The slight imbalance of these forces creates subsonic, inward motion that gradually accelerates. Previous researchers performing direct simulations have also documented slow motion prior to runaway collapse. However, the characteristics of this phase, and indeed whether it led to collapse or rebound, depended on the unavoidable artifices of a direct simulation: the imposed overdensity in the cloud, the numerical accuracy of the initial state, and the precise treatment of the central few zones (Hunter 1977;Boss & Black 1982;Foster & Chevalier 1993;Ogino et al. 1999). We have demonstrated physically how an initially static, marginally stable cloud inevitably contracts. This contraction is essentially the inward phase of a slow oscillation that smoothly leads to free-fall collapse.
Motivated by the observed infall signature of many starless dense cores, others have offered different pictures. Keto & Field (2005) hypothesized that some cores are born in a gravitationally unstable state. They followed the collapse of such an object numerically, accounting for internal cooling by molecular lines and thermal dust emission. This work extends that of Zhou et al. (1990), who modeled the (starred) core B335 as undergoing collapse from a singular isothermal sphere. 3 As already noted, the early history of unstable objects is problematic; Keto & Field (2005) speculate that the cloud fragmented from a larger, turbulent velocity field. Myers & Lazarian (1998) attributed this localized fragmentation to enhanced dissipation via ion-neutral friction, leading to a pressure-driven, inward flow.
Our general concern about using unstable or actively collapsing states to match observations is their brevity. To illustrate the point more quantitatively, suppose the infall signatures reflect first-order motion, i.e., that the cloud's fundamental eigenfrequency is nonzero. The nondimensional perturbation r 1 (m, t) would then be expressed as The function ξ 1 (m) obeys equation (37), but with ω 2 replaced by −|ω • | 2 . If we again normalize the coefficient c • in equation (39) to unity, then ξ 1 (m) resembles the curve ξ 10 in Figure 2. Notice the overall minus sign in equation (65) that ensures contraction prior to free-fall collapse. Over what time does this transition occur?
The physical velocity, in units of the sound speed, is ǫ ∂r 1 /∂t = ǫ |ω • | r 1 . If ǫ is still about 0.2, and if we are to match velocities of 0.2 times the sound speed at the present epoch (t = 0), then the growth rate |ω • | −1 must be about unity. The dimensional time t for the velocity to increase by a factor e is, from equations (9) and (12), The numerical solution to the altered equation (37) tells that the nondimensional cloud massM corresponding to |ω • | = 1 is 4.02. The center-to-edge density contrast of this object is 30, a figure that is marginally consistent with observations (recall Fig. 1). Returning to L1544, its gas temperature is 10 K (Benson & Myers 1989), while the most recent mass estimate is 2 M ⊙ (Ohashi et al 1999). Equation (66) then gives an e-folding time of 2 × 10 5 yr. The statistical prevalence of starless dense cores makes it unlikely that they are evolving over such a brief, dynamical interval. Note finally that if the initial perturbation amplitude ǫ were smaller than 0.2, then |ω • | would be correspondingly larger, driving down the evolutionary time even more.
We emphasize again the general nature of our theoretical finding. A more realistic model for a starless core should certainly include the anisotropic supporting force from the interstellar magnetic field. Such magnetostatic configurations are themselves subject to global oscillations. The lowest eigenfrequency vanishes in the marginally stable state (Tomisaka et al. 1988). These states, analogues of the Bonnor-Ebert configuration used here, will also undergo slow contraction prior to collapse, even under the approximation of flux freezing. They are close to being magnetically supercritical, so that ambipolar diffusion may enhance the contraction process (Ciolek & Basu 2001). In any event, it will be interesting to track the evolution through the transonic phase into full collapse, both in our spherical model and its magnetized generalization.
Returning to the observational motivation of this study, it will also be instructive to calculate, within our spherical model, the predicted profiles for molecular emission lines of varying optical depth. For our representative case of ǫ = +0.2, we find a maximum infall speed of 0.2 times the sound speed after about 4 free-fall times. This speed is comparable to the 0.08 km s −1 inferred for the well-studied starless core L1544 through N 2 H + observations (Williams et al. 1999) We are encouraged by this finding, but stress the need for more comprehensive and detailed comparisons. Note especially that the amplitude ǫ is not readily observable with any precision; nor is the evolutionary time t. And yet contraction models of lower ǫ and larger t broadly mimic, in their velocity profiles, those with higher ǫ and shorter t. Hopefully, calculated line profiles will differ enough to resolve this ambiguity. If some of the profiles successfully match observations, we will not only have gained new insight into the mechanism of dense core contraction, but also a new measure for their longevity prior to collapse.
We are grateful for stimulating conversations with Phil Chang and Steve Shore during the inception of this project. We also thank the referee, Tom Hartquist, for comments that improved the original manuscript. S. S. was partially supported by NSF Grant AST-0639743.
We may also replace any column in M by the vector y. For example, consider the matrix M ′ given by There is no solution at all for the second derivatives unless det M ′ = 0. This condition gives us the further differential relation ∆ ∂r 2 ∂r 0 − ∆ ∂r 2 ∂t = g ∆r 0 , along the + characteristic, and − ∆ ∂r 2 ∂r 0 − ∆ ∂r 2 ∂t = g ∆r 0 , along the -characteristic. Figure 7 is a schematic spacetime diagram that illustrates the practical procedure. We set up uniform grids along the r 0 -and t-axes. At t = 0, we have ∂r 2 /∂t = 0, while the solution to equation (51) gives the initial values of ∂r 2 /∂r 0 for the fiducial initial state. (An analogous equation is solved for the alternative state; see §4.4.) Starting from any two adjacent points on the r 0 -axis, such as A and B, we use equations (A11) and (A12) to solve simultaneously for the two first partial derivatives at point D, where the two characteristics intersect. Similarly, we find the two derivatives at E starting from the pair B and C. The derivatives at D and E then yield those at point G, and so on.
We also need to propagate the information contained in the central and surface boundary conditions. Since r 2 (0, t) = 0, it is also true that ∂r 2 /∂t = 0 at points such as F in the figure. This condition, along with equation (A12) for the -characteristic joined to point D, allows us to solve for ∂r 2 /∂r 0 at F . Since both partial derivatives are now known at F and G, this information can then be used to find the derivatives at points further advanced in time. Near the cloud's outer edge, we establish the two partial derivatives at point L in the usual manner. Knowing this information at L, we find one relation between the two derivatives at the boundary point N by using equation (A11) for the + characteristic. The outer boundary condition, equation (58) then supplies a second relation between the derivatives. To find r 2 itself at N , we further use ∆r 2 = ∂r 2 ∂t ∆t + ∂r 2 ∂r 0 ∆r 0 .  Shown is the first-order displacement ξ 1 as a function of radius, for the primary (ξ 10 ), first harmonic (ξ 11 ), and second harmonic (ξ 12 ). All curves have been normalized to have the same initial slope. -Evolution of the displacement r 2 with two initial conditions. The left panel shows that, at t = 0, the r 2 -profile for the altered initial state (solid curve) differs substantially from that associated with the fiducial initial cloud that had a purely first-order density perturbation (dashed curve). Nevertheless, as the right panel shows, the profile is close to the previous one at t = 2. Fig. 7.-Schematic spacetime diagram illustrating the solution procedure for the fundamental partial differential equation (55). The straight, diagonal line segment from A to D lies along a + characteristic, while the line segment from B to D lies along a -characteristic.