text:
# On the use of a non-extensive distribution in the solar plasma

## Abstract

We describe the physical motivations why we are urged to look for an ionic equilibrium distribution function slightly different from the Maxwellian to calculate nuclear reaction rates in the solar core plasma answering, at the same time, to a recent paper by Bahcall et al. (astro-ph/0010055).

In a recent paper, Bahcall et al.  discuss motivations why Salpeter formula for screening of nuclear reactions in the Sun is correct (in contrast with conclusions contained in several other papers) and show five different derivations of the Salpeter screening factor $`f_0`$.

A paper of us  is also included in this critical discussion (although our work is not concerned directly on the factorizable screening factor $`f_0`$ but, rather, on the evaluation of the reaction rates) because of the use we make of an ionic equilibrium distribution which differs very slightly from the Maxwellian distribution. The reason of using this equilibrium distribution is related to the need of including many-body effects. It was already proposed ad hoc by Clayton  and we can see, a posteriori, that coincides to the equilibrium distribution of the non-extensive thermostatistics developed by Tsallis  and used recently in many different physical applications , as for instance, gravitational  and high energy problems .

In our calculations the screening factor $`f_0`$ is taken to be unity because its effect on the reaction rates is negligible if compared to the important depletions (or enhancements) of the standard rates due to other effects (electromagnetic fluctuations in a plasma) responsible of the non-extensive equilibrium distribution function. Also, by following the derivation of $`f_0`$ in Ref. and using non-extensive distribution, we can calculate for the solar core about the same value of the screening factor. For comments based on Molecular Dynamics Model (MDM) on the validity of Salpeter screening formula see Ref. . Several authors are actually working on MDM within non-extensive Tsallis statistics . It seems to us obvious that, if all the derivations of the Salpeter formula are based on the validity of the assumptions needed for the Boltzmann-Gibbs statistics, the equilibrium distribution must be, of course, the Maxwellian distribution. However the motivations that deviations from these assumptions are not correct should be explicitly given.

Results within non-extensive statistics go toward a better agreement with the experimental results of the calculated neutrino fluxes. In spite of this achievement, we never excluded the possibility of neutrino oscillation or new physics because only experiments will give the correct answer on this problem. At the same time, we have also shown that modifications of the reaction rates do not affect bulk properties of the gravitationally stabilized solar core such as sound speed or hydrostatic equilibrium, depending on mean values obtained averaging over Maxwellian distribution function .

We confirm that Tsallis distribution is an equilibrium distribution and that the argument that there is enough time to arrive at equilibrium after a reaction is fully accomplished by the Tsallis equilibrium distribution that can be derived in, at least, three different ways containing not alternative reasons.

In one approach we exploit the knowledge that the distribution in the solar interior cannot be too much different from the Maxwellian one and add small corrections (higher-order terms in a derivative expansion) to the coefficients of the standard Fokker-Plank equation. Tsallis’ distribution is immediately generated .

A second approach focuses on the electric microfields that have been shown to exists in the solar plasma: the time-spatial fluctuations in the particles positions produce specific fluctuations of microscopic electric field in a given point of the plasma. We have shown that such electric microfields in the solar plasma (with energy density of the order of $`10^{16}`$ MeV/fm3) implies a deviation from the Maxwellian distribution; the entity depending on the value of the plasma parameter and on the correlation among ions .

The third approach has been just started and aims to connect the distribution of collective variables to memory effects and long-time correlations between velocities . There should exist solutions compatible with the Tsallis’ distribution and/or other non-Maxwellian distributions.

These three approaches are not exhaustive and not necessary alternative. Nevertheless, it is suggestive that all of them point in the same direction: the Maxwell-Boltzmann distribution of velocity should have small but nonnegligible corrections in the solar plasma and the Tsallis’ distribution could provide a better description.

Normal stellar matter, such as the one in the Sun, is non-degenerate, i.e., quantum effects are small (in fact, they are small for electrons and completely negligible for ions), it is non-relativistic, and it is in good thermodynamical equilibrium. On this ground, the particle velocity distribution is almost universally taken to be a Maxwell-Boltzmann distribution, without much questioning.

However, derivations of the ubiquitous Maxwell-Boltzmann distribution are based on several assumptions . In a kinetical approach, one assumes (1) that the collision time be much smaller than the mean time between collisions, (2) that the interaction be sufficiently local, (3) that the velocities of two particles at the same point are not correlated (Boltzmann’s Stosszahlansatz), and (4) that energy is locally conserved when using only the degrees of freedom of the colliding particles (no significant amount of energy is transferred to collective variables and fields). In the equilibrium statistical mechanics approach, one uses the assumption that the velocity probabilities of different particles are independent, corresponding to (3), and that the total energy of the system could be expressed as the sum of a term quadratic in the momentum of the particle and independent of the other variables, and a term independent of momentum, but if (1) and (2) are not valid the resulting effective two-body interaction is not local and depends on the momentum and energy of the particles. Finally, even when the one-particle distribution is Maxwellian, additional assumptions about correlations between particles are necessary to deduce that the relative-velocity distribution, which is the relevant quantity for rate calculations, is also Maxwellian.

At least in one limit the Maxwell-Boltzmann distribution can be rigorously derived: systems that are dilute in the appropriate variables, whose residual interaction is small compared to the one-body energies. In spite of the fact that the effects of the residual interaction cannot be neglected, as a good first approximation the solar interior can be studied in this dilute limit; therefore, it is reasonable to suppose that the velocity distribution in the Sun is not too far from the Maxwellian one.

Likewise, one often assumes that the solar core could be treated as an ideal (Debye) plasma. However, there are physical conditions and/or specific applications that needs higher accuracy for which becomes necessary to take into account modifications of the standard plasma theory.

When $`\mathrm{\Gamma }0.05÷0.1`$, the mean Coulomb energy potential is not much smaller of the thermal kinetic energy and the screening length $`R_Da`$. It is not possible to clearly separate individual and collective degrees of freedom. The presence of at least two different scales of energies of the same rough size produces deviations from the standard statistics which describe the system in terms of a single scale, $`kT`$.

The reaction time necessary to build up screening after a hard collision can be estimated from the inverse solar plasma frequency $`t_{pl}=\omega _{pl}^1=\sqrt{m/(4\pi ne^2)}10^{17}`$ sec, and it is comparable to the collision time $`t_{coll}=\sigma vn^110^{17}`$ sec. Therefore, several collisions are likely necessary before the particle looses memory of the initial state and the scattering process can not be considered Markovian. In addition, screening starts to become dynamical: the time necessary to build up again the screening after hard collisions is not negligible any more.

Many remarks we have just made are reported in well known books on statistical mechanics, kinetic and plasma theory and in works by astrophysicist as, for instance, Paquette et al.  and Cox et al. . Let us report, specifically, from Cox et al.: “one usually assumes (1) negligible radiative forces, (2) complete ionization, (3) Maxwellian velocity distributions and the same kinetic temperature for all ions and electrons, (4) diffusion velocities which are much smaller than mean thermal velocities, (5) no magnetic fields, (6) collisions dominated by “classical” interactions between two point particles, and (7) plasma which can be considered a diluite gas, i.e., the ideal gas equation of state applies. \[..\] When the number of particles in a Debye sphere around an ion is not $`1`$, assumption (6) and (7) above are inappropriate. Ions are no longer effectively screened from their surrounding, and multiple particle collisions and collective effects become important. The “classical” approach to calculating transport properties of ions using Boltzmann’s equation becomes invalid. In the Sun, at number for particles per Debye sphere is only about 1.4 at the base of the convection zone, and 2.2 at the core, so the “diluite gas” approximation fails.”

Finally, let us recall that the introduced equilibrium distribution contains the factor $`e^{(1q)/2(E/kT)^2}`$ where $`q`$ is the Tsallis parameter. In the solar core plasma we have found for it the value $`|(1q)/2|0.003`$ ($`q=1`$ corresponds to the Boltzmann-Gibbs statistics and Maxwellian distribution). In spite of the very small deviation of the behavior at high energy of the Tsallis distribution respect to the standard curve the effect on the reaction rates (depletion or enhancement) is quite important .


format: text

text:
# Perfect Gauge Actions on Anisotropic LatticesWork supported in part by Schweizerischer Nationalfonds Presented by P. Rüfenacht

## 1 INTRODUCTION

Measuring heavy states (e.g. glueballs) on moderately large lattices is feasible using anisotropic lattices whose resolution in temporal direction is finer than the one in spatial directions . Having rather coarse (spatial) lattices calls for the use of improved actions; a radical choice there is the classically perfect FP action  which has been constructed for (isotropic) SU(3) gauge theory . Our new parametrization  can serve as a basis for the construction of the corresponding anisotropic action as will be shown later. Anisotropic classically perfect gauge actions are an independent approach for measuring heavy states in pure gauge theory and thus probe universality.

## 2 CREATING PERFECT ACTIONS

In this paper we consider SU(N) pure gauge theory in $`d=4`$ Euclidean space on a periodic lattice, the numerical work is done for SU(3).

For asymptotically free theories, the classically perfect FP action is determined by a saddle point equation :

$$𝒜^{\text{FP}}(V)=\underset{\{U\}}{\mathrm{min}}\left\{𝒜^{\text{FP}}(U)+T(U,V)\right\},$$
(1)
where $`T(U,V)`$ is the blocking kernel that defines the block transformation connecting the initial (fine) configuration $`U`$ and the blocked (coarse) configuration $`V`$.

In order to get a FP action for coarse configurations occuring in MC simulations, eq. (1) is used as a recursion relation connecting an action on a coarse configuration $`V`$ to one on a finer configuration $`U`$. The recursion step is iterated $`34`$ times using some appropriate simple starting action $`𝒜(U)`$ on the smoothest configuration. In principle the final FP action contains infinitely many couplings and thus has to be parametrized such that it may be used in a numerical simulation.

## 3 GETTING ANISOTROPIC

We pursue two different possibilities of modifying the method sketched in the previous section to obtain an *anisotropic* action. Firstly, the procedure may be repeated with an anisotropic starting action for smooth configurations, which may be a simple anisotropic discretization of the continuum action. Due to the isotropic nature of the blocking the anisotropy is relayed up to the level of coarse configurations.

Secondly, we can start with the parametrized isotropic FP action on coarse configurations and perform one or several additional *purely spatial* blocking steps each doubling the anisotropy $`\xi =a_s/a_t`$ (if one is working with scale-2 blocking steps) such that one obtains actions with anisotropies $`\xi `$=2, 4, 8, $`\mathrm{}`$. These actions are not *Fixed Point* actions in the narrower sense as the action on the l.h.s. of eq. (1) has got an anisotropy twice as large as the one on the r.h.s. However, the actions are still classically perfect. In particular, they possess exactly scale invariant lattice instanton solutions.

## 4 SCALAR FIELDS, PERTURBATIVE TESTS

For scalar fields it can be proven analytically that both methods mentioned in the last section are completely equivalent and yield the same result. Furthermore it can be shown that the couplings for anisotropic actions on scalar fields can be made short-ranged in all directions, this has already been studied earlier in .

In addition, we study gauge fields in the quadratic approximation, where the action is

$$𝒜=\frac{1}{2N_c}\underset{n,r}{}\rho _{\mu \nu }(r)\text{Tr}[A_\mu (n+r)A_\nu (n)]$$
(2)
neglecting terms of order 3 and higher in the gauge potential. Fig. 1 displays the nice spectra obtained for different anisotropies. Also, in the perturbative potential the anisotropic classically perfect actions perform as well as the isotropic FP action.

## 5 ANISOTROPIC STARTING ACTION

Repeating the isotropic FP action program  with an anisotropic starting action has got the advantage that the recursion steps are exactly the same as for the isotropic action, i.e. exactly the same RG transformation may be used. Furthermore, one may construct perfect actions of any desired anisotropy $`\xi `$. However, the procedure is costly as one has to perform again several RGT and fitting steps for each value of $`\xi `$. Besides, in each step it has to be carefully checked that the action parametrizations used are accurate on the respective fluctuations and that the renormalization of the anisotropy is small. That is why we prefer to perform a small number of purely spatial blocking steps starting with the isotropic action parametrized in .

## 6 SPATIAL BLOCKING STEP

To perform a purely spatial blocking step we modify the symmetric blocking kernel :

$$\begin{array}{c}T(U,V)=\hfill \\ \hfill \frac{\kappa }{N_c}\underset{n_B,\mu }{}\left(\text{Re}\text{Tr}[V_\mu (n_B)Q_\mu ^{}(n_B)]𝒩_\mu ^\beta \right),\end{array}$$
(3)
where $`Q_\mu (n_B)`$ is a mean of products of fine links that connect two coarse lattice sites (we use $`Q_\mu `$ of type III blocking, defined in ). Thus, the fine links are first smeared and subsequently blocked to the coarse lattice. The normalization term $`𝒩_\mu ^\beta `$ assures that the partition function remains unchanged under the corresponding RG transformation. The parameter $`\kappa `$ determines the weight of the blocking kernel relative to the action value $`𝒜(U)`$ and thus the freedom of the fine link products to fluctuate with respect to the coarse links.

We split the sum in $`T(U,V)`$ into a spatial and a temporal part with different $`Q_\mu (n_B)`$ and values of $`\kappa `$. Whereas $`Q_i(n_B)`$ is kept the same for the spatial directions $`i=1,2,3`$, $`Q_4(n_B)`$ now connects two coarse lattice sites whose (temporal) lattice distance is the same as on the initial fine lattice, i.e. the links are only smeared but not blocked. The parameters $`\kappa _t`$ and $`\kappa _s`$ are independent of each other and can be chosen differently to make the resulting action accurately parametrizable with a small number of parameters.

## 7 ANISOTROPIC PARAMETRIZATION

To parametrize the anisotropic actions we generalize our “fat link” parametrization described in :

$$𝒜[U]=\frac{1}{N_c}\underset{x}{}\underset{\mu <\nu }{}\underset{k,l}{}p_{kl}u_{\mu \nu }(x)^kw_{\mu \nu }(x)^l,$$
(4)
$`u_{\mu \nu }`$ $`=`$ $`\text{Re Tr}(1U_{\mu \nu }^{\text{pl}}),`$
$`w_{\mu \nu }`$ $`=`$ $`\text{Re Tr}(1W_{\mu \nu }^{\text{pl}}),`$ (5)
where $`U_{\mu \nu }^{\text{pl}}`$, $`W_{\mu \nu }^{\text{pl}}`$ denote the plaquettes built of standard and smeared gauge links, respectively.

Firstly, the coefficients $`p_{kl}`$ are chosen differently for plaquette orientations parallel ($`p_{kl}^{\text{tm}}`$) and perpendicular ($`p_{kl}^{\text{sp}}`$) to the temporal direction. Secondly, the construction of the smeared plaquettes is altered: There are three different kinds of smeared links contributing to a smeared plaquette $`W_{\mu \nu }^{\text{pl}}`$:

temporal links (always contributing to a temporal plaquette), spatial links that contribute to a spatial plaquette and spatial links that contribute to a temporal plaquette.

In addition, there are four kinds of staples contributing to a smeared link which are “asymmetric” either because they lie in the same plane as the plaquette that is going to be built or because they go in time direction:

Because of this we introduce four different asymmetry parameters $`\eta `$ (weights of the staples) and for each iteration level $`i`$ three different parameters $`c_i`$ (specifying the expansion of the smeared links) . These additional parameters do not influence the overhead of the action significantly, however the time used for parametrizing increases considerably.

Tests with reduced sets of actions (smaller number of non-zero $`\eta `$ parameters or only one $`c_i`$ value for all plaquette orientations) show that in order to fit well the perfect action the sophisticated parametrization is indeed needed.

## 8 THE $`\xi =\mathrm{𝟐}`$ ACTION

The quadratic approximation suggests a ratio of $`\kappa _t/\kappa _s=\xi ^2`$ for the $`\kappa `$-parameters in the anisotropic blocking kernel. As the blocking is unchanged in spatial direction, $`\kappa _s`$ is kept at the optimal value of $`\kappa _s=8.8`$ obtained for the isotropic case; thus $`\kappa _t=35.2`$. Fits on 40 configurations are performed and yield an action with 3 levels of $`c_i`$, $`0<k+l4`$ for $`p_{kl}^{\text{sp}}`$ and $`0<k+l3`$ for $`p_{kl}^{\text{tm}}`$. The fit is stable as the result does not change considerably when adding further independent configurations, and extending the parameter set by increasing the number of $`c`$-levels or the number of linear parameters does not further decrease $`\chi ^2`$ significantly.

## 9 THE RENORMALIZATION OF THE ANISOTROPY

The anisotropy of the above action is measured using the torelon dispersion relation as described in . On a $`12^2\times 6\times 30`$ lattice at $`\beta =3.3`$ with periodic boundary conditions we measure the Polyakov line around the short spatial direction on APE-smeared configurations of smearing levels 3, 6, 9, 12, 15 with smearing parameter $`\lambda _s=0.1`$ (see ). 23600 alternating Metropolis and over-relaxation sweeps and 4720 measurements are performed. A fully correlated fit using variational techniques yields for the renormalized anisotropy $`\xi _R=1.949(28)`$; the renormalization of the input anisotropy $`\xi _0=2`$ is thus very small. A rough estimate of the scale yields $`a_s=0.153(8)`$ fm and thus $`a_t=0.079(5)`$ fm. Fig. 2 shows the lower part of the dispersion relation including the result of the correlated fit in the range $`p^2=\mathrm{1..4}`$. Taking into account larger fitting ranges up to $`p^2=\mathrm{0..8}`$, we obtain results that are all consistent with the quoted value of $`\xi _R`$ within statistical errors.

## 10 CONCLUSIONS AND PROSPECTS

The $`\xi =2`$ classically perfect gauge action has been constructed using one spatial blocking step starting with the isotropic FP action. It shows a very small renormalization of the anisotropy at $`\beta =3.3`$ corresponding to $`a_s0.15`$ fm. Scaling tests will be performed on this action measuring the static quark–antiquark potential (including off–axis separations to probe rotational invariance) and glueball masses. Furthermore, the spatial blocking procedure sketched above will be repeated using the $`\xi =2`$ parametrized classically perfect action to build a perfect $`\xi =4`$ action.

Acknowledgements: We would like to thank Ferenc Niedermayer for useful suggestions and discussions.


text:
# PROBING PORES USING ELEMENTARY QUANTUM MECHANICS

## Abstract

The relaxation of polarized spins in a porous medium has been utilized as a probe of its structure. We note that the governing diffusion problem has a close parallel to that of a particle in a box, an elementary Quantum mechanics toy model. Following the spirits of “free electron” model, we use generic properties of the eigen spectrum to understand features common to a wide variety of pore geometry, consistent with large scale numerical simulations and experimental data.

Keywords: Magnetic Resonance; Pores; Diffusion; Quantum Mechanics.

INTRODUCTION

Relaxation of polarized spins by paramagnetic impurities embedded in pore matrix provides a valuable probe into its structure, especially when elaborate techniques such as MRI imaging are not easily applicable. Exploration for hydrocarbon inside a bore hole is one such case where it is critical to extract limited, but pertinent information in a timely manner from measurements often done with a hand tied. Over the past years, a lot of progress has been made both in development of magnetic resonance tools1 as well as theoretical aspects2 geared specifically toward extracting surface-to-volume ratio from the observed magnetization evolution. There are several more sophisticated MR techniques now being routinely applied to rocks in the laboratories and thus potentially much more detailed information on the structure of the pore space is available.3 In addition to the pore structure which controls practically important properties such as ite permeability, one is often interested in the details of pore filling fluids or gas such as how is is made up out of water and oil, diffusion constants of its constituents, and wetting condition. As our scope begins to embrace this degree of complexity, it becomes essential to bring some concepts up to date since sorting out the complexity of the pore matrix from that of its content is a minimum requirement in such endeavor.

Although we are specifically motivated by MR applied to porous media, the underlying diffusion phenemenon under complex boundary conditions occurs in a variety of real life situations. To illustrate a useful point, let us consider an extreme case: zombies wandering on top of a shopping mall whose fence has occasional openings at the edge. Depending on the size of the space, the speed of a zombie, and the frequency of the openings, the depletion of its population occurs either in a homogeneous way or with more interesting profile in space. The conventional approach sans imaging relied on deducing these parameters from observing how the overall population decays.

QUANTUM MECHANICS

The analysis is simple if we take the continuum approach, treating the population as a scalar variable, subject to the classic diffusion equation.4,5 The diffusive current $`j_D`$ arises from inhomogeneous density distribution $`m`$ via $`j_D=Dm`$. Conservation of local flux $`j_D=\frac{d}{dt}m`$ then leads to $`D^2m=\frac{d}{dt}m`$ with the boundary condition $`\widehat{n}j_d=\rho m`$ on the wall. Treating it as a boundary value problem, the problem is solved in terms of eigenmode expansion, as was studied in detail by Brownstein and Taar for simple geometries years earlier. In this work, we employ an analogy to Quantum mechanics by observing that $`D\frac{\mathrm{}}{2m}`$ and $`tit`$ turns the equation for $`m`$ into the Schrödinger eqaution for a particle confined to a leaky container with $`=\frac{\mathrm{}^2}{2m}^2`$. $`m`$ is identified with the wave function $`\mathrm{\Psi }`$ of the QM particle in an imaginary time, but we note that the $`m`$ is directly observed in our classical system, rather than $`|\mathrm{\Psi }|^2`$. The Hamiltonian operator is simply that of kinetic energy of the particle and the boundary condition is of mixed type, $`[\widehat{n}D+\rho ]m=0`$. The existence, orthonormality of eigenmodes are guarranteed with $`\mathrm{\Psi }_n`$ satisfying $`\mathrm{\Psi }_n=E_n\mathrm{\Psi }_n.`$ By multiplying both sides with $`\mathrm{\Psi }_n^{}`$ and integrating over the pore volume, with the help of the boundary condition, one can show

$$\left(\frac{\mathrm{}^2}{2m}d^3r|\mathrm{\Psi }_n|^2+\rho \mathrm{}|d𝐬||\mathrm{\Psi }_n|^2\right)=E_n,$$
(1)
which shows that $`E_n0`$ with equality allowed only for $`\rho =0`$ with a uniform mode $`\mathrm{\Psi }_0=\frac{1}{\sqrt{V}}`$. The energy (relaxation rate) is broken down into two components: the kinetic energy part and the potential energy part which is localized in the surface region. For each eigenmode, we have a length scale $`\mathrm{}_n=\frac{|d𝐬||\mathrm{\Psi }_n|^2}{{\scriptscriptstyle d^3r|\mathrm{\Psi }_n|^2}}`$ with which we have

$$E_n=\frac{\mathrm{}^2}{2m}d^3r|\mathrm{\Psi }|^2(1+\frac{2m\rho }{\mathrm{}}\mathrm{}_n).$$
(2)
Since $`\mathrm{}_n\frac{𝒮}{V}\frac{1}{k_n^2}`$ decreases with increasing $`n`$, this clearly shows how the high (excited) modes become progressively independent of $`\rho .`$

From now on, we introduce the bra-ket notation, $`\mathrm{\Psi }_n(𝐫)=<𝐫|n>`$ for compactness. The eigenmodes form a complete set with $`|n>_{n=0}^{\mathrm{}}<n|=1`$ and the time evolution of an arbitrary state $`|\mathrm{\Psi }_i>`$ is given by

$$|\mathrm{\Psi }_t>=\underset{n}{}|n>e^{\frac{E_n}{\mathrm{}}t}<n|\mathrm{\Psi }_i>.$$
(3)
This is all very formal. The hard part lies in the specific spectrum for a given boundary shape and the strength of the sink. The distribution of eigenvalues encodes much information of the boundary shape, i.e. pore structure.6 It is possible that there exists significant degeneracy, since different shapes may lead to an identical spectrum. Despite this, we anticipate that we can relate some characteristics of the eigenspectrum to a given shape of the pore wall. Borrowing the concepts from random matrix theory and Quantum chaos, we may study the correlation among the levels via evaluating $`𝒫(ϵ_{nm})`$ where $`ϵ=E_nE_m.`$ This function generally displays a universal characterisitc for large $`E_n`$’s, while at small values (correlation among slow diffusion modes) unique fingerprint of the boundary structure manifests itself. Note that $`\rho `$ acts as a knob taking the system from a Neunmann ($`\rho 0`$) to a Dirchlet ($`\rho \mathrm{}`$) boundary condition. This variation in sink strength largely affect the lowest eigen mode in a non-trivial manner while for higher levels, amounts only to a uniform phase shift by $`\frac{\pi }{2}`$. The effect of a large scale pore shape on the spectrum can be most easily appreciated in the model system of particle in a rectangle. There the spectrum follows ( in the $`\rho 0`$ limit for simplicity):

$$\frac{E_𝐧}{\mathrm{}}=\frac{1}{T_𝐧}=\frac{\mathrm{}\pi ^2}{2m}[(\frac{n_x}{L_x})^2+(\frac{n_y}{L_y})^2+(\frac{n_z}{L_z})^2].$$
(4)
One can easilty obtain the density of states from the surface area of an ellisoid in the $`𝐤`$ space as in the free electron model. Some representative spectra are shown in Figure 2 in which we varied the aspect ratio of a rectangle while keeping its volume fixed at $`L^3=1`$ and $`\rho =0.`$ Qualitatively, the shape of the bouding boxes represents a 1D-, 2D- and 3D space respectively and the corresponding spectrum displays a characteristic distribution. As mentioned earlier, what is probably more meaningful in respect of the boundary shape is the auto-correlation function of the level positions. Practically, this requires that a given measurement picks out the full spectrum with extremely high resolution. Unfortunately, this step is highly compromised by the inherent limitations of inverse Laplace transform with which one transforms the time-domain data into the spectrum. In addition to this, the preparation and detection states usually have high degree of symmetry, thus picking out only a subset of the full spectrum.

APPLICATION

Usually, the final state is detected by measuring the net magnetization. This is equivalent to projecting the final state onto a uniform state. If one also prepared a uniform initial state, $`|\mathrm{\Psi }_i>=|\mathrm{\Psi }_0>,`$ then we obtain

$$(t)=<\mathrm{\Psi }_0|\mathrm{\Psi }>_t=\underset{n}{}e^{\frac{E_n}{\mathrm{}}t}|<n|\mathrm{\Psi }_0>|^2.$$
(5)
For small $`\rho `$, we know from Equation 1 that the lowest eigenmode becomes close to the uniform mode $`\mathrm{\Psi }_0\frac{1}{\sqrt{V}}`$, with $`E_0\rho S/V`$ following directly from the equation. Furthermore, from orthogonality condition, we expect $`<n|0>=\delta _{n0}`$ and therefore, $`(t)e^{\frac{E_0}{\mathrm{}}t}.`$ This immediately shows the limitations of detecting uniform magnetization as a probe for multiple length scales in a given pore. Even for determination of the average surface-to-volume ratio, $`\frac{𝒮}{V}`$ from the observed $`E_0\rho \frac{𝒮}{V}`$ is problematic unless the value of $`\rho `$ is independently and accurately known since the often dominant, lowest eigenmode is extremely sensitive to $`\rho `$ and the all higher modes which remain robust as $`\rho `$ is varied, unfortunately do not show up in measurements. To illustrate this quantitatively, we show in Figure 3 partial spectrum for a spherical pore of radius $`a`$ as we vary the control parameter $`\rho a/D`$ from $`0.05`$ to $`5`$. The general eigenmodes for a sphere is given by $`\mathrm{\Psi }_{k,l,m}=j_l(kr)Y_l^m`$ with $`l=0,1,\mathrm{}`$ and $`m=l,\mathrm{}0,\mathrm{}l.`$ For given $`l`$ and $`m`$, there is a manifold of eigenmodes with $`k`$ determied from

$$ka\frac{lj_{l1}(ka)(l+1)j_{l+1}(ka)}{2l+1}+\frac{\rho a}{D}j_l(ka)=0.$$
(6)
Usually, only the subset with $`l=m=0`$ is considered on symmetrical grounds with isotropic initial and detected states $`|\mathrm{\Psi }_0>`$, but for more general cases (See discussion below), full angular variation should be retained even with the isotropic boundary condition. The most salient feature in the figure is the extreme sensitivity of the lowest eigenmode on $`\rho `$. The higher modes, including the lowest modes in each $`l0`$ manifolds, are less sensitive and their positions give direct information on the size of the sphere $`a`$. Another important characteristics of the higher modes is that they generally lead to a universal density of states, a la free electron model, which reflects the fact that the most pore walls look very similar on microscopic length scales, and the fast decay modes are insensitive to the large scale variation of the pore wall.

There are ways to get around the barrier to the excited states. One way is to create and detect a localized state. This extreme case begins with a localized polarization $`<𝐫|\mathrm{\Psi }_i>=\delta (𝐫𝐫_0)`$, and detecting again at the same spot, $`<\mathrm{\Psi }_i|`$. The observed signal will follow:

$$_{loc}(𝐫_0)=\underset{n}{}e^{t/T_n}|<n|𝐫_0>|^2$$
(7)
with the spectral weight $`|<𝐫_0|n>|^2`$ for each mode democratically distributed over a broad range. The two cases considered so far are two extreme limits where the prepared and the detected states are oriented in the Hilbert space with some overlap with the (uniform) ground state of the given boundary value problem, as schematically depicted in Figure 4. The most general situation may be described in terms of a unitary rotation of the uniform state by the application of a local phase encoding-decoding operation, $`e^{i\alpha (𝐫)}`$:

$$_\alpha =\underset{n}{}e^{t/T_n}|<n|e^{i\alpha (𝐫)}|\mathrm{\Psi }_0>|^2.$$
(8)
This operation can indeed be realized in the laboratory via application of $`90`$\- and $`180`$-degree rf pulses, a common MR techque. The critical ingredient for this to work is the presence of internal field inhomogeneity on sub-pore length scales. In the presence of a uniform applied field along the z-axis, the small magnetic susceptibility contrast between the matrix and its filling fluid induces a local field proportional to $`\mathrm{}\chi .`$ Therefore, one may rotate the spins onto the plane transverse to the polarizing field, and let them precess under the local field $`h(𝐫)`$ for a given time $`\tau `$, accumulating a phase $`\varphi (𝐫)=\gamma h(𝐫)\tau `$.

The spins are then rotated back, orthogonal to the transverse plane and allowed to relax over time $`t`$ before one apply the reversing process. This winding-unwinding procedure would work perfectly if the molecules do not diffuse. However, the wound (or encoded) state now evolves with $`|n>e^{t/T_n}<n|`$ with significantly enhanced overlap with all excited states $`|n>`$ with $`n>0.`$ One can immediately see a potential benefit if the modulating field $`h(𝐫)`$ closely follows the underlying geometry: The rotation of $`|\mathrm{\Psi }_i>`$ will occur in such a way that the maximal overalp with the class of eigenmodes with length scales matching the variation of the boundary. One can qualitatively understand this by noting that the local field can be described as arising from an inhomogeneous distribution of surface dipoles whose strength is proportional to $`𝐡\widehat{n}`$. Therefore, the local field tends to get enhanced in the vicinity of the pore wall, and is sensitive to its local orientation. The encoded pattern created by flashing of the local field (in practice, the spins are momentarily rotated to feel its presence) therefore reflects the contour of the pore wall in a convoluted manner. Notwithstanding some detailed issues, this is quite sufficient in bringing out the full eigenspectrum into light.

Finally, we note that the widely studied $`D(t)|𝐫(t)𝐫(0)|^2`$ is analogous to the polarizability $`_n<\mathrm{\Psi }_f|\widehat{𝐫}|n>e^t<n|\widehat{𝐫}|\mathrm{\Psi }_i>`$ which for a rectangle leads to $`\frac{1}{12}(L_x^2+L_y^2+L_z^2)+\frac{1}{\pi ^2}_{n>0}_ie^{t\frac{D_0\pi ^24n^2}{L_i^2}}\frac{L_i^2}{n^2}(1)^n.`$ This, along with generic properties of the excited states, leads to an analytic for $`D(t)`$ which covers a wide range of $`t`$ turning it into an effeftive probe into the pore shape as will be reported elsewhere.

To be concrete, we are carrying out large scale numerical simulations of spin diffusion in various pore structures such as random glass bead pack and real rocks obtained with 3D microtomographic technique.8 The complex 3D pores are represented as coupled cluster of micron length scale cubes in which a large number of random walkers (typically $`10^410^5`$ in numbers) are deployed to simulate a variety of MR measurements. For $`_\alpha ,`$ we obtain the internal field by numerically evaluating dipole contributions from all surrounding chambers. Looking at the results, an example of which is shown in Figure 5, one can indeed verify the claim made earlier on the connection between the pore shape and the local field configuration.

Figure 6 shows a typical simulation results performed on the 3D microtomographic data of Fontainebleu sample. The result for $``$ shows excellent agreement with experimental data on the same class of rocks. Furthermore, the simulated $`_\alpha `$ excavates the length scales hidden in $``$ as shown in the inset for $`|<n|e^{i\alpha (𝐫)}|\mathrm{\Psi }_0>|^2`$. The presence of multiple length scales manifests itself in apprearence of the peaks, and their spectral weight depends on the strength of $`\alpha \mathrm{}\chi B\tau `$. More detailed numerical studies on the variety of pore structure is under progress and will be reported elsewhere.


text:
# Infall models of Class 0 protostars

## 1 Introduction

An important aim of the study of young stellar objects (YSOs) has been to classify them along an evolutionary sequence. According to the widely used classification scheme of Lada & Wilking (1984) and Lada (1987), YSOs are divided into three classes based on the slopes of their infrared spectral energy distribution (SED). The reddest objects, dubbed Class I sources, are thought to be protostars surrounded by infalling envelopes. Their SEDs are consistent with dust emission from envelopes with infall rates near those predicted by the standard theory of protostellar formation, which invokes free-fall collapse from an initially thermally-supported, isothermal, singular molecular cloud core in hydrostatic equilibrium (Adams & Shu 1986; Adams, Lada & Shu 1987). Kenyon, Calvet & Hartmann (1993, hereafter KCH), in particular, were able to match the observed SEDs of Class I sources in Taurus-Auriga with rotating infall models of the type described by Terebey, Shu & Cassen (1984). This agreement between theory and observation is generally seen as testimony to the success of the standard model of protostellar collapse.

Recently, Andr$`\stackrel{´}{\mathrm{e}}`$, Ward-Thompson & Barsony (1993, hereafter AWB) identified YSOs characterized by blackbody-like SEDs that peak at submillimeter wavelengths. The low dust temperatures ($`T_d`$ $``$ 20K) and large dust masses inferred for these sources led AWB to suggest that these objects constitute a new class of objects, called “Class 0”, which might be even younger than typical Class I objects. Observationally, the defining characteristics of Class 0 sources include a high ratio of submillimeter to bolometric luminosity (L$`{}_{}{}^{\lambda >350\mu m}{}_{submm}{}^{}`$/Lbol $`>`$ 5 $`\times `$ 10-3), invisibility at near-infrared wavelengths, and the presence of molecular outflows (Barsony 1994). Andr$`\stackrel{´}{\mathrm{e}}`$ & Montmerle (1994) argued that Class 0 sources have more circumstellar mass than stellar mass. In this picture, Class 0 sources have not yet assembled the bulk of their final stellar mass and are still in the main accretion phase. Consequently, the authors argue that Class 0 sources are the youngest protostars, and subsequently evolve into Class I sources (see André, Ward-Thompson, & Barsony 2000 $`=`$ AWB2).

Following AWB’s identification of the prototype Class 0 source VLA 1623 in the rich star-forming core $`\rho `$ Ophiuchi A, several other objects of the class have been found in Orion B (Ward-Thompson et al. 1995), Serpens (Hurt & Barsony 1996), Cep E (Lefloch et al. 1996), and Perseus (O’Linger et al. 1997) molecular clouds. In a spectroscopic survey of 47 candidate protostars, Mardones et al. (1997) report that Class 0 sources have stronger infall signatures –i.e., larger negative velocity shifts– than Class I sources, and attribute that result to a difference in the physical conditions in the circumstellar envelopes of the two source types.

In order to understand better the observed characteristics of this ‘youngest’ population of YSOs, we have performed detailed radiative transfer calculations of dust envelope emission in Class 0 sources, following the methods described by KCH and Calvet et al. (1994). We also study the observed radial intensity profiles of the spatially-resolved prototype Class 0 source VLA 1623 at three submillimeter wavelengths. In both cases, we compare the observations to the predictions of the standard model, and explore the possibility that environments exterior to the infall region can affect the spectral energy distribution and the inferred density gradients.

## 2 Standard Model

### 2.1 Methods and Parameters

Following KCH, we adopt the Terebey, Shu & Cassen (1984, hereafter TSC) solution for a uniformly rotating protostellar cloud with a radial density distribution initially approaching that of a singular isothermal sphere, $`\rho (r)r^2`$. However, unlike TSC, we assume an outer boundary to the envelope, $`r_{out}`$. Inside the collapse region, the infall velocity approaches the free fall value, $`vr^{1/2}`$, and the density distribution converges toward $`\rho r^{3/2}`$, because the mass infall rate, $`\dot{M}_i`$ is constant. Inside the centrifugal radius, $`R_c`$, material falls onto a circumstellar disk rather than radially onto the central object.

We use the angle-averaged density distribution to construct the spherically symmetric source function, $`S_\nu =\sigma _\nu J_\nu +(1\sigma _\nu )B_\nu `$, where $`\sigma _\nu `$ is the albedo and $`B_\nu `$ is the Planck function. Scattering increases the flux at short wavelengths, $`\lambda 3\mu `$m. We then derive the emergent flux from the envelope using $`S_\nu `$ for the source function and the angle-dependent density, $`\rho (r,\theta )`$, for the opacity.

The angular momentum of the infalling material modifies the density distribution from the purely-radial infall solution at $`rR_c`$. However, for the models considered here, the observable emission originates from envelope radii much larger than typical disk sizes $`100`$ AU, and so the value of $`R_c`$ is unimportant. Similarly, we do not include envelope cavities such as those produced by bipolar outflows. KCH pointed out that envelope cavities were needed to explain the near-infrared emission of some Taurus Class I sources; but the very large optical depths characteristic of the models considered here prevent the escape of any significant amount of near-infrared scattered light from the envelope.

Thus, our models consist of a central stellar source, with bolometric luminosity *L*, surrounded by an envelope that is essentially spherically-symmetric, with inner radius $`r_{in}`$ and outer radius $`r_{out}`$. The inner radius of the envelope is set by the dust destruction radius, which is assumed to occur at $`T`$ 1500K. For $`r_{out}`$, we adopt values roughly equivalent to the observed sizes of the sources being modeled. The models have then three free parameters: the source luminosity *L*, $`r_{out}`$, and the density of the envelope (see below). We used Draine & Lee (1984) opacities for $`\lambda 200\mu m`$ and a shallower $`\beta =`$1.5 (where $`\kappa _\nu \nu ^\beta `$) opacity law for longer wavelengths.

### 2.2 Model Results

#### 2.2.1 Spectral energy distribution

We began by constructing a grid of models with a wide range of envelope densities for source luminosities consistent with the observed values for Class 0 sources of interest. These models were defined in terms of $`\rho _1`$, the density at *r*= 1 AU in the limit that $`R_c`$ 0. The mass infall rate is related to this density parameter and the central mass by

$$\dot{M}=1.9\times 10^6M_{\mathrm{}}\mathrm{yr}^1\left(\frac{\rho _1}{10^{14}\mathrm{g}\mathrm{cm}^3}\right)\left(\frac{M_{}}{1M_{\mathrm{}}}\right)^{1/2}$$
(1)
The envelope temperature distribution is not very sensitive to the SED of the central source because its radiation is immediately absorbed by the innermost (i.e., large optical depth) regions of the envelope and reradiated at longer wavelengths.

These models show general agreement with the analytic approximation of KCH for the peak wavelength of the SED:

$$\lambda _m60\left(\frac{L}{10^{33}\mathrm{erg}\mathrm{s}^1}\right)^{1/12}\left(\frac{\rho _1}{10^{13}\mathrm{g}\mathrm{cm}^3}\right)^{1/3}\mu \mathrm{m}.$$
(2)
For typical source luminosities and expected central masses, in order to shift the peak of the SED to $`\lambda 200\mu m`$, as is observed in Class 0 sources, the models require infall rates $`10^4M_{\mathrm{}}yr^1`$.

Unfortunately, SEDs of most known Class 0 sources are not well defined observationally. Many of them have been observed only at a couple of sub-millimeter wavelengths, and the IRAS upper limits often do not provide strong constraints in the mid- and far-infrared. Therefore, we restrict our detailed modeling to two sources with relatively well-defined SEDs:

VLA 1623 Located in the $`\rho `$ Oph A molecular cloud, VLA 1623 is a centrally-peaked source of millimeter emission (Andr$`\stackrel{´}{\mathrm{e}}`$ et al. 1990b), a weak 6 cm radio source (Leous et al. 1991), and the driving source of a strong CO bipolar outflow (Andr$`\stackrel{´}{\mathrm{e}}`$ et al. 1990a). The highly collimated outflow has a mechanical luminosity of $`L_{flow}0.52L_{\mathrm{}}`$ and a short dynamical timescale $`t_D30006000`$ yrs, making it one of the youngest known flows. VLA 1623 is undetected at near- and mid-infrared wavelengths to very low limits. Based on these properties as well as its remarkably cold temperature ($`T`$ 20K), high internal obscuration, apparently massive circumstellar structure and high $`L_{submm}/L_{bol}`$ ratio, AWB proposed it as the prototype of the newly-defined Class 0 sources. Observations with the JCMT-CSO interferometer show an unresolved central source –presumably a circumstellar disk– in VLA 1623. However, its contribution to the total flux at 1.3mm and 800$`\mu `$m is $``$ 10% (Pudritz et al. 1996).

In our infall model for VLA 1623, we set $`L=1L_{\mathrm{}}`$ and $`r_{out}`$=2000AU, which corresponds to 12” at the distance of Ophiuchus and is roughly the observed source size in the sub-millimeter. (We adopt $`d=160`$ pc to Oph, the same as AWB.) For the model that best fits the observed SED of VLA 1623 (Fig. 1a), the mass infall rate is $`2\times 10^4M_{\mathrm{}}\mathrm{yr}^1`$ for a central stellar mass of 0.5$`M_{\mathrm{}}`$. A lower infall rate model, with $`\dot{M}10^5M_{\mathrm{}}\mathrm{yr}^1`$, shown as a dotted line in Fig. 1a, does not match the data. Using the same parameterization of opacity as AWB – $`\kappa _\nu =0.1\mathrm{cm}^2\mathrm{g}^1(\nu /1000GHz)^{1.5}`$– we also obtain a total envelope mass of $`0.6M_{\mathrm{}}`$ for the best-fit model.

HH24MMS This bright source of millimeter emission near the Herbig-Haro object HH24 in the Orion B molecular cloud at a distance of 460 pc was discovered by Chini et al. (1993). It is not detected in the near- or far-infrared, although IRAS limits are rather poor because of source confusion. VLA observations have revealed a 3.6cm source coincident with the millimeter peak (Bontemps, Andr$`\stackrel{´}{\mathrm{e}}`$ & Ward-Thompson 1995). Using JCMT observations at 350, 450, 800 and 1100$`\mu `$m, Ward-Thompson et al. (1995) have produced a reasonably good SED for HH24MMS, even though its peak is not tightly constrained due to a lack of shorter wavelength detections.

We modeled HH24MMS with $`L=2.5L_{\mathrm{}}`$ and $`r_{out}`$=4000AU, corresponding to an angular cloud diameter of 18”. Our best fit model requires a mass infall rate of $`3\times 10^4M_{\mathrm{}}\mathrm{yr}^1`$, assuming $`M_{}=0.5M_{\mathrm{}}`$ (Fig. 1b).

#### 2.2.2 Radial intensity profiles

Since VLA 1623 is slightly extended at sub-millimeter wavelengths, we attempted to compare our “best-fit” model profiles to the observed azimuthally-averaged radial intensity profiles at 350$`\mu `$m, 450$`\mu `$m and 800$`\mu `$m of AWB. The model output was convolved with the appropriate beam size at each wavelength and scaled to a distance of 160 pc. All the profiles, shown as dotted lines in Figure 2, have been normalized to the peak flux values. As Andr$`\stackrel{´}{\mathrm{e}}`$ et al. noted, a modified TSC infall model with $`\rho r^{3/2}`$ cannot reproduce the observed profiles, which are significantly shallower.

In our basic model, the temperature at the outer edge of the infalling cloud is not constrained and in fact drops to $``$ 7 K. Given that external heating is likely to keep the cloud surface at $``$ 20K, we explored the effect of adopting a minimum temperature in the envelope of $`T_{min}=`$ 20K. The result, shown as dashed lines in Figure 2, is to produce somewhat shallower intensity distributions, but still not shallow enough to fit the observations well.

Our calculations differ from those of Terebey, Chandler & André (1993) in that we assume an outer boundary to the envelope. The inability to determine the zero flux level in the observed intensity profiles is an additional source of uncertainty in our fits.

## 3 Contribution of exterior regions

The presence of a compact VLA source and a collimated CO bipolar outflow as well as the centrally peaked submillimeter emission strongly suggest that a central source –i.e., a protostar– has already formed in VLA 1623. Moreover, the presence of a compact ($`<70`$ AU) mm and submm source consistent with a dusty circumstellar disk (Pudritz et al. 1996) emphasizes the likelihood that substantial infall has already occurred. Thus, one would expect to find an infall region characterized by a density distribution of roughly $`\rho r^{3/2}`$ (see §4.1). In this context the $`\rho r^{1/2}`$ distribution inferred by AWB from the radial intensity profiles is surprising. 111A $`\rho r^{1/2}`$ density distribution corresponds to $`\mathrm{\Sigma }r^{1/2}`$ surface density distribution. Thus, for the slowly-varying temperature distributions characteristic of the outer envelope $`Tr^{0.4}`$, the surface brightness distribution would actually increase with increasing radius without including an outer cutoff radius; this parameter has a crucial effect on the resulting intensity distribution.

One possibility is that material outside the infall region contributes to the observed intensity profile. As clearly seen in the submillimeter continuum maps of $`\rho `$ Ophiuchi (AWB; Motte et al. 1998), VLA 1623 is indeed surrounded by –and likely embedded in– larger-scale cloud emission. We have attempted to model the exterior cloud material as a uniformly emitting “slab” at 20K. The resulting radial intensity profiles, which include contributions from the infall zone as well as the exterior region, are shown as solid lines in Figure 2. We find that in order to match the observed intensity distributions at 350$`\mu `$m, 450$`\mu `$m and 800$`\mu `$m, the “slab” needs to have a column density $`<N_{H2}>6\times 10^{22}`$ cm-2 or $`A_V60`$, implying $`M_{slab}0.3M_{\mathrm{}}`$ in the 2000 AU radius region. This compares well with typical column densities of $`10^{23}`$ in the larger $`\rho `$ Oph A cloud derived by Motte et al. (1998).

Could extinction by this outer region make a Class I source appear as a Class 0? A “slab” of the kind we propose, with $`A_V60`$, would significantly extinct emission from the inner regions at short wavelengths. However, because the extinction of such an outer region is $`A_\lambda <1`$ at wavelengths $`\lambda 30\mu m`$, it would not shift the peak of a Class I SED to $`\lambda >100\mu m`$. Therefore, it is unlikely that a Class I source with a relatively flat SED in $`\lambda `$ Fλ would appear as a Class 0 simply due to external extinction.

Our best-fit, composite model for the SED of VLA 1623 is shown in Figure 3. Using the same opacities as AWB, this object has a total mass of $`0.6M_{\mathrm{}}`$, divided roughly equally between the central infall zone and the outer, nearly constant-density region. Our composite model still requires a high infall rate of $`10^4M_{\mathrm{}}\mathrm{yr}^1(M_{}/0.5M_{\mathrm{}})^{1/2}`$, i.e., the necessary infall rate is only reduced by a factor of two by the inclusion of the constant density outer region.

Interestingly, model fits to the ‘youngest’ (i.e., most embedded) protostars in Chandler & Richer’s (2000) sample imply that the power-law index of the envelope density distribution $`p`$ (defined by $`\rho r^p`$) is lower than that predicted by the standard collapse model. They suggest that a shallow density profile may be maintained in the outer envelope by magnetic fields and/or turbulence. On the other hand, Hogerheijde & Sandell (2000) find that the density distribution in envelopes around four embedded YSOs in Taurus are fit equally well by a radial power-law index $`p=`$1.0–2.0 or by a Shu collapse model.

## 4 Discussion

### 4.1 Environment and high infall rates

In the previous section we showed that models of a collapsing isolated cloud can fit the SEDs of Class 0 source if a high infall rate is adopted, but that the resulting surface brightness distributions are inconsistent with observations, as originally found by AWB. We have extended this result to include the effects of external heating, and still find a discrepancy between theory and observation. The surface brightness distributions can be explained with simple infall models only by including an additional emitting component, which could be material in the near environment of the collapsing cloud.

Infall models generally exhibit surface density distributions $`\mathrm{\Sigma }r^{1/2}`$ (implying $`\rho r^{3/2}`$ in spherical geometry), whether from spherical or flattened clouds, similarity solutions for collapse or numerical time-dependent calculations, with or without cloud magnetization (e.g., Larson 1969; Penston 1969; Shu 1977; Foster & Chevalier 1993; Galli & Shu 1993a,b; Hartmann et al. 1994; Hartmann, Calvet, & Boss 1996; Li & Shu 1996). In isothermal clouds at the point of collapse, many calculations show density distributions which tend to steepen more towards $`r^2`$ (e.g., Mouschovias 1991, and references therein; Foster & Chevalier 1993). These models strongly suggest that the observations of Class 0 sources encompass not only an infall region, but external material not necessarily in collapse.

The question then arises: if the external environment contributes a large fraction of the sub-mm flux, can we in fact say anything about infall rates? A static cloud, with density distribution $`\rho r^{1/2}`$ as considered by AWB, could, if physically realizable, fit the observations in the absence of any infall. (Note that such a static cloud is inconsistent with the standard model of isothermal clouds, in which a $`\rho r^2`$ density distribution is expected). It seems likely that some fraction of the sub-mm emission arises from an infall region, because some infall must have occurred to produce the central source, and given that pre-stellar cores appear to be much less centrally condensed than Class 0 sources (Ward-Thompson et al. 1994; André et al. 1996; Ward-Thompson, Motte, & André 1999; AWB2).

Nevertheless, let us take the extreme limit by supposing that all of the sub-mm emission arises from a static cloud. We may then estimate the rate at which this cloud would collapse. For a cloud mass of $`0.6M_{\mathrm{}}`$ concentrated within a radius of 2000 AU, the free-fall timescale is $`2\times 10^4`$ years. Thus, if this cloud were to collapse, it would exhibit an average infall rate of $`3\times 10^5M_{\mathrm{}}yr^1`$. This estimated infall rate is not much smaller than the $`10^4M_{\mathrm{}}yr^1`$ estimated for the pure infalling envelope model based on the SED fitting. Moreover, if the region that fell in to form the outflow source was originally denser than the outer static region, as seems probable, the infall rate could be higher than this estimate. Thus, one is lead to the same order of magnitude infall rate from consideration of the spatial concentration of the local gas, independent of the fraction of emission currently arising from an external cloud rather than the infall region.

High infall rates are also suggested by observations of outflows. In a survey of nearby embedded YSOs, Bontemps et al. (1996) found that Class 0 sources have unusually powerful outflows. The Class 0 objects in their sample lie an order of magnitude above the well-known correlation between the outflow momentum flux (integrated over time) and bolometric luminosity that holds for Class I sources. Since models of bipolar outflows (Shu et al. 1994, Ferreira & Pelletier 1995) predict that the mass-loss rate is some fraction of the mass-accretion rate, the stronger outflows in Class 0s also call for significantly higher infall rates than in Class Is. Henriksen, André, & Bontemps (1997) attempt to derive accretion rates as a function of time for Class 0 sources which are consistent with infall at non-constant rates, again finding very rapid mass accumulation is required to explain outflow properties.

In addition, from a spectroscopic survey of 47 candidate protostars, Mardones et al. (1997) report that Class 0 sources have stronger infall signatures –i.e., larger negative velocity shifts– than Class I sources. One possible interpretation of their results is that Class 0s have higher density envelopes than Class Is, again implying higher infall rates in Class 0 objects.

Thus, our high infall rates, which imply that Class 0 sources are a short-lived phenomenon, are consistent with previous conclusions in the literature (e.g., AWB). A comparison of source statistics in Taurus-Auriga and Ophiuchus support this conclusion. With an infall rate of $`10^4\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$, a 0.5 $`M_{\mathrm{}}`$ star would be accreted in $``$5000 years, a timescale consistent with the dynamical lifetime of VLA 1623. André & Montmerle (1994) and Barsony (1994) estimate a lifetime of $`10^4`$ yr for Class 0 sources in Oph on statistical grounds. In addition, as there are around 200 known pre-main-sequence objects in Taurus, one would expect 5000 yrs$`/10^6`$ yrs $`\times `$ 200 or $``$1 Class 0 source; and Taurus-Auriga contains only two “border-line” Class 0 sources, L1527 and L1551 IRS 5 (Chen et al. 1995). (AWB2 suggest that two Class 0 sources are present in Taurus, L1527 and the extremely low-luminosity source IRAM 04191+1522 \[André, Motte, & Bacmann 1999\]).

One problem with the picture of high infall rates for Class 0 sources is that much higher accretion rates onto the central protostar should lead to high accretion luminosities. However, Class 0 objects are not very much more luminous than their Class I counterparts. Henriksen, André, & Bontemps (1997) and Andr$`\stackrel{´}{\mathrm{e}}`$ (1997) have suggested several factors that may reduce the theoretical accretion luminosity: in Class 0 sources, the central stellar mass might be lower, the stellar radius probably larger, and the amount of accretion energy dissipated in the wind is also likely to be larger. Another natural explanation is that material falls onto a disk, instead of falling directly onto the central star. However, the material must still somehow flow from the disk to the star at some point in the protostellar phase to build up the central star. It may be that Class 0 sources also undergo “FU Orionis-type” outbursts once enough material piles up in the disk to make it unstable (Kenyon et al. 1990, 1994).

Note also that an “external region” is not necessarily physically unconnected with the Class 0 source; i.e., this material may eventually collapse onto and be incorporated into the central star(s). Whether this is necessarily the case, however, is not clear.

### 4.2 Class 0 vs. Class I: the evolutionary model vs. initial density

AWB suggested that the Class 0 sources represent the earliest phase of protostellar collapse, based on the large amount of circumstellar material in the near environment (compared with the likely mass of the central protostellar object). Henriksen, André, & Bontemps (1997) elucidated a picture which in principle provides a theoretical basis for this evolutionary scheme.

Henriksen et al. pointed out that theoretical models of clouds initially in hydrostatic equilibrium, such as the isothermal Bonnor-Ebert sphere (e.g., Foster & Chevalier 1993), tend to have flat inner density distributions surrounded by an outer envelope of roughly $`\rho r^2`$. In addition, Henriksen et al. emphasized the observational evidence for flattened central regions in pre-stellar cores mentioned in §4.1. When the initially flat inner regions of such a cloud fall in, high mass infall rates are produced (see also Foster & Chevalier 1993); material from the outer envelope falls in later at slower rates. This leads to time-dependent infall rates, with the highest rates occuring during the earliest phases of collapse, identified as the Class 0 phase.

This picture, however, does not explain why Class 0 sources tend to be found in dense star forming regions. Moreover, while the general picture of decreasing infall rate with increasing age is quite plausible, it is not clear whether such a Class 0 phase would dominate protostellar infall, because the amount of mass in the flattened inner portions of the initial cloud is generally smaller than the mass in the envelope, implying that most of the protostellar mass must be added during the envelope infall (Class I) phase. Henriksen et al. (1997) show that, for a flat region of a cloud to contain half of the total cloud mass, the $`\rho r^2`$ envelope can only extend out a factor of 1.6 in radius beyond the flat region – essentially, the cloud is almost entirely at constant density, with only a small outer envelope region, inconsistent with observations of typical pre-stellar cores. Henriksen et al. suggest that clouds with steeper edges, $`\rho r^3`$, could alleviate this problem while maintaining the required large factor $`10`$ between Class I and Class 0 lifetimes; but they do not show that the very slow infall implied by these models for the Class I phase is consistent with SED models.

Alternatively, one might suppose that the high environmental densities characterizing the birthsites of Class 0 sources guarantee that envelope masses will be very large on small scales. In other words, the Class 0 sources might be the Class I sources of dense star-forming regions.

As a specific example, consider two protostellar clouds of the same mass but formed in regions of differing average density. The characteristic infall rates will scale roughly as the cloud mass divided by the free fall time, which in turn scales roughly as the inverse square root of the mean density; thus, $`\dot{M}M/<\rho >^{1/2}<\rho >^{1/2}R^{3/2}`$, where $`R`$ is the outer radius of the cloud. (Let us fix the central star mass in order to discuss the difference in infall rates.) If we assume that the Class 0 source has an infall rate of $`10^4M_{\mathrm{}}yr^1`$, and a Class I source has an infall rate of $`4\times 10^6M_{\mathrm{}}yr^1`$, then if the two protostellar clouds originally had the same mass, the factor of 25 in $`\dot{M}`$ implies that the Class I source would have to be larger in radius by a factor of 8.5. In other words, the mass within a 2000 AU (12 arcsec at 160 pc) radius in a Class 0 source should be compared to the mass within a radius of 16,000 AU ($`100`$ arcsec) in a Class I source to properly assess the relative evolutionary states of these objects.

The molecular surveys of Myers and collaborators (e.g., Benson & Myers 1989; Zhou et al. 1989) have indicated that many cores, such as those in Taurus, have substantial spatial extensions. For instance, even the lowest-mass ammonia cores in the survey of Benson & Myers (1989) have $`0.80.6M_{\mathrm{}}`$ of material, comparable to that of the Class 0 sources, but distributed within structures of (half-power diameters) $`10,000`$ AU $`0.05`$ pc rather than 2000 AU $`=0.01`$ pc. In many other cases the ammonia core mass is considerably larger than $`1M_{\mathrm{}}`$ and spread over larger scales. Whether a substantial portion of the material in a typical Taurus molecular cloud core actually is incorporated into a particular protostellar source is not clear –or whether it makes only one star. However, the same could be said of Class 0 sources. Neglect of the possible importance of extended material could bias the comparison between the two types of sources.

Bontemps et al. (1996) attempted a systematic assessment of Class 0 and Class I source masses from measurements of 1.3mm continuum emission. In computing dust masses they adopted a “typical” source diameter of 1 arcmin over which to integrate the emission. In doing so, Bontemps et al. recognized that they might be measuring only the “inner” envelope mass of “isolated” protostars such as those in Taurus. Moreover, if our suggestion of the presence of significant “external” material in VLA 1623 is correct, then Class 0 mass estimates might be biased by the inclusion of mass that may not end up on the protostar, especially if integrated over a radius of 30 arcsec $`4800`$ AU, rather than the 2000 AU radius adopted in the previous sections. Given the large differences in densities observed in star-forming regions, the systematic biases in mass estimates from adopting uniform source sizes can be of an order of magnitude, which could easily affect assessments of the relative evolutionary states of Class 0 and Class I sources. Unfortunately, the observational problem of deciding where a protostellar cloud “stops” is likely to be quite difficult. Yet, without consideration of this potential bias, it seems premature to conclude that all Class 0 sources must be in a radically different phase of evolution than Class I sources.

In Taurus, the source statistics of Class I sources suggest a lifetime of $`12\times 10^5`$ yr (Kenyon et al. 1990, 1994). If one adopts the infall rates of $`4\times 10^6M_{\mathrm{}}yr^1`$ typically estimated from SED fitting (ALS; KCH), then one would expect that approximately $`0.40.8M_{\mathrm{}}`$ is added to protostars during the Class I phase - which is comparable to the total mass accreted onto the average Taurus T Tauri star. It is difficult to see why the Class I period is not an important part of the protostellar phase in Taurus, unless somehow infall rates have been overestimated by an order of magnitude. If, on the other hand, Class I sources are the protostars of low-density regions like Taurus, then there is no discrepancy to explain. The two borderline Class 0 objects in Taurus (L1551 IRS 5 and L1527; Chen et al. 1995; or L1527 and IRAM 04191+1522; André et al. 1999) could be among the youngest protostellar sources in Taurus (Ladd et al. 1998); however, they might also be objects formed from particularly dense subregions within Taurus.

To summarize, we agree with AWB and AWB2 that the Class 0 sources are very young in an absolute sense because their envelopes must collapse so rapidly. However, we suggest that the high infall rates of Class 0 sources can result also from the high densities of their initial pre-stellar cloud cores rather than from always being the very first phases of collapse. Whether Class 0 sources are systematically younger than the Class I sources in a relative sense, i.e. have a much smaller fraction of their total mass accumulated into a central protostar, is far from clear given the observational biases and uncertainties in defining envelope limits. Some Class 0 sources may also be young in this relative sense, but it is difficult to support the conclusion that they are all the “youngest protostars”.

### 4.3 Dense environments and the standard picture

Infall rates on the order of $`10^4M_{\mathrm{}}\mathrm{yr}^1`$ seem to be needed to account for Class 0 sources. This infall rate is much higher than the typical infall rates $`4\times 10^6M_{\mathrm{}}\mathrm{yr}^1`$ for Class I sources in Taurus-Auriga (KCH). In the standard singular isothermal sphere model, the infall rate is $`\dot{M}c_s^3/G`$; thus the mean Taurus infall rate requires a sound speed $`c_s0.25\mathrm{km}\mathrm{s}^1`$ or a gas temperature of $`18`$ K assuming a mean molecular weight of 2.3, only slightly larger than typical temperatures found in the region. In contrast, to obtain infall rates of $`10^4M_{\mathrm{}}\mathrm{yr}^1`$ in a cloud supported only by thermal gas pressure would require unrealistically high cloud temperatures of about 150K. Thus, some other source of support must be involved if the cloud is initially in hydrostatic equilibrium.

These considerations emphasize the importance of the density of the protostellar environment to an understanding of Class 0 sources. The related question is how do such dense regions of material arise? If such regions are in hydrostatic equilibrium, turbulent motions undoubtedly play a role in supporting denser clouds (e.g., Lizano & Shu 1989; Fuller & Myers 1993). On the other hand, since the sources of turbulence are presumably non-isotropic and inhomogeneously distributed (i.e., winds, magnetic waves excited by moving dense gas clumps, etc.) arranging static situations seems very difficult. In addition, supersonic magnetic turbulence may decay rapidly (Stone, Ostriker, & Gammie 1998). An alternative view is that dense streams of gas in the ISM, however produced, collide and dissipate energy in a shock (e.g., Ballesteros-Paredes, Vázquez-Semadeni & Scalo 1999), resulting in dense regions which then undergo gravitational collapse. In such a general picture, there is no obvious reason why densities of protostellar regions need be restricted in any way to levels supportable by thermal gas pressure. In this way protostars with dense envelopes and rapid infall rates could be produced without reference to the usual constraints of thermal support (ALS). The Class 0 sources, or at least some of them, could be telling us more about protostellar formation in extreme environments than about the very earliest phases of protostellar collapse.

A related question is, why did the standard model of singular isothermal sphere collapse work so well in Taurus? We suggest that this is simply due to initial conditions. In the absence of any complicating factors –magnetic field support, turbulent motions, etc. – thermal gas pressure provides an absolute lower limit to the forces resisting gravitational collapse, especially as cosmic ray heating seems to maintain cloud temperatures at a minimum value of $`10`$ K in the absence of any other heat sources (Goldsmith 1988). Absent substantial turbulent support, or high-density flows, Taurus may only represent a limiting case of star formation.

## 5 Summary

We have described radiative transfer calculations for the spectral energy distributions of recently identified Class 0 protostars. Our models require infall rates on the order of $`10^4\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$ to account for the observed SEDs. The sub-millimeter intensity profiles of the Class 0 source VLA 1623 apparently imply a shallower density distribution than that produced by infall. But we show that material exterior to the infall region may contribute to the observed intensity profiles and thus infall models (with exterior emission) can be consistent with observations. The high infall rates in Class 0s, as required by our models and implied by a variety of observations, may be inevitable given their dense environments. Thus, we suggest that protostars in dense regions are more likely to appear as Class 0 sources while those in lower-density regions manifest Class I characteristics. The Class 0/Class I division may thus represent initial conditions more than evolutionary status, at least for some sources.

We are grateful to the referee, Dr. S. Terebey, for a thorough review. This work was supported by the Smithsonian Institution.

Figure Captions

Figure 1 - (a) Spectral energy distribution of VLA 1623 using data from AWB. (The expected contribution from free-free emission at long wavelengths has been subtracted.) The solid curve is an infall model with $`\dot{M}2\times 10^4M_{\mathrm{}}`$, and the dotted curve is a model with $`\dot{M}10^5M_{\mathrm{}}`$. (b) Spectral energy distribution of HH24MMS using data from Ward-Thompson et al. (1995). The solid curve is an infall model with $`\dot{M}3\times 10^4M_{\mathrm{}}`$.

Figure 2 - Comparison of the radial intensity profiles of VLA 1623 at 800$`\mu `$m, 450$`\mu `$m, and 350$`\mu `$m with the basic infall model in which $`\rho r^{3/2}`$ (dotted line), the infall model where $`T_{min}`$=20K (dashed line), and the composite model described in Section 3 which includes a contribution from regions exterior to the infalling core (solid line). The model output was convolved with the appropriate beam size at each wavelength and scaled to a distance of 160 pc (13.5”=2160AU, 7.5”=1200AU, and 6.3”=1000AU, respectively).

Figure 3 - Comparison of the spectral energy distribution of VLA 1623 to the best-fit composite model (solid line). The dotted line represents emission of the infall region while the dashed line shows the contribution of the exterior “slab.”


text:
# References

UNILE-CBR1

Borel Resummation of the Perturbative Free Energy of Hot Yang-Mills Theory

Rajesh R. Parwani

Departimento di Fisica, Universita’ di Lecce,

and

Istituto Nazionale di Fisica Nucleare

Sezione di Lecce

Via Arnesano, 73100 Lecce, Italy.

## Abstract

The divergent perturbative expansion of the free-energy density of thermal $`SU(3)`$ gauge theory is resummed into a rapidly convergent series using a variational implementation of the method of conformal mapping of the corresponding Borel series. The resummed result differs significantly from non-perturbative lattice simulations and the discrepancy is attributed to the presence of a pole on the positive axis of the Borel plane. The position of that pole is determined numerically and the difference between the lattice data and the resummed series is related to a phenomenological bag ‘constant’.

It is generally believed that the collision of heavy-ions at sufficiently high energies will lead to the formation of a new phase of matter, the quark-gluon plasma, and experiments to produce such a plasma are underway at Brookhaven and CERN. As the effective coupling, $`\alpha `$, of quantum chromodynamics (QCD) decreases with an increase in energy, theorists have used perturbative methods to study properties of the plasma at high temperature ($`T`$). For example, a completely analytical calculation of the free-energy density of QCD to order $`\alpha ^{5/2}`$ was performed a few years ago . The purely gluonic contribution is given by

$`{\displaystyle \frac{F}{F_0}}`$ $`=`$ $`1{\displaystyle \frac{15}{4}}\left({\displaystyle \frac{\alpha }{\pi }}\right)+30\left({\displaystyle \frac{\alpha }{\pi }}\right)^{3/2}+\left(67.5\mathrm{ln}\left({\displaystyle \frac{\alpha }{\pi }}\right)+237.220.63\mathrm{ln}\left({\displaystyle \frac{\overline{\mu }}{2\pi T}}\right)\right)\left({\displaystyle \frac{\alpha }{\pi }}\right)^2`$ (1)
$`\left(799.2247.5\mathrm{ln}\left({\displaystyle \frac{\overline{\mu }}{2\pi T}}\right)\right)\left({\displaystyle \frac{\alpha }{\pi }}\right)^{5/2},`$
where $`F_0=\frac{8\pi ^2T^4}{45}`$ is the contribution of non-interacting gluons and $`\overline{\mu }`$ is the renormalization scale in the $`\overline{MS}`$ scheme. Unfortunately (1) is an oscillatory, non-convergent, series even for $`\alpha `$ as small as $`0.2`$, which is close to the value of physical interest.

Pade Approximants (P.A.’s) were used in Ref. to resum the series (1). It was found that the dependence on the arbitrary scale $`\overline{\mu }`$ was reduced and the convergence of the series somewhat improved. However P.A.’s have a number of well-known drawbacks. For those and other reasons, in Refs. the authors abandoned the expansion of the free-energy density with respect to the coupling constant and considered instead selective resummations of gauge-invariant diagrams. Though the results of  compare favorably with lattice simulations , calculations beyond leading order are complicated and thus it seems that the issue of convergence is left open. In Ref., yet another procedure was used to study the free-energy of hot QCD: short distance perturbative effects were handled analytically while long-distance effects were described  by an effective three-dimensional theory and studied numerically.

It appears that a number of separate issues concerning Eq.(1) have become confused in the literature. The first issue is whether the given divergent series can be resummed into a convergent series, preferably in a systematic and well-motivated way. The second issue is whether such a resummed series accurately represents the physical quantity. The third issue is whether one obtains any new physical insight in the process. It will be the attempt of this paper to shed some light on these and related questions.

Recall that the divergence of perturbative expansions in quantum field theory is is a generic phenomena . Given a series

$$\widehat{S}_N(\lambda )=1+\underset{n=1}{\overset{N}{}}f_n\lambda ^n,$$
(2)
where $`\lambda `$ is the coupling constant, one expects the coefficients $`f_n`$ to grow as $`n!`$ for large $`n`$. It is natural then to introduce the Borel transform

$$B(z)=1+\underset{n=1}{\overset{N}{}}\frac{f_nz^n}{n!}$$
(3)
which has better convergence properties than (2). The series (2) may be recovered from (3) through a Laplace transform

$$\widehat{S}_N(\lambda )=\frac{1}{\lambda }_0^{\mathrm{}}𝑑ze^{z/\lambda }B(z).$$
(4)
In order to proceed non-trivially, one first performs an approximate summation of the series (3) so that Eq.(4) then gives the resummed version of (2). Now, suppose that the only singularity of $`B(z)`$ is at $`z=1/p`$, with $`p`$ real and positive. Then the radius of convergence of the Borel series is $`1/p`$. In order to perform the integral in (4), one needs to extend the domain of convergence of the Borel series. One way to do this is by the method of conformal mapping. Define

$$w(z)=\frac{\sqrt{1+pz}1}{\sqrt{1+pz}+1}$$
(5)
which maps the Borel plane to a unit circle. The inverse of (5) is given by

$$z=\frac{4w}{p}\frac{1}{(1w)^2}$$
(6)
The idea is to rewrite (3) in terms of the variable $`w`$. Therefore, using (6), $`z^n`$ is expanded to order $`N`$ in $`w`$ and substituted into (3,4). The result is

$$S_N(\lambda )=1+\frac{1}{\lambda }\underset{n=1}{\overset{N}{}}\frac{f_n}{n!}\left(\frac{4}{p}\right)^n\underset{k=0}{\overset{Nn}{}}\frac{(2n+k1)!}{k!(2n1)!}_0^{\mathrm{}}e^{z/\lambda }w(z)^{(k+n)}𝑑z,$$
(7)
where $`w(z)`$ is given by (5). Equation (7) represents a resummation of the original series (2). This technique has been used in determining critical exponents in statistical systems where the singularity at $`z=1/p`$ is due to instantons , and in QCD at zero temperature where the singularity is due to renormalons .

Currently no information is available about the exact location of singularities in the Borel plane of thermal QCD though undoubtedly there is at least one on the negative semi-axis. Therefore in order to apply the resummation (7) to Eq.(1), a new idea is introduced in this paper: It is first assumed that the only singularity is at $`z=1/p,p>0`$, with the value of $`p`$ determined by the condition that it be the position of an extremum of (7). That is, $`p`$ is chosen to be a solution of

$$\left(\frac{S_N(\lambda ,p)}{p}\right)_{\lambda =\lambda _0}=0.$$
(8)
Since $`S_N(\lambda ,p)`$ depends on the coupling $`\lambda `$, one first fixes $`\lambda `$ at some reference value $`\lambda =\lambda _0`$ (say, at the mid-point of the range of interest) in order to solve Eq.(8). Fortunately, it turns out that the solution of (8), and hence the convergence of $`S_N`$, is not very sensitive to the exact value of $`\lambda _0`$.

Let me illustrate the technique by applying it to two cases where exact results are known. Consider first the integral

$$I(\lambda )=_0^{\mathrm{}}𝑑z\frac{e^z}{1+z\lambda }.$$
(9)
If the right-hand-side of (9) is expanded as a power series in $`\lambda `$ one obtains at $`N`$-th order

$$I_N(\lambda )=1+\underset{n=1}{\overset{N}{}}\lambda ^n(1)^nn!.$$
(10)
Clearly, from (9), the exact Borel transform of this series is $`B(z)=1/(1+z)`$ with a singularity at $`z=1`$. Ignoring this information, let us resum the divergent series (10) using (7) with $`f_n=(1)^nn!`$ and values of $`p`$ determined for each $`N`$ through equation (8) at the reference value $`\lambda _0=0.5`$. The results for (8) are as follows: There is no extremum for $`N=1`$. For $`N=2`$, there is a minimum at $`p=2.65`$. For $`N=3`$ there is a local maximum at $`p=1.5`$ and a minimum at $`p=5.1`$. For $`N=4`$ there is a local minimum at $`p=1.3`$, a local maximum at $`p=2.3`$ and a global minimum at $`p=8.4`$. Jumping ahead to $`N=8`$, there is a local minimum at $`p=1.075`$, a local maximum at $`p=1.5`$ and a global minimum at $`p=3.6`$.

Thus in general (8) has more than one solution for a given $`N`$ and $`\lambda _0`$. In Fig.(1), Eq.(7) is plotted for $`N=2,3,4`$ and $`8`$, at the respective minimum. Notice how the curves converge rapidly to the exact value given by (9). Alternatively, one might choose for each $`N`$, the value of $`p`$ (from the multiple solutions of Eq.(8)) which seems to be part of a converging sequence. In this case the values, $`p=2.65(N=2),p=1.5(N=3),p=1.3(N=4)`$, and $`p=1.075(N=8)`$, appear to converge to the exact value $`p=1`$. The curves are shown in Fig.(2). Clearly the curves in Fig.(2) converge faster to the exact value than those of Fig.(1), but unless one is interested in very high numerical accuracy, the difference is not significant. For example, at $`\lambda =0.5`$, the exact value of (9) is $`0.722657`$, while the resummed value for $`N=8`$ is given at the global minimum $`p=1.075`$ by $`S_8(\lambda =0.5,p=1.075)=0.722652`$, and at the local minimum $`p=3.6`$ by $`S_8(\lambda =0.5,p=3.6)=0.722524`$.

The main points illustrated by this example, which seem to be common to the other cases studied, are, (i) the rapid convergence of the resummed series represented by (7) compared to the original wildly oscillating series (10), (ii) the relative insensitivity of the convergence and numerical value of the resummed series to the particular extremum chosen among the possible multiple solutions of (8) (for a given $`N`$ and $`\lambda _0`$), even if the chosen value of $`p`$ is quite different from the exact value, (iii) the relative insensitivity of Eq.(8), and hence the convergence of (7) to the precise value of $`\lambda _0`$.

For another example, consider thermal $`O(M)`$ $`\lambda ^2\varphi _4^4`$ field theory in the limit $`M\mathrm{}`$. The exact free-energy density in this case has been determined in Ref.. In Eq.(5.8) of that paper111Note that the definition of the coupling constant used here is different from that in Ref.. the perturbative expansion, in $`\lambda `$, of the free-energy density is also given up to $`\lambda ^6`$. Defining $`S=(F(T)F(0))/F_{ideal}`$ and choosing $`\overline{\mu }=T`$ for simplicity, the values of $`f_n`$ for $`2n6`$ can be read off from Eq.(5.8) of Ref. and (8) solved at some reference value, say $`\lambda _0=4`$. The solutions of (8) in this case are: $`N=3,p=0.1`$(min); $`N=4,p=0.05`$(local max), $`p=0.2`$(min); $`N=5,p=0.025`$ (local min), $`p=0.1`$ (local max), $`p=0.3`$(global min); $`N=6,p=0.1`$ (local max), $`p=0.45`$ (global min). As in the first example, the convergence of the resummed series is found to be rapid even for large coupling, in contrast to the oscillatory behaviour of the ordinary pertubation expansion observed in . Fig.(3) shows the curves for $`S_N(\lambda ),3N6`$, for the value $`p=0.1`$, which seems to be the value $`p`$ converges to as $`N`$ increases. The exact value of the free-energy density at $`\lambda =8`$ taken from Fig.(6) of Ref. is about $`0.875`$. By comparison the resummed value predicted here is given at sixth order by $`S_6(\lambda =8,p=0.1)=0.889`$. On the other hand, if one evaluates $`S_6`$ not at $`p=0.1`$ but rather at its global minimum $`p=0.45`$, one gets $`S_6(\lambda =8,p=0.45)=0.849`$, a difference of less than $`5\%`$.

Actually, no information is available about the singularities in the Borel plane for the $`O(M\mathrm{})`$ scalar field studied in . The good agreement of the results obtained here with the exact results of Ref. leads one to conjecture that for the free-energy of this model, the singularity closest to the origin in the Borel plane might be near $`p=0.1`$, that is, $`z=10`$.

Note that although the coupling constant in the scalar field theory model is $`\lambda ^2`$, the perturbative expansion of the free-energy density contains the odd powers $`\lambda ^3`$ and $`\lambda ^5`$ which is typical of thermal field theories with massless particles (or at very high temperatures) . Physically these are due to collective effects such as Debye screening and it is sometimes suggested in the literature that these terms should be treated on a different footing. However, as the analysis above shows, from a mathematical point of view these odd powers are no different from the other terms in the expansion of the free-energy density and can be resummed as part of a single series.

Finally, the resummation technique of Eqs.(7,8) is applied to the free-energy density of $`SU(3)`$ gauge theory given in Eq.(1). As in Refs., I replace $`(\frac{\alpha }{\pi })^{\frac{1}{2}}`$ by the approximate two-loop running coupling constant defined by

$$\lambda (c,x)=\frac{2}{\sqrt{11L(c,x)}}\left(1\frac{51}{121}\frac{\mathrm{ln}(L(c,x))}{L(c,x)}\right)$$
(11)
where $`L(c,x)=\mathrm{ln}((2.28\pi cx)^2)`$, $`c=\overline{\mu }/2\pi T`$ and $`x=T/T_c`$, with $`T_c270MeV`$ the critical temperature which separates the low and high temperature phases . Furthermore, as in Ref., I have absorbed the $`\mathrm{ln}(\alpha )`$ term which appears at three-loop order into the coefficient of the $`\alpha ^2`$ term in (1).

Fixing first the reference values $`c_0=1,x_0=3`$ (which fixes the reference value of $`\lambda _0`$), the results of (8) are: $`N=2`$, no extremum; $`N=3,p=3.2`$(min); $`N=4,p=7.6`$(min); $`N=5,p=13.1`$(min). Since $`c`$ and $`x`$ appear in (1) and (11) only logarithmically, changing these values in the range of say, $`0.5<c<2`$, $`2<x<5`$, has almost no impact on the solution of Eq.(8) and hence on the optimal values of $`p`$.

The curves for the resummed series are plotted in Fig.(4) for $`c=1`$, that is at the renormalization scale $`\overline{\mu }=2\pi T`$. Again the rapid and monotonic convergence is manifest, the result approaching the ideal gas value even at moderate temperatures $`2T_c`$. One can estimate the effect of the unknown higher order, $`\lambda ^6`$, contribution. It turns out to be negligible . Therefore one feels confident that the $`N=5`$ curve in Fig.(4) is numerically close to the (unknown) exact sum of the perturbation series. In Fig.(5) the curve for $`S_5(x,c,p=13.1)`$ is plotted for three values of $`c`$ to indicate its mild dependence (less than $`1\%`$) to the arbitrary renormalization scale $`\overline{\mu }`$.

Lattice results for the free-energy density of pure $`SU(3)`$ theory are shown in Fig(6). The lattice community has indicated that their errors are under control (less than $`5\%`$). In that case, I am left with the task of explaining the significant difference (e.g. $`15\%`$ at $`T=3T_c`$) between the best analytically resummed result represented by the $`N=5`$ curve in Fig.(4) and the lattice data. At zero temperature, it is known that non-abelian gauge theories are not Borel summable . That is, $`B(z)`$ contains singularities for positive $`z`$, rendering the integral in (4) ambiguous. The situation is not expected to be different at non-zero temperature. Usually , the presence of such singularities is taken to indicate the existence of non-perturbative corrections. One can estimate the ambiguity, $`\delta S`$, and hence the non-perturbative correction, as the residue of the integrand in (4) at the location of the singularity . If the singularity of $`B(z)`$, on the positive semi-axis, closest to the origin is a pole at $`z=q`$, then from (4)

$$\delta S=\frac{A}{\lambda }e^{q/\lambda }$$
(12)
where $`A`$ is a constant. Assuming that the difference between the lattice data and Fig.(5), is due to (12), the constants $`A`$ and $`q`$ can be determined by rewriting (12) as

$$\mathrm{ln}(\lambda \delta S)=\mathrm{ln}(A)q/\lambda $$
(13)
and using for $`\delta S`$ the difference between Fig.(6) and the median curve in Fig.(5) (i.e. $`c=1`$). Fig.(7) shows the left-hand-side of (13) plotted against $`1/\lambda `$. This gives $`A=e^{8.7}`$ and $`q=2.62`$, and so, with $`\lambda (x)\lambda (c=1,x)`$,

$$S_{latt}=S_{pert}\frac{1}{\lambda (x)}e^{8.72.62/\lambda (x)},$$
(14)
where $`S_{latt}`$ represents the lattice data for the free-energy, and $`S_{pert}`$ the Borel resummed perturbative result, both normalized with respect to the ideal gas value.

It is extremely reassuring that both the sign and magnitude of $`q`$ determined in this way are self-consistent with the assumptions made. In particular, the singularity at $`z=q=2.62`$ is more than 30 times away from the origin than the singularity at $`z=1/p=1/13.1`$ and justifies a posteriori the resummation procedure (7) which considered only the nearest singularity. Furthermore since $`S_{pert}`$ is extremely close to the ideal gas value, Eq.(14) may be interpreted as a generalisation of phenomenological equations of state  for the free-energy where the second term on the right-hand-side of (14) is called a ’bag constant’. In our case the ‘constant’ is really temperature dependent and represents a non-perturbative contribution to the free-energy that vanishes at infinite temperature. Note that this interpretation of the second term on the right-hand-side of (14) is consistent with the usual one only because the resummed perturbative result lies above the lattice data.

Let me now summarize the main results of this paper. A new procedure, a variational version of the well-known conformal-mapping of Borel series, was introduced and its practicability illustrated. It was shown that the badly divergent series for the free-energy density of thermal $`SU(3)`$ gauge theory (1) could be resummed in a systematic and relatively well-motivated way into a rapidly convergent series. However the final result differed significantly from non-perturbative lattice data, suggesting that the discrepancy is due to the non-Borel-summability of the theory. Numerically, the difference was parametrised in terms of two constants (see (14)), and it is suggested that the ambiguity in the Borel integral is the bag ‘constant’ used in phenomenological models for the free-energy density.

In physical terms, the situation may be described as follows. If the Borel resummed perturbative expansion had agreed with the lattice data, then one would have argued that the high-temperature phase of $`SU(3)`$ theory is appropriately described by weakly coupled gluons. However, if the deviation of lattice data from the resummed and convergent perturbation expansion found here is taken at face value, then one is led to conclude that even at temperatures as high as three times $`T_c`$ ($`700MeV`$), the phase of thermal $`SU(3)`$ is not accurately described by weakly coupled gluons but rather by more complicated structures which are responsible for the non-perturbative bag contribution. If one accepts this conclusion, then one must also be open to the possibility that what will be produced at CERN and Brookhaven might be better characterised as something other than a quark-gluon plasma. For some alternative descriptions of the high-temperature phase of QCD see, for example, .

An extension of the analysis presented here to QCD and other thermal gauge theories, a further development of the methodology itself and its applications to other physical problems, will be presented in an accompanying series of papers .

Acknowledgements: I thank Claudio Coriano for stimulating discussions and hospitality at Martignano, and the Department of Physics at Lecce for financial support. I also thank the YITP at Stony Brook for hospitality during the final stages of this work.

Figure Captions

Figure 1: Plot of Eq.(7) for the model in Eq.(10), for $`N=2,3,4`$ and $`8`$ at the respective global minimum of Eq.(8). The curves move upwards with increasing $`N`$. The curves for $`N=3,4`$ practically coincide, while the curve for the exact expression (9) is indistinguisable from that for $`N=8`$.

Figure 2: Plot of Eq.(7) for the model in Eq.(10), for $`N=2,3,4`$ and $`8`$ at the respective values of $`p`$ from Eq.(8) which converge to $`1`$. The values of the curves at $`\lambda =0.5`$ are, $`N=2(0.704),N=3(0.726),N=4(0.7219),N=8(0.7227)`$. The curve for the exact expression (9) is hardly distinguisable from that for $`N=3,4,8`$.

Figure 3: Plot of Eq.(7) for the model in Ref. for $`N=3,4,5`$ and $`N=6`$ at the value $`p=0.1`$. The curves move upwards as $`N`$ increases. For $`N=5,6`$ they are practically identical. The exact curve in Ref. lies slightly below that for $`N=6`$: At $`\lambda =8`$ the exact curve has the value $`0.875`$.

Figure 4: The resummed (7) free-energy density of hot $`SU(3)`$ gauge theory for $`N=3,4`$ and $`5`$, at the renormalization point $`\overline{\mu }=2\pi T`$. The curves move upwards as $`N`$ increases.

Figure 5: The fifth order resummed free-energy density (7) of hot $`SU(3)`$ gauge theory at $`p=13.1`$, for the renormalization scale values $`\overline{\mu }=0.5,1`$ and $`2`$. The free-energy density increases with increasing $`\overline{\mu }`$.

Figure 6: Mean lattice results for the free-energy density of hot $`SU(3)`$ gauge theory from Ref.. Here $`S_{latt}`$ refers to the free-energy divided by the free-energy of an ideal gas of gluons.

Figure 7: Plot of the left-hand-side of Eq.(13) against $`1/\lambda `$ at several points (temperatures).