text:
# The impact of mass loss on star cluster formation. I. Analytic results

## 1 Introduction

Most if not all stars are formed in clusters or associations, and therefore this is the predominant mode of star formation contributing to the Galactic-field population (see reviews in Grebel & Brandner 2002). To understand how such aggregates form and evolve remains a severe challenge however, since no theory exists yet that accounts in detail for the (presumably) simpler case of the formation of individual stars. On the observational front, surveys of star-forming regions suggest that the mass-fraction of gas used up at birth does not exceed 30% of the total when averaged over typical cluster scales (Lada 1999). This low star-formation efficiency (SFE) implies that young stars or clusters of stars should be embedded in gas (see Lada & Lada 1991). Yet populous clusters with ages $`\begin{array}{c}<\hfill \\ \hfill \end{array}`$ a few Myrs (e.g. the Orion Nebula Cluster, R136 in 30 Doradus) are already void of gas. The mass-fraction of gas left after the epoch of star formation depends crucially on the star formation history of the cluster. The formation of proto-stars by hydrodynamical collapse points to the rapid formation of opaque stellar cores followed by time-dependent accretion, with the more massive cores accreting at higher rates (e.g. Boily & Lynden-Bell 1995). Stellar masses are then accrued over a period of typically $`10^510^6`$ years (Low & Lynden-Bell 1976; Shu 1977; Wuchterl & Klessen 2001). Observational evidence in support of rapid accretion as opposed to a slow star formation process sustained by ambipolar diffusion of magnetic field is discussed by Hartmann et al. (2001), who point out that the age spreads between cluster members are too small compared with that expected from ambipolar diffusion (see Jones et al. 2001 and references therein). Because opaque stellar cores will not all have the same mass (or indeed chemical composition), it is unlikely that stars will end their accretion phase at the same universal time. Consequently, the star cluster will not be entirely depleted of gas when the more massive stars reach the Hydrogen-burning phase.

OB stars provide an efficient mechanism by which gas is driven out. Stellar winds from OB stars yield a momentum flux $`8\times 10^3M_{\mathrm{}}\mathrm{km}\mathrm{sec}^1\mathrm{yr}^1`$ (Churchwell 1999). This is sufficient to unbound gas from an $`10^3M_{}`$, 1 pc radius cluster in $`10^5`$ years, which is comparable to its crossing time $`t_{\mathrm{cr}}=2R/\sigma 1/\sqrt{G\rho (0)}`$. Here $`\rho (r)`$ is the mass density and $`\sigma `$ the velocity dispersion. In addition, the ionising radiation heats the gas to $`10^4`$ K, causing an over-pressure and expansion at the sound velocity $`10`$ km/s. Furthermore, the star-formation timescale of $`10^6`$ yrs far exceeds a single cluster crossing time $`t_{\mathrm{cr}}`$: the magnitude of stellar velocities should therefore reflect the depth of the cluster gravitational potential since protostars would have time to explore the entire cluster volume. This suggests that a rich cluster is already close to virial equilibrium when OB winds set in.

The question whether or not the stars will re-establish a bound gas-free equilibrium has bearings on the spatial distribution and statistics of stellar populations in the Galaxy (Kroupa 2002, Grebel & Brandner 2002). A long-standing goal has been to understand in what way gas-ejection and the low SFE inferred from observations set constraints on the star cluster properties at birth (see e.g. Lada 1999). Hills (1980) has argued from this standpoint and using the virial theorem that if the SFE is less or equal to 50% the entire cluster will dissolve. His classic argument holds for clusters that have achieved virial equilibrium and which then undergo a sudden loss of mass. Theoretical and numerical work support this conclusion when gas evacuation proceeds rapidly (Tenorio-Tagle et al. 1986; Geyer & Burkert 2001), while the numerical N-body study of Lada et al. (1984) already found bound stellar cores regardless of the gas evacuation time-scales. This was also found to be the case in recent, larger-N studies (Goodwin 1997; Kroupa, Aarseth & Hurley 2001). Adams (2000) has recently pointed out that if the SFE peaks with the central cluster density, a small core of stars will always remain bound, even for an SFE well below the reference 50% obtained from the virial. This effectively solves the problem of forming bound clusters but at the expense of simplicity, by assuming a radial variation for the SFE, the nature of which is yet poorly constrained. Our goal is to show that one may solve for the bound fraction of stars analytically without introducing further complications.

The numerical N-body results and the analysis of Adams (2000) hint that Hills’ (1980) original argument is lacking an ingredient to account for bound clusters at low SFE. We note that the derivation based on the virial theorem is a statement about the global properties of the cluster with regards to the impact of the SFE. A refinement is therefore possible if we reformulate the problem from the point of view of the equilibrium stellar distribution function $`(`$DF), $`F(𝒙,𝒗)`$, as follows : First, select stars by binding energy $`(E)`$ to identify bound ($`E<0`$) and candidate escaping ($`E>0`$) stars. The DF may then be modified to take account of the cluster mass loss due to unbound gas and stars. As a second step, we setup an iterative scheme to identify bound stars self-consistently once gas and escapers are subtracted from the gravitational potential. This scheme allows us to solve analytically for the surviving fraction of bound stars. We show that the result is function of both SFE and the chosen DF by contrasting results obtained for polytropes, and in particular the Plummer sphere. The time-evolution of clusters, and applications to other DF’s, are deferred to a companion paper (Paper II).

## 2 Hills’ problem

A cluster of stars forms converting a fraction $`ϵ`$ of the gas into stars in the process. Stellar winds or a supernova event blow/s out the remainder on a time-scale $`\tau t_{\mathrm{cr}}`$. What fraction of the stars remain to form a bound cluster? Since the star formation epoch will have lasted as long as or longer than a cluster dynamical time $`t_{\mathrm{cr}}`$, the system as a whole is close to virial equilibrium before gas is expelled. In the following, we consider the limit $`\tau 0`$, which defines an initial value problem for an open cluster’s mass profile and velocity field.

Hills (1980) argued that equilibrium, self-gravitating stellar systems would dissolve if more than half the mass is lost instantaneously: stars then preserve their kinetic energy established under a deeper potential well, hence may escape if their binding energy $`E`$ becomes positive. We write the total energy

$$E=\frac{GM^2}{R}+\frac{1}{2}Mv^2<0,$$
(1)
where $`M,R`$ are the mass and gravitational radius of a spherical distribution of mean square velocity $`v^2`$. Before gas-expulsion the total mass in gas and stars is $`M=M_{\mathrm{init}}=M_{\mathrm{gas}}+M_{}`$, while after the gas is expelled $`M=M_{}`$ but the velocity dispersion, $`v^2=GM_{\mathrm{init}}/R`$, remains unchanged. The new, total system energy is now

$$E_{}=\frac{GM_{}^2}{R}+\frac{1}{2}M_{}v^2\frac{GM_{}^2}{2R_{}}$$
(2)
where the last equality follows from applying the virial theorem to the stellar system of mass $`M_{}`$ after it has reached a new equilibrium, of energy $`E_{}`$ but with radius $`R_{}>R`$. We may rearrange (2) to obtain

$$\frac{R_{}}{R}=\frac{1}{2}\frac{M_{}}{M_{}\frac{1}{2}M}=\frac{1}{2}\frac{MM_{\mathrm{gas}}}{\frac{1}{2}MM_{\mathrm{gas}}}.$$
(3)
The system radius $`R_{}\mathrm{}`$ when $`M_{\mathrm{gas}}=M/2`$ and hence the stellar system has zero or positive total binding energy when the mass is reduced by 50% or more $`(ϵ=M_{}/M1/2)`$.

Hills’ derivation applies to any mass distribution. It does not, however, say anything specific about the evolution of the mass profile except to relate the final virial radius $`R_{}`$ to the initial system radius $`R`$. Our intuition is that a fraction of the stars with larger velocities will escape before we reach the reference 50% gas mass loss and hence the residual stellar potential should not be evaluated from the global stellar mass fraction alone. In other words, the stellar velocity DF should determine which stellar orbit remains bound according to its binding energy: this is confirmed when considering orbits in a harmonic potential as shown in Appendix A. The example of harmonic motion will serve us later in interpreting results obtained for systems with realistic mass profiles and DF’s. Below we introduce a scheme to pin down more accurately the relation of bound mass to star formation efficiency, $`ϵ`$.

## 3 Distribution functions

### 3.1 Iterative scheme

Consider an isotropic mixture of stars and gas within a spherical volume. The stellar DF $`F(𝒓,𝒗)`$ is then a function of energy alone. As a simplifying assumption we take the SFE to be independent of position in the cluster. Hence the SFE parameter $`ϵ=`$ constant. Thus gas and stars are initially mixed in the same proportion throughout, however note that $`ϵ`$ does not fix the profiling of stellar masses with radius, rather the total mass of gas made into stars :

$$\frac{_\mathrm{i}(\mathrm{number}\mathrm{of}\mathrm{stars}\mathrm{in}\mathrm{mass}\mathrm{bin}i)\times (\mathrm{stellar}\mathrm{mass}m_i)}{m_{\mathrm{gas}}}=\mathrm{constant}$$

at any radius $`r`$. Both gas and stars are distributed according to the same DF. The total system mass is then

$$M=M_{}+M_{\mathrm{gas}}=ϵ^1M_{}=ϵ^1_VF(𝒓,𝒗)\text{d}^3𝒓\text{d}^3𝒗=ϵ^1_{E_c}^0F(E)𝒲(E)dE,$$
(4)
where $`|E_c|`$ is the maximum binding energy and $`𝒲`$ the density of states of energy $`E`$. The cluster potential $`\varphi (r)`$ follows from integrating twice Poisson’s equation

$$\varphi (r;M,R)=_{\mathrm{}}^r\frac{\mathrm{d}u}{u^2}_0^u4\pi G\rho (w)w^2dwy(r/R)\frac{GM}{R},$$
(5)
where $`y(x)`$ is a dimensionless function. Keeping $`R`$ constant while removing the gas instantly so the potential includes stars only, we have

$$\varphi (r;M,R)\varphi _{}=ϵ\varphi ,$$
(6)
and hence the fraction of stars at radius $`r`$ which now have positive energy is

$$1\lambda _e=\frac{{\displaystyle _{v_{e,}}^{v_e}}F(𝒓,𝒗)\text{d}^3𝒓v^2dv}{{\displaystyle _0^{v_e}}F(𝒓,𝒗)\text{d}^3𝒓v^2dv}=\frac{{\displaystyle _{v_{e,}}^{v_e}}f(v)v^2dv}{{\displaystyle _0^{v_e}}f(v)v^2dv},$$
(7)
with $`v_e`$ the (local) escape velocity computed from $`\varphi `$ (before gas-expulsion), and similarly for $`v_{e,}`$ obtained from $`\varphi _{}`$. For convenience of notation, we write the stellar velocity distribution function $`f(v)`$ defined at constant radius $`𝒓`$ as

$$f(𝒗)\text{d}^3𝒗F(𝒓,𝒗)\text{d}^3𝒓\text{d}^3𝒗.$$

For a given DF and $`ϵ`$, we may compute the quantity $`1\lambda _e`$ at all radii. Note that should the SFE be a function of radius, the re-normalisation (6) leading to $`\varphi _{}`$ would not apply, however (7) may still be computed if the stars’ potential is given and $`v_e`$ known from the initial gas + stars mixture. The key step is to adjust the potential $`\varphi _{}`$ itself, since potential and velocity field do not match any longer. This is normally computed numerically as an initial value problem, for instance using N-body integration which takes into account all star-star interactions and the redistribution of energy between them in the time-dependent potential. We take a different approach. If we thought that stars acquire only little extra energy during the time that they escape, then the fraction of stars remaining might be computed as follows.

Since the fraction (7) of positive-energy stars is known at all $`r`$, we re-compute the gravitational potential counting only stars with $`E0`$ at each radius. Neglecting dynamical evolution, the cluster radius $`R`$ and structure function $`y(r/R)`$ are unchanged and hence the potential can be re-computed by integration, from $`\mathrm{}`$, inwards. Once the new potential is known, a fraction of the remaining stars will again be unbound by virtue of the stars lost during the previous iteration. We therefore repeat the procedure until finally the cluster mass converges to a finite quantity, in which case no more stars escape and the original distribution function is depleted from all escapers in a self-consistent manner.

We consider the case of scale-free DF in detail, then turn to more physically-motivated cases in the following sections.

### 3.2 Example : scale free DF’s

The case of a power-law isotropic velocity distribution,

$$f(𝒗)\text{d}^3𝒗v^{\beta 2}\text{d}^3𝕧v^\beta \text{d}v,$$

where $`\beta `$ is constant, provides a tractable starting point. The DF is truncated at the local escape velocity and we take the potential $`\varphi (r)`$ to be a power of the radius: the system is therefore scale-free. The one-dimensional velocity dispersion obtained from Jeans’ equation which for power-law density and potential yields

$$\sigma ^2(r)=\frac{1}{\rho (r)}_r^{\mathrm{}}\rho (x)\varphi (x)\text{d}x=2K\varphi (r)=Kv_e^2$$
(8)
at each radius for some constant of proportionality, $`K`$. Substituting $`\rho (r)ϵ\rho (r)`$ in (5), we find the new potential from (6) and escape velocity

$$v_{e,}(r)=\sqrt{2\varphi _{}(r)}=ϵ^{1/2}v_e(r).$$
(9)
Equation (7) becomes (with $`uv/v_e`$)

$$1\lambda _e=\frac{{\displaystyle _{v_{e,}}^{v_e}}f(v)v^2dv}{{\displaystyle _0^{v_e}}f(v)v^2dv}=\frac{{\displaystyle _{ϵ^{1/2}}^1}f(u)u^2du}{{\displaystyle _0^1}f(u)u^2du}=1ϵ^{(\beta +1)/2},$$
(10)
where we took $`\beta >1`$. Thus $`\lambda _e=ϵ^{(\beta +1)/2}`$ is independent of radius. Repeating the procedure to take account of the positive-energy stars, we substitute $`ϵ\lambda _eϵ`$, etc, so that after $`j`$ iterations

$$\lambda _e[j]=(\lambda _e[j1]ϵ)^{(\beta +1)/2}$$

with $`\lambda _e[0]=1`$. Assuming convergence of the series with $`j\mathrm{}`$ we may solve for $`\lambda _e`$ to obtain

$$\lambda _e=ϵ^{{\displaystyle \frac{1+\beta }{1\beta }}}.$$
(11)
(See Boily & Kroupa 2002 for an alternative derivation.) In terms of the stellar mass fraction, the predicted mass of bound stars, $`M_{}^b`$, is then simply

$$M_{}^b=\lambda _eϵM=\lambda _eM_{}$$
(12)
(and similarly for the stellar density $`\rho _{}[r]`$ wrt $`\rho [r]`$). The eigenvalue $`\lambda _e`$ in (12) converges to a non-zero (positive) value as $`j\mathrm{}`$ only for $`\beta <1`$, since $`ϵ1`$. All DF’s with $`\beta >1`$ lead to cluster disruption, because the high-velocity range of the DF is too densely populated, leading to catastrophic stellar loss after the expulsion of any amount of gas.

Fig. 1 graphs the fraction of bound stars $`M_{}^b/M_{}`$ obtained from (12) as function of the SFE $`ϵ`$ for five values of the power index $`\beta `$. All curves meet at $`ϵ=1`$ and $`M_{}^b/M_{}=1`$. Note how a flat DF ($`\beta =0)`$ predicts a linear decline of bound star fraction.

The volume density $`\rho (r)`$ is recovered from integrating $`F(𝒓,𝒗)`$ over all velocities at constant radius $`r`$ :

$$\rho (r)\text{d}^3𝒓=_0^{v_e}F(𝒓,𝒗)\text{d}^3𝒓\text{d}^3𝒗=_0^{v_e}f(v)\mathrm{\hspace{0.17em}4}\pi v^2\text{d}v_0^{v_e}v^\beta \text{d}v=\frac{v_e^{\beta +1}}{\beta +1}=(2\varphi [r])^{(\beta +1)/2}(2\varphi )^n$$
(13)
where the last relation emphasises the similarity with polytropes of index $`n`$, to which we return. Since (12) was derived assuming a power-law potential, it follows from (13) that $`\rho (r)`$ is also a power of the radius. Power-law density profiles may be of application to young clusters such as NGC2282 (Horner et al. 1997). Letting $`\rho (r)r^\alpha `$, we deduce from (5) that $`\varphi (r)r^{2\alpha }`$; a centrally condensed cluster requires $`\alpha >0`$. The power indices $`\alpha `$ and $`\beta `$ are then related to each other using (13)

$$\alpha =2\frac{\beta +1}{\beta 1}$$

and hence $`\alpha >0`$ for $`\beta >1`$ or $`\beta <1`$. We recall that $`\beta >1`$ from (10), hence we must choose $`\beta >1`$ for a self-consistent solution. NGC2282 shows an axially-averaged surface density $`\mathrm{\Sigma }R^1`$ which can be obtained from projecting a spherical cluster of volume density $`\rho r^2`$, or $`\alpha =2`$ (singular isothermal sphere). The index $`\beta =(\alpha +2)/(\alpha 2)\mathrm{}`$ when $`\alpha =2`$ which is unrealistic. Therefore no self-consistent solution of this sort exists with $`\beta >1`$ applicable to NGC2282. The large value of $`\beta `$ obtained for this case is at best indicative of rapid dissolution for any $`ϵ<1`$.

The example of power-law DF’s suggests a strong dependence of the effective SFE when cluster disruption occurs, and the shape of the DF itself. Naive application to NGC2282 of a power-law DF based on the Jeans equations of equilibria leads to a nonphysical solution. In the following we apply our iterative scheme to other known DF’s.

## 4 Plummer model and the family of polytropes

We wish to extend our basic result (11) to a class of DF’s with properties similar to those of observed star clusters. The power-law relation of the scale-free DF’s suggests that we consider the family of polytropes defined as

$$F(E)=(E)^{n\frac{3}{2}}$$
(14)
where the power index $`n>1/2`$ may take non-integer values (BT+87). The case $`n=5`$ has been studied in detail by Plummer (1911). Then

$$F(E)=(E)^{7/2}=(|\varphi (r)|\frac{1}{2}v^2)^{7/2}.$$
(15)
The total system mass of a Plummer sphere is finite but infinite in extent; the density profile

$$\rho (r)=\frac{3M}{4\pi R_p^3}\left[1+r^2/R_p^2\right]^{5/2}\varphi ^5(r)=\left(2v^2[r]\right)^5$$
(16)
where $`M`$ is the total system mass, $`R_p`$ the Plummer radius and both the velocity dispersion and density maximise at the centre (see e.g. Spitzer 1987). At large distances, the density profile $`\rho (r)r^5`$ decreases like a power of the radius. Remembering (8), the Plummer velocity DF $`f(v)(1[v/v_e]^2)^{7/2}`$. Inserting this into (7) and with

$$p(ϵ)1\frac{1210}{105}ϵ+\frac{2104}{105}ϵ^2\frac{1488}{105}ϵ^3+\frac{384}{105}ϵ^4$$

we find

$$1\lambda _e=\frac{{\displaystyle _{ϵ^{1/2}}^1}(1u^2)^{7/2}u^2du}{{\displaystyle _0^1}(1u^2)^{7/2}u^2du}=1\frac{2}{\pi }\left(\mathrm{sin}^1ϵ^{1/2}p(ϵ)ϵ^{1/2}\sqrt{1ϵ}\right).$$
(17)
We note that, as for the case of power-laws, the solution is independent of radius. Therefore the potential $`\varphi _{}`$ and escape velocity $`v_e`$ follow from (6). The fraction of stars labeled for escape are removed from the potential and the integral re-evaluated with the substitution of $`ϵ^{1/2}`$ by $`(ϵ\lambda _e)^{1/2}`$ as the lower bound of integration in (17). This follows from the same transformation leading to (9). We may repeat the procedure until $`\delta (1\lambda _e)0`$ on successive iterations so no additional stars are lost. The converged fraction $`\lambda _e`$ is then solution of

$$\lambda _e=\frac{2}{\pi }(\mathrm{sin}^1(\lambda _eϵ)^{1/2}p(\lambda _eϵ)(\lambda _eϵ)^{1/2}\sqrt{1\lambda _eϵ}).$$
(18)
This may be expressed in parameterised form. Defining

$$x\lambda _eϵ;x[0,1],$$
(19)
we have the solution $`\lambda _e(x)`$ from (18). The full solution follows from solving for $`ϵ`$ as a function of $`x`$,

$$ϵ(x)=\frac{x}{\lambda _e(x)}$$
(20)
which is known for all $`x`$ in the allowed range. The relation $`\lambda _e(x)`$ to $`ϵ(x,\lambda _e)`$ saves us from iterating on (18), which otherwise yields the same result. In Appendix B we extend the solution (18) to configurations out of dynamical equilibrium.

### 4.1 Analytic results for Plummer spheres

From the definition (19) we must have $`ϵ=\lambda _e=1`$ when $`x=1`$. Inserting (18) into (20) we find

$$\underset{x0}{lim}ϵ(x)x^{1/2}\mathrm{}$$

while $`\lambda _e(x0)=0`$. Since $`ϵ(1)=1`$ and $`ϵ(x<1)1`$ for a physically meaningful solution, it must reach a minimum somewhere in the interval $`x:[0,1]`$. We find $`\text{d}ϵ/\text{d}xϵ^{}(x)=0`$ from basic calculus at

$$x=0.225259;\mathrm{minimum}\mathrm{of}ϵ=0.442201$$

when $`\lambda _e=0.509404`$. Solutions exist only for an SFE in the range $`ϵ:[0.442201,1]`$. The result (20) is displayed on Fig. 2 where we graph the same quantity $`\lambda _e`$ as in Fig. 1 but for a Plummer DF. The curves are best understood if one starts at $`(ϵ,\lambda _e)=(1,1)`$, upper right-hand corner of the plot, shifting to the left as the SFE $`ϵ`$ is reduced. The most striking feature is the sudden drop in bound fraction $`\lambda _e`$ for an SFE $`44.22`$% (solid curve, Fig. 2). The lower branch of the parameter solution $`\lambda _e(x)`$ is shown as dash but corresponds to no physical state of the system. For comparison, we have also displayed as open circle the results of a series of N-body calculations, which show remarkable agreement with (18). Details of the N-body calculations will be found in Paper II.

Around the critical value, the predicted bound fraction varies from $`\lambda _e0.63`$ upward of an SFE $`ϵ=0.443`$, to 0 for an SFE $`<0.442`$. For $`ϵ=0.50`$, we compute a fraction of bound stars $`\lambda _e86\%`$. This counts as bound nearly all stars initially in the system, when we expected complete dissolution from (3).

To illustrate the significance of repeated iterations on (17), we also plot on the figure the result when no iterations are performed and the solution $`\lambda _e(ϵ)`$ taken directly from (17) (dotted line, Fig. 2). This graphs a selection of stars by energy but without taking into account the shallower gravitational potential as a result of stars leaving the system. As the bound fraction remains finite for any positive SFE, and no sudden dissolution is predicted, we conclude that feedback on the gravitational potential from stars lost changes the character of the solution very significantly. This should not come as a surprise since cluster dissolution implies a stream of stars moving outwards. What this shows is that the impact of this stream must be taken into account when quantifying the fraction of stars which ultimately will form a bound equilibrium after a phase of dynamical evolution. By contrast, Adams (2000) allows $`ϵ`$ to vary with position but does not iterate to establish final membership.

### 4.2 Analytic results for all polytropes

The results of the previous section may be extended to all polytropes of index $`n>1/2`$ defined by (14). The general solution $`\lambda _e(x)`$ may be written for any index $`n`$ according to the definition (19). We find

$$\lambda _e(x)=\frac{x_2^{3/2}F_1(\frac{3}{2},\frac{3}{2}n,\frac{5}{2},x)}{{}_{2}{}^{}F_{1}^{}(\frac{3}{2},\frac{3}{2}n,\frac{5}{2},1)}=\frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(n\frac{1}{2})}\frac{4x^{3/2}}{3\sqrt{\pi }}_2F_1(\frac{3}{2},\frac{3}{2}n,\frac{5}{2},x)$$
(21)
where $`\mathrm{\Gamma }(u)`$ is the $`\mathrm{\Gamma }`$function, and the hyper-geometric function $`{}_{2}{}^{}F_{1}^{}(k,l,m,z)`$ assumes polynomial expressions for integer values of $`n`$ (see e.g. Weisstein 1999 for details). We may extend the steps for the Plummer sphere to any polytrope, identifying the minimum of $`ϵ(x)`$ in each case. This is done easily with software such as mathematica or maple.

Results for $`x`$ where $`ϵ^{}=0`$ are given in Table 1 for several polytropes of index $`n`$ ranging from $`n=\frac{3}{2}`$ to $`40`$. The minimum of $`ϵ`$ is also listed along with the stellar mass fraction $`\lambda _e`$. The case of $`n=3/2`$ is special. Then $`ϵ=x^{{\scriptscriptstyle \frac{1}{2}}}`$ everywhere in the allowed interval, and so $`ϵ^{}(x)<0`$ for all $`x`$’s. Indeed no solutions exist with $`ϵ(x)<1`$ for all indices $`n`$ in the range $`\frac{1}{2}<n\frac{3}{2}`$. The $`n=\frac{3}{2}`$ polytrope admits $`ϵ=\lambda _e=1`$ as a special solution satisfying our constraints on bound mass, however $`ϵ^{}(1)0`$ does not satisfy the mathematical requirement of a local minimum.

All polytropes of index $`n5`$ stretch to infinity, and those with $`n>5`$ all have infinite mass. We note that $`ϵ`$ decreases smoothly with increasing $`n`$ and in particular $`ϵ0`$ as $`1/n`$. The same is true of the mass fraction $`\lambda _e`$, though for $`n=40`$ we find $`\lambda _e0.40`$ for a relatively low SFE of $`ϵ5\%`$. The bound fraction $`\lambda _e`$ at low SFE recovered from the polytropic DF overlaps with observational estimates if the index $`n`$ is sufficiently large. However, all polytropes with $`n<5`$ dissolve at a larger SFE than the Plummer case ($`n=5`$). Since $`n<5`$ polytropes show a broad region of near-uniform volume density, these data suggest that the core-halo structure plays a role in the survival of clusters, in the sense that more concentrated polytropes (larger index $`n`$) resist dissolution better.

## 5 Time-evolution and other considerations

The analytic approach of Section 3 did not aim to relate final equilibrium to initial conditions, as we have completely ignored details of the time-evolution. We turn here to basic considerations of cluster dynamics that were left out until now.

Adams (2000) has recently revisited the problem of forming bound clusters through dynamical evolution. His analytical approach leads to the conclusion that in order to form bound clusters, of low SFE averaged over the entire cluster, the SFE must peak with the central stellar density. In his problem, therefore, the large concentration of stars results from efficient star formation in a small volume. Our own results would support this view, in as much as we find the more centrally concentrated polytropes to dissolves at lower SFE. We did not, however, consider an SFE varying radially with the density. The problem would seemingly be solved from the freedom in profiling the SFE with the stellar potential. Tenorio-Tagle et al. (1986) had noted that gas density gradients have a strong impact on the outcome, since more bound gas (steep gradients) requires more mechanical energy for removal, which demands more time under a constant UV flux: this in effect reduces the ratio of dynamical time to gas-evacuation time and shifts the dynamics toward a slow evolution regime, which we briefly outline below.

In condensed star clusters, the dynamical time $`t_{\mathrm{cr}}`$ at the half-mass radius may be long compared to a finite gas-evacuation timescale $`\tau `$, while being shorter than $`\tau `$ at the centre of the cluster. When $`t_{\mathrm{cr}}/\tau 1`$, the integral over an orbital path $`l`$

$$J_\varphi _l𝒗\text{d}𝒙vr(E/M)^{1/2}rM^{1/2}r^{1/2}$$
(22)
is an adiabat so $`J_\varphi =`$ constant during evolution (Weinberg 1993; BT+87, §3.6). Here $`M`$ is the cluster mass, $`E`$ the binding energy and $`v,r`$ the orbital velocity and radius of a star, respectively. The regime $`t_{\mathrm{cr}}/\tau 1`$ reduces to an adiabatic, quasi-static transformation of the initial mass profile, with mass loss acting as a positive source of energy. Notice that in order for (22) to apply, no limit is set on the fraction of mass lost through gas expulsion, only that it proceeds slowly: then no stars are lost through adiabatic evolution. Indeed under such condition clusters prove much more robust to total disruption (Lada et al. 1984; Goodwin 1997). From (22) we deduce $`t_{\mathrm{cr}}\sqrt{r^3/M}r^2/J_\varphi `$. Since the system expands, it will move out of the adiabatic regime $`t_{\mathrm{cr}}/\tau 1`$ (i.e., $`t_{\mathrm{cr}}\tau `$) if the expansion factor

$$\frac{R_{final}}{R_{initial}}=\sqrt{\frac{\tau }{t_{\mathrm{cr},initial}}}.$$
(23)
We can draw constraints from observations of Scorpius clusters which suggest an expansion factor of up to five between clusters of different ages (e.g. Brown 2001). If we take $`\tau 10^5`$ years as reference (see the introduction), and inserting a ratio = 5 on the left-hand side of (23), the central region of a cluster would evolve adiabatically throughout the gas expulsion phase if the core $`t_{\mathrm{cr}}`$ is initially on the order of $`4\times 10^3`$ years. For a condensed cluster such that initially $`t_{\mathrm{cr}}\tau `$ at the half-mass radius, we may assume that $`t_{\mathrm{cr}}1/\sqrt{G\rho }\mathrm{}>\tau `$ from that point on. The ratio of central to half-mass densities $`\rho (0)/\rho _{1/2}=t_{\mathrm{cr}}^2(r_{1/2})/t_{\mathrm{cr}}^2(0)>\tau ^2/t_{\mathrm{cr}}^2(0)`$ would then be $`>(10^5/4\times 10^3)^2=62510^3`$. A Michie-King cluster model (cf. King 1966) with $`\mathrm{\Psi }/\sigma ^2\mathrm{}>9`$ meets this condition, however a Plummer sphere does not. This points to yet greater dependence on the DF than was obtained for the family of polytropes. Results obtained for the King family of DF will be presented in Paper II.

This analysis reinforces our view that instantaneous gas removal provides the most severe criterion for the survival of bound cluster, since any realistic situations will be a compromise between rapid gas removal and the adiabatic ($`t_{\mathrm{cr}}\tau `$) regime. Hills (1980) had pointed out that explosive mass loss may yet occur at a time when the cluster has not fully reached dynamical equilibrium. If the stars’ velocities were small initially, the cluster would contract radially until $`T/|W|1`$ so that $`Q=2T/|W|=2`$ at most (e.g., Boily, Athanassoula & Kroupa 2002). Two possibilities may occur: in the first one, the gas is drawn in with the stars so that the mass ratio of gas to stars is unchanged along the path of the stars. In this case we may compute the effect of instantaneous gas-expulsion by removing stars with positive energy and setting $`Q>1`$ in the DF (see Appendix B). We find for Plummer spheres that total disruption would occur if the SFE were less than about 48.8% when the virial ratio assumes its maximum value of $`Q=2`$. Clearly this has the contrary effect to the one sought. A second possibility is that the gas is expelled when the stars have not acquired large velocities and $`Q<1`$, a more likely situation given that a cluster initially at rest would spend more time in this dynamical phase as can be deduced from Kepler’s third law. We find in this case for a Plummer sphere and for $`Q1/4`$ bound stars at an SFE $`24\%`$, as in observations of star forming regions. This provides a solution for bound clusters where stellar ages are less than a mean dynamical cluster crossing time and virial equilibrium has not yet been established. This view has to be weighted against the likelihood of observing clusters out of equilibrium and the reduced time available to form stars.

## 6 General conclusions

We have presented a simple method for determining the fraction of stars that will survive a phase of rapid gas loss following the birth of the first heavy stars to form a bound cluster. This approach, based on an iterative selection of stars by binding energy, is a modified version of the problem discussed by Hills (1980). The star formation efficiency $`ϵ`$ is taken as free parameter which is held constant throughout the cluster. Some highlights from this study are:

1) The choice of distribution function determines the fraction of bound stars that remain to form a bound cluster following the rapid loss of gas (cf. Eq. 7). Taking account of the DF allowed us to estimate more accurately the bound mass than the predicted disruption for SFE $`<50\%`$ based on the virial theorem alone (cf. Eq. 3).

In Appendix A, we show using a harmonic oscillator model that a large fraction of a self-gravitating stellar system in virial equilibrium will survive gas expulsion despite a low SFE ($`<0.4`$) if the stellar velocity distribution function favours stars with low velocities. It is, therefore, important to measure stellar velocities in a cluster-forming molecular cloud core to obtain observational constraints on the shape of the DF.

2) Analytic solutions were obtained for all polytropes of index $`n\frac{3}{2}`$. Polytropes with index $`n1`$ yield a bound fraction for very low SFE, easily matching the range suggested by observations.

These findings are backed up by a series of N-body calculations, where the orbital evolution of individual stars is fully taken into account: such calculations provide a good match to the solution curve (18) for the Plummer sphere (Fig. 2). This encourages us to apply (12) to other DF and to consider the dynamics fully in order to identify which, if any, initial conditions DF may lead to bound clusters despite very low SFE. We propose to cover these topics in greater details in a companion paper (Paper II).

## acknowledgments

We are very grateful to D.C. Heggie for a thorough read of a draft version of this paper. His insight led to a much-improved section 4. We thank him and A. Lançon for comments that led to clarifications of the arguments. CMB wishes to thank A.G.A. Brown and E.K. Grebel who took him gently through the observational literature. We are grateful to them for their help and suggestions. We extend our gratitude to an anonymous referee for constructive comments. CMB was funded under the SFB 439 program in Heidelberg and wishes to thank his host R. Spurzem at the ARI for unwavering support.

## Appendix A: Bound orbits in a harmonic potential

To see how stars on different orbits may escape at different rates, as argued in Section 2, consider a self-gravitating homogeneous mass distribution of density $`\rho `$. The stellar orbital time $`t_{\mathrm{cr}}`$ is then independent of position. Stars are often distributed uniformly in space over a small volume so this situation has bearing on actual systems. Our aim is to relate the membership of a star to its orbital parameters prior to gas expulsion.

The equation of motion for this case reduces to that of a harmonic oscillator which in one dimension reads

$$\ddot{r}(t)=\varphi (r)=\frac{GM(<r)}{r^2}=\omega ^2r(t)$$
(24)
where we set $`\omega ^2=4\pi G\rho /3`$. Integrating (24) we have

$$r(t)=r_o\mathrm{cos}(\omega t+\theta _o)$$
(25)
where the amplitude $`r_o`$ and angle $`\theta _o`$ are set by the initial conditions. When the homogeneous region is cut off at radius $`r=R`$, the potential everywhere in space is

$$\varphi (r)=\frac{3}{2}\frac{GM}{R}+\frac{\omega ^2}{2}r^2rR,$$
(26)
and $`\varphi (r)=GM/r`$ for $`r>R`$. The derivative of $`\varphi `$ is discontinuous at $`r=R`$ but this has no bearing on the argument presented here since we consider motion with $`r(t)<R`$ only. We imagine an ensemble of stars sharing the same orbit of amplitude $`r_o`$ but otherwise phase-mixed. The energy per unit mass for that orbit

$$=\varphi (r)+\frac{1}{2}v^2(r)=\varphi (r_o)$$

is modified at time $`t=t_o`$ when a fraction of the mass is lost and only stars remain: $`\varphi _{}(r)=ϵ\varphi (r)`$ is now the gravitational potential. The new energy per unit mass at radius $`r`$

$$_{}=\varphi _{}(r)+\frac{1}{2}v^2(r)=ϵ\varphi (r)+\varphi (r)$$

is positive if $`v^2>2\varphi _{}(r)=2ϵ\varphi (r)`$ which is possible provided (regrouping terms and remembering  )

$$\omega ^2r_o^2+ϵ\omega ^2r^2ϵ\frac{3GM}{R}=ϵ3\omega ^2R^2.$$
(27)
Since the original orbit (25) has $`r(t)r_o`$, we find from (27) a minimal amplitude $`r_o`$ for unbound orbits as function of the SFE $`ϵ`$ :

$$\sqrt{\frac{3ϵ}{1+ϵ}}\frac{r_o}{R}<1,$$
(28)
where the upper limit follows from our definition of $`R`$. Equation (28) states that all stars on orbits with amplitude $`r_o`$ sufficiently small will remain bound, while at least a fraction of those on large-amplitude orbits will escape. Otherwise said, stars on low-binding energy orbits, ie which achieve high velocities at some point on their orbit, are more likely to unbind following gas expulsion. All orbits with $`r_oR`$ would lose none or a smaller fraction of stars than those that visit the edge of the system.

## Appendix B: Deviations from dynamical equilibrium

The virial ratio $`Q2T/|W|`$ = 1 for a system in dynamical equilibrium. In practice the DF may not satisfy this condition exactly and this requires modifications to (17). We consider the case of the Plummer sphere only but a similar treatment would apply to all polytropes.

Substituting absolute values for $`E`$ in (15) allows to treat sub- as well as super-virial dynamics. Then the integrals (17) read

$$1\lambda _e=\frac{{\displaystyle _{ϵ^{1/2}}^{Q^{1/2}}}(|1u^2|)^{7/2}u^2du}{{\displaystyle _0^{Q^{1/2}}}(|1u^2|)^{7/2}u^2du}=1\frac{(ϵ)}{(Q)+𝒢(Q)}$$
(29)
where the function $``$ is defined by

$$(x)\left(\mathrm{sin}^1x^{1/2}p(x)x^{1/2}\sqrt{1x}\right)$$
(30)
with $`p(x)`$ defined by (17), and

$$𝒢(Q>1)\frac{7}{256}\left[\sqrt{Q}\sqrt{Q1}p(Q)+\mathrm{log}(\sqrt{Q}+\sqrt{Q1})\right]$$
(31)
and $`𝒢(Q1)=0.`$ Clearly the numerator in (29) evaluates to zero if $`ϵQ`$ when the solution $`\lambda _e=1`$. Since the SFE $`ϵ`$ is bounded to $`[0,1]`$, this is possible only for sub-virial conditions $`(Q<1)`$. We recover (17) by setting $`Q=1`$ in (29) since $`(1)=\pi /2`$.

Sub-virial conditions are relevant to the problem of forming bound clusters (e.g. Klessen & Burkert 2000; Clarke et al. 2000), with the caveat that sub-virial conditions at the gas-expulsion time require synchronous star formation throughout the cluster to within a fraction of the dynamical time, when the spread in colours suggest at least a formation time of a few million years (Hartmann et al. 2001). We wish to show that a point of the parameterised solution (18) and (20) shifts to lower $`ϵ`$ as $`Q`$ is reduced from 1. Treating $`Q`$ as the variable in (29) and differencing we have

$$\frac{\text{d}ϵ(x_o)}{\text{d}Q}=\frac{x_o}{(x_o)}^{}(Q);\frac{\text{d}\lambda _e(x_o)}{\text{d}Q}=\frac{\lambda _e(x_o)}{^2(x_o)}^{}(Q),$$
(32)
where $`x_o`$, the parameter defined in (19), identifies a point of the solution curve $`(\lambda _e,ϵ)`$ when $`Q=1`$. The derivative

$$^{}(Q)=\frac{128}{7}\sqrt{Q}\left(1Q\right)^{7/2}>0Q<1.$$

It is easy to see from the above equation and (32) that reducing the virial ratio $`(\text{d}Q<0)`$ shifts the abscissa $`ϵ`$ to smaller values while increasing the bound fraction $`\lambda _e`$, for any point $`x_o`$. Therefore ‘cooler’ initial conditions lead to bound clusters for lower SFE $`ϵ`$. The gap from finite $`\lambda _e`$ to zero (total dissolution) is more pronounced but occurs at lower values of $`ϵ`$ as $`Q`$ is reduced.

Since $`^{}(Q)`$ is near zero around $`Q=1`$, the net effect of reducing the virial ratio $`Q`$ from unity is at first marginal and only becomes significant once $`Q\frac{1}{2}`$. Similar conclusions apply to super-virial ($`Q>1`$) conditions. Total dissolution occurs for larger SFE $`ϵ`$ in these cases. A cluster collapsing from rest would reach $`Q=2`$ at most during violent relaxation. Dissolution occurs for an $`ϵ48.4\%`$, only marginally larger than the 44.2% value obtained for equilibrium configurations. The solution for $`Q=2`$ and $`ϵ=1`$ (SFE of 100%) yields a bound fraction $`\lambda _e0.91`$, such that 9% of the stars escape. This is remarkably similar to what is observed in N-body realisations of collapsing cold spheres.


format: text

text:
# Renormalization Group Summation and the Free Energy of Hot QCD

## I Introduction and Summary

The renormalization procedure in quantum field theory inevitably introduces a renormalization scale parameter $`\mu ^2`$ into perturbative results; as $`\mu ^2`$ is unphysical, its value is in principle arbitrary. Consequently, the calculation of any physical quantity would necessarily result in the dependence on $`\mu ^2`$ disappearing. However, at any finite order of perturbation theory, residual dependence on $`\mu ^2`$ renders the result ambiguous, for changing the value of $`\mu ^2`$ changes the predicted value of the physical quantity that has been computed Stevenson:1981vj .

This problem has proved to be particularly acute in the calculation of the thermodynamic free energy in thermal field theory Kap:FTFT . In quantum chromodynamics (QCD) at temperatures much larger than the deconfinement temperature the free energy has been calculated Arnold:1995eb ; Zhai:1995ac ; Braaten:1996jr  to order $`\alpha _s^{5/2}`$, four terms in the perturbative series beyond the leading ideal-gas term. Whereas the first few approximations turn out to show little sign of convergence for any temperature of practical interest, the result to order $`\alpha _s^{5/2}`$ happens to be centered about the results obtained in lattice gauge theory, but the dependence on the renormalization scale parameter $`\mu ^2`$ is so large that it has little predictive power.

Recently, in Ref. Ahmady:2002fd  in the context of standard model calculations, it was shown how the renormalization group (RG) equation can be used to sum in a systematic manner the leading log (LL), next-to-leading log (NLL), …effects given the perturbative results to one loop, two loop, …order. This so-called “renormalization group summation” (RG$`\mathrm{\Sigma }`$) has been found to lead to a considerable reduction of the dependence on the parameter $`\mu ^2`$ within a given renormalization scheme.

In the case of thermal field theory this procedure requires generalization because the perturbative series also involves half-integer powers and logarithms of $`\alpha _s`$. This is carried out in Sec. II and applied to the available three-loop result for the free energy in hot QCD in the modified minimal subtraction ($`\overline{\mathrm{MS}}`$) scheme. The RG$`\mathrm{\Sigma }`$ result which includes all LL and NLL effects turns out to be only weakly dependent on $`\mu ^2`$. Unfortunately, as discussed in Sec. III, this does not increase substantially the predictive power of the three-loop calculation, because even within a fixed renormalization scheme (here the $`\overline{\mathrm{MS}}`$ scheme) there arises a new ambiguity in the form of the initial conditions for the differential equations, whose integration carries out the RG-summation. This ambiguity is somewhat larger than at zero temperature because the RG-summation of the thermal perturbative series leads to two uncoupled sets of differential equations rather than one, reflecting the particular difficulties that thermal perturbation theory present.

In the Appendix, we briefly discuss the difference between the strictly perturbative solution for the running coupling $`\alpha _s`$ to two-loop order and the exact solution at two-loop order, which can be given in closed form in terms of Lambert’s W function.

## II Renormalization Group Summation in Hot QCD

In QCD with $`n_f`$ flavors of quarks, the thermodynamic free energy at high temperature has been computed to be Arnold:1995eb ; Zhai:1995ac ; Braaten:1996jr

$``$ $`=`$ $`{\displaystyle \frac{8\pi ^2}{45}}T^4\{(1+{\displaystyle \frac{21}{32}}n_f)+{\displaystyle \frac{15}{4}}(1+{\displaystyle \frac{5}{12}}n_f){\displaystyle \frac{\alpha _s}{\pi }}+30\left[(1+{\displaystyle \frac{n_f}{6}})\left({\displaystyle \frac{\alpha _s}{\pi }}\right)\right]^{3/2}`$ (1)
$`+\{237.2+15.97n_f0.413n_f^2+{\displaystyle \frac{135}{2}}(1+{\displaystyle \frac{n_f}{6}})\mathrm{ln}\left[{\displaystyle \frac{\alpha _s}{\pi }}(1+{\displaystyle \frac{n_f}{6}})\right]`$
$`{\displaystyle \frac{165}{8}}(1+{\displaystyle \frac{5}{12}}n_f)(1{\displaystyle \frac{2}{33}}n_f)\mathrm{ln}{\displaystyle \frac{\overline{\mu }}{2\pi T}}\left\}\right({\displaystyle \frac{\alpha _s}{\pi }})^2`$
$`+(1+{\displaystyle \frac{n_f}{6}})^{1/2}[799.221.96n_f1.926n_f^2`$
$`+{\displaystyle \frac{495}{2}}(1+{\displaystyle \frac{n_f}{6}})(1{\displaystyle \frac{2}{33}}n_f)\mathrm{ln}{\displaystyle \frac{\overline{\mu }}{2\pi T}}]\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^{5/2}+𝒪(\alpha _s^3\mathrm{ln}\alpha _s)\}`$
where $`\overline{\mu }`$ is the renormalization scale parameter of the $`\overline{\mathrm{MS}}`$ scheme and $`\alpha _s(\overline{\mu })`$ is the running coupling in this scheme whose form to three-loop order is given in Groom:2000in , though we shall restrict ourselves to its two-loop version in the following (see the Appendix for more discussion).

Changing the renormalization scale parameter $`\overline{\mu }`$ in principle should not alter the value of $`F`$, as changes in the explicit $`\overline{\mu }`$ are compensated for by changes of $`\alpha _s(\overline{\mu })`$, and indeed, the result (1) is independent of $`\overline{\mu }`$ to order $`\alpha _s^{5/2}`$. However, numerically the dependence on $`\overline{\mu }`$ is large, in fact larger than that of the result to order $`\alpha _s^1`$ unless $`\overline{\mu }`$, which has to be of the same order as $`T`$ to avoid large logarithms, is much larger than even the electroweak scale.

Extrapolating from Eq. (1), we assume the complete all-order result of $``$ can be represented by a series of the form

$$/_0=1+\underset{n=0}{\overset{\mathrm{}}{}}\left(R_n(u)x^{n+1}+S_n(u)x^{n+\frac{3}{2}}+T_n(u)x^{n+2}\mathrm{ln}x\right)$$
(2)
where $`_0`$ is the ideal-gas value, $`x=\alpha _s(\overline{\mu })/\pi `$, $`u=xL`$, $`L=\mathrm{ln}(\overline{\mu }^2/(2\pi T)^2)`$ and

$$R_n(u)=\underset{m=0}{\overset{\mathrm{}}{}}A_{n+m,m}u^m,S_n(u)=\underset{m=0}{\overset{\mathrm{}}{}}B_{n+m,m}u^m,T_n(u)=\underset{m=0}{\overset{\mathrm{}}{}}C_{n+m,m}u^m,$$
(3)
although only the coefficients with $`n1`$ are accessible by thermal perturbation theory. In fact, all of the perturbatively accessible cofficients with $`m=0`$ have been calculated already, with the exception of $`C_{1,0}`$, which is forthcoming111Y. Schröder, private communication.

In Ref. Ahmady:2002fd  it has been shown how to sum all RG-accessible logarithms when the lowest-order coefficients to the sums in Eq. (3), $`A_{n,0}`$, $`B_{n,0}`$, and $`C_{n,0}`$, and the $`\beta `$-function coefficients (see Eq. (9) below) are known. In this paper we extend this to a perturbative series of the form (2).

While the perturbative expression (1) is in powers of $`x^{1/2}`$, successive RG$`\mathrm{\Sigma }`$-perturbative expressions are given by

$`_{\mathrm{RG}\mathrm{\Sigma }}^{(1)}/_0`$ $`=`$ $`1+xR_0(xL),`$ (4)
$`_{\mathrm{RG}\mathrm{\Sigma }}^{(2)}/_0`$ $`=`$ $`1+xR_0(xL)+x^{3/2}S_0(xL),`$ (5)
$`_{\mathrm{RG}\mathrm{\Sigma }}^{(3)}/_0`$ $`=`$ $`1+xR_0(xL)+x^{3/2}S_0(xL)+x^2(R_1(xL)+T_0(xL)\mathrm{ln}x),`$ (6)
$`_{\mathrm{RG}\mathrm{\Sigma }}^{(4)}/_0`$ $`=`$ $`1+xR_0(xL)+x^{3/2}S_0(xL)+x^2(R_1(xL)+T_0(xL)\mathrm{ln}x)+x^{5/2}S_1(xL)`$ (7)
which all are perturbatively accessible in hot QCD, as we shall see.

The explicit dependence of $``$ on $`\overline{\mu }^2`$ and its implicit dependence through $`x(\overline{\mu }^2)`$ are such that

$$\overline{\mu }^2\frac{d}{d\overline{\mu }^2}=0=\left(\overline{\mu }^2\frac{}{\overline{\mu }^2}+\beta (x)\frac{}{x}\right)$$
(8)
where

$$\beta (x)=\overline{\mu }^2\frac{x}{\overline{\mu }^2}=(b_2x^2+b_3x^3+b_4x^4+\mathrm{}).$$
(9)
Substitution of (2) into (8) yields

$`0`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\{[R_n^{}x^{n+2}+(b_2x^2+\mathrm{})(uR_n^{}+(n+1)R_n+xT_n)x^n]`$ (10)
$`+\left[S_n^{}x^{n+\frac{5}{2}}+(b_2x^2+\mathrm{})(uS_n^{}+[n+{\displaystyle \frac{3}{2}}]S_n)x^{n+\frac{1}{2}}\right]`$
$`+[T_n^{}x^{n+3}+(b_2x^2+\mathrm{})(uT_n^{}+(n+2)T_n)x^{n+1}]\mathrm{ln}x\}`$
Lowest order terms of the form $`x^n`$, $`x^{n+\frac{1}{2}}`$ and $`x^n\mathrm{ln}x`$ give rise to differential equations for $`R_0(u)`$, $`S_0(u)`$ and $`T_0(u)`$, respectively. The boundary conditions on these equations are the computed values of $`A_{0,0}`$, $`B_{0,0}`$ and $`C_{0,0}`$ respectively, which can be read off Eq. (1). This gives

$$R_0(u)=A_{0,0}w^1,S_0(u)=B_{0,0}w^{3/2},T_0(u)=C_{0,0}w^2,w(1+b_2u).$$
(11)
These functions incorporate the LL contributions to $``$ to all orders.

To next order in the coupling $`x`$, we find differential equations for $`R_1`$, $`S_1`$, and $`T_1`$ (which rely on knowning the above solutions for $`R_0`$, $`S_0`$, and $`T_0`$). Solving these equations gives

$`R_1`$ $`=`$ $`w^2\left[A_{1,0}\left({\displaystyle \frac{b_3}{b_2}}A_{0,0}+C_{0,0}\right)\mathrm{ln}w\right]`$ (12)
$`S_1`$ $`=`$ $`w^{5/2}\left[B_{1,0}{\displaystyle \frac{3}{2}}{\displaystyle \frac{b_3}{b_2}}B_{0,0}\mathrm{ln}w\right]`$ (13)
$`T_1`$ $`=`$ $`w^3\left[C_{1,0}2{\displaystyle \frac{b_3}{b_2}}C_{0,0}\mathrm{ln}w\right]`$ (14)
which incorporate the NLL contributions to $``$ to all orders.

Continuing in this way, we obtain

$`R_2`$ $`=`$ $`w^3[A_{2,0}(2{\displaystyle \frac{b_3}{b_2}}A_{1,0}+{\displaystyle \frac{b_3^2}{b_2^2}}A_{0,0}+C_{1,0}+{\displaystyle \frac{b_3}{b_2}}C_{0,0})\mathrm{ln}w`$ (15)
$`+{\displaystyle \frac{b_3}{b_2}}({\displaystyle \frac{b_3}{b_2}}A_{0,0}+2C_{0,0})\mathrm{ln}^2w+({\displaystyle \frac{b_3^2}{b_2}}b_4)A_{0,0}u]`$
$`S_2`$ $`=`$ $`w^{7/2}\left[B_{2,0}{\displaystyle \frac{b_3}{b_2}}\left({\displaystyle \frac{5}{2}}B_{1,0}+{\displaystyle \frac{3}{2}}{\displaystyle \frac{b_3}{b_2}}B_{0,0}\right)\mathrm{ln}w+{\displaystyle \frac{15}{8}}{\displaystyle \frac{b_3^2}{b_2^2}}B_{0,0}\mathrm{ln}^2w+{\displaystyle \frac{3}{2}}\left({\displaystyle \frac{b_3^2}{b_2}}b_4\right)B_{0,0}u\right]`$ (16)
$`T_2`$ $`=`$ $`w^4\left[C_{2,0}{\displaystyle \frac{b_3}{b_2}}\left(3C_{1,0}+2{\displaystyle \frac{b_3}{b_2}}C_{0,0}\right)\mathrm{ln}w+3{\displaystyle \frac{b_3^2}{b_2^2}}C_{0,0}\mathrm{ln}^2w+2\left({\displaystyle \frac{b_3^2}{b_2}}b_4\right)C_{0,0}u\right]`$ (17)
these give the NNLL contributions to $``$.

The coefficients $`A_{0,0}`$, $`B_{0,0}`$, $`C_{0,0}`$, $`A_{1,0}`$, and $`B_{1,0}`$ are determined in the $`\overline{\mathrm{MS}}`$ scheme by Eq. (1) and this allows us to construct the RG$`\mathrm{\Sigma }`$ approximation $`_{\mathrm{RG}\mathrm{\Sigma }}^{(4)}`$ given in Eq. (7). The next approximation would also involve $`T_1`$, $`R_2`$, $`S_2`$, and $`T_2`$. For these, one would need results for $`C_{1,0}`$ and $`A_{2,0}`$, $`B_{2,0}`$, $`C_{2,0}`$, but only $`C_{1,0}`$ is computable in perturbation theory because the order $`\alpha _s^3`$-contribution to $``$ is inherently nonperturbative Linde:1980ts ; Gross:1981br . The remaining coefficients would have to be derived from a nonperturbative framework such as lattice gauge theory, or they might be estimated by nonperturbative resummation techniques such as Padé approximations Kastening:1997rg ; Hatsuda:1997wf . Given that, the above formulae would allow the construction of the next RG$`\mathrm{\Sigma }`$ approximation $`_{\mathrm{RG}\mathrm{\Sigma }}^{(5)}`$.

## III Discussion

In Fig. 1 the renormalization scale dependence of the purely perturbative result (1) for $`/_0`$ at $`T=3T_c`$ for $`n_f=0`$ is displayed using a two-loop running coupling constant $`\alpha _s(\overline{\mu }=e^\lambda 2\pi T)`$ with $`\lambda `$ varying between $`1`$ and 1 (long-dashed line). The two dots on the vertical axis give two recent lattice results from Refs. Boyd:1996bx ; Okamoto:1999hi , and the perturbative result is seen to agree well with these for $`\overline{\mu }=2\pi T`$ (i.e. $`\lambda =0`$), but deviates strongly for different choices of the renormalization scale.

The RG$`\mathrm{\Sigma }`$ result as obtained above is given by the solid line, and it shows a rather weak dependence on the renormalization scale. The systematic RG summation as described above is thus able to absorb almost all of the scale dependence.

However, this does not mean that the RG$`\mathrm{\Sigma }`$ approach predicts the perturbative result at $`\overline{\mu }=2\pi T`$. Rather, the latter has been used as the intial condition for the differential equations determining the functions $`R_0`$, $`S_0`$, $`T_0`$, $`R_1`$, and $`S_1`$. Organizing the RG summation in terms of a variable $`L_\kappa =\mathrm{ln}[\overline{\mu }^2/(\kappa 2\pi T)^2]`$ instead of the variable $`L`$ introduced in Eq. (2) would have led to a different definition of the constants $`A_{0,0}`$, $`B_{0,0}`$, $`C_{0,0}`$ that provide the initial conditions. Moreover, the series involving half-integer powers of $`x`$ in Eq. (2) leads to differential equations that decouple from those responsible for the series involving integer powers and logarithms of $`x`$. The former series could therefore have been derived by introducing a different variable $`L_\nu =\mathrm{ln}[\overline{\mu }^2/(\nu 2\pi T)^2]`$ for $`S_n`$.

In identifying the ambiguity in $`L_\kappa `$ and $`L_\nu `$ (as parametrized by $`\kappa `$ and $`\nu `$ respectively) we are keeping the value of $`x`$ unaltered—both $`\overline{\mu }`$ is held fixed and the form of $`x(\overline{\mu })`$ dictated by having chosen to work in the $`\overline{\mathrm{MS}}`$ scheme is not changed. Rather, we note that the invariance of the perturbative result (1) under the change

$$A_{1,0}+A_{1,1}\mathrm{ln}\frac{\overline{\mu }}{2\pi T}(A_{1,0}+A_{1,1}\mathrm{ln}\kappa )+A_{1,1}\mathrm{ln}\frac{\overline{\mu }}{\kappa 2\pi T}A_{1,0}^{(\kappa )}+A_{1,1}L_\kappa $$
(18)
(with analogous equations for $`B_{1,0}^{(\nu )}`$ and $`C_{1,0}^{(\kappa )}`$) is lost when one deals with the RG$`\mathrm{\Sigma }`$ expressions $`R_1`$, $`S_1`$ and $`T_1`$. The boundary condition for the differential equations for $`R_1`$, $`S_1`$ and $`T_1`$, and the logarithm in the solution to these equations have a dependence on $`\kappa `$ (or $`\nu `$) that no longer automatically compensates as in the perturbative result.

In Fig. 1 the short-dashed lines forming a big Z show the result of varying $`\kappa `$ and $`\nu `$ around the value 1 by a factor $`e^\lambda `$ with $`\lambda `$ between $`1`$ and $`1`$. The upper bar in the big Z is formed by $`\nu =e^1`$, $`\kappa =e^\lambda `$, the diagonal by $`\nu =\kappa =e^\lambda `$, and the lower bar by $`\nu =e^1`$, $`\kappa =e^\lambda `$. Evidently, when $`\kappa `$ and $`\nu `$ are identified, the ambiguity is a little bit smaller than the one given by the renormalization scale dependence of the purely perturbative result, but varying $`\nu `$ independently of $`\kappa `$ leads to even larger variations, at their extremes even exceeding the ideal gas result for $`||`$. Unfortunately, when varying $`\kappa `$ and $`\nu `$, either together or independently, there is no saddle point for $`_{\mathrm{RG}\mathrm{\Sigma }}^{(4)}`$ that would allow one to eliminate these ambiguities by a principle of minimal sensitivity Stevenson:1981vj .

The renormalization group equation (8) states that the thermodynamic potential is independent of the renormalization scale $`\overline{\mu }`$. When one uses this equation to incorporate all the logarithmic contributions coming from higher order perturbation theory whose form is implied by this equation, one indeed finds the dependence of a perturbative approximation to $``$ on $`\overline{\mu }`$ being diminished, as expected. However, we have identified another source of ambiguity (characterized by $`\kappa `$ and $`\nu `$) which leads to large variations of the perturbative RG$`\mathrm{\Sigma }`$ result for $``$. This highlights the numerically strong scheme dependence of the perturbative result which still remains after having eliminated the strong dependence on the renormalization scale $`\overline{\mu }`$ by the RG$`\mathrm{\Sigma }`$ method. (There are, of course, other ambiguities in the RG$`\mathrm{\Sigma }`$ result that could arise due to changing the renormalization scheme from $`\overline{\mathrm{MS}}`$; these we have not addressed.)

As can be seen from Fig. 1, the extra dependence on the parameter $`\nu `$ associated with the part of the perturbative series involving half-integer powers in $`\alpha _s`$ is in fact the one which dominates the uncertainties of the perturbative result (the short-dashed lines are rather flat when only $`\kappa `$ is varied, but depend strongly on $`\nu `$). This part of the perturbative series is exclusively associated with “soft” collective phenomena such as screening and Landau damping and calls for a more complete treatment than conventional perturbation theory is able to achieve. Recent attempts in this direction that have been put forward include separate Padé approximations to soft and hard contributions Cvetic:2002ju , optimization of perturbation theory using the hard-thermal-loop effective action ABS , and approximately self-consistent propagator resummation in the so-called $`\mathrm{\Phi }`$-derivable approach BIR .

###### Acknowledgements.

Financial support has been provided by NSERC. D. G. C. McK. would like to thank the Technical University of Vienna where this work was performed. R. and D. MacKenzie provided helpful suggestions. Y. Schröder and G. Moore pointed out typographical errors in the first version of this paper.

*

## Appendix A

The beta function to two-loop order,

$$\mu ^2\frac{dx}{d\mu ^2}=\frac{dx}{dL}=b_2x^2+b_3x^3,$$
(19)
can be integrated in closed form

$$_{L_0}^L𝑑L^{}=LL_0=_{x_0}^x\frac{dx^{}}{b_2x^2+b_3x^3}=\frac{1}{b_2}\left[\rho \mathrm{ln}\frac{1+\rho x}{x}\frac{1}{x}\right]_{x_0}^x(\rho =b_3/b_2)$$
(20)
which becomes, if $`x_0\mathrm{}`$ as $`L_0\mathrm{ln}\mathrm{\Lambda }^2`$,

$$\stackrel{~}{L}\mathrm{ln}\frac{\mu ^2}{\mathrm{\Lambda }^2}=\frac{\rho }{b_2}\mathrm{ln}\left[We^W(\rho e)\right]\frac{\rho }{b_2}\mathrm{ln}\rho $$
(21)
where

$$W=1\frac{1}{\rho x},\text{or}x=[\rho (W+1)]^1.$$
(22)
Hence, the solution to Eq. (19) can be written Gardi:1998qr  as a Lambert W function Corless:1996  with

$$W(z)e^{W(z)}=z=\frac{1}{e}\left(\frac{\mu ^2}{\mathrm{\Lambda }^2}\right)^{b_2/\rho },$$
(23)
where the real branch $`W_1`$ with $`W<1`$ has to be taken.

On the other hand, the standard perturbative two-loop result is usually given as

$$x=\frac{1}{b_2\stackrel{~}{L}}\left(1+\frac{\rho }{b_2}\frac{\mathrm{ln}\stackrel{~}{L}}{\stackrel{~}{L}}\right)$$
(24)
which assumes small $`x`$ and correspondingly large $`L`$.

In Fig. 2 we compare the “exact” two-loop result (23) for $`x\alpha _s/\pi `$ (solid line) with the strictly perturbative one from Eq. (24) (dashed line) for the case of QCD with $`n_f=0`$ ($`b_2=11/4`$, $`\rho =51/22`$) as a function of $`\stackrel{~}{L}`$. The divergence at $`\stackrel{~}{L}=0`$ makes itself felt significantly earlier (as $`\stackrel{~}{L}`$ approaches zero) in the standard perturbative result (24) than in the exact two-loop result (23). However, for $`\stackrel{~}{L}4`$, the difference between the two coupling is less than 1.5%. In the perturbative treatment of hot QCD, if one chooses a renormalization scale $`\overline{\mu }=2\pi T`$ this is indeed the case for all $`T>T_c\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ (in pure glue QCD a typical value, which we have adopted in this paper, is $`T_c=1.14\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$); a noticeable difference thus arises only for smaller choices of $`\overline{\mu }/T`$ in combination with $`T`$ sufficiently close to $`T_c`$.


text:
# The 𝑠𝑢(1,1)-Algebraic Boson Model in the Deformed Boson Scheme

## Acknowledgements

The main part of this work was performed when two of the authors (Y.T. and M.Y.) stayed at Coimbra in August and September of 2002. They express their sincere thanks to Professor J. da Providência, co-author of this paper.