text:
# Stable Models of Elliptical Galaxies

## 1 Introduction

Consider a large ensemble of mass points (stars) which interact only by the gravitational field which they create collectively. Such a collisionless, self-gravitating gas is used to model galaxies or globular clusters, cf.  and the references there. The time evolution of the density $`f=f(t,x,v)0`$ of the stars in phase space is governed by the following nonlinear system of partial differential equations:

$$_tf+v_xf_xU_vf=0,$$

$$U=4\pi \rho ,\underset{|x|\mathrm{}}{lim}U(t,x)=0,$$

$$\rho (t,x)=f(t,x,v)𝑑v.$$

Here $`t\mathrm{IR}`$ denotes time, $`x,v\mathrm{IR}^3`$ denote position and velocity respectively, $`\rho `$ is the spatial mass density of the stars, and $`U`$ the gravitational potential which the stars induce collectively. Collisional or relativistic effects are neglected. This system, which Sir J. H. Jeans introduced at the beginning of the last century for the above modeling purposes, is sometimes in the astrophysics literature referred to as the collisionless Boltzmann-Poisson system. Following the mathematics convention we call it the Vlasov-Poisson system.

In the present paper we are interested in the steady states of this system. Besides their very existence, questions of interest are the possible shapes of such steady states, in particular their symmetry type, and whether or not they are dynamically stable. We approach these questions from a mathematics point of view, in particular, only such information as can be extracted from the system stated above is to enter our arguments. Our “ideal” reader is an astrophysicist interested in these same questions. We find it deplorable that there is little communication between astrophysicists and mathematicians investigating these problems, and the present paper is an attempt to change this. Axially symmetric steady states of the the Vlasov-Poisson system were obtained in  as perturbations of spherically symmetric ones via the implicit function theorem. From an astrophysics point of view axisymmetric models of elliptical galaxies have been investigated in , but as to their stability no rigorous results are known to us. In the present paper we follow the variational approach developed in : We obtain axially symmetric steady states with vanishing and with non-vanishing velocity field as minimizers of appropriately defined energy-Casimir functionals. These steady states are no longer restricted to small perturbations of spherically symmetric ones, and, most importantly, they are nonlinearly stable. The question we would like to raise is: Do our mathematical constructions yield suitable models for real-world galaxies, or, if not, how can the mathematical approach be modified to do so?

If $`U_0=U_0(x)`$ is a stationary potential then any function $`f_0=\varphi (E)`$ of the particle energy

$$E=E(x,v)=\frac{1}{2}|v|^2+U_0(x)$$
(1.1)
solves the Vlasov equation with the potential $`U_0`$, since $`E`$ is a conserved quantity along particle trajectories. So if $`U_0`$ is the potential induced by $`f_0`$ then this is a steady state of the Vlasov-Poisson system—at this point it must be emphasized that to obtain a self-consistent model by this approach, i. e., to make sure that $`U_0`$ is indeed the potential induced by $`f_0`$, is mathematically non-trivial. Any steady state obtained by this ansatz is necessarily spherically symmetric, as follows from a result of Gidas, Ni & Nirenberg . In the present paper we make $`f_0`$ depend on an additional invariant of the particle trajectories. If $`U_0`$ is axially symmetric, i. e., invariant under all rotations about, say, the $`x_3`$-axis, then the corresponding component of angular momentum,

$$P=P(x,v)=x_1v_2x_2v_1$$
(1.2)
is such an additional invariant. To obtain steady states with a number density $`f_0`$ which depends on $`E`$ and $`P`$ we minimize the energy-Casimir functional

$$\frac{1}{2}|v|^2f(x,v)𝑑v𝑑x\frac{1}{2}\frac{\rho (x)\rho (y)}{|xy|}𝑑x𝑑y+Q(f(x,v),P(x,v))𝑑v𝑑x$$

subject to the constraint that $`f`$ is axially symmetric and has prescribed mass $`M>0`$:

$$f(x,v)𝑑v𝑑x=M.$$

The corresponding Euler-Lagrange equation then shows that any minimizer is a steady state of the desired form. In order to conclude that this steady state is nonlinearly stable it is essential that the energy-Casimir functional is conserved along solutions. This is always true for the sum of the first two terms which are the kinetic and the potential energy of the system. However, the third term is conserved only if $`P`$ is constant along particle orbits which is true in general only if the solution is axially symmetric. Thus, the stability result which we derive is restricted to axially symmetric perturbations of the steady state. The fact that the third term in our functional is not conserved along general solutions makes the terminology “Casimir functional” questionable, but we stick to it for lack of a better alternative. The situation is similar to the one of anisotropic, spherically symmetric steady states with $`f_0`$ depending on the particle energy and the modulus of the angular momentum, $`|x\times v|`$, where stability against spherically symmetric perturbations holds, cf. , but the general case is open. For isotropic steady states where $`f_0`$ depends only on the particle energy no such restriction is necessary, cf. .

We want to stress the fact that in this article every conclusion is rigorously derived in the mathematical sense. Our method stems from our long term research project in which both the existence and the stability of galaxy configurations has been mathematically investigated. Since the Vlasov-Poisson system can be viewed as an infinite dimensional Hamiltonian system, even the celebrated Antonov stability theorems for polytropes are, strictly speaking, only valid in the linearized sense. This is like the case of center equilibria in a $`2\times 2`$ Hamiltonian system: in general no dynamical stability assertion can be deduced from a linearized analysis, and some Liapunov type methods are needed. Recently, we were able to verify the validity of Antonov’s theorems in a dynamical, nonlinear sense . Our variational method also yields stability of certain Camm models .

The paper proceeds as follows: In the next section we formulate the basic set up and state our main results. In order to put these into perspective we include a short discussion of stability concepts from a mathematics point of view. In the third section we discuss various examples and some further properties of the steady states which we obtain. For interested readers we collect the details of our mathematical proofs in an appendix.

We conclude the introduction with some references to the mathematical literature. In addition to our work mentioned above the stability of spherically symmetric steady states of the Vlasov-Poisson system in the present stellar dynamics case is also investigated in . Global classical solutions to the initial value problem for the Vlasov-Poisson system were first established in , cf. also  and . For the plasma physics case where the sign in the Poisson equation is reversed, the stability problem is better understood, and we refer to . A rather general condition which guarantees finite mass and compact support of steady states, but not their stability, is established in .

## 2 Main results

Let us first fix some notation. For a measurable function $`f=f(x,v)0`$ we define its induced spatial density

$$\rho _f(x):=f(x,v)𝑑v,x\mathrm{IR}^3,$$

and its induced gravitational potential

$$U_f:=\rho _f\frac{1}{||}.$$

Next we define

$`E_{\mathrm{kin}}(f)`$ $`:=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle |v|^2f(x,v)𝑑v𝑑x},`$
$`E_{\mathrm{pot}}(f)`$ $`:=`$ $`{\displaystyle \frac{1}{8\pi }}{\displaystyle |U_f(x)|^2𝑑x}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\rho _f(x)\rho _f(y)}{|xy|}𝑑x𝑑y},`$
$`𝒞(f)`$ $`:=`$ $`{\displaystyle Q(f(x,v),P(x,v))𝑑v𝑑x}`$
where $`Q`$ is a given function satisfying certain assumptions specified below and $`P`$ is defined in (1.2). We will minimize the energy-Casimir functional

$$_C:=E_{\mathrm{kin}}+E_{\mathrm{pot}}+𝒞$$

under a mass constraint, i. e., over the set

$$_M:=\{fL^1(\mathrm{IR}^6)f0,f(x,v)𝑑v𝑑x=M,𝒞(f)+E_{\mathrm{kin}}(f)<\mathrm{}\},$$

where $`M>0`$ is prescribed. Eventually, we need to restrict ourselves to axially symmetric functions in this set:

$`_M^S:=\{f_M`$ $``$ $`f(Ax,Av)=f(x,v),x,v\mathrm{IR}^3,`$
$`A\text{any rotation about the }x_3\text{-axis}\}.`$
For the function $`Q`$ determining the Casimir functional we make the following

Assumptions on $`Q`$: $`Q,_fQ,_f^2QC([0,\mathrm{}[\times \mathrm{IR})`$, $`Q0`$, $`Q(0,)=0=_fQ(0,)`$, and with constants $`C,C^{},F,F^{}>0`$, and $`0<k,k^{}<3/2`$,

* $`Q(f,P)Cf^{1+1/k},fF`$, $`P\mathrm{IR}`$,
* $`_f^2Q(f,P)>0,f>0,P\mathrm{IR}`$,
* $`Q(f,)`$ is increasing on $`]\mathrm{},0[`$ and decreasing on $`]0,\mathrm{}[`$, $`f0`$,
* $`Q(f,0)C^{}f^{1+1/k^{}},fF^{}`$.

The first two assumptions are essential whereas the assumptions (Q3) and (Q4) can be modified in various ways, cf. the next section. Note that under the assumptions above $`_fQ(,P):[0,\mathrm{}[[0,\mathrm{}[`$ is one-to-one and onto for any $`P\mathrm{IR}`$. The steady states which we shall obtain will be of the form

$$f_0(x,v)=\varphi (E_0E,P),$$

where the particle energy $`E`$ is defined in terms of the induced potential $`U_0`$ as in (1.1),

$$\varphi (E,P):=\{\begin{array}{ccc}(_fQ(,P))^1(E)& ,& E0,P\mathrm{IR},\hfill \\ 0& ,& E<0,P\mathrm{IR},\hfill \end{array}$$
(2.1)
and $`E_00`$ is some cut-off energy. We should point out that if $`Q`$ does not depend on $`P`$—a case which is included in our assumptions—then no symmetry assumptions need to be made in the choice of the set $`_M`$ nor anywhere else. We now state our first main result:

###### Theorem 1

For every $`M>0`$ there is a minimizer $`f_0_M^S`$ of the energy-Casimir functional $`_C`$. Let $`U_0`$ denote the potential induced by $`f_0`$. Then $`f_0`$ is a function of the particle energy $`E`$ and angular momentum $`P`$ as defined in (1.1) and (1.2),

$$f_0(x,v)=\varphi (E_0E,P)$$

where $`\varphi `$ is defined in terms of $`Q`$ by (2.1). The parameter $`E_0`$ plays the role of a Lagrange multiplier and is given by

$$E_0:=\frac{1}{M}\left[E+_fQ(f_0,P)\right]f_0𝑑v𝑑x.$$

In particular, $`f_0`$ is an axially symmetric steady state of the Vlasov-Poisson system with total mass $`M`$.

The crucial part here is the existence assertion a detailed proof of which is given in the appendix. If one computes the Euler-Lagrange condition corresponding to our variational problem one obtains the assertions on the form of the minimizer. For the details of this tedious but straight forward argument we refer to  or .

Remark 1. If one is interested only in the existence of a steady state the usual approach is to prescribe some suitable function $`f_0(x,v)=\varphi (E_0E,P)`$, compute the corresponding spatial density $`\rho _0`$ which becomes, via $`E`$, a functional of the potential $`U_0`$, and one is then left with the problem of proving that the semilinear Poisson equation

$$U_0=4\pi \varphi (E_0\frac{1}{2}|v|^2+U_0,P(x,v))𝑑v$$
(2.2)
has a suitable solution. The non-trivial problem here is to decide which solutions have finite total mass and compact support in space. The first requirement for $`\varphi `$ from the point of view of our theorem is that $`\varphi `$ is a strictly increasing function of $`E_0E`$ so that a corresponding function $`Q`$ and the Casimir functional can be defined. Of course, for the existence question alone this monotonicity condition which in astrophysics textbooks appears as a stability condition is not necessary. A mathematical approach which makes no such monotonicity assumption can be found in . The point is that the minimizer which we obtain is not just any old steady state, but a nonlinearly stable one.

Some remarks on stability concepts. It may be useful to briefly review various concepts of stability before we state and discuss our stability theorem; for a more detailed such discussion we refer to . Consider a nonlinear dynamical system which for short we write as

$$\dot{x}=A(x),$$

$`A`$ being some nonlinear operator on the possibly infinite dimensional state space where solutions $`tx(t)`$ take their values. Let $`x_0`$ be a steady state, i. e., $`A(x_0)=0`$.

* The “appropriate” stability approach is certainly one where we can deal with our problem directly, i. e., without simplifying or linearizing it in any way. The steady state is called nonlinearly stable or Liapunov stable if one can guarantee that sufficiently small perturbations of the steady state launch solutions which are arbitrarily small perturbations of the steady state for future times. More formally: One must find two norms $`,||`$ on the state space such that for any (small) number $`ϵ>0`$ there always exists another (small) number $`\delta >0`$ such that every solution $`x(t)`$ with $`x(0)x_0<\delta `$ satisfies $`|x(t)x_0|<ϵ`$ for all $`t0`$. Some comments are in order:

+ It is desirable but not always possible to use the same norms for measuring the initial perturbation and the one at time $`t`$; only on finite dimensional spaces are all norms equivalent.
    + Sometimes one may not even obtain norms but more general tools to measure the deviations with.
    + From a physics point of view it is unsatisfactory if no rule is given on how to obtain the $`\delta `$ if $`ϵ`$ is prescribed. From this point of view a better result would be an estimate of the form

    $$|x(t)x_0|Cx(0)x_0,t0$$

    possibly with some explicit constant $`C`$. Unfortunately, mathematicians are not always (not often?) able to provide this.
    + A global existence result which guarantees that solutions exist for all time $`t0`$ at least for initial data close to the steady state is a necessary prerequisite and integral part for all the above.
* Since assessing nonlinear stability in the sense above is difficult, the problem is often approached via linearization. By Taylor expansion one can compute the linearization of $`A`$ at $`x_0`$; formally $`A_{\mathrm{lin}}=D_xA(x_0)`$. There are now at least two sub-concepts of linearized stability: One is to repeat the definition for Liapunov stability given above, but replace $`A`$ by its linearization $`A_{\mathrm{lin}}`$. After all, saying that $`A`$ is nonlinear does not forbid $`A`$ happening to be linear. This is what we want to call linearized stability. A somewhat different concept is spectral stability: One considers solutions of the linearized equation of the form $`e^{\lambda t}x`$ and tries to find out what the possible eigenvalues $`\lambda `$ are. If the real parts of all these eigenvalues are strictly negative then one calls the steady state spectrally stable. Again, some comments are in order:

+ In general there is no guarantee that the linearized problem actually has solutions of the form $`e^{\lambda t}x`$. If it does, then in order to draw any conclusions one must know that there are sufficiently many such solutions to get the general solution of the linearized problem by superposition (completeness of the eigenfunctions). And even if this goes well, there is simply no general result which will allow any conclusion on the behavior of the original, nonlinear system, unless the system is finite dimensional!
    + If one is dealing with a conservative system with some sort of Hamiltonian structure (such as the Vlasov-Poisson system), then the best to expect as far as spectral stability is concerned may be that all the eigenvalues are purely imaginary (due to inherent symmetries of the spectrum). In this case no stability follows for the nonlinear system, not even in the finite dimensional case.
    + Assume one can establish the existence of a growing mode, i. e., of a solution of the form $`e^{\lambda t}x`$ where $`\lambda `$ has positive real part. Then again there is no general result saying that $`x_0`$ is now also nonlinearly unstable, unless the system is finite dimensional. However, such a growing mode is a valuable first step toward proving a nonlinear instability result. How difficult it still may be to get to a mathematically rigorous nonlinear instability result from there is illustrated by

The upshot from all this is: One should be very careful to draw conclusions about stability for nonlinear, infinite dimensional systems from linearization, in particular, for Hamiltonian ones.

All this said let us now return to our minimizer $`f_0`$ from Theorem 1. To investigate its dynamical stability we note first that if $`f=M`$,

$$_C(f)_C(f_0)=d(f,f_0)\frac{1}{8\pi }U_fU_0_2^2,$$
(2.3)
where $`_2`$ denotes the usual norm on $`L^2(\mathrm{IR}^3)`$ and

$$d(f,f_0):=\left[Q(f,P)Q(f_0,P)+(EE_0)(ff_0)\right]𝑑v𝑑x;$$

we are allowed to subtract the term $`E_0(ff_0)`$ from the integrand since its integral vanishes. Since $`Q(,P)`$ is convex, the integrand can be estimated from below by

$$\left[_fQ(f_0,P)+EE_0\right](ff_0).$$

According to Theorem 1, this quantity is zero where $`f_0>0`$, while it equals $`(EE_0)f0`$ where $`f_0=0`$. Thus we see that

$$d(f,f_0)0,f_M.$$

Now we note that the left hand side of (2.3) is constant along axially symmetric solutions. Moreover, and this is the crucial point, the term $`U_fU_0_2^2`$ which seems to have the wrong sign in (2.3) vanishes along minimizing sequences by Theorem A1 stated in the appendix.

Before we state our stability result we point out that if we shift a minimizer in the $`x_3`$ direction we obtain another minimizer. Moreover, we do in general not know whether the minimizers are unique up to spatial shifts. As discussed in the next section a minimizer must a-posteriori have the following additional symmetry property

$$f_0(x_1,x_2,x_3^{}x_3,v_1,v_2,v_3)=f_0(x_1,x_2,x_3^{}+x_3,v_1,v_2,v_3),x,v\mathrm{IR}^3,$$

for some $`x_3^{}\mathrm{IR}`$. If we take without loss of generality $`x_3^{}=0`$ we refer to this as reflexion symmetry. This symmetry propagates along solutions of the time dependent problem, and, by restricting ourselves to data with this symmetry, we can at least ignore the non-uniqueness due to shifts:

###### Theorem 2

For every $`ϵ>0`$ there is a $`\delta >0`$ such that for any solution $`tf(t)`$ of the Vlasov-Poisson system with $`f(0)C_c^1(\mathrm{IR}^6)_M^S`$ and reflexion symmetric,

$$d(f(0),f_0)+\frac{1}{8\pi }U_{f(0)}U_0_2^2<\delta $$

implies that for all $`t0`$,

$$d(f(t),f_0)+\frac{1}{8\pi }U_{f(t)}U_{f_0}_2^2<ϵ,$$

provided $`f_0`$ is unique or at least isolated in the reflexion symmetric subset of $`_M^S`$.

Since after what was said above the proof is simple and instructive we include it here:

Proof. Assume the assertion were false. Then there exist $`ϵ>0`$, $`t_n>0`$, and $`f_n(0)C_c^1(\mathrm{IR}^6)_M^S`$ reflexion symmetric and such that for all $`n\mathrm{IN}`$,

$$d(f_n(0),f_0)+\frac{1}{8\pi }U_{f_n(0)}U_0_2^2<\frac{1}{n}$$
(2.4)
but

$$d(f_n(t_n),f_0)+\frac{1}{8\pi }U_{f_n(t_n)}U_0_2^2ϵ.$$
(2.5)
By (2.4) and (2.3),

$$\underset{n\mathrm{}}{lim}_C(f_n(0))=h_M^S,$$

where $`h_M^S`$ is the finite, negative minimum of $`_C`$ on $`_M^S`$, cf. Lemma A 3. Since $`_C`$ is conserved along classical, axially symmetric solutions as launched by $`f_n(0)`$, since mass is conserved, and since the assumed symmetry propagates,

$$\underset{n\mathrm{}}{lim}_C(f_n(t_n))=h_M^S\text{and}f_n(t_n)_M^S,n\mathrm{IN},$$

i. e., $`(f_n(t_n))`$ is a reflexion symmetric minimizing sequence for $`_C`$ in $`_M^S`$. Up to a subsequence we may therefore assume by Theorem A 1 that

$$U_{f_n(t_n)}U_{f_0}_2^20;$$
(2.6)
note that due to the reflexion symmetry of the minimizing sequence and any minimizer the shifts along the $`x_3`$-axis in the statement of Theorem A 1 must be zero in the present situation. It is at this point that the uniqueness or isolation of the minimizer $`f_0`$ is used. Since $`lim_n\mathrm{}_C(f_n(t_n))=h_M^S=_C(f_0)`$ we conclude by (2.6) and (2.3) that

$$d(f_n(t_n),f_0)0,n\mathrm{},$$

and we arrive at a contradiction to (2.5). $`\mathrm{}`$

Remark 2. If the minimizer is not isolated, i. e., if arbitrarily close to it with respect to the “distance” used in the theorem there are other minimizers which do not result from a simple translation in space, then we obtain a stability result of the following form:

Let $`_M_M^S`$ denote the set of all minimizers of $`_C`$ in $`_M^S`$. Then for every $`ϵ>0`$ there is a $`\delta >0`$ such that for any solution $`tf(t)`$ of the Vlasov-Poisson system with $`f(0)C_c^1(\mathrm{IR}^6)_M^S`$,

$$\underset{f_0_M}{inf}\left[d(f(0),f_0)+\frac{1}{8\pi }U_{f(0)}U_0_2^2\right]<\delta $$

implies that for all $`t0`$,

$$\underset{f_0_M}{inf}\left[d(f(t),f_0)+\frac{1}{8\pi }U_{f(t)}U_{f_0}_2^2\right]<ϵ.$$

Remark 3. Although we only showed that $`d(f,f_0)0`$ for $`f_M`$, one may think of this term as a weighted $`L^2`$-difference of $`f`$ and $`f_0`$, and by a Taylor expansion a stronger and more explicit estimate can be obtained if $`_f^2Q`$ is bounded away from zero.

Remark 4. The restriction $`f(0)_M`$ for the perturbed initial data would be acceptable from a physics point of view: A perturbation of a given galaxy, say by the gravitational pull of some outside object, results in a perturbed state which is just a rearrangement of the original state, in particular, its mass remains unchanged. But in general symmetry properties of the steady state will be destroyed by such a perturbation. However, we need that $`f(0)`$ is axially symmetric so that $`𝒞`$ is conserved along the resulting axially symmetric solution $`f(t)`$. So while this symmetry restriction for the perturbation seems necessary within the present mathematical framework it is quite undesirable from a physics point of view. If the steady state does not depend on $`P`$ and thus is isotropic, then no symmetry restrictions are necessary anywhere in our results, and the stability is with respect to quite general perturbations.

Remark 5. Our stability result is of the nonlinear Liapunov type discussed above in general, but it suffers from the defect that given $`ϵ`$ it does not say how $`\delta `$ needs to be chosen. The reason for this is that our proof is by contradiction and relies on a compactness result.

## 3 Examples and further properties

In this section we investigate some additional properties of the resulting steady states, keeping the discussion somewhat informal. First we present some examples for functions $`Q`$ which satisfy our assumptions or possible variations of them:

Examples. A simple class of examples for functions $`Q`$ which satisfy the assumptions (Q1)–(Q4) is given by

$$Q(f,P)=f^{1+1/k}g(P),f0,P\mathrm{IR},$$
(3.1)
where $`0<k<3/2`$ and $`g`$ is positive, continuous, bounded, bounded away from zero, increasing on $`]\mathrm{},0[`$, and decreasing on $`]0,\mathrm{}[`$. Note that in this case,

$$f_0(x,v)=C(E_0E)^kg(P)^k$$

on its support. Examining the proof of Lemma A 3 in the appendix, which is the only place where the assumptions (Q3) and (Q4) enter, one can see that these assumptions can for example be replaced by

* $`Q(\lambda f,P)\lambda ^{1+1/k^{}}Q(f,P),f0,0\lambda 1,P\mathrm{IR}`$, with $`1/2k^{}<3/2`$,
* $`Q(f,P)C^{}f^{1+1/k^{\prime \prime }},fF^{},|P|P_0`$, with constants $`C^{}>0,F^{}>0,P_0>0`$, and $`0<k^{\prime \prime }<3/2`$.

Examples which satisfy the assumptions (Q1), (Q2), (Q3’), (Q4’) but not the original ones are provided by (3.1) with $`1/2k<3/2`$ and $`g`$ which has all the properties stated above except the monotonicity.

Regularity and boundary condition. Since by Theorem 1 $`f_0`$ is a function of the quantities $`E`$ and $`P`$ which are constant along particle trajectories we are justified to call $`f_0`$ a steady state provided $`U_0`$ is sufficiently regular to allow for the definition of particle trajectories to begin with. Now by construction $`U_0=4\pi \rho _0`$ on $`\mathrm{IR}^3`$, at least in the sense of distributions, and from the very construction of $`f_0`$ certain integrability properties of $`\rho _0`$ follow, cf. the appendix. Applying the usual Sobolev space arguments one can eventually conclude that $`U_0`$ is twice continuously differentiable and $`lim_{|x|\mathrm{}}U_0(x)=0`$, cf. \[11, Thm. 3\] for the technical details.

Finite mass and compact support. By construction, the steady states have finite mass $`M>0`$ which we prescribe by our mass constraint. To continue we note that by Theorem 1 and a change of variables,

$`\rho _0(x)`$ $`=`$ $`{\displaystyle f_0(x,v)𝑑v}`$ (3.2)
$`=`$ $`2\pi {\displaystyle _{U_0(x)}^{E_0}}{\displaystyle _{\sqrt{2(EU_0(x))}}^{\sqrt{2(EU_0(x))}}}\varphi (E_0E,r(x)p)𝑑p𝑑E`$
where $`U_0(x)<E_0`$, and $`\rho _0`$ is zero else; here $`r(x):=\sqrt{x_1^2+x_2^2}`$. Assume that $`lim_{|x|\mathrm{}}U_0(x)=0`$, see the discussion above. Then clearly the cut-off energy cannot be positive, i. e., $`E_00`$, since otherwise we get infinite mass “at spatial infinity”. Since $`U_0<0`$ everywhere, (3.2) shows that $`\rho _0`$ will have compact support if and only if $`E_0<0`$. An additional, sufficient condition which implies this is

$$f_fQ(f,P)3Q(f,P),f0,P\mathrm{IR},$$
(3.3)
which holds for example for $`Q(f,P)=f^{1+1/k}g(P)`$ with $`1/2k<3/2`$ and some function $`g`$. To see this we rewrite the formula for $`E_0`$ from Theorem 1:

$$E_0=\frac{1}{M}\left[E_{\mathrm{kin}}(f_0)+2E_{\mathrm{pot}}(f_0)+f_0_fQ(f_0,P)dvdx\right].$$

Now we note that for solutions of the Vlasov-Poisson system,

$$\frac{d}{dt}xvf=2E_{\mathrm{kin}}(f)+E_{\mathrm{pot}}(f)$$

so that for a steady state the right hand side is zero. Thus by (3.3),

$`E_{\mathrm{kin}}(f_0)+2E_{\mathrm{pot}}(f_0)+{\displaystyle f_0_fQ(f_0,P)dvdx}`$ $``$
$`E_{\mathrm{kin}}(f_0)+2E_{\mathrm{pot}}(f_0)+3{\displaystyle Q(f_0,P)𝑑v𝑑x}`$ $`=`$
$`3E_{\mathrm{kin}}(f_0)+3E_{\mathrm{pot}}(f_0)+3𝒞(f_0)`$ $`=`$ $`3_C(f_0)<0`$
by Lemma A 3 (a). An alternative condition which also guarantees compact support is the following:

$$\begin{array}{c}\varphi (E_0E,P)C_1(E_0E)^{k_1},E\mathrm{},P\mathrm{IR},\\ \varphi (E_0E,P)C_2(E_0E)^{k_2},EE_0,P\mathrm{IR}\end{array}\}$$
(3.4)
for positive constants $`C_1,C_2`$ and $`0<k_1,k_2<3/2`$. This assumption also implies that $`E_0<0`$, cf. \[24, Thm. 3\].

In the spherically symmetric case our approach provides an explicit bound on the radius of the steady state, cf. . Note that in this case one can under appropriate assumptions decide whether $`U_0`$ crosses the cut-off energy level $`E_0`$ by direct examination of the semilinear Poisson equation (2.2), since this equation reduces to an ordinary differential equation with respect to the radial variable, cf. . However, in the axially symmetric case this is no longer true, and the corresponding analysis of the genuine partial differential equation (2.2) would be much more difficult.

Symmetry. By construction the steady states which we obtain are axially symmetric, but so are the spherically symmetric ones whose existence was known already. In order to make sure that we have found qualitatively new steady states we must show that they are in general not spherically symmetric. To see this we take $`Q`$ such that $`\varphi (E_0E,)`$ is not constant on any neighborhood of $`P=0`$. We claim that neither $`\rho _0`$ nor $`U_0`$ are spherically symmetric in this case. Indeed, if $`\rho _0`$ were spherically symmetric the same were true for $`U_0`$. If $`U_0`$ were spherically symmetric we take some $`x\mathrm{IR}^3`$ with $`r(x)=\sqrt{x_1^2+x_2^2}0`$ small and take $`A\mathrm{SO}(3)`$ such that $`Ax=(0,0,|x|)`$, hence $`r(Ax)=0r(x)`$ but $`U_0(Ax)=U_0(x)`$. Inserting this into the formula (3.2) for $`\rho _0`$ we see that in general $`\rho _0(Ax)\rho _0(x)`$ so $`\rho _0`$ is not spherically symmetric.

Indeed, this is a bit more than just saying that $`f_0`$ is not spherically symmetric. To see this, take a spherically symmetric steady state $`f_0`$ with induced density $`\rho _0`$, potential $`U_0`$, and particle energy $`E=|v|^2/2+U_0(x)`$. Let $`g_0(x,v)=\psi (E,P)`$ be any function which is odd in $`P`$, and such that $`f_0+g_00`$. Since $`P`$ is not invariant under general rotations neither is $`f_0+g_0`$, but the fact that $`\psi `$ is odd in $`P`$ implies that $`\rho _{f_0+g_0}=\rho _0`$ and the same holds for the potential. So this trivial construction gives an axially symmetric steady state where the phase space density is not spherically symmetric but the macroscopic quantities $`\rho _0`$ and $`U_0`$ are. Our steady states are in general not of this trivial type.

Another symmetry issue refers to the dependence of the minimizer on the variable $`x_3`$. Given a minimizer $`f_0`$ of which we do not yet know any symmetry with respect to $`x_3`$ we can do a symmetric decreasing rearrangement of $`f_0`$ with respect to $`x_3`$ while keeping all other variables fixed. Denote this rearrangement by $`f_0^{}`$; as to the rearrangement concept we refer the reader to \[16, Ch. 3\]. Obviously, the rearrangement does not change the kinetic energy, and neither does it change the Casimir functional, since $`P`$ does not depend on $`x_3`$. By \[16, Thms. 3.7, 3.9\] it can at most decrease the potential energy, with equality iff $`f_0`$ is already symmetric and decreasing in $`x_3`$ up to a shift in $`x_3`$. But since $`f_0`$ minimizes the energy-Casimir functional it must posses this symmetry, and without loss of generality we can assume that $`f_0`$ is reflexion symmetric in $`x_3`$. The corresponding symmetry in the variable $`v_3`$ follows from the form which $`f_0`$ must have due to Theorem 1.

Stationary versus static solutions. If instead of the Vlasov-Poisson system we consider a self-gravitating fluid as described by the Euler-Poisson system then every so-called static solution, i. e., every steady state with vanishing velocity field, is spherically symmetric, cf. . This turns out to be false for the steady states of the Vlasov-Poisson system which we obtain. By definition the velocity field equals $`j_0/\rho _0`$ on the support of $`\rho _0`$ where the mass current density $`j_0`$ is given by

$$j_0(x)=vf_0(x,v)𝑑v=2\pi _{U_0(x)}^{E_0}_{\sqrt{2(EU_0(x))}}^{\sqrt{2(EU_0(x))}}s\varphi (E_0E,r(x)p)𝑑p𝑑Ee_t(x)$$

with $`e_t(x):=(x_2,x_1,0)/r(x)`$; on the $`x_3`$-axis the velocity field vanishes. Obviously, if $`\varphi `$ is even in $`P`$ then $`j_0`$ vanishes identically, so there exist static solutions which are not spherically symmetric among the ones we obtain. On the other hand, if, say, $`\varphi (E,P)>\varphi (E,P)`$ for all $`E0`$ and $`P>0`$ the velocity field is non-trivial and corresponds to an average rotation about the axis of symmetry in the counterclockwise direction.

Less regular dependence on $`P`$. When examining the appendix the reader may notice that at no place we really make use of the fact that all the functions under consideration can be taken to be axially symmetric. This is because as opposed to spherical symmetry, exploiting axial symmetry for example in the estimates for the potential is technically quite unpleasant. However, it is possible and even necessary if one wishes to study examples of the form

$$Q(f,P)=f^{1+1/k}P^{2l/k},f0,P\mathrm{IR}$$
(3.5)
which for $`l0`$ do not satisfy the assumptions on $`Q`$ which we stated above. The corresponding steady states would take the form

$$f_0(x,v)=C(E_0E)^kP^{2l}$$

which is analogous to the classical, spherically symmetric polytropes, except that the square of the third component of the angular momentum replaces the square of the modulus of the angular momentum. Clearly, (3.5) satisfies the crucial convexity condition (Q2), provided $`k>0`$. When we examine the scaling arguments in Lemma A 3 we find that they go through, provided

$$l>1,0<k<l+\frac{3}{2}.$$

The place where significant changes become necessary are the estimates in Lemma A 1 and Lemma A 6. Clearly, one can no longer control $`f_{1+1/k}`$ and $`\rho _f_{1+1/n}`$ in terms of $`E_{\mathrm{kin}}(f)`$, $`𝒞(f)`$, and $`M`$. Instead, one gains control of the weighted norms

$$P^{l/k}(x,v)f^{1+1/k}(x,v)𝑑v𝑑x,r^{2l/n}(x)\rho _f^{1+1/n}(x)𝑑x$$

where $`r(x)=\sqrt{x_1^2+x_2^2}`$. We conjecture that the weight factor can be dealt with in the estimates for the potential energy etc. if one makes proper use of the axial symmetry of the potential.

## 4 Appendix

We now give the proof of the existence part in Theorem 1 and the additional results used in the proof of Theorem 2. First we collect some estimates for $`\rho _f`$ and $`U_f`$ induced by an element $`f_M`$. These estimates make no use of symmetry and rely on the assumption (Q1). Their main point is to see that $`_C`$ is bounded from below on $`_M`$ and to establish certain bounds along minimizing sequences of $`_C`$. Constants denoted by $`C`$ are positive, may depend only on $`Q`$ and $`M`$, and their value may change from line to line. By $`_p`$ we denote the usual $`L^p`$-norm of functions over $`\mathrm{IR}^3`$ or $`\mathrm{IR}^6`$ as the case may be.

###### Lemma A 1

Let $`n:=k+3/2<3`$. Then for any $`f_M`$ the following holds:

* $`{\displaystyle f^{1+1/k}(x,v)𝑑v𝑑x}C(1+𝒞(f))`$
* $`{\displaystyle \rho _f^{1+1/n}(x)𝑑x}C\left(1+E_{\mathrm{kin}}(f)+𝒞(f)\right)`$
* $`U_fL^6(\mathrm{IR}^3)`$ with $`U_fL^2(\mathrm{IR}^3)`$, the two forms of $`E_{\mathrm{pot}}(f)`$ stated in Section 2 are equal, and

$$|U_f|^2𝑑xC\rho _f_{6/5}^2C\left(1+E_{\mathrm{kin}}(f)+𝒞(f)\right)^{n/3}.$$

Proof. Part (a) is a direct consequence of (Q1), splitting the $`v`$-integral in the definition of $`\rho _f`$ into small and large $`v`$’s and optimizing the split yields (b), while the extended Young’s inequality and interpolation together with (b) imply (c). For details cf.  or . $`\mathrm{}`$

The estimates above have an immediate but important consequence:

###### Lemma A 2

The energy-Casimir functional $`_C`$ is bounded from below on $`_M`$, i. e.,

$$h_M:=\underset{f_M}{inf}_C(f)>\mathrm{},$$

and along any minimizing sequence of $`_C`$ in $`_M`$ the quantities $`E_{\mathrm{kin}}(f)`$, $`𝒞(f)`$, $`f_{1+1/k}`$, and $`\rho _f_{1+1/n}`$ are bounded.

Proof. Lemma A 1 yields the estimate

$$_C(f)E_{\mathrm{kin}}(f)+𝒞(f)C\left(1+E_{\mathrm{kin}}(f)+𝒞(f)\right)^{n/3},f_M,$$

and since $`n<3`$ the assertions follow. $`\mathrm{}`$

Lemma A 1 and Lemma A 2 remain valid if we replace $`_M`$ by its axially symmetric subset $`_M^S`$. Of course the infimum of $`_C`$ on this smaller set may be larger; we denote

$$h_M^S:=\underset{f_M^S}{inf}_C(f).$$

The assumptions (Q3) and (Q4) determine the behavior of $`_C`$ under scaling transformations which we use to show that $`h_M`$ is negative and to relate the $`h_M`$’s for different values of $`M`$:

###### Lemma A 3

* Let $`M>0`$. Then $`\mathrm{}<h_M<0`$.
* For $`0<\overline{M}M`$,

$$h_{\overline{M}}\left(\overline{M}/M\right)^{5/3}h_M.$$

Proof. Given any function $`f`$, we define a rescaled function $`\overline{f}(x,v)=af(bx,cv)`$, where $`a,b,c>0`$. Then

$$\overline{f}𝑑v𝑑x=a(bc)^3f𝑑v𝑑x,$$
(4.1)
i. e. $`f_M`$ iff $`\overline{f}_{\overline{M}}`$ where $`\overline{M}=a(bc)^3M`$. Next

$`_C(\overline{f})`$ $`=`$ $`ab^3c^5E_{\mathrm{kin}}(f)+a^2b^5c^6E_{\mathrm{pot}}(f)`$ (4.2)
$`+(bc)^3{\displaystyle Q(af(x,v),(bc)^1P(x,v))𝑑v𝑑x}.`$
To prove (a) we fix a function $`f_1`$ with $`fF^{}`$ and let $`a=M(bc)^3`$ so that $`\overline{f}_M`$. Then by (Q3) and (Q4) and with positive constants $`C_1,C_2,C_3`$ which depend on $`f`$,

$`_C(\overline{f})`$ $``$ $`Mc^2E_{\mathrm{kin}}(f)+M^2bE_{\mathrm{pot}}(f)+Ma^1{\displaystyle Q(af(x,v),0)𝑑v𝑑x}`$
$``$ $`C_1c^2C_2b+C_3a^{1/k^{}},`$
provided $`a1`$ so that (Q4) can be applied; note that $`E_{\mathrm{pot}}(f)<0`$. We want the negative term to dominate as $`b0`$, so we let $`c=b^{\gamma /2}`$ for some $`\gamma >0`$. Then $`a=Mb^{3(1\gamma /2)}`$, and

$$_C(\overline{f})C_1b^\gamma C_2b+C_3M^{1/k^{}}b^{3(1\gamma /2)/k^{}}.$$

Since $`0<k^{}<3/2`$ we can fix $`\gamma ]1,2[`$ such that $`3(1\gamma /2)/k^{}>1`$. For $`b>0`$ sufficiently small $`_C(\overline{f})`$ will then be negative and $`a=Mb^{3(1\gamma /2)}<1`$ as required. This proves part (a) of the lemma. To prove (b) we choose $`a=c=1`$ and $`b=(\overline{M}/M)^{1/3}`$ so that the mapping $`_M_{\overline{M}}`$, $`f\overline{f}`$ is one-to-one and onto and $`b^11`$. By (4.2),

$`_C(\overline{f})`$ $`=`$ $`b^3E_{\mathrm{kin}}(f)+b^5E_{\mathrm{pot}}(f)+b^3{\displaystyle Q(f(x,v),b^1P(x,v))𝑑v𝑑x}`$
$``$ $`b^5E_{\mathrm{kin}}(f)+b^5E_{\mathrm{pot}}(f)+b^5{\displaystyle Q(f(x,v),P(x,v))𝑑v𝑑x}`$
$`=`$ $`b^5_C(f);`$
we multiplied the two positive terms by $`b^21`$ and used the monotonicity of $`Q`$ in the $`P`$-variable which we required in (Q3). By the choice of $`b`$ and the definitions of $`h_M`$ and $`h_{\overline{M}}`$ the proof is complete. $`\mathrm{}`$

It is again obvious that the assertions of the lemma remain valid if we restrict ourselves to the axially symmetric functions in $`_M`$. The reader might also check that the scaling arguments work under the assumptions (Q3’) and (Q4’) from Section 3 as well.

Next we provide a splitting estimate which will be used to show that along a minimizing sequence the mass cannot escape to infinity; here and in the following we denote for $`0<R<S\mathrm{}`$,

$`B_R`$ $`:=`$ $`\{x\mathrm{IR}^3||x|R\},`$
$`B_{R,S}`$ $`:=`$ $`\{x\mathrm{IR}^3|R|x|<S\}.`$
###### Lemma A 4

Let $`f_M`$. Then

$$\underset{a\mathrm{IR}^3}{sup}_{a+B_R}f(x,v)𝑑v𝑑x\frac{1}{RM}\left(2E_{\mathrm{pot}}(f)\frac{M^2}{R}\frac{C\rho _f_{1+1/n}}{R^{(5n)/(n+1)}}\right),R>1.$$

The proof follows from splitting the potential energy into three parts according to $`|xy|1/R`$, $`1/R<|xy|<R`$, and $`R|xy|`$; for the details we refer to  or \[23, Lemma 3\] The splitting estimate above has the following important consequence for minimizing sequences:

###### Lemma A 5

Let $`(f_i)_M`$ be a minimizing sequence of $`_C`$. Then there exist $`\delta _0>0`$, $`R_0>0`$, $`i_0\mathrm{IN}`$, and a sequence of shift vectors $`(a_i)\mathrm{IR}^3`$ such that

$$_{a_i+B_R}f_i(x,v)𝑑v𝑑x\delta _0,ii_0,RR_0.$$

If $`(f_i)_M^S`$, i. e., the minimizing sequence consists of axially symmetric functions then one can choose $`a_i=(0,0,z_i)`$ with $`z_i\mathrm{IR}`$, i. e., only shifts along the axis of symmetry need to be admitted.

Proof. By Lemma A 2, $`(\rho _{f_i}_{1+1/n})`$ is bounded. By Lemma A 3 (a) we have

$$E_{\mathrm{pot}}(f_i)_C(f_i)\frac{1}{2}h_M<0,ii_0,$$

for a suitable $`i_0\mathrm{IN}`$. Thus by Lemma A 4 there exist $`\delta _0>0`$, $`R_0>0`$, and a sequence of shift vectors $`(a_i)\mathrm{IR}^3`$ as required.

Now assume that the functions $`f_i`$ are axially symmetric, and let $`a_i=(\overline{a}_i,z_i)`$ with $`\overline{a}_i\mathrm{IR}^2`$ and $`z_i\mathrm{IR}`$. Suppose we can show that the sequence $`(\overline{a}_i)`$ is bounded. Then we can pick some $`\overline{R}>0`$ such that

$$(\overline{a}_i,0)+B_{R_0}B_{\overline{R}},i\mathrm{IN}$$

and thus

$$a_i+B_{R_0}(0,0,z_i)+B_{\overline{R}},i\mathrm{IN},$$

which implies our assertion in the symmetric case: For $`R\overline{R}`$,

$$_{(0,0,z_i)+B_R}f_i(x,v)𝑑v𝑑x_{a_i+B_{R_0}}f_i(x,v)𝑑v𝑑x\delta _0,ii_0.$$

So it suffices to show that $`(\overline{a}_i)`$ is bounded. Assume the contrary and choose some integer $`N\mathrm{IN}`$ such that $`N\delta _0>M`$. If $`|\overline{a}_i|`$ is sufficiently large then simple geometry tells us that there exist rotations $`A_1,\mathrm{},A_N`$ about the $`x_3`$-axis such that the balls $`A_ka_i+B_{R_0},k=1,\mathrm{},N`$, are pairwise disjoint. By axial symmetry of $`f_i`$,

$$_{A_ka_i+B_{R_0}}f_i(x,v)𝑑v𝑑x=_{a_i+B_{R_0}}f_i(x,v)𝑑v𝑑x\delta _0,k=1,\mathrm{},N,$$

and since the balls on the left hand side are pairwise disjoint,

$$f_i(x,v)𝑑v𝑑x\underset{k=1}{\overset{N}{}}_{A_ka_i+B_{R_0}}f_i(x,v)𝑑v𝑑xN\delta _0>M$$

which violates the mass constraint and gives the desired contradiction. $`\mathrm{}`$

We will also need to exploit the well known compactness properties of the solution operator of the Poisson equation:

###### Lemma A 6

Let $`(\rho _i)L^{1+1/n}(\mathrm{IR}^3)`$ be bounded and

$$\rho _i\rho _0\text{weakly in}L^{1+1/n}(\mathrm{IR}^3).$$

* For any $`R>0`$,

$$U_{\mathrm{𝟏}_{B_R}\rho _i}U_{\mathrm{𝟏}_{B_R}\rho _0}\text{strongly in}L^2(\mathrm{IR}^3).$$
* If in addition $`(\rho _i)`$ is bounded in $`L^1(\mathrm{IR}^3)`$, $`\rho _0Ł^1(\mathrm{IR}^3)`$, and for any $`ϵ>0`$ there exists $`R>0`$ and $`i_0\mathrm{IN}`$ such that

$$_{|x|R}|\rho _i(x)|𝑑x<ϵ,ii_0$$

then

$$U_{\rho _i}U_{\rho _0}\text{strongly in}L^2(\mathrm{IR}^3).$$

For the details of the proof we refer to  or \[23, Lemma 5\]; it relies on the fact that over bounded sets the solution operator to the Poisson equation is compact in the appropriate spaces and the fact that in the context of the lemma we can make the contribution to the potential energy coming from outside sufficiently large balls small.

We are now ready to prove the existence of a minimizer, more precisely we prove:

###### Theorem A 1

Let $`(f_i)_M^S`$ be a minimizing sequence of $`_C`$. Then there exists a sequence of shifts $`(z_i)\mathrm{IR}`$ and a subsequence, again denoted by $`(f_i)`$, such that for any $`ϵ>0`$ there exists $`R>0`$ with

$$_{(0,0,z_i)+B_R}f_i(x,v)𝑑v𝑑xMϵ,i\mathrm{IN},$$

$$Tf_i:=f_i(+(0,0,z_i),)f_0\text{weakly in}L^{1+1/k}(\mathrm{IR}^6),i\mathrm{},$$

and

$$_{B_R}f_0Mϵ.$$

Finally,

$$U_{Tf_i}U_0\text{strongly in}L^2(\mathrm{IR}^3),i\mathrm{},$$

and $`f_0_M^S`$ is a minimizer of $`_C`$.

Proof. We split $`f_M^S`$ into three different parts:

$$f=\mathrm{𝟏}_{B_{R_1}\times \mathrm{IR}^3}f+\mathrm{𝟏}_{B_{R_1,R_2}\times \mathrm{IR}^3}f+\mathrm{𝟏}_{B_{R_2,\mathrm{}}\times \mathrm{IR}^3}f=:f_1+f_2+f_3;$$

the parameters $`R_1<R_2`$ of the split are yet to be determined. The induced spacial densities are denoted by $`\rho _k`$, their masses by $`M_k`$, and the induced potentials by $`U_k`$, $`k=1,2,3`$. With

$$I_{lm}:=\frac{\rho _l(x)\rho _m(y)}{|xy|},l,m=1,2,3,$$

we have

$$_C(f)=_C(f_1)+_C(f_2)+_C(f_3)I_{12}I_{13}I_{23}.$$

If we choose $`R_2>2R_1`$ then

$$I_{13}\frac{C}{R_2}.$$

Next, we use the Cauchy-Schwarz inequality, the extended Young’s inequality, and interpolation to get

$`I_{12}+I_{23}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}\left|{\displaystyle (U_1+U_3)U_2dx}\right|C\rho _1+\rho _3_{6/5}U_2_2`$
$``$ $`C\rho _{1+1/n}^{(n+1)/6}U_2_2.`$
Using the estimates above and Lemma A 3 (b) we find

$`h_M^S_C(f)`$ $``$ $`\left(1\left({\displaystyle \frac{M_1}{M}}\right)^{5/3}\left({\displaystyle \frac{M_2}{M}}\right)^{5/3}\left({\displaystyle \frac{M_3}{M}}\right)^{5/3}\right)h_M^S`$ (4.3)
$`+C\left(R_2^1+\rho _{1+1/n}^{(n+1)/6}U_2_2\right)`$
$``$ $`{\displaystyle \frac{C}{M^2}}\left(M_1M_2+M_1M_3+M_2M_3\right)h_M^S`$
$`+C\left(R_2^1+\rho _{1+1/n}^{(n+1)/6}U_2_2\right)`$
$``$ $`Ch_M^SM_1M_3+C\left(R_2^1+\rho _{1+1/n}^{(n+1)/6}U_2_2\right);`$
observe that by Lemma A 3 (a) $`h_M^S<0`$ and that constants denoted by $`C`$ are positive and depend on $`M`$ and $`Q`$, but not on $`R_1`$ or $`R_2`$. We want to use (4.3) to show that up to a subsequence and a shift $`M_3`$ becomes small along any minimizing sequence for $`i`$ large provided the splitting parameters are suitably chosen.

So let $`(f_i)_M^S`$ be a minimizing sequence and define $`R_0`$, $`\delta _0`$, and $`(z_i)`$ according to Lemma A 5. The shifted sequence $`(Tf_i)`$ is again minimizing so by Lemma A 2 we can pick a subsequence, again denoted by $`(f_i)`$, with a weak limit as stated in the theorem. Now choose $`R_1>R_0`$ so that by Lemma A 5, $`M_{i,1}\delta _0`$ for $`i`$ large. By (4.3),

$$Ch_M^S\delta _0M_{i,3}\frac{C}{R_2}+CU_{0,2}_2+CU_{i,2}U_{0,2}_2+_C(Tf_i)h_M^S$$
(4.4)
where $`U_{i,l}`$ is the potential induced by $`f_{i,l}`$ which in turn has mass $`M_{i,l}`$, $`i\mathrm{IN}\{0\}`$, and the index $`l=1,2,3`$ refers to the splitting. Given any $`ϵ>0`$ we increase $`R_1>R_0`$ such that the second term on the right hand side of (4.4) is small, say less than $`ϵ/4`$. Next choose $`R_2>2R_1`$ such that the first term is small. Now that $`R_1`$ and $`R_2`$ are fixed, the third term in (4.4) converges to zero by Lemma A 6 (a). Since $`(Tf_i)`$ is minimizing the remainder in (4.4) follows suit. Therefore, for $`i`$ sufficiently large,

$$_{(0,0,z_i)+B_{R_2}}Tf_i=MM_{i,3}M(Ch_M^S\delta _0)^1ϵ.$$
(4.5)
Clearly, $`f_00`$ a. e., and $`f_0`$ is axially symmetric. By weak convergence we have that for any $`ϵ>0`$ there exists $`R>0`$ such that

$$M_{B_R}f_0𝑑v𝑑xMϵ$$

which in particular implies that $`f_0L^1(\mathrm{IR}^6)`$ with $`f_0𝑑v𝑑x=M`$. The functional $`𝒞`$ is convex, so by Mazur’s Lemma and Fatou’s Lemma

$$𝒞(f_0)\underset{i\mathrm{}}{lim\; sup}𝒞(Tf_i).$$

The strong convergence of the gravitational fields now follows by Lemma A 6 (b), and in particular,

$$_C(f_0)\underset{i\mathrm{}}{lim\; sup}_C(Tf_i)=h_M^S$$

so that $`f_0_M^S`$ is a minimizer of $`_C`$. $`\mathrm{}`$

Acknowledgment: The research of YG is supported in part by an NSF grant and a Sloan Fellowship. GR acknowledges support by the Wittgenstein 2000 Award of P. A. Markowich.


format: text

text:
# The Formation of Massive Stars and Star Clusters

## 1. Introduction

If there were no massive ($`m_{}>8M_{}`$) stars, the heavens would appear quite dull. Deprived of many of the heavy elements, it is not clear that Earth-like astronomers could even exist to contemplate the obscurity of such a universe. Fortunately, back in our own cosmos, we can be thankful for the existence of high-mass stars and perhaps devote some effort to understand them and how they come about. All massive stars appear to form in clusters, where they are greatly outnumbered by their lower-mass cousins, and so the questions of how an individual massive star forms and how the cluster forms are intimately entwined.

In this article we first outline a theory for the formation of individual stars, including massive ones, by accretion (§2, taken from McKee & Tan 2002; hereafter MT02). High pressures in the star-forming region result in rapid accretion rates, which allow massive stars to form in spite of the fierce glare of their radiation (Wolfire & Cassinelli 1987), and to do so swiftly before stellar evolution has time to age the newborn stars. In certain regions of the Galaxy, the required high pressures are known to result from the weight of hundreds to thousands of solar masses of strongly self-gravitating and turbulent gas, condensed into clumps of order a parsec across. A large fraction of this material inevitably collapses into a stellar cluster, with most of the mass partitioned into relatively low-mass stars. The massive stars are expected to form preferentially at the center of the clump, where the pressures are highest. We apply our theory for individual star formation to model the birth of the cluster (§§3 & 4, taken from Tan & McKee 2002a; hereafter TM02). Extrapolating models of outflows from low-mass protostars to those from higher masses, we use observations of protocluster outflow intensities to estimate the star formation rates and cluster formation times.

## 2. Collapse of Turbulent Cloud Cores and Massive Star Formation

Our basic assumption is that a star forms from the gravitational collapse of a molecular cloud core that begins in approximate hydrostatic equilibrium; in particular, it does not form by the coalescence of low-mass stars (Bonnell, Bate, & Zinnecker 1998) or by the sweeping up of ambient clump material (Bonnell et al. 1997). On dimensional grounds we then expect that the protostellar accretion rate is given by

$$\dot{m}_{}=\varphi _{}\frac{m_{}}{t_{\mathrm{ff}}},$$
(1)
where $`t_{\mathrm{ff}}=(3\pi /32G\rho )^{1/2}`$ is the free-fall time and $`\varphi _{}`$ is a dimensionless constant of order unity. This equation has the same dependence on dimensional parameters as that for isothermal collapse, $`\dot{m}_{}c_{\mathrm{th}}^3/G`$ (Shu 1977), if the thermal sound speed, $`c_{\mathrm{th}}`$, is replaced by the virial velocity $`(Gm_{}/R)^{1/2}`$, since $`(Gm_{}/R)^{3/2}/Gm_{}(Gm_{}/R^3)^{1/2}m_{}/t_{\mathrm{ff}}`$. We have shown (MT02) that if the collapse is spherical and self-similar, then $`\varphi _{}`$ is quite close to unity provided that the value of $`\rho `$ entering $`t_{\mathrm{ff}}`$ is evaluated at the radius in the initial cloud that just encloses the gas that goes into the star when its mass is $`m_{}`$. We summarize these results below. Collapse to the star will naturally proceed via a disk, where we assume instabilities drive high inflow rates if the disk mass grows to be a significant fraction of the star’s (Shu et al. 1990).

First we allow for the disruptive effect of protostellar outflows on the gas core so that only a fraction $`ϵ_{\mathrm{core}}=m_{}/M`$ of the core mass can accrete onto the star. The value of $`ϵ_{\mathrm{core}}`$ is uncertain for massive stars, since it depends on the unknown properties of their protostellar winds. We assume that it is constant over the accretion history of the star, and for numerical estimates set $`ϵ_{\mathrm{core}}=0.5`$, similar to the value calculated for low-mass stars (Matzner & McKee 2000).

For the core, we assume a spherical power-law ambient medium, which ensures accretion is self-similar. With $`\rho r^{k_\rho }`$ and $`Pr^{k_P}`$, it follows that the core is a polytrope with $`P\rho ^{\gamma _p}`$. In hydrostatic equilibrium (HSE) $`k_\rho =2/(2\gamma _p)`$ and $`k_P=\gamma _pk_\rho =2\gamma _p/(2\gamma _p)`$ (McLaughlin & Pudritz 1996). Let $`c(P/\rho )^{1/2}`$ be the effective sound speed. The equation of HSE then gives $`M=(k_Pc^2r)/G`$ and $`\rho =Ac^2/(2\pi Gr^2)`$, with $`A=\gamma _p(43\gamma _p)/(2\gamma _p)^2`$. Equation (1) then implies

$$\dot{m}_{}=\varphi _{}ϵ_{\mathrm{core}}\frac{4}{\pi \sqrt{3}}k_PA^{1/2}\frac{c^3}{G},$$
(2)
which is a generalization of the isothermal accretion rate to the nonthermal case.

For the polytropic sphere, $`\rho r^{k_\rho }`$ and $`Mr^{3k_\rho }`$, which implies $`\rho M^{k_\rho /(3k_\rho )}=M^{2/(43\gamma _p)}`$. Thus we have $`\dot{m}_{}m_{}\rho ^{1/2}m_{}^{11/(43\gamma _p)}`$ for the mass dependence of the accretion rate. Integration yields

$$m_{}=m_f\left(\frac{t}{t_f}\right)^{43\gamma _p},$$
(3)
where $`m_f`$ is the final stellar mass, which is attained at a time $`t_f`$ (McLaughlin & Pudritz 1997, hereafter MP97). Note that for $`\gamma _p<1`$ the accretion rate accelerates (MP97). Equation (3) implies

$$\dot{m}_{}=(43\gamma _p)\frac{m_{}}{t}=(43\gamma _p)\frac{m_f}{t_f}\left(\frac{t}{t_f}\right)^{33\gamma _p}.$$
(4)
As discussed by MP97, termination of the accretion breaks the self-similarity once the expansion wave reaches $`m_f`$. This occurs at a time they label $`t_{\mathrm{ew}}`$, which is about $`0.45t_f`$. Thereafter, equation (4) becomes approximate, but they argue that the approximation should be reasonably good. From equations (1) and (4), the star-formation time is

$$t_f=\frac{(43\gamma _p)}{\varphi _{}}t_{\mathrm{ff}}.$$
(5)
Using the results of MP97, we can evaluate $`\varphi _{}`$ in the non-magnetic case,

$$\varphi _{\mathrm{non}}=(43\gamma _p)\frac{t_{\mathrm{ff}}}{t_f}=\pi \sqrt{3}\left[\frac{(2\gamma _p)^2(43\gamma _p)^{(76\gamma _p)/2}m_0}{8^{(53\gamma _p)/2}}\right]^{1/(43\gamma _p)},$$
(6)
where $`m_0`$ is tabulated by MP97. For example, for the singular isothermal sphere the parameter $`m_0=0.975`$ and $`\varphi _{\mathrm{non}}=0.975\pi \sqrt{3}/8=0.663`$. For other values of $`\gamma _p`$ in the range $`0\gamma _p1`$, $`\varphi _{\mathrm{non}}1.13/(1+0.7\gamma _p^2)`$. We conclude that $`\dot{m}_{}m_{}/t_{\mathrm{ff}}`$ to within a factor 1.5 for spherical cores in which the effective sound speed increases outward.

We can estimate the effect of magnetic fields on the accretion rate from the work of Li & Shu (1997), who considered collapse of self-similar, isothermal, magnetized, toroidal clouds. The equilibrium surface density is $`\mathrm{\Sigma }=(1+H_0)c_{\mathrm{th}}^2/(\pi G\varpi )`$, where $`\varpi `$ is the cylindrical radius and $`H_0`$ is a parameter that increases from zero as the magnetic field is increased. They show that the accretion rate is $`\dot{m}_{}=1.05(1+H_0)c_{\mathrm{th}}^3/G`$, which is larger than the isothermal case by about a factor $`(1+H_0)`$. However, equation (1) predicts $`\dot{m}_{}M\rho ^{1/2}M^{3/2}/\varpi ^{3/2}\mathrm{\Sigma }^{3/2}\varpi ^{3/2}(1+H_0)^{3/2}`$. To reconcile this result with the correct value, we require $`\varphi _{}\varphi _{\mathrm{non}}/(1+H_0)^{1/2}`$. They regard $`H_01`$ as typical for low-mass star formation, which would lead to a reduction in the accretion rate (as expressed in eqs.  and ) by a factor of about 1.4.

The power-law structure of cloud cores is truncated by the ambient pressure in the star-forming clump, which is thus equivalent to the core surface pressure, $`P_s`$. We parameterize the rate of individual core collapse in terms of this pressure, which we then relate to the clump’s surface density. A higher ambient pressure means that a core of a given mass is truncated at higher density and thus its free-fall and star formation timescales are shorter.

Since $`\rho _s=P_s/c_s^2`$, $`M=(k_Pc^2r)/G`$ and $`\rho =Ac^2/(2\pi Gr^2)`$, the density at the surface of the core is

$$\rho _s=\left(\frac{Ak_P^2ϵ_{\mathrm{core}}^2P_s^3}{2\pi G^3m_f^2}\right)^{1/4}.$$
(7)
The general expression for the accretion rate in terms of the surface pressure and the final stellar mass can now be inferred from equations (1) and (7):

$$\dot{m}_{}=4.02\times 10^4\varphi _{}(Ak_P^2ϵ_{\mathrm{core}}^2)^{\frac{1}{8}}\left(\frac{m_f}{30\mathrm{M}_{}}\right)^{\frac{3}{4}}\left(\frac{P_s}{10^8\mathrm{Kcm}^3}\right)^{\frac{3}{8}}\left(\frac{m_{}}{m_f}\right)^{\frac{3(2k_\rho )}{2(3k_\rho )}}\mathrm{M}_{}\mathrm{yr}^1$$
(8)
and the corresponding value of the star-formation time is

$$t_f=7.47\times 10^4\left(\frac{43\gamma _p}{\varphi _{}A^{1/8}k_P^{1/4}ϵ_{\mathrm{core}}^{1/4}}\right)\left(\frac{m_f}{30\mathrm{M}_{}}\right)^{1/4}\left(\frac{10^8\mathrm{K}\mathrm{cm}^3}{P_s}\right)^{3/8}\mathrm{yr}.$$
(9)
The accretion rate and star formation time depend to some extent on the value of $`k_\rho `$ in the protostellar cores. No data are available on the structure of cores that are forming very massive stars. For “high-mass” cores in Orion, Caselli & Myers (1995) find $`k_\rho 1.45`$ with a dispersion of $`\pm 0.2`$. According to van der Tak et al. (2000), the clumps in which high-mass cores are embedded have values of $`k_\rho `$ ranging from 1 to 2, centered around 1.5. For numerical evaluation of the accretion rate, we shall adopt $`k_\rho =1.5`$, which corresponds to $`\gamma _p=2/3`$, $`k_P=1`$, $`A=3/4`$, and $`\varphi _{\mathrm{non}}=0.90`$. In this case, the protostellar accretion rate varies linearly with time, $`\dot{m}_{}t`$. So long as $`A\gamma _p=2(k_\rho 1)/k_\rho `$ is not near zero, the star-formation time is relatively insensitive to $`k_\rho `$. For $`k_\rho =1.5`$, the numerical coefficient in equation (8) is $`3.49\times 10^4/(1+H_0)^{1/2}\mathrm{M}_{}\mathrm{yr}^1`$ and $`\dot{m}_{}(m_{}/m_f)^{0.5}`$. Similarly, the numerical coefficient in equation (9) becomes $`1.72\times 10^5(1+H_0)^{1/2}\mathrm{yr}`$.

At present, direct observations of the total pressure at the surface of a core that will form a massive star are not available, so we must estimate this pressure theoretically. Since the core is embedded in a turbulent clump, its pressure will fluctuate, but the time-averaged surface pressure for the typical core should be approximately the same as the spatial average of the pressure in the clump, $`P_s=\varphi _P\overline{P}`$, where $`\varphi _P`$ is a parameter of order unity. The mean pressure in a clump is (MT02):

$$\overline{P}=\left(\frac{3\pi f_{\mathrm{gas}}\varphi _B\alpha _{\mathrm{vir}}}{20}\right)G\mathrm{\Sigma }^2,$$
(10)
where $`f_{\mathrm{gas}}`$ is the fraction of the clump mass in gas (as opposed to stars), $`\alpha _{\mathrm{vir}}5\sigma ^2R_{\mathrm{cl}}/GM_{\mathrm{cl}}M_{\mathrm{vir}}/M_{\mathrm{cl}}`$ is the virial parameter, $`\sigma `$ is the one-dimensional velocity dispersion, $`x`$ indicates a mass-averaged quantity, $`\varphi _B=c^2/\sigma ^21.3+1.5/m_A^2`$ represents the contribution of the magnetic pressure, $`m_A`$ is an average of the Alfven Mach number, and $`\mathrm{\Sigma }=M_{\mathrm{cl}}/(\pi R_{\mathrm{cl}}^2)`$ is the surface density of the clump, including the mass of any embedded stars. In terms of $`\mathrm{\Sigma }_{\mathrm{vir}}`$, the surface density determined from the virial theorem, we have $`\mathrm{\Sigma }=\mathrm{\Sigma }_{\mathrm{vir}}/\alpha _{\mathrm{vir}}`$.

Regions of high-mass star formation studied by Plume et al. (1997, hereafter P97) are characterized by virial masses $`M_{\mathrm{vir}}3800\mathrm{M}_{}`$ and radii $`0.5\mathrm{pc}`$. The corresponding mean column density is $`\mathrm{\Sigma }_{\mathrm{vir}}1\mathrm{g}\mathrm{cm}^2`$; the corresponding visual extinction is $`A_V=(N_\mathrm{H}/2\times 10^{21}`$ cm$`{}_{}{}^{2})`$ mag $`=(\mathrm{\Sigma }_{\mathrm{vir}}\times 214)`$ mag. These column densities are far greater than those of giant molecular clouds (GMCs) (0.035 g cm-2—Solomon et al. 1987) or of regions of low-mass star formation (the average $`\mathrm{\Sigma }`$ in the C18O cores in Taurus is 0.032 g cm-2—Onishi et al. 1996). On the other hand, the central stellar surface density in the Orion Nebula Cluster is about 1 g cm-2 (Hillenbrand & Hartmann 1999). The mean density in the Plume et al. clumps is $`n_\mathrm{H}2\times 10^5`$ cm-3. As discussed by MT02, the virial parameter $`\alpha _{\mathrm{vir}}1.31.4`$ for GMCs, whereas it is quite close to unity for cores that are actively forming stars. Although there is no direct determination of $`\alpha _{\mathrm{vir}}`$ for massive star-forming clumps, P97 regard the virial mass as the most accurate, and we adopt $`\alpha _{\mathrm{vir}}=1`$ for our numerical estimates.

MT02 discuss the estimates for the other parameters, $`\varphi _B2.8`$, $`f_{\mathrm{gas}}2/3`$ and $`\varphi _P2`$ needed to determine $`P_s`$. Thus from equation (10)

$$P_s=2.28\times 10^8\left(f_{\mathrm{gas}}\varphi _P\varphi _B\alpha _{\mathrm{vir}}\right)\mathrm{\Sigma }^2\mathrm{K}\mathrm{cm}^38.5\times 10^8\mathrm{\Sigma }^2\mathrm{K}\mathrm{cm}^3.$$
(11)
With $`ϵ_{\mathrm{core}}=0.5`$ and $`H_0=1`$, the corresponding protostellar accretion rate is

$`\dot{m}_{}`$ $`=`$ $`4.75\times 10^4ϵ_{\mathrm{core}}^{1/4}{\displaystyle \frac{(f_{\mathrm{gas}}\varphi _P\alpha _{\mathrm{vir}}\varphi _B)^{3/8}}{(1+H_0)^{1/2}}}\left({\displaystyle \frac{m_f}{30\mathrm{M}_{}}}\right)^{3/4}\mathrm{\Sigma }^{3/4}\left({\displaystyle \frac{m_{}}{m_f}}\right)^{0.5}\mathrm{M}_{}\mathrm{yr}^1`$ (12)
$``$ $`4.6\times 10^4\left({\displaystyle \frac{m_f}{30\mathrm{M}_{}}}\right)^{3/4}\mathrm{\Sigma }^{3/4}\left({\displaystyle \frac{m_{}}{m_f}}\right)^{0.5}\mathrm{M}_{}\mathrm{yr}^1,`$
and the star formation time is

$`t_f`$ $`=`$ $`{\displaystyle \frac{1.26\times 10^5}{ϵ_{\mathrm{core}}^{1/4}(f_{\mathrm{gas}}\varphi _P\alpha _{\mathrm{vir}})^{3/8}}}\left[{\displaystyle \frac{(1+H_0)^{1/2}}{\varphi _B^{3/8}}}\right]\left({\displaystyle \frac{m_f}{30\mathrm{M}_{}}}\right)^{1/4}\mathrm{\Sigma }^{3/4}\mathrm{yr}`$ (13)
$``$ $`1.29\times 10^5\left({\displaystyle \frac{m_f}{30\mathrm{M}_{}}}\right)^{1/4}\mathrm{\Sigma }^{3/4}\mathrm{yr}.`$ (14)
For example, with the fiducial values of our parameters, a $`100\mathrm{M}_{}`$ star forming in a clump with $`\mathrm{\Sigma }=1`$ g cm-2 has a final accretion rate of $`1.1\times 10^3\mathrm{M}_{}\mathrm{yr}^1`$ and a star-formation time of $`1.75\times 10^5\mathrm{yr}`$. Our result for the accretion rate is somewhat lower than that of Osorio et al. (1999), who found accretion rates by matching the infrared spectra of embedded protostars. Our result for the star formation time is comparable to that of Behrend & Maeder (2001), who performed an analysis based on protostellar outflows, although our result for the time evolution of the accretion rate is quite different.

The above solutions apply to purely nonthermal cores. We have also allowed for an additional thermal ($`T100\mathrm{K}`$) component of pressure support (MT02). For typical high pressure star-forming regions cores with $`M>1\mathrm{M}_{}`$ are dominated by the nonthermal component. This implies that cores are supersonically turbulent and therefore clumpy, which means we expect the accretion rate to exhibit large fluctuations about the mean value. Radiation pressure can reduce the accretion flow to the most massive stars from the values predicted above (Wolfire & Cassinelli 1987; Jijina & Adams 1996; Tan & McKee 2002b).

## 3. Protostellar Evolution, Luminosities and Outflows

The properties of accreting protostars depend on their accretion rates (e.g., Stahler, Shu & Taam 1980; Stahler 1988; Palla & Stahler 1992; Nakano et al. 2000). In particular, the accretion luminosity $`L_{\mathrm{acc}}=f_{\mathrm{acc}}Gm_{}\dot{m}_{}/r_{}`$, where $`f_{\mathrm{acc}}`$ is a factor of order unity accounting for energy radiated by an accretion disk or used to drive protostellar outflows, and the stellar radius $`r_{}`$ may depend sensitively on $`\dot{m}_{}`$. Massive stars join the main sequence during their accretion phase at a mass that also depends on the accretion rate (Palla & Stahler 1992). The intensity of protostellar outflows, at least from low-mass stars, is thought to depend on the accretion rate and the Keplerian velocity, $`v_K`$, at the stellar surface (e.g. Shu et al. 1994). To predict protostellar properties and thus compare theory with observation, we require a model of protostellar evolution.

As previous investigations have mostly considered constant accretion rates, we have developed a simple model for protostellar evolution based on that of Nakano et al. (2000), which allows easy implementation of time dependent accretion. The model accounts for the total energy of the protostar as it accretes and dissociates matter, expands and contracts, and, if the central temperature $`T_c>10^6\mathrm{K}`$, burns deuterium. We have modified Nakano et al.’s model to include additional processes, such as deuterium shell burning, and calibrated these modifications against the more detailed calculations of Stahler (1988) and Palla & Stahler (1992). Protostellar evolution is followed from $`m_{}=0.1\mathrm{M}_{}`$ to the end of accretion ($`m_{}=m_f`$). Massive protostars reach the main sequence (Schaller et al. 1992) before they have finished accreting. The evolution of protostellar radius as a star grows in mass is shown in Figure 1 for several different final stellar masses, all forming in a clump with $`\mathrm{\Sigma }=1\mathrm{g}\mathrm{cm}^2`$ and $`T=100\mathrm{K}`$. Also shown is the evolution in bolometric luminosities (including the contribution from an accretion disk), and our use of these models to predict the masses and accretion rates of protostars thought to be illuminating several nearby hot molecular cores (HMCs). MT02 consider the effect of varying $`k_\rho `$ and $`\mathrm{\Sigma }`$ on these estimates and also compare results with Osorio et al. (1999), who have modeled the same sources.

Powerful outflows are ubiquitous from regions of both low and high-mass star formation (e.g., Richer et al. 2000). Models of magnetically driven outflows from a disk (Blandford & Payne 1982) and from the inner edge of a disk (Shu et al. 1994) have been presented. A common feature of these models is the production of a bipolar outflow with momentum distribution $`p_w(\mathrm{sin}\theta )^2`$ for $`\theta >\theta _0`$, where $`\theta `$ is measured from the outflow axis and $`\theta _010^2`$ (Matzner & McKee 1999). On scales large compared to the source (Matzner & McKee 1999)

$$\frac{d\dot{p}_w}{d\mathrm{\Omega }}=\frac{\dot{p}_w}{4\pi \mathrm{ln}(2/\theta _0)(1+\theta _0^2\mathrm{cos}^2\theta )}.$$
(15)
We parameterize the total momentum flux escaping in the outflow via

$$\dot{p}_w=f_{w,\mathrm{esc}}\dot{m}_wv_w=f_{w,\mathrm{esc}}f_w\dot{m}_{}v_w=\varphi _w\dot{m}_{}v_K,$$
(16)
where $`f_{w,\mathrm{esc}}`$ is the fraction of outflow momentum that escapes from the core, $`f_w=\dot{m}_w/\dot{m}_{}`$, $`\varphi _w=f_{w,\mathrm{esc}}f_wv_w/v_K`$ and $`v_K`$ is the Keplerian velocity at the equatorial radius of the star. Najita & Shu (1994) considered the acceleration of winds from an accreting $`0.5\mathrm{M}_{}`$ protostar. For several different X-wind model boundary conditions their results gave quite a large variation in $`f_w0.10.8`$, but only a small spread in $`f_wv_w/v_K0.6`$ within $`30`$%. Najita & Shu (1994) also showed that $`v_w`$ is approximately independent of $`\theta `$ so that $`\dot{m}_w`$ has the same distribution as $`\dot{p}_w`$, and thus $`f_{\mathrm{w},\mathrm{esc}}`$ applies to wind mass as well as momentum. If all the protostellar wind material escapes in two cones of solid angle $`\mathrm{\Omega }_{w,\mathrm{esc}}`$ (opening angle $`\theta _{w,\mathrm{esc}}`$), and if core material originally inside these cones is also swept up, then

$$ϵ_{\mathrm{core}}=\frac{12\mathrm{\Omega }_{w,\mathrm{esc}}/(4\pi )}{1+f_{w,\mathrm{esc}}f_w}.$$
(17)
For $`ϵ_{\mathrm{core}}=0.5`$ (Matzner & McKee 2000) and $`f_w=0.1,0.3,0.8`$ we solve equations (15) (with $`\dot{p}_w\dot{m}_w`$) and (17) for $`f_{w,\mathrm{esc}}`$ and $`\theta _{w,\mathrm{esc}}`$. We find $`f_{w,\mathrm{esc}}=0.89,0.86,0.79`$ and $`\theta _{w,\mathrm{esc}}=57^{},51^{},36^{}`$, respectively. The escape efficiencies are quite high because of the concentration of the wind near the outflow axis (eq. ). Assuming the wind is atomic, the swept-up gas is molecular and they escape together in a momentum conserving shell, the fraction of momentum leaving the core that is atomic in these three cases is 0.49, 0.62, 0.78, respectively.

We adopt $`f_w=0.3`$ and $`f_{w,\mathrm{esc}}=0.86`$ as our fiducial values so that $`\varphi _w0.5`$, since $`f_wv_w/v_K0.6`$ (Najita & Shu 1994). This estimate is broadly consistent with observations of outflows from low-mass stars (Bontemps et al. 1996). We apply the above model for low-mass protostellar outflows to stars of all masses as they accrete. This application is uncertain: for the X-wind model (Shu et al. 1988) to apply, the stellar surface field strength required by a $`10\mathrm{M}_{}`$ protostar with radius $`10\mathrm{R}_{}`$ accreting at a few $`\times 10^4\mathrm{M}_{}\mathrm{yr}^1`$ is $`10^3\mathrm{G}`$. The momentum injection rates and total momenta from protostars forming in a $`\mathrm{\Sigma }=1\mathrm{g}\mathrm{cm}^2`$ environment are shown in Figure 2. A fit to these results is

$$p_wp_1m_f^{\alpha _p}=68\left(\frac{\varphi _w}{0.5}\right)\left(\frac{m_f}{\mathrm{M}_{}}\right)^{1.4}\mathrm{M}_{}\mathrm{km}\mathrm{s}^1.$$
(18)
This relation depends only weakly on $`\mathrm{\Sigma }`$ (TM02).

## 4. Probing the Rate of Star Formation in Clusters with Outflows

We model the formation of star clusters with $`𝒩_f=M_f/\overline{m}_f`$ stars, where $`M_f`$ is the total stellar mass and $`\overline{m}_f`$ is the mean individual stellar mass. Stars are drawn from a power-law initial mass function (IMF) truncated above $`m_u`$ and below $`m_l`$ and of the form $`dF/d\mathrm{ln}m_{}=A_{}m_{}^\alpha _{}`$, where $`F`$ is the fraction of all stars that have a mass greater than $`m_{}`$. The Salpeter (1955) IMF has $`\alpha _{}=1.35`$, $`A_{}=0.0603`$ and $`\overline{m}_{}=0.353\mathrm{M}_{}`$ for $`m_l=0.1\mathrm{M}_{}`$ and $`m_u=120\mathrm{M}_{}`$. We adopt this as our fiducial IMF. While observational determinations of the IMF indicate deviations from this simple power-law form at lower masses, these have little effect on most of the feedback processes operating in a protocluster, which are controlled by the massive stars. The total protostellar wind momentum, $`p_{\mathrm{w},\mathrm{tot}}`$, from the forming star cluster does receive important, approximately equal, contributions from all decades of the IMF, since $`\alpha _{}\alpha _p`$. However, even in this case, the variation caused by uncertainties in the low-mass IMF is minor. The low-mass stars do dominate the total stellar mass, and this may be parameterized via $`\mu _h`$, the total stellar mass per high-mass ($`m_{}>8\mathrm{M}_{}`$) star, most of which end their life in a core-collapse supernova explosion. For our adopted IMF, $`\mu _h=134\mathrm{M}_{}`$, which is close to the value for the Miller & Scalo (1979) IMF.

We use Monte Carlo simulations to find the median values of the post-accretion luminosity, $`L_f`$, and total protostellar outflow momentum released during formation of a cluster of mass $`M_f`$. Typical lower-mass clusters fail to fully sample the IMF. For $`M_f/\mathrm{M}_{}>1600`$, $`L_f/M_f630\mathrm{L}_{}\mathrm{M}_{}^1`$, while for $`100<M_f/\mathrm{M}_{}<1600`$, $`L_f/M_f436(M_f/1000\mathrm{M}_{})^{0.8}\mathrm{L}_{}\mathrm{M}_{}^1`$. For the outflows, with the fiducial values of the parameters $`\mathrm{\Sigma }=1\mathrm{g}\mathrm{cm}^2`$, $`T=100\mathrm{K}`$, $`k_\rho =1.5`$ and $`\varphi _w=0.5`$, we find

$$\frac{p_{w,\mathrm{tot}}}{M_f}87\left(\frac{\varphi _w}{0.5}\right)\left(\frac{M_f}{1000\mathrm{M}_{}}\right)^{0.14}\mathrm{km}\mathrm{s}^1(M_f<1000\mathrm{M}_{}).$$
(19)
For $`M_f>1000\mathrm{M}_{}`$, $`p_{w,\mathrm{tot}}/M_f87(\varphi _w/0.5)\mathrm{km}\mathrm{s}^1`$. The loss of momentum in outflow-outflow interactions would lower these estimates. However, because of the collimation of momentum in jets (eq.), we expect this effect to be relatively minor, particularly for smaller clusters with fewer simultaneous independent outflows.

Unless triggered by an external influence such as a supernova blast wave, the overall rate at which stars form within a clump is not expected to exceed $`M_f/\overline{t}_{\mathrm{ff}}`$, where $`\overline{t}_{\mathrm{ff}}=(3\pi /32G\overline{\rho })^{1/2}`$ is the free-fall time of the clump evaluated at the mean density (stars + gas) at a typical stage in the evolution (i.e. when half the stars have formed). We write the star formation rate in terms of the numerical parameter $`\eta 1`$ as

$$\dot{M}_{}=\frac{M_f}{\eta \overline{t}_{\mathrm{ff}}}=\frac{2M_{,\mathrm{typ}}}{\eta \overline{t}_{\mathrm{ff}}}=\frac{2M_{\mathrm{cl},\mathrm{typ}}}{3\eta \overline{t}_{\mathrm{ff}}},$$
(20)
where $`M_{,\mathrm{typ}}=0.5M_f`$ is the typical current mass of stars, $`M_{\mathrm{cl},\mathrm{typ}}`$ is the typical current total mass of the clump, and, for $`ϵ_{\mathrm{core}}=0.5`$, $`M_{,\mathrm{typ}}=M_{\mathrm{cl},\mathrm{typ}}/3`$, if the protostellar outflows are efficiently ejected from the clump and most of the clump mass is in cores (see TM02). The star formation rate is then

$$\dot{M}_{}=9.7\times 10^3\left(\frac{1}{\eta }\right)\left(\frac{M_{\mathrm{cl},\mathrm{typ}}}{1000\mathrm{M}_{}}\right)^{3/4}\mathrm{\Sigma }^{3/4}\mathrm{M}_{}\mathrm{yr}^1$$
(21)
and the cluster formation time is

$$t_{\mathrm{cluster}}=\eta \overline{t}_{\mathrm{ff}}=6.86\times 10^4\eta \left(\frac{M_{\mathrm{cl},\mathrm{typ}}}{1000\mathrm{M}_{}}\right)^{1/4}\mathrm{\Sigma }^{3/4}\mathrm{yr}.$$
(22)
We wish to estimate the cluster star formation rate, and thus $`\eta `$, from observations. One method is to relate the summed total of many individual stellar accretion events to the production of a protocluster wind, which is the superposition of many individual bipolar protostellar outflows. From equations (19) and (22) the expected momentum injection rate from a forming star cluster is

$`\dot{p}_{w,\mathrm{tot}}`$ $`=`$ $`{\displaystyle \frac{p_{w,\mathrm{tot}}}{t_{\mathrm{cluster}}}}`$ (23)
$`=`$ $`\{\begin{array}{cc}\frac{0.80}{\eta }\left(\frac{\varphi _w}{0.5}\right)\left(\frac{M_{\mathrm{cl},\mathrm{typ}}}{10^3\mathrm{M}_{}}\right)^{0.89}\mathrm{\Sigma }^{3/4}\mathrm{M}_{}\mathrm{km}\mathrm{s}^1\mathrm{yr}^1\hfill & \frac{M_{\mathrm{cl},\mathrm{typ}}}{\mathrm{M}_{}}<1500\text{,}\hfill \\ \frac{0.85}{\eta }\left(\frac{\varphi _w}{0.5}\right)\left(\frac{M_{\mathrm{cl},\mathrm{typ}}}{10^3\mathrm{M}_{}}\right)^{3/4}\mathrm{\Sigma }^{3/4}\mathrm{M}_{}\mathrm{km}\mathrm{s}^1\mathrm{yr}^1\hfill & \frac{M_{\mathrm{cl},\mathrm{typ}}}{\mathrm{M}_{}}>1500\text{.}\hfill \end{array}`$
We compare equation (23) against observations of outflows from protoclusters (Choi et al. 1993, hereafter C93). The observational determination of $`\dot{p}_{w,\mathrm{tot}}`$ involves measuring the total outflow momentum over some region and identifying a characteristic timescale of the flow over the scale of this region. The CO observations of C93 cover regions with sizes similar to $`R_{\mathrm{cl}}`$ for which virial mass estimates have been made from CS transitions, tracing dense gas. The spatial extent of the outflow observations is limited by sensitivity and the flows almost certainly continue beyond $`R_{\mathrm{cl}}`$. This is also evident from the fact that the flow timescales are much shorter than $`\overline{t}_{\mathrm{ff}}`$. These considerations allow us to neglect the effects of evolution in determining an estimate of the instantaneous star formation rate from the observed flows. Though poorly resolved, the morphologies of the flows appear quasi-spherical, as would be expected if there were multiple driving sources.

From CO line strengths, C93 derive the mass of the flow, correcting for optical depth and excitation temperature effects and assuming standard abundance ratios. The line of sight momentum is then summed over velocity channels. We have further corrected these estimates for projection effects (a factor of two increase for a spherical flow). The total momentum (and mass) of the flow may be greater by the amount contained in atomic gas. For example we have seen in our fiducial model for protostellar winds from an individual core that 62% of the total momentum is in this form (i.e. resides in gas that was processed through the accretion disk). As this flow sweeps up gas beyond the core, the atomic fraction is reduced. TM02 estimate that about 25% of the total momentum is in atomic form, so to recover the total momentum we therefore increase the observed value by 33%. This question may be examined observationally, though current surveys have placed only upper limits on the fraction of atomic momentum in the flows (Choi et al. 1994).

Given estimates of mass and momentum, we determine a characteristic flow velocity, $`\overline{v}_p`$, and thus an outflow timescale $`t_p=R_p/\overline{v}_p`$ (here we use the data uncorrected for projection effects, which should approximately cancel in this ratio). The value of $`\dot{p}_{w,\mathrm{tot}}`$ is then calculated, and with knowledge of the clump’s mass and surface density, used to estimate $`\eta `$. The properties of our sample of protoclusters and their outflows are listed in Tables 1 and 2.

Excluding the cases of Cep A and NGC6334 I(N) (below), we find the logarithmic mean value of $`\eta `$ is 10.4. The standard deviation is about a factor of two, consistent with the observational uncertainties. This gives a cluster formation time (eq.) of $`7.1\times 10^5(M_{\mathrm{cl},\mathrm{typ}}/1000\mathrm{M}_{})^{1/4}\mathrm{\Sigma }^{3/4}\mathrm{yr}`$, which is consistent with estimates of $`<10^6\mathrm{yr}`$ from the spread in pre-main-sequence ages in the Orion Nebula Cluster (Palla & Stahler 1999), which has $`M_{\mathrm{cl}}4500\mathrm{M}_{}`$ (Hillenbrand & Hartmann 1998). We note that the effect of wind-wind interactions, which we have ignored in our analysis, would shorten this estimate by the fraction of wind momentum that is dissipated.

Cep A is the most extreme of the C93 protoclusters in terms of its surface density. The estimates of virial mass and surface density are made uncertain by the fact that 2D maps of the dense gas show moderate elongation (Y. Shirley, private communication). Furthermore, the outflow observations do not symmetrically probe the central region of Cep A (M. Choi, private communication) and so the total $`\dot{p}_{w,\mathrm{tot}}`$ may be underestimated, thus artificially raising $`\eta `$. For these reasons we have excluded Cep A from our statistical sample. NGC6334 I(N)’s outflow was observed by Megeath & Tieftrunk (1999) and the outflow parameters determined using different techniques to those of C93. This source is thought to be in an extremely early stage of evolution; for example the observed outflow extends over only a small fraction of the clump, which has a relatively low luminosity. Again we exclude this source from the statistical sample, but consider it in the context of an example of the earliest stage of cluster formation.

A common observational diagram used in the study of outflows from protostars and protoclusters plots $`\dot{p}_{w,\mathrm{tot}}`$ versus $`L_{\mathrm{bol}}`$. We consider a simple description of the evolution of protoclusters in this diagram. Assuming a constant star formation rate (eq.) implies that, averaged over a population of clusters, $`\dot{p}_{w,\mathrm{tot}}`$ is a constant over most of the evolution. In the earliest stages the luminosity is small and dominated by that resulting from accretion. The minimum luminosity is $`G\overline{m}_{}\dot{M}_{}/r_{}(\overline{m}_{})`$, where $`\overline{m}_{}<\overline{m}_f`$ is the mean protostellar mass. We have $`\dot{M}_{}=\dot{p}_{w,\mathrm{tot}}/(\varphi _w\overline{v}_K)`$, where $`\overline{v}_K`$ is the Keplerian velocity at the stellar surface of an average protostar. The mean radius is about $`3\mathrm{R}_{}`$ (Figure 1) and thus

$$L_{\mathrm{min}}=4.3\times 10^4\left(\frac{0.5}{\varphi _w}\right)\left(\frac{\overline{m}_{}}{0.3\mathrm{M}_{}}\frac{3\mathrm{R}_{}}{r_{}(\overline{m}_{})}\right)^{1/2}\left(\frac{\dot{p}_{w,\mathrm{tot}}}{\mathrm{M}_{}\mathrm{km}\mathrm{s}^1\mathrm{yr}^1}\right)\mathrm{L}_{}.$$
(24)
Observed protoclusters, selected for showing signs of massive star formation, will tend to be over-luminous compared to equation (24). Those in the earliest stages of evolution should lie systematically closer to the birthline.

The maximum luminosity is approximately the post-accretion luminosity, $`L_f`$, for all but the smallest clusters (where accretion can make a significant contribution). We have determined $`L_f(M_f)`$ from Monte Carlo simulations (above). Relating $`M_f=2M_{\mathrm{cl},\mathrm{typ}}/3`$ for a typical forming cluster (TM02), we then eliminate $`M_{\mathrm{cl},\mathrm{typ}}`$ in equation (23). The resulting expression defines the “cluster birthline” (the 3 cases correspond to values of $`M_f/\mathrm{M}_{}`$ of $`10^210^3`$, $`10^31.6\times 10^3`$ and $`>1.6\times 10^3`$, respectively):

$$L_f=\{\begin{array}{cc}3.5\times 10^7\varphi _L^{2.0}\mathrm{\Sigma }^{1.5}\left(\frac{\dot{p}_{w,\mathrm{tot}}}{\mathrm{M}_{}\mathrm{km}\mathrm{s}^1\mathrm{yr}^1}\right)^{2.0}\mathrm{L}_{}\hfill & 6.9\times 10^3<\frac{L_f}{\mathrm{L}_{}}<4.4\times 10^5\hfill \\ 7.9\times 10^7\varphi _L^{2.4}\mathrm{\Sigma }^{1.8}\left(\frac{\dot{p}_{w,\mathrm{tot}}}{\mathrm{M}_{}\mathrm{km}\mathrm{s}^1\mathrm{yr}^1}\right)^{2.4}\mathrm{L}_{}\hfill & 4.4\times 10^5<\frac{L_f}{\mathrm{L}_{}}<1.0\times 10^6\hfill \\ 1.1\times 10^7\varphi _L^{4/3}\mathrm{\Sigma }^1\left(\frac{\dot{p}_{w,\mathrm{tot}}}{\mathrm{M}_{}\mathrm{km}\mathrm{s}^1\mathrm{yr}^1}\right)^{4/3}\mathrm{L}_{}\hfill & 1.0\times 10^6<\frac{L_f}{\mathrm{L}_{}}\hfill \end{array}$$
(25)
where $`\varphi _L=(\varphi _w/0.5)(10/\eta )`$.

The protocluster and cluster birthlines are shown in Figure 3. In Figure 3a we also plot the results of Monte Carlo simulations of many clusters for two cases: $`M_f=200,1000\mathrm{M}_{}`$. For each cluster, stars are drawn randomly from the IMF and stochastically in time and then followed through our protostellar evolution model for pressures corresponding to $`\mathrm{\Sigma }=1\mathrm{g}\mathrm{cm}^2`$. The central/outer lines show the median/68-percentile values of each cluster distribution. Points mark equal time intervals in the evolution. In Figure 3b we plot the observed protoclusters, giving $`\dot{p}_{w,\mathrm{tot}}`$ a factor of four uncertainty. Sources in an early stage of formation (NGC2071 and NGC6334 I(N)) are close to the protocluster birthline, while those in a late stage (W3($`\mathrm{H}_2\mathrm{O}`$)—as evidenced by its multiple H II regions) are close to or below the cluster birthline.

## 5. Conclusions

We have presented a theoretical model for high-mass star formation that determines the accretion rate in terms of the final and instantaneous masses of the star, the ambient pressure surrounding the gas core, and the polytropic index of pressure support. We related the pressure due to self-gravity in the centers of typical Galactic star-forming clumps to their surface density; $`\overline{P}/k`$ approaches $`10^9\mathrm{K}\mathrm{cm}^3`$ for typical values of $`\mathrm{\Sigma }1\mathrm{g}\mathrm{cm}^2`$. These pressures allow a $`100\mathrm{M}_{}`$ star to form in $`10^5\mathrm{yr}`$ with a final accretion rate $`10^3\mathrm{M}_{}\mathrm{yr}^1`$. Modeling protostellar evolution, we predicted the properties of several nearby massive protostars. Under the hypothesis that the outflows from massive protostars are driven magneto-centrifugally in the same manner as those from low-mass protostars, we found that each decade of the IMF contributes approximately equally to the total outflow momentum. This helps to explain the generally poor overall collimation of outflows from high-mass star-forming regions. We utilized these models and observations of outflow intensities from several nearby protoclusters to estimate their star formation rates, finding that these clusters are typically forming over at least several free-fall times. We have proposed theoretical birthlines in the $`\dot{p}_{w,\mathrm{tot}}`$ versus $`L_{\mathrm{bol}}`$ diagram to delimit the phase of cluster formation.

### Acknowledgments.

We thank Henrik Beuther, Minho Choi, Neal Evans, Ralph Pudritz, Yancy Shirley, Steve Stahler and Malcolm Walmsley for helpful discussions. Our research is supported by NSF grants AST-9530480 and AST-0098365, by a NASA grant supporting the Center for Star Formation Studies, and (for JCT) by a Spitzer-Cotsen fellowship from Princeton University.

## References

Behrend, R., & Maeder, A. 2001, A&A, 373, 190

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 1997, MNRAS,285,201

Bonnell, I. A., Bate, M. R., & Zinnecker, H. 1998, MNRAS, 298, 93

Bontemps, S., Andre, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858

Caselli, P., & Myers, P. C. 1995, ApJ, 446, 665

Choi, M., Evans, N. J. (II), & Jaffe, D. T. 1993, ApJ, 417, 624 (C93)

Choi, M., Evans, N. J. (II), Jaffe, D. T., & Walker, C. K. 1994, ApJ, 435, 734

Hillenbrand, L.A., & Hartmann, L.W. 1998, ApJ, 492, 540

Jijina, J., & Adams, F. C. 1996, ApJ, 462, 874

Li, Z.-Y., & Shu, F. H. 1997, ApJ, 475, 237

Matzner, C. D., & McKee, C. F. 1999, ApJ, 526, L109

Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364

McKee, C. F., & Tan, J. C. 2002, in preparation (MT02)

McLaughlin, D. E., & Pudritz, R. E. 1996, ApJ, 469, 194

McLaughlin, D. E., & Pudritz, R. E. 1997, ApJ, 476, 750 (MP97)

Megeath, S. T., & Tieftrunk, A. R. 1999, ApJ, 526, L113

Miller, G. E., & Scalo, J. M. 1979, ApJS, 41, 513

Najita, J. R., & Shu, F. H. 1994, ApJ, 429, 808

Nakano, T., Hasegawa, T., Morino, J.-I., & Yamashita, T. 2000, ApJ, 534, 976

Onishi, T., Mizuno, A., Kawamura, A., Ogawa, H., & Fukui, Y. 1996, ApJ,465,815

Osorio, M., Lizano, S., & D’Alessio, P. 1999, ApJ, 525, 808

Palla, F., & Stahler, S. W. 1992, ApJ, 392, 667

Palla, F., & Stahler, S. W. 1999, ApJ, 525, 722

Plume, R., Jaffe, D. T., Evans, N. J. (II), Martin-Pintado, J., & Gomez-Gonzalez, J. 1997, ApJ, 476, 730 (P97)

Richer, J., Shepherd, D., Cabrit, S., Bachiller, R., & Churchwell, E. 2000, in Protostars & Planets IV, eds. V. Mannings, A. P. Boss, & S. S. Russell (Tucson: The University of Arizona Press), 867

Salpeter, E. E. 1955, ApJ, 121, 161

Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269

Shu, F. H. 1977, ApJ, 214, 488

Shu, F. H., Lizano, S., Ruden, S. P., & Najita, J. 1988, ApJ, 328, L19

Shu, F. H., Tremaine, S., Adams, F. C., & Ruden, S. 1990, ApJ, 358, 495

Shu, F. H., Najita, J., Ostriker, E., Wilkin, F., Ruden, S., & Lizano, S. 1994, ApJ, 429, 781

Solomon, P.M., Rivolo, A.R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730

Stahler, S. W. 1988, ApJ, 332, 804

Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637

Tan, J. C., & McKee, C. F. 2002a, in preparation (TM02)

Tan, J. C., & McKee, C. F. 2002b, in preparation

van der Tak, F. F. S., van Dishoeck, E. F., Evans, N. J. (II), & Blake, G. A. 2000, ApJ, 537, 283

Wolfire, M. G., & Cassinelli, J. 1987, ApJ, 319, 850