text:
# Dynamical turbulent flow on the Galton board with friction and magnetic field

\[

## Abstract

We study numerically and analytically the dynamics of charged particles on the Galton board, a regular lattice of disc scatters, in the presence of constant external force, magnetic field and friction. It is shown that under certain conditions friction leads to the appearance of a strange chaotic attractor. In this regime the average velocity and direction of particle flow can be effectively affected by electric and magnetic fields. We discuss the applications of these results to the charge transport in antidot superlattices and stream of suspended particles in a viscous flow through scatters.

\]

It is well known that dissipation can lead to the appearance of strange chaotic attractors in nonlinear nonautonomous dynamical systems . In this case the energy dissipation is compensated by an external energy flow so that stationary chaotic oscillations set in on the attractor. Such an energy flow is absent in the Hamiltonian conservative systems and therefore the introduction of dissipation or friction is expected to drive the system to simple fixed points in the phase space. This rather general expectation is surly true if the system phase space is bounded. However a much richer situation appears in the case of unbounded space, where unexpectedly a strange attractor can be induced by dissipation in an originally conservative system. To investigate this situation we study the dynamics of particles on the Galton board in the presence of constant external fields and friction. This board, introduced by Galton in 1889 , represents a triangular lattice of rigid discs with which particles collide elastically. For the case of free particle motion, the collisions with discs make the dynamics completely chaotic on the energy surface as it was shown by Sinai . In this paper we study how the dynamics of a charged particle in a presence of electric and magnetic fields is affected by a friction force $`𝐅_f=\gamma 𝐯`$ directed against particle velocity $`𝐯`$ (see Fig. 1). Without discs an external in-plane electric field $`𝐄`$ and a perpendicular magnetic field $`𝐁`$ create a stationary particle flow with the velocity $`𝐯_f=(𝐟+F_L)/\gamma `$. Here $`𝐟=e𝐄`$ is the effective force, $`𝐅_L=e𝐯_𝐟\times 𝐁`$ is the Lorentz force and $`e,m`$ are the particle charge and mass. All perturbations decay to this flow with a rate proportional to $`\gamma `$ so that this laminar flow can be considered as a simple attractor. The effects of friction inside one cell of the Galton board at $`𝐁=0`$ have been studied in  and it has been found that friction leads to the appearance of a nontrivial strange attractor. At present the effects of energy dissipation are actively investigated with the aim to construct equilibrium and nonequilibrium steady states in a deterministic way (see  and Refs. therein). Here the Nosé-Hoover thermostat with a momentum-dependent friction coefficient allows to reach a number of interesting results with applications to molecular dynamics and nonequilibrium liquids . In our studies, contrary to , we concentrate mainly on the spatial structure of the turbulent chaotic particle flow appearing in the presence of friction. We show that the flow direction can be efficiently affected by a magnetic field. The obtained results describe the electron dynamics in antidot superlattice which has been experimentally realized in semiconductor heterostructures . In such structures the effects of classical chaos play an important role  and the effects of friction we discuss here can appear for relatively strong electric fields.

To study numerically the dynamics of this model, we fix the disc radius $`a=1`$, and $`e=m=1`$ so that the system is characterized only by the distance $`R`$ between discs packed into equilateral triangles all over the $`(x,y)`$ plane, the friction coefficient $`\gamma `$, the external force of strength $`f=|𝐟|`$ and the cyclotron frequency $`\omega _c=B`$ . For $`RR_c=4/\sqrt{3}`$ there is no straight line, by which a particle can cross the whole plane without collisions (at $`f=0,B=0,\gamma =0`$), and we start the discussion from this case. The particle dynamics is simulated numerically by using the exact solution of Newton equations between collisions, and by determining the collision points with the Newton algorithm. This way the trajectories are computed with high precision, and for example for $`\gamma =0`$ the total energy is conserved with a relative precision of $`10^{14}`$ . A typical example of the chaotic flow formed by an ensemble of particle trajectories is shown in Fig. 1. To illustrate the mixing properties of the flow we attributed a color to each trajectory that allows to follow their interpenetration and spreading. This figure shows that there is a certain penetration depth of one color into another, however this depth is finite since on average the initial color repartition is still visible.

The properties of color penetration can be understood from the analysis of single trajectory dynamics. Such a typical example is presented in Fig. 2. It shows that the particle moves with an average constant velocity $`v_f`$ under some angle $`\alpha `$ to the external force $`𝐟`$, except for $`B=0`$ where $`\alpha =0`$. This drift velocity is constant only on average since on a smaller scale the particle moves chaotically between scatters following a strange chaotic attractor. In Fig. 2 for $`\gamma =0.1`$ the drift velocity is relatively large and the particle does not have enough time to move around many scatters in the direction perpendicular to the flow. As a result the penetration depth for color mixing is not very large. Surprisingly, for smaller friction, the drift velocity becomes smaller and the penetration depth increases so that the particle makes many turns around discs as it is shown in Fig. 2 for $`\gamma =0.004`$. This dependence is opposite to the case without scatters where $`v_f`$ drops with the $`\gamma `$ growth.

This result can be understood on the basis of the following physical arguments (for $`B=0`$ see also ). In the regime with weak friction the particles start to diffuse among discs in a chaotic manner with the diffusion rate $`D=vl/2`$ where $`v`$ is the particle velocity and $`l`$ is the mean free path . For $`R1`$ we have $`lR1`$ while the dependence of $`l`$ on $`R`$ will be discussed in more detail latter. During the dissipative time scale $`\tau _\gamma =m/\gamma `$ this diffusion leads to the particle displacement $`\mathrm{\Delta }r\sqrt{Dm/\gamma }`$ along the direction of the drift velocity $`𝐯_f`$. This gives the change in the potential energy $`Uf\mathrm{\Delta }r\mathrm{cos}\alpha f_{ef}\sqrt{Dm/\gamma }`$ where $`f_{ef}=f\mathrm{cos}\alpha `$. In the stationary regime at time $`t\tau _\gamma `$, this potential energy should be comparable with the kinetic energy of the particle so that $`Umv^2`$, where $`v^2=<v_x^2+v_y^2>`$ is the average velocity square. Hence

$$v^2(f\mathrm{cos}\alpha )^{4/3}(l/\gamma m)^{2/3}.$$
(1)
This relation allows to determine the drift velocity of the flow $`𝐯_f`$. Indeed during the time $`\tau _c`$ between collisions the particle is accelerated by average forces $`𝐟`$ and $`𝐅_L=e𝐯_f\times B`$ that gives the average drift velocity $`𝐯_f=(𝐟+F_L)\tau _c/m`$. Since the dynamics is chaotic the direction of velocity is changed randomly after each collision so that $`𝐯_f`$ is accumulated only between collisions. The time $`\tau _c`$ is determined by the mean free path $`l`$ and the average velocity $`v`$: $`\tau _c=l/v=2D/v^2`$. Thus for the angle $`\alpha `$ between $`𝐯_f`$ and $`𝐟`$ we obtain the relation:

$$\mathrm{tan}\alpha =eB\tau _c/m\frac{eBl^{2/3}\gamma ^{1/3}}{(mf\mathrm{cos}\alpha )^{2/3}}.$$
(2)
The amplitude of fluctuations around this direction is $`\mathrm{\Delta }r\sqrt{Dm/\gamma }`$ which also determines the color mixing depth (see Fig. 1).

From (1) and (2) we obtain the drift velocity amplitude

$$v_f=f_{ef}\tau _c/ml^{2/3}(\gamma f_{ef}/m^2)^{1/3},$$
(3)
with $`f_{ef}=f\mathrm{cos}\alpha `$. For $`B=0`$ the particles flow in the $`𝐟`$ direction. In this case their mobility is $`\mu =v_f/f=\tau _c/m=D/(mv^2/2)`$. This is in fact the Einstein relation according to which the mobility is given by ratio of the diffusion rate to the average kinetic energy (temperature) . At $`B=0`$ the relations (1), (2) and (3) are in the agreement with those in  and with the numerical data shown in Fig. 3. The values of $`v_f`$ and $`v^2`$ are obtained from one very long trajectory (with a length of $`10^310^4R`$ ) or 10 shorter trajectories. Within statistical fluctuations this gives the same $`v_f`$ and $`v^2`$ independent of their initial values.

The relations (1), (2), (3) allow to estimate the value of the Lyapunov exponent $`\lambda `$. Indeed, the particle moves with a typical velocity $`v`$ and as in the case of the Sinai billiard with $`lR1`$ we have $`\lambda v/l`$. Therefore in the regime when:

$$\gamma <\gamma _c\sqrt{mf_{ef}/l},$$
(4)
the value of $`\lambda `$ is much larger than the dissipation rate $`\gamma /m`$. As a result for $`\gamma \gamma _c`$ the strange attractor is fat and its fractal dimension is close to the maximal dimension $`4`$, which is determined by the number of degrees of freedom (we remind that contrary to the nondissipative case the energy is not conserved). For $`\gamma \gamma _c`$ the dissipation time $`\tau _\gamma `$ becomes much shorter than the time between collisions $`\tau _c`$. In this case the dissipation dominates chaos and the strange attractor degenerates into a simple attractor. For $`\gamma <\gamma _c`$ our numerical simulations performed with high computer accuracy show that trajectories remain chaotic for displacements from the origin being larger than $`10^5R`$. Also at $`B=0`$ it can be shown that $`\gamma _cmv/l\sqrt{fma}/R`$ for $`Ra`$ and $`\gamma _c\sqrt{mf/a^3}\mathrm{\Delta }R`$ for $`\mathrm{\Delta }R=R2aa`$.

The above changes in the mean free path $`l`$ at $`B=0`$ also affect the drift velocity of the flow through the relation (3). Indeed, for $`\mathrm{\Delta }R1`$ we expect $`l\mathrm{\Delta }R`$ that gives $`v_f\mathrm{\Delta }R^{2/3}`$. This dependence is close to the numerical data shown in Fig. 4 even if the numerical value of the exponent is approximated better by $`0.5`$. In the other limit $`R1`$, we have $`lR^2/a`$ that gives $`v_fR^{4/3}`$. This power dependence is in a satisfactory agreement with the data in Fig. 4 although the numerical value of the exponent is closer to $`1`$. We attribute these small deviations in the exponent values to a restricted interval of variation in $`R`$. Actually we can not use very large/small values of $`\mathrm{\Delta }R`$ since in these limits the value of $`\gamma `$ becomes comparable with $`\gamma _c`$ and chaotic attractor disappears. We note that according to our data (see Fig. 4) the strange attractor exists even in the case $`R>R_c=4/\sqrt{3}`$ when at $`f=0,\gamma =0`$ there are straight trajectories crossing the whole plane without collision. Apparently the contribution of these orbits is not significant if $`\gamma >0`$ and if $`𝐟`$ is not directed along these lines.

The introduction of magnetic field allows to change efficiently the direction of the flow. The numerical data for the variation of $`\mathrm{tan}\alpha `$ with the magnetic field and other system parameters are presented in Fig. 5. The average dependence is in a good agreement with the equation (2) for a large region of parameter variation where $`\mathrm{tan}\alpha `$ changes by two orders of magnitude. At the same time for moderate angles $`\alpha <1`$ the flow velocity $`v_f`$ is weakly affected by $`B`$. For example, for the case of Fig. 2, $`v_f`$ remains practically the same for $`B=2`$ and $`B=0`$ ($`v_f0.15`$). For $`\alpha >1`$ magnetic field starts to change significantly $`v^2`$ and $`v_f`$ in agreement with (1) and (3). While the data in Fig. 5 in average follow the dependence (2) the fluctuations around the average are rather large. Their origin becomes clear from Fig. 6 where only the magnetic field is changed. In this case, $`\mathrm{tan}\alpha `$ has a pronounced peak which is located near the value of $`B`$ where the cyclotron radius $`r_c=v/B`$ is equal to the disc radius. Indeed for $`r_c>a`$ a trajectory can make a full turn around a disc that allows to increase $`\alpha `$ and to reach a strong deviation of the global flow from the direction of the electric field. The growth of $`\alpha `$ leads to a drop of currant (conductivity) in the $`𝐟`$ direction and hence to the increase of resistivity. In fact, the peaks in resistivity near $`r_ca`$ had been observed experimentally  and explained theoretically  in the linear response regime. Our data show that the peaks should also exist in the regime with strange chaotic attractor and relatively strong electric field where the $`IV`$ characteristics become nonlinear. For a denser package of discs the above peak is still present ($`R=2.24`$ in Fig. 6) but it is less pronounced. It is interesting to note that in this case, $`\alpha `$ can even be negative so that the particles flow against the average Lorentz force. We note that the possibility of such a flow had been discussed for systems with Hamiltonian dynamical chaos in . Generally for a dense disc package the contribution of resonant orbits with $`r_ca`$ starts to be significant and deviations from the average dependence (2) become rather strong.

The above dynamical turbulent flow has rather interesting and unusual properties and it would be interesting to study it in experiments with antidot superlattices like in . The regime we discussed should appear when the steady state velocity $`v`$ from (1) becomes larger than the Fermi velocity $`v_F`$ of the charge carriers. Thus, in addition to (4), the condition $`v>v_F`$ determines the threshold for appearance of a strange attractor. The experimental investigation of this phenomenon should also shade light on the role of quantum effects in such regime. We note that in the derivation of (1) - (4) we didn’t use any specific properties of the disc distribution in the plane and therefore the results remain valid for randomly distributed discs. Our study represents also certain interest for transport properties of neutral/charged particles suspended in a viscous flow streaming through a system of scatters. Indeed, a laminar stream with the velocity $`v_s`$ creates an effective force $`f=v_s\gamma _{ef}`$ where $`\gamma _{ef}`$ is the effective friction created by the viscosity of the liquid. Such kind of transport can be studied experimentally with viscous liquids and its investigation can contribute to a better understanding of the interplay between dissipation, turbulence and chaos.

We thank R. Klages for constructive critical remarks and for pointing to us Refs. , and W. Hoover for his stimulating interest to our work.


format: text

text:
# Radiative Equilibrium and Temperature Correction in Monte Carlo Radiation Transfer

## 1 Introduction

There is an ever increasing wealth of observational evidence indicating the non-sphericity of almost every type of astronomical object (e.g., extended circumstellar environments, novae shells, planetary nebulae, galaxies, and AGNs). To accurately interpret this data, detailed two- and three-dimensional radiation transfer techniques are required. With the availability of fast workstations, many researchers are turning to Monte Carlo techniques to produce model images and spectra for the asymmetric objects they are investigating. In Monte Carlo radiation transfer simulations, packets of energy or “photons” are followed as they are scattered and absorbed within a prescribed medium. One of the features of this technique is that the locations of the packets are known when they are absorbed, so we can determine where their energy is deposited. This energy heats the medium, and to conserve radiative equilibrium, the absorbed energy must be reradiated at other wavelengths, depending on the opacity sources present. Tracking these photon packets, while enforcing radiative equilibrium, permits the calculation of both the temperature structure and emergent spectral energy distribution (SED) of the envelope. The ability of Monte Carlo techniques to easily follow the transfer of radiation through complex geometries makes them very attractive methods for determining the temperature structure within non-spherical environments — a task which is very difficult with traditional ray tracing techniques.

Previous work on this problem for spherical geometries includes the approximate solutions by Scoville & Kwan (1976), who ignored scattering, Leung (1976), and diffusion approximations by Yorke (1980). The spherically symmetric problem has been solved exactly by Rowan-Robinson (1980), Wolfire & Cassinelli (1986), and Ivezić & Elitzur (1997), who used a scaling technique. Extensions of the exact solution to two dimensions have been performed by Efstathiou & Rowan-Robinson (1990, 1991), while approximate two-dimensional models have been presented by Sonnhalter, Preibisch, & Yorke (1995) and Men’shchikov & Henning (1997).

Radiative equilibrium calculations using the Monte Carlo technique have been presented by Lefevre, Bergeat, & Daniel (1982); Lefevre, Daniel, & Bergeat (1983); Wolf, Henning, & Secklum (1999); and Lucy (1999). Most of these authors (Lucy being exceptional) use a technique in which stellar and envelope photon packets are emitted separately. The number of envelope packets to be emitted is determined by the envelope temperature, while the envelope temperature is determined by the number of absorbed packets. Consequently these techniques require iteration, usually using the absorbed stellar photons to provide an initial guess for the envelope temperature. The iteration proceeds until the envelope temperatures converge. Note that the stellar luminosity is not automatically conserved during the simulation; only after the temperatures converge is the luminosity approximately conserved.

In contrast, Lucy adopts a strategy in which absorbed photon packets are immediately re-emitted, using a frequency distribution set by the current envelope temperature. Although the frequency distribution of the reprocessed photons is incorrect (until the temperatures have converged), his method automatically enforces local radiative equilibrium and implicitly conserves the stellar luminosity. The insight of Lucy’s method is that conservation of the stellar luminosity is more important than the spectral energy distribution when calculating the radiative equilibrium temperatures. Nonetheless, this method requires iteration.

The primary problem faced by Lucy’s method is the incorrect frequency distribution of the re-emitted photons. In this paper we develop an adaptive Monte Carlo technique that corrects the frequency distribution of the re-emitted photons. Essentially, our method relaxes to the correct frequency and temperature distribution. Furthermore it requires no iteration as long as the opacity is independent of temperature. Such is the case for astrophysical dust. In Section 2, we describe the temperature correction algorithm. We compare the results of our code with a spherically symmetric code in Section 3, and in Section 4 we present results for two dimensional axisymmetric density structures.

## 2 Monte Carlo Radiative Equilibrium

We wish to develop a method to calculate the temperature distribution throughout an extended dusty environment for use with Monte Carlo simulations of the radiation transfer. The radiation transfer technique we employ has been described in detail in other papers: Code & Whitney (1995); Whitney & Hartmann (1992, 1993); Wood et al. (1996), so we only summarize it here. The basic idea is to divide the luminosity of the radiation source into equal-energy, monochromatic “photon packets” that are emitted stochastically by the source. These packets are followed to random interaction locations, determined by the optical depth, where they are either scattered or absorbed with a probability given by the albedo. If the packet is scattered, a random scattering angle is obtained from the scattering phase function (differential cross section). If instead the packet is absorbed, its energy is added to the envelope, raising the local temperature. To conserve energy and enforce radiative equilibrium, the packet is re-emitted immediately at a new frequency determined by the envelope temperature. These re-emitted photons comprise the diffuse radiation field. After either scattering or absorption plus reemission, the photon packet continues to a new interaction location. This process is repeated until all the packets escape the dusty environment, whereupon they are placed into frequency and direction-of-observation bins that provide the emergent spectral energy distribution. Since all the injected packets eventually escape (either by scattering or absorption followed by reemission), this method implicitly conserves total energy. Furthermore it automatically includes the diffuse radiation field when calculating both the temperature structure and the emergent spectral energy distribution.

We now describe in detail how we calculate the temperature structure and SEDs of dusty environments illuminated by a radiation source. This radiation can come from any astrophysical source, either internal or external, point-like or extended.

### 2.1 Radiative Equilibrium Temperature

Initially we divide the source luminosity, $`L`$, into $`N_\gamma `$ photon packets emitted over a time interval $`\mathrm{\Delta }t`$. Each photon packet has the same energy, $`E_\gamma `$, so

$$E_\gamma =L\mathrm{\Delta }t/N_\gamma .$$
(1)
Note that the number of physical photons in each packet is frequency-dependent.

When the monochromatic photon packet is injected into the envelope, it will be assigned a random frequency chosen from the spectral energy distribution of the source. This frequency determines the dust absorptive opacity, $`\kappa _\nu `$, and scattering opacity, $`\sigma _\nu `$ (both per mass), as well as the scattering parameters for the ensuing random walk of the packet through the envelope. The envelope is divided into spatial grid cells with volume $`V_i`$, where $`i`$ is the cell index. As we inject source photon packets, we maintain a running total, $`N_i`$, of how many packets are absorbed in each grid cell. Whenever a packet is absorbed in a grid cell, we deposit its energy in the cell and recalculate the cell’s temperature. The total energy absorbed in the cell

$$E_i^{\mathrm{abs}}=N_iE_\gamma .$$
(2)
We assume that the dust particles are in local thermodynamic equilibrium (LTE), and for simplicity we adopt a single temperature for the dust grains. Note that although we use dust for the continuous opacity source, we could replace the dust by any continuous LTE opacity source that is independent of temperature. In radiative equilibrium, the absorbed energy, $`E_i^{\mathrm{abs}}`$, must be reradiated. The thermal emissivity of the dust $`j_\nu =\kappa _\nu \rho B_\nu (T)`$, where $`B_\nu (T)`$ is the Planck function at temperature $`T`$, so the emitted energy is

$`E_i^{\mathrm{em}}`$ $`=`$ $`4\pi \mathrm{\Delta }t{\displaystyle 𝑑V_i\rho \kappa _\nu B_\nu (T)𝑑\nu }`$ (3)
$`=`$ $`4\pi \mathrm{\Delta }t{\displaystyle \kappa _P(T)B(T)\rho 𝑑V_i},`$ (4)
where $`\kappa _P=\kappa _\nu B_\nu 𝑑\nu /B`$ is the Planck mean opacity, and $`B=\sigma T^4/\pi `$ is the frequency integrated Planck function. If we adopt a temperature that is constant throughout the grid cell, $`T_i`$, then

$$E_i^{\mathrm{em}}=4\pi \mathrm{\Delta }t\kappa _P(T_i)B(T_i)m_i,$$
(5)
where $`m_i`$ is the mass of the cell.

Equating the absorbed (2) and emitted (5) energies, we find that after absorbing $`N_i`$ packets, the dust temperature is given by

$$\sigma T_i^4=\frac{N_iL}{4N_\gamma \kappa _P(T_i)m_i}.$$
(6)
Note that the Planck mean opacity, $`\kappa _P`$, is a function of temperature, so equation (6) is actually an implicit equation for the temperature, which must be solved every time a packet is absorbed. Since this equation is solved so many times, an efficient algorithm is desirable. Fortunately $`\kappa _P`$ is a slowly varying function of temperature, which implies a simple iterative algorithm may be used to solve equation (6). To do so efficiently, we pre-tabulate the Planck mean opacities for a large range of temperatures and evaluate $`\kappa _P(T_i)`$ by interpolation, using the temperature from the previous guess. After a few steps, we have the solution for $`T_i`$. Note that because the dust opacity is temperature-independent, the product $`\kappa _P\sigma T_i^4`$, which is proportional to $`\kappa _\nu B_\nu 𝑑\nu `$, increases monotonically with temperature. Consequently $`T_i`$ always increases when the cell absorbs an additional packet.

### 2.2 Temperature Correction

Now that we know the temperature after absorbing an additional packet within the cell, we must reradiate this energy so that the heating always balances the cooling. Prior to absorbing this packet, the cell previously has emitted packets that carried away an energy corresponding to the cell’s previous emissivity $`j_\nu ^{}=\kappa _\nu B_\nu (T_i\mathrm{\Delta }T)`$, where $`\mathrm{\Delta }T`$ is the temperature increase arising from the packet absorption. Note that these previous packets were emitted with an incorrect frequency distribution corresponding to the previous temperature, $`T_i\mathrm{\Delta }T`$. The total energy that should be radiated now corresponds to $`j_\nu `$ at the new temperature, $`T_i`$. Thus the additional energy to be carried away is given by

$$\mathrm{\Delta }j_\nu =j_\nu j_\nu ^{}=\kappa _\nu \left[B_\nu (T_i)B_\nu (T_i\mathrm{\Delta }T)\right],$$
(7)
which is the shaded area shown in Figure 1. As long as the packet energy $`E_\gamma `$ is not too large (this may be assured by choosing a large enough number of photon packets, $`N_\gamma `$, to use for the simulation), the temperature change $`\mathrm{\Delta }T`$ is small, so the temperature correction spectrum

$$\mathrm{\Delta }j_\nu \kappa _\nu \mathrm{\Delta }T\frac{dB_\nu }{dT}.$$
(8)
Note that $`\mathrm{\Delta }j_\nu `$ is everywhere positive because $`\mathrm{\Delta }T>0`$, and the Planck function is a monotonically increasing function of temperature. Therefore to correct the previously emitted spectrum, we immediately re-emit the packet (to conserve energy), and we choose its frequency using the shape of $`\mathrm{\Delta }j_\nu `$. This procedure statistically reproduces $`\mathrm{\Delta }j_\nu `$ for the distribution of the re-emitted packets. Normalizing this distribution, we find the temperature correction probability distribution

$$\frac{dP_i}{d\nu }=\frac{\kappa _\nu }{K}\left(\frac{dB_\nu }{dT}\right)_{T=T_i},$$
(9)
where $`dP_i/d\nu `$ is the probability of re-emitting the packet between frequencies $`\nu `$ and $`\nu +d\nu `$, and the normalization constant $`K=_0^{\mathrm{}}\kappa _\nu (dB_\nu /dT)𝑑\nu `$.

Now that the packet’s frequency has changed, we change the opacity and scattering parameters accordingly and continue with scattering, absorption, temperature correction, and re-emission until all the photon packets finally escape from the system. In principle we could also account for back-warming of the source. Whenever a packet hits the source, the source must reradiate this new energy. This will change the temperature of the source, and the new source photons can be emitted using a difference spectrum similar to equation (9).

When we begin our calculation, no packets have been absorbed, so the initial temperature is zero throughout the envelope. This means that the initial temperature change is not small as is required by equation (8). One could use equation (7) to re-emit the first packet that is absorbed; however, this is not necessary. The number of packets producing this initial temperature change is small; it is of order the number of spatial grid cells. Furthermore, these packets generally are re-emitted at such long wavelengths that they are not observable. Consequently the error arising from using equation (8) to re-emit every packet is too small to be of importance.

As the simulation runs, the envelope starts at zero temperature. It then heats up, and the radiation field “relaxes” to its final spectral shape. The temperature correction procedure is simply a way of re-ordering the calculation (which frequencies are being used at a given moment) so that in the end, all the frequencies have been properly sampled. Consequently, after all the stellar photon packets have been transported, the final envelope temperature is the correct radiative equilibrium temperature, and the emergent spectral energy distribution has the correct frequency shape. Furthermore, the photon re-emission automatically accounts for the diffuse envelope emission. Note that energy is necessarily conserved, and there is no time-consuming iteration in our scheme; calculating the radiative equilibrium temperature requires no more computational time than an equivalent pure scattering calculation in which the temperature structure is held fixed. Similarly, there is no issue of convergence in our method. Unlike $`\mathrm{\Lambda }`$-iteration, the photon packets carry energy over large distances throughout the envelope because they are never destroyed, and of course there is no iteration at all. After running $`N_\gamma `$ packets, we have the final answer. The simulation does not continue running until some convergence criterion is met, and the only source of error in the calculation is the statistical sampling error inherent in Monte Carlo simulations.

## 3 Benchmark Verification

To validate our method for determining the radiative equilibrium temperatures and emergent fluxes, we compared our results against a set of benchmark calculations recently developed by Ivezić et al. (1997) for testing spherically symmetric dust radiative equilibrium codes. The parameters listed by Ivezić et al. enable us to exactly reproduce the same set of physical conditions (i.e, input spectrum, dust destruction radius, optical depth, opacity frequency distribution, and radial density structure).

For all benchmark tests, Ivezić et al. used a point source star, radiating with a black body spectrum whose temperature $`T_{}=2500`$K. The dust density distribution was a power law with radius,

$$\rho =\rho _0\left(\frac{R_{\mathrm{dust}}}{r}\right)^p,$$
(10)
where

$$\rho _0=\frac{\tau _\lambda }{(\kappa _\lambda +\sigma _\lambda )(1R_{\mathrm{dust}}/R_{\mathrm{max}})}\{\begin{array}{cc}1/R_{\mathrm{max}}\hfill & \text{(}p=0\text{),}\hfill \\ 1/R_{\mathrm{dust}}\hfill & \text{(}p=2\text{).}\hfill \end{array}$$
(11)
The inner radius of the envelope is the dust destruction radius, $`R_{\mathrm{dust}}`$, the outer radius is $`R_{\mathrm{max}}`$, and the total radial optical depth is $`\tau _\lambda `$, specified at $`\lambda =1\mu \mathrm{m}`$. The dust absorptive opacity, $`\kappa _\nu `$, and scattering opacity, $`\sigma _\nu `$, were taken to be

$`{\displaystyle \frac{\kappa _\nu }{(\kappa +\sigma )_{1\mu \mathrm{m}}}}`$ $`=`$ $`0.5\{\begin{array}{cc}1\hfill & \text{(}\lambda <1\mu \mathrm{m}\text{),}\hfill \\ (1\mu \mathrm{m}/\lambda )\hfill & \text{(}\lambda >1\mu \mathrm{m}\text{),}\hfill \end{array}`$ (12)
$`{\displaystyle \frac{\sigma _\nu }{(\kappa +\sigma )_{1\mu \mathrm{m}}}}`$ $`=`$ $`0.5\{\begin{array}{cc}1\hfill & \text{(}\lambda <1\mu \mathrm{m}\text{),}\hfill \\ (1\mu \mathrm{m}/\lambda )^4\hfill & \text{(}\lambda >1\mu \mathrm{m}\text{).}\hfill \end{array}`$ (13)
Since the total optical depth at $`1\mu \mathrm{m}`$ is independently specified, these opacities have been normalized to that at $`1\mu \mathrm{m}`$ for convenience. The wavelength-dependent scattering albedo is given by $`a=\sigma _\nu /(\kappa _\nu +\sigma _\nu )`$, and the scattering was assumed to be isotropic (note that in dust simulations, we would normally use a non-isotropic phase function for the scattering).

In principle, the inner radius of the dust shell, $`R_{\mathrm{dust}}`$, is determined by the dust condensation temperature, chosen by Ivezić et al. to be $`T_{\mathrm{cond}}=800\mathrm{K}`$. However, we are only testing the temperature correction procedure, so we have not implemented a scheme to solve self consistently for the dust destruction radius. The values for $`R_{\mathrm{dust}}`$ were calculated instead using Ivezić et al.’s eq. (4) and data from their Table 1. Finally, the outer radius of the dust shell was set to be $`R_{\mathrm{max}}=10^3R_{\mathrm{dust}}`$. The parameters describing the various test simulations are summarized in Table 1.

To begin the simulation, we release stellar photon packets with a black body frequency distribution, given by the normalized Planck function

$$b_\nu (x)=\frac{15}{\pi ^4}\frac{x^3}{\mathrm{e}^x1},$$
(14)
where $`x=h\nu /kT_{}`$. A particularly simple method for sampling the black body distribution is given by Carter & Cashwell (1975). Since this reference is somewhat obscure and difficult to obtain, we summarize the method here. First, generate a uniform random number, $`\xi _0`$, in the range 0 to 1, and determine the minimum value of $`l`$, $`l_{\mathrm{min}}`$, that satisfies the condition

$$\underset{i=1}{\overset{l}{}}i^4\frac{\pi ^4}{90}\xi _0.$$
(15)
Next obtain four additional uniform random numbers (in the range 0 to 1), $`\xi _1`$, $`\xi _2`$, $`\xi _3`$, and $`\xi _4`$. Finally, the packet frequency is given by

$$x=\mathrm{ln}(\xi _1\xi _2\xi _3\xi _4)/l_{\mathrm{min}}.$$
(16)
After emitting these packets from the star, we track them through the envelope. To determine the envelope temperature, we must count how many packets are absorbed in each grid cell. Since the envelope is spherically symmetric, we employ a set of spherical shells for our grid. To obtain the best Poisson error statistics, we should ideally construct the grid positions so that equal numbers of packets are absorbed in each cell. Since the probability of absorbing a photon packet is proportional to the optical depth, we choose equal radial optical depth grid locations,

$$\frac{r_i}{R_{\mathrm{dust}}}=\{\begin{array}{cc}i\mathrm{\Delta }r+1\hfill & \text{(}p=0\text{),}\hfill \\ N_f/(N_fi)\hfill & \text{(}p=2\text{),}\hfill \end{array}$$
(17)
where $`N_r=200`$ is the total number of cells we used, $`\mathrm{\Delta }r=(R_{\mathrm{max}}/R_{\mathrm{dust}}1)/N_r`$, and $`N_f=N_r/(1R_{\mathrm{dust}}/R_{\mathrm{max}})`$. Integrating the density, eq. (10), over the cell volume to obtain the mass, we find from eq. (6) that the temperature in each grid cell is given by

$$T_i^4=T_{}^4\{\begin{array}{cc}\frac{N_iN_r(R_{}/R_{\mathrm{dust}})^2[(\kappa _{1\mu \mathrm{m}}+\sigma _{1\mu \mathrm{m}})/\kappa _P(T_i)]}{4N_\gamma \tau _{1\mu \mathrm{m}}[(i^2i+1/3)\mathrm{\Delta }r^2+(2i1)\mathrm{\Delta }r+1]}\hfill & \text{(}p=0\text{),}\hfill \\ \frac{N_i(N_fi)(N_fi+1)(R_{}/R_{\mathrm{dust}})^2(1R_{\mathrm{dust}}/R_{\mathrm{max}})}{4N_\gamma \tau _{1\mu \mathrm{m}}N_f[\kappa _P(T_i)/(\kappa _{1\mu \mathrm{m}}+\sigma _{1\mu \mathrm{m}})]}\hfill & \text{(}p=2\text{),}\hfill \end{array}$$
(18)
where we have used $`L=4\pi R_{}^2\sigma T_{}^4`$ for the stellar luminosity.

We then proceed with the radiation transfer, temperature calculation, and reemission as described in Section 2 until all packets exit the envelope. When the packets escape, they are placed into $`N_\nu =1024`$ uniform frequency bins,

$$\nu _k=\nu _{\mathrm{max}}\frac{k}{N_\nu },$$
(19)
where $`h\nu _{\mathrm{max}}=16kT_{}`$. The width of each bin $`\mathrm{\Delta }\nu =\nu _{\mathrm{max}}/N_\nu `$. Since the envelope is spherically symmetric, the observed flux $`F_{\nu _k}=N_kE_\gamma /4\pi d^2\mathrm{\Delta }t\mathrm{\Delta }\nu `$, where $`N_k`$ is the number of packets in the $`k^{\mathrm{th}}`$ frequency bin, and $`d`$ is the observer’s distance from the star. Normalizing to the total flux, $`F=L/4\pi d^2`$, the SED is given by

$$\left(\frac{\nu F_\nu }{F}\right)_k=(k1/2)\frac{N_k}{N_\gamma }.$$
(20)
The factor $`(k1/2)`$ arises from using the frequency at the center of the bin.

In Figures 2a and 2b we show the results of our simulation compared with the output of one of the codes tested by Ivezić et al. This code, called DUSTY, is publicly available and is described in Ivezić & Elitzur (1997). We see that our Monte Carlo calculations reproduce both the correct temperature structure and SED. Note the large errors at the longest and shortest wavelengths in the Monte Carlo calculations. At these wavelengths, the number of emerging packets is small, resulting in large errors (the relative error in each flux bin equals $`1/\sqrt{N_k}`$, owing to the Poisson sampling statistics inherent in Monte Carlo simulations). The excellent agreement of the comparisons shown in Figures 2a and 2b (the differences are smaller than the numerical error of the DUSTY calculations) demonstrates the validity of our temperature correction procedure described in Section 2.

Now that we have verified our our basic radiative equilibrium algorithm, we can proceed to investigate the temperature structure and SEDs of other geometries. Owing to the inherently three-dimension nature of Monte Carlo simulations (even our 1-D spherically symmetric code internally tracks the photon packets in three dimensions), our code is readily modified for arbitrary geometries. We now show the results of an application to axisymmetric circumstellar environments.

## 4 Axisymmetric 2-D Calculations

For the purpose of this illustrative simulation, we adopt a stellar blackbody with $`T_{}=3500\mathrm{K}`$, envelope inner radius $`R_{\mathrm{dust}}=10R_{}`$, and a simple ellipsoidal parameterization for the circumstellar density. The isodensity contours are elliptical with $`a/b`$ being the ratio of the semi-major to semi-minor axis. The density is given by

$$\rho =\rho _0\frac{(R_{\mathrm{dust}}/r)^2}{1+f^2\mathrm{cos}^2\theta },$$
(21)
where $`\rho _0`$ is given by eq. (11; $`p=2`$ case), $`\theta `$ is the polar angle, and the “flattening factor”

$$f=\sqrt{(a/b)^21}.$$
(22)
Note that the equatorial to polar density ratio at a given radius is $`\rho _{\mathrm{eq}}/\rho _{\mathrm{pole}}=(a/b)^2`$.

As before, we divide the circumstellar environment into cells with $`N_r=200`$ radial and $`N_\mu =20`$ latitudinal grid points. Note that the envelope is symmetric about the equator, so we combine the cells below the equator with their counterparts above the equator. This is automatically accomplished by using $`\mu _j=\mathrm{cos}(\theta _j)`$ for the latitudinal grid point coordinate. Spacing the grid so that the radial and latitudinal optical depths are the same for each cell, we find

$`r_i`$ $`=`$ $`{\displaystyle \frac{N_f}{N_fi}}R_{\mathrm{dust}},`$ (23)
$`\mu _j`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1+(a/b)^2\mathrm{tan}^2(\pi j/2N_\mu )}}}.`$ (24)
With these cell coordinates, equation (6) for the temperature in cell $`(i,j)`$ becomes

$$T_{i,j}^4=T_{}^4\frac{N_{i,j}(N_fi)(N_fi+1)(R_{}/R_{\mathrm{dust}})^2(1R_{\mathrm{dust}}/R_{\mathrm{max}})f}{4N_\gamma \tau _V^{\mathrm{eq}}N_f[\kappa _P(T_{i,j})/(\kappa _V+\sigma _V)][\mathrm{tan}^1(f\mu _{j1})\mathrm{tan}^1(f\mu _j)]},$$
(25)
where $`N_{i,j}`$ is the number of packets absorbed in the cell, and we have chosen $`\lambda =5500`$Å ($`V`$-band) for our equatorial radial optical depth parameter, $`\tau _V^{\mathrm{eq}}`$.

For the dust opacity, we adopt a standard MRN interstellar grain mixture (Mathis, Rumpl, & Nordsieck 1977), using optical constants from Draine & Lee (1984). Figure 3 shows the opacity and albedo in graphical form. The data for this figure is available in tabular form<sup>1</sup><sup>1</sup>1ftp://gradj.pa.uky.edu/dusty/distribution/ism-stnd.dat from the DUSTY Web site, http://www.pa.uky.edu/$``$moshe/dusty. Note the prominent silicate absorption features at $`10\mu `$m and $`18\mu `$m. For the current demonstration, we assume that the scattering is isotropic, but we can easily accommodate any phase function, analytic or tabulated when calculating the emergent SED.

Unlike the spherically symmetric benchmarks, the SED now depends on viewing angle, so we must bin the escaping packets both in direction and frequency. Since the envelope is axisymmetric, the observed flux only depends on the inclination angle, $`i`$, of the envelope symmetry axis (i.e., the stellar rotation axis). To obtain approximately equal numbers of escaping packets, we choose $`N_{\mathrm{inc}}=10`$ direction bins, with equal solid angles, $`\mathrm{\Delta }\mathrm{\Omega }=4\pi /N_{\mathrm{inc}}`$. The escaping directions (inclination angles), $`i_l`$, for these bins are given by

$$\mu _l^{\mathrm{esc}}=l/N_{\mathrm{inc}},$$
(26)
where $`\mu _l^{\mathrm{esc}}=\mathrm{cos}(i_l)`$. In addition to the direction bins, we use the frequency bins given by eq. (19).

Now that we have specified both the direction and frequency bins, the observed flux is $`F_{\nu _{k,l}}=N_{k,l}E_\gamma /\mathrm{\Delta }\mathrm{\Omega }d^2\mathrm{\Delta }t\mathrm{\Delta }\nu `$, where $`N_{k,l}`$ is the number of escaping packets in the $`(k,l)`$ bin. Normalizing to the bolometric flux F gives the emergent SED,

$$\left(\frac{\nu F_\nu }{F}\right)_{k,l}=\frac{(k1/2)N_{k,l}N_{\mathrm{inc}}}{N_\gamma }.$$
(27)
We now choose to investigate the SEDs produced by two density structures with different degrees of flattening. The first has a density ratio $`\rho _{\mathrm{eq}}/\rho _{\mathrm{pole}}=1000`$ to represent a disk-like structure. The second has a density ratio $`\rho _{\mathrm{eq}}/\rho _{\mathrm{pole}}=10`$, which is mildly oblate, representative of an infalling protostellar envelope. For each density structure, we perform optically thick (in the mid-IR) and optically thin calculations. The optically thick calculations have an equatorial $`V`$-band optical depth $`\tau _V^{\mathrm{eq}}=200`$, while the optically thin calculations have $`\tau _V^{\mathrm{eq}}=20`$. Table 2 summarizes these density structures. For comparison, we also have performed calculations for $`p=2`$ spherically symmetric models containing the same total mass as our disk and envelope densities. Keeping the same inner and outer radii, the optical depth for the equivalent spherical model is

$$\tau _V^{\mathrm{sp}}=\tau _V^{\mathrm{eq}}\frac{\mathrm{tan}^1f}{f}.$$
(28)
### 4.1 Disk Model

Figure 4a shows the incident stellar spectrum and the emergent SED as a function of viewing angle for the disk-like model, as well as the result of the equivalent spherically symmetric calculation. For both disk simulations, when viewing the system pole-on, we are looking through optically thin circumstellar dust (see Table 2) and can see the star at optical wavelengths. In the IR, there is an excess arising from the circumstellar disk, which reprocesses the stellar radiation. Note that, for pole-on viewing, we see the silicate features in emission, since the disk is optically thin in the vertical direction. As we go to higher viewing angles, the optical depth to the central star increases, and as a result, the star becomes more extincted in the optical region. Note that at almost edge-on viewing, a “shoulder” appears around $`1\mu `$m for the optically thick case. This arises due to the dominant effects of scattering the stellar radiation at these wavelengths. These scattering shoulders are also present in the axisymmetric calculations presented by Efstathiou & Rowan-Robinson (1991), Sonnhalter et al. (1995), Men’schikov & Henning (1997), and D’Alessio et al. (1999). At wavelengths longer than about $`1\mu `$m, the albedo begins to drop rapidly (see Fig. 3), and the disk thermal emission begins to dominate, so the shoulder terminates. Beyond $`30\mu `$m, the envelope becomes optically thin, so the spectrum is independent of inclination and is dominated by the dust emission.

The two dimensional temperature structure for the disk-like models is displayed in Figure 4b. We see that at the inner edge of the envelope ($`R_{\mathrm{dust}}`$) there is little variation of the temperature with latitude, while at large radii there is a clear latitudinal temperature gradient, with the dust in the denser equatorial regions being cooler than dust at high latitudes.

In the polar region, the material is optically thin to the stellar radiation, so it heats up to the optically thin radiative equilibrium temperature. This temperature has a power law behavior, $`Tr^{0.4}`$ for dust opacity ($`\kappa \lambda ^1`$). As can be seen in Figure 4b, the polar temperature does indeed have a power law decrease with a slope of approximately $`0.4`$.

In contrast, the disk only displays this power law behavior at large radii. At the inner edge of the disk, the disk sees the same mean (stellar) intensity as is present in the polar region ($`J=WB`$, where $`W=0.5\{1[1(R/r)^2]^{1/2}\}`$ is the dilution factor). Consequently, the inner edge of the disk heats up to the optically thin radiative equilibrium temperature — the same as the polar temperature. However at larger radii, the opaque material in the inner regions of the disk shields the outer regions from direct heating by the stellar radiation. This shielding reduces the mean intensity ($`J<WB`$). As a result, the outer disk is only heated to a fraction of the optically thin radiative equilibrium temperature. Thus the equatorial region is cooler than the polar region. Eventually at large enough radii, the disk becomes optically thin to the heating radiation. At that point, it sees a radially streaming radiation field from an effective photosphere that has a much lower temperature than the star. From that point outward, the disk temperature displays an optically thin power law decrease with a slope that parallels the polar temperature.

The spherically symmetric calculation overestimates both the emergent flux and the disk temperature. In the optically thin limit, the IR continuum is proportional to the mass of the circumstellar material, so one would expect that a spherically equivalent mass would reproduce the long wavelength spectrum. Recall however that we can see the star at pole-on viewing angles. This implies that the disk does not reprocess the entire bolometric luminosity of the star. Consequently, the IR excess is less than that in the spherical model. Note that one can nonetheless reproduce the edge-on SED using a spherical model if we allow the density power law to depart from $`1/r^2`$ and change the size of the spherical envelope. This also was noted by Sonnhalter et al. (1997), who found that they could fit the IR continuum of their disk models by changing the radial dependence of the circumstellar density and the outer radius for the same total envelope mass.

### 4.2 Envelope Model

The spectral energy distribution and temperature structure for the oblate envelope are shown in Figures 5a and 5b. The optically thin envelope displays significant extinction of the star at all viewing angles. Close to edge-on, scattering shoulders appear around $`1\mu `$m that are similar to those seen in the disk model (Fig. 4a). Finally, because the model is optically thin in the mid-IR, the silicate feature is always in emission, and the SED is independent of viewing angle for wavelengths longer than a few microns.

For the denser envelope, the star is extremely faint at optical wavelengths. Along with the increased extinction in the optical, the scattering shoulders are less prominent, due to the thermal emission by the envelope. In the mid-IR, the envelope is optically thick edge-on and optically thin pole-on. Consequently, the silicate features go from absorption to emission as the viewing angle changes from edge-on to pole-on. The general shape of the SEDs, which now peak in the mid-IR for all inclinations, are reminiscent of spectra from embedded (Class I) T Tauri stars ,which are commonly modeled using flattened axisymmetric dusty envelopes (e.g., Adams, Lada, & Shu 1987; Kenyon, Calvet, & Hartmann 1993; Men’shchikov & Henning 1997; D’Alessio, Calvet, & Hartmann 1997). Of these T Tauri simulations, only Efstathiou & Rowan-Robinson (1991) have performed an exact calculation for the SED and circumstellar temperature for the Terebey, Shu, & Cassen (1984) collapse model.

The temperature structure for our envelope models (Fig. 5b) is qualitatively similar to the temperature structure of the disk models (Fig. 4b). The temperature at the inner edge of the envelope is independent of latitude, while at larger radii the equatorial regions are cooler than the polar regions. The primary difference is that the latitudinal temperature gradient is not as extreme.

Finally, we note that the equivalent spherically symmetric model better reproduces the far IR spectrum of the envelope models than the disk models. This is because the envelope models reprocess the entire bolometric luminosity, while in disk models, some of the stellar luminosity escapes through the polar region.

## 5 Discussion

We have developed a temperature correction procedure for use in Monte Carlo radiation transfer codes. We have tested our method against other spherically symmetric benchmark codes and successfully matched their results. After verifying our method, we applied it to obtain sample temperature distributions and SEDs for 2-D axisymmetric disk-like models and mildly oblate envelopes. These simulations illustrate the important role envelope geometry can play when interpreting the SEDs of embedded sources.

The primary limitation of our temperature correction procedure is that it only applies if the opacity is independent of temperature. This is not true for free-free opacity or hydrogen bound-free opacity, but it is true for dust opacity, which is the dominant opacity source in many astrophysical situations. The reason our method is limited to temperature-independent opacities is that two associated problems occur if the opacity varies when the cell’s temperature changes. The first is that the cell will have absorbed either too many or too few of the previous packets passing through the cell. The second is the associated change in the interaction locations of the previous packets, which implies that the paths of the previous photon packets should have been different. These problems do not occur if the opacity is independent of temperature.

Lucy (1999) has proposed a slightly different method to calculate the equilibrium temperature. Instead of sampling photon absorption, he directly samples the photon density (equivalent to the mean intensity) by summing the path length of all packets that pass through a cell. Typically more packets pass through a cell than are absorbed in the cell, so this method potentially produces a more accurate measurement of the temperature for a given total number of packets, especially when the envelope is very optically thin. The disadvantage of his method is that it requires iteration to determine the envelope temperature. In principle one could partially combine both methods. First, run our simulation, adding pathlength counters to each cell. After running all packets, use the pathlength information to calculate a final temperature. This will provide a more accurate temperature that can be used to calculate the source function. After obtaining the source function for each cell, it is a simple matter to integrate the transfer equation to obtain the SED.

Another limitation of the Monte Carlo method is that it is not well suited to envelopes with very high optical depths ($`\tau >100`$ – 1000), unless there is an escape channel for the photons. For example, geometrically thin disks can be extremely optically thick in the radial direction but optically thin in the polar direction, providing an escape channel for the photons. Similarly, dense envelopes can be very optically thick to the illumination source, but optically thin to the reprocessed radiation, which is a another escape channel. In the event no such escape channels exist, one must turn to other methods. We are currently investigating how to couple Monte Carlo simulation in the optically thinner regions with other methods in the optically thick interior.

The focus of the present paper has been on the development and implementation of the temperature correction procedure. For this reason, we have made several simplifying assumptions regarding the circumstellar opacity and geometry. For example, we have assumed a single temperature for all dust grains regardless of their size and composition. This differs from other investigations where different types of grains can have different temperatures at the same spatial location (e.g., the spherically symmetric radiation transfer code developed by Wolfire & Cassinelli 1986). Similarly, we have not implemented a procedure to solve for the location of the dust destruction radius, which will be different for grains of differing size and composition. These issues are currently under investigation.

With the speed of today’s computers, Monte Carlo radiation transfer simulations can be performed in a reasonably short time. The 2-D simulations presented in this paper employed $`10^8`$ packets, requiring about two hours of CPU time; the spherically symmetric cases used $`10^6`$ packets, requiring about one minute. So for continuum transfer, Monte Carlo simulation is proving to be a very powerful technique for investigating arbitrary density structures and illuminations.

We would like to thank Barbara Whitney and Mário Magalhães for many discussions relating to this work. An anonymous referee also provided thoughtful comments which helped improved the clarity of our presentation. This work has been funded by NASA grants NAG5-3248, NAG5-6039, and NSF grant AST-9819928.


text:
# 1 Introduction

## 1 Introduction

The well known secular shift of the pericenter $`\omega `$ of a test body induced by the Schwarzschild metric generated by a static, spherically symmetric distribution of mass \[Misner et al., 1973; Ciufolini and Wheeler, 1995\] is one of the classical tests of Einstein’s General Relativity.

### 1.1 The Solar System tests

The advance of pericenter, which can be expressed in terms of $`\frac{2+2\gamma \beta }{3}`$, where $`\beta `$ and $`\gamma `$ are the Eddington–Robertson–Schiff PPN parameters \[Will, 1993\], has been measured in the Solar System, on the Mercury perihelion, by radar signals transmitted from Earth to Mercury and back to Earth \[Shapiro et al., 1972\], yielding $`\frac{2+2\gamma \beta }{3}=1.005\pm 7\times 10^3`$. The contribution of possible systematic errors raises the uncertainty to about $`2\times 10^2`$. From the analysis of ten years of data the relative error published in \[Shapiro et al., 1976\] amounts to 5$`\times 10^3`$. The major sources of systematic error in this measurement are the poorly-known variations in the topography of the planet’s surface and the uncertainties in the radar scattering law \[Shapiro, 1990; Pitjeva, 1993\]. For a review of periastron advance measurements on binary pulsars see \[Will, 2001\] and references therein.

Several space missions have been proposed in order to measure, among other effects, this phenomenon by using a variety of techniques: the most recent and promising are SORT, IPLR and ASTROD \[Ni, 2001\].

More accurate measurements of $`\frac{2+2\gamma \beta }{3}`$ might be performed in the future by means of the ESA BepiColombo Mercury orbiter \[Balogh et al., 2000; Milani et al., 2001a; Milani et al., 2001b\]. In \[Milani et al., 2001b\] it is claimed that the BepiColombo mission, which should be launched in 2009, should allow to measure $`\gamma `$ and $`\beta `$ with a relative accuracy of the order of $`2\times 10^6`$.

### 1.2 The proposed experiment

In this paper we propose to perform a measurement of the general relativistic “gravitoelectric” pericenter advance in the gravitational field of the Earth by using some suitable combinations of the orbital residuals of existing or proposed, spherical, passive, geodetic laser-ranged satellites \[Iorio, 2002a\] with particular emphasis to LAGEOS and LAGEOS II in order to exploit the relevant experience obtained with the gravitomagnetic Lense-Thirring measurements in \[Ciufolini et al., 1996; 1997; 1998\]. In \[Ciufolini and Matzner, 1992\] a determination of the LAGEOS general relativistic gravitoelectric perigee shift in the field of the Earth is reported, but the accuracy amounts only to $`2\times 10^1`$.

This paper is organized as follows: in section 2 we compare the present experimental accuracy in satellite laser-ranging measurements with the relativistic expression of the perigee shift of the two LAGEOS satellites. In section 3 we analyze some of the most important sources of systematic errors. The accuracy achievable in measuring $`\beta `$ and $`\gamma `$ is also discussed. By detecting the relativistic perigee rate of the proposed LARES \[Ciufolini, 1998\] laser-ranged satellite, in conjunction with the node and the perigee of LAGEOS II, it might be possible to further improve the accuracy in the measurement of these parameters over a long enough time span. This is the subject of section 4 together with an analysis of the potential role of the existing laser-ranged satellite Starlette. Section 5 is devoted to the conclusions.

## 2 The relativistic perigee precession of LAGEOS type satellites

As known, in the slow-motion and weak-field approximation, the Schwarzschild metric generated by a static, spherically symmetric distribution of mass-energy induces an additional post-Newtonian “gravitoelectric force” which acts on the orbit of a test body by shifting its pericenter; in the PPN formalism the pericenter rate can be written as:

$$\dot{\omega }_{\mathrm{GE}}=\frac{3nGM}{c^2a(1e^2)}\times \frac{2+2\gamma \beta }{3},$$
($`1`$)
in which $`G`$ is the Newtonian gravitational constant, $`c`$ is the speed of light in vacuum, $`M`$ is the mass of the central object, $`a`$ and $`e`$ are semimajor axis and eccentricity, respectively, of the orbit of the test body and $`n=\sqrt{GM/a^3}`$ is its mean motion. In the following we define $`\nu \frac{2+2\gamma \beta }{3}`$. General Relativity predicts that the perigee shifts for LAGEOS and LAGEOS II amount to 3,312.35 milliarcseconds per year (mas/y in the following) and 3,387.46 mas/y respectively.

The actual experimental precision allows for detecting such rates for both the LAGEOS satellites \[Ciufolini, 1996\]. Indeed, for the perigee the observable quantity is $`rea\dot{\omega }`$; at present the RMS of the residuals of the best orbital fits amounts to $`\delta r_{\mathrm{exp}}1`$ cm, or less, for the two LAGEOS satellites, over several orbits and with a set of force models (i.e. not including modeling errors). Since the LAGEOS eccentricity is $`e_\mathrm{I}=4.5\times 10^3`$ the accuracy in detecting the perigee is $`\delta \omega _{\mathrm{exp}}^\mathrm{I}\delta r_{\mathrm{exp}}/(e_\mathrm{I}a_\mathrm{I})37`$ mas. So that, over 1 year the relative accuracy in the measurement of the relativistic perigee shift would be of the order of $`1\times 10^2`$. For LAGEOS II this measurement accuracy is better than for LAGEOS. Indeed, the LAGEOS II eccentricity is $`e_{\mathrm{II}}=1.4\times 10^2`$, thus $`\delta \omega _{\mathrm{exp}}^{\mathrm{II}}12`$ mas; this may yield an accuracy of the order of $`3\times 10^3`$ over 1 year. Over 8 years, by assuming the same fit error as before, it would amount to $`4\times 10^4`$. These considerations rule out the possibility of directly using the perigee of the other existing, spherical, geodetic laser-ranged satellites Etalon-1, Etalon-2, Ajisai, Stella, Westpac-1 because their eccentricities are even smaller than that of LAGEOS. On the contrary, Starlette has an eccentricity of the order of $`2\times 10^2`$; however, since it orbits at a lower altitude is more sensitive than the LAGEOS satellites to atmospheric drag and to Earth’s zonal harmonics, so that would be difficult to process its data at an acceptable level of accuracy (See also section 4). Accordingly, in order to detect the gravitoelectric relativistic shift in the gravitational field of Earth, we will focus on the perigee of LAGEOS and especially of LAGEOS II .

## 3 The systematic errors

The perigee of an Earth satellite is a ”non-clean” orbital element in the sense that it is affected by a large number of aliasing classical forces inducing systematic errors which must be carefully identified and analyzed.

### 3.1 The static geopotential error

The most important source of systematic errors in such a measurement is represented by the mismodeling of the classical perigee precession induced by the errors in even zonal harmonics of the Earth’s gravitational field \[Kaula, 1966\]. This error is really critical because the resulting aliasing trend cannot be removed from the data and nothing can be done about it apart from assessing the related error as reliably as possible.

By using the covariance matrix of the EGM96 Earth gravity model \[Lemoine et al., 1998\] and adding the correlated terms in a root-sum-square fashion up to degree $`l=20`$, we obtain for LAGEOS a systematic error $`\delta \nu /\nu _{\mathrm{zonals}}=8.1\times 10^3`$ whereas for LAGEOS II we have $`\delta \nu /\nu _{\mathrm{zonals}}=1.5\times 10^2`$. Since the major source of error lies in the uncertainty of the first two mismodeled even zonal harmonics $`\delta J_2`$ and $`\delta J_4`$, in order to eliminate most of the static and dynamical even zonal terms of the geopotential, following the method in \[Ciufolini, 1996\] for the measurement of the gravitomagnetic Lense-Thirring effect, we search for suitable combinations of the orbital residuals of the existing laser-ranged satellites. In Tab.1 we report the most promising combinations: their general form is

$$\delta \dot{\omega }^{\mathrm{II}}+\underset{\mathrm{i}=1}{\overset{\mathrm{N}}{}}c_\mathrm{i}\delta \dot{\mathrm{\Omega }}^\mathrm{i}+c_{\mathrm{N}+1}\delta \dot{\omega }^\mathrm{I}=x_{\mathrm{GR}}\nu ,$$
($`2`$)
in where $`N`$ is the number of the nodes of different laser-ranged satellites employed, $`x_{\mathrm{GR}}`$ is the slope, in mas/y, of the relativistic trend to be measured, in Tab.1 $`\delta \nu /\nu _{\mathrm{zonals}}`$ is the systematic error induced by the even zonal harmonics up to degree $`l=20`$ calculated with the EGM96 covariance matrix. The coefficients $`c_i`$ are determined in terms of the orbital parameters of the satellites entering the combinations.

It is important to stress that the use of the LAGEOS satellites, due to their altitude, makes our measurement substantially insensitive to the errors in the zonal harmonics of degree $`l>20`$, so that our estimates of $`\delta \nu /\nu _{\mathrm{zonals}}`$ presented here are valid even in the case that the EGM96 covariance matrix for degrees higher than $`l>20`$ would not be accurate enough.

In addition to LAGEOS and LAGEOS II we have only considered Ajisai \[Iorio, 2002a\] since it is well tracked, contrary for example to the Etalon satellites, and it would be less demanding than for the other satellites to reduce its laser-ranged data to a level of accuracy comparable to that of LAGEOS and LAGEOS II. Moreover, the other laser-ranged satellites orbit at lower altitudes, therefore they are more sensitive to terms of the geopotential of higher degree. Consequently, as confirmed by numerical calculations, the inclusion of their data in the combined residuals would increase $`\delta \nu /\nu _{\mathrm{zonals}}`$111About combination 3, it should be noted that, in order to obtain more reliable and accurate estimates of the systematic error due to the even zonal harmonics of the geopotential, the calculation of the error should include terms with degree higher than $`l=20`$ due to the sensitivity of Ajisai to such higher degree terms. .

It is important to point out that the uncertainties quoted in Tab.1 for $`\delta \nu /\nu _{\mathrm{zonals}}`$ will be reduced in the near future when the new gravity models from the CHAMP and GRACE missions will be released \[Pavlis, 2002a\].

### 3.2 The time–dependent systematic errors

In regard to the evaluation of the impact of the other sources of systematic errors, the role of the coefficients $`c_i`$ entering the combinations has to be taken into account, indeed, according to their magnitude, they can reduce or emphasize the effects of the perturbing forces acting on the orbital elements weighted by them.

For example, for the combination 2 of Table 1, over 1 year, the impact of the error in measuring the perigee rate of LAGEOS amounts to about $`2.5\times (3.1\times 10^3)=7.7\times 10^3`$ while for combination 3 it is about $`1.37\times (4.7\times 10^3)=6.4\times 10^3`$.

In this paper we only consider in detail the combination 1 of Table 1

$$\delta \dot{\omega }^{\mathrm{II}}0.87\times \delta \dot{\mathrm{\Omega }}^{\mathrm{II}}2.86\times \delta \dot{\mathrm{\Omega }}^\mathrm{I}=3,387.46\times \nu ,$$
($`3`$)
so that we exploit the insight acquired with the Lense-Thirring LAGEOS experiment \[Ciufolini et al., 1996; 1997; 1998\].

The long-period harmonic perturbations, according to their periods $`P`$ and to the adopted observational time span $`T_{\mathrm{obs}}`$, may turn out to be less insidious than the mismodeled secular perturbations due to the zonal harmonics of the geopotential since, if $`P<T_{\mathrm{obs}}`$ and $`T_{\mathrm{obs}}=nP,n=1,2,\mathrm{}`$ they average out; if their periods are shorter than the time span they can be viewed as empirically fitted quantities which can be subsequently removed from the signal.

The results and estimates recently obtained for the LAGEOS gravitomagnetic experiment can be applied to our proposed measurement. Indeed, the solid Earth and ocean tidal perturbations on LAGEOS and LAGEOS II and their impact on the Lense–Thirring effect have been studied in \[Iorio, 2001; Iorio and Pavlis, 2001; Pavlis and Iorio, 2002\], while the role played by the non–gravitational perturbations has been analyzed in \[Lucchesi, 2001; 2002\]. Regarding the atmospheric drag, whose impact might be thought of as a serious drawback, especially for the perigee of LAGEOS II, in \[Ciufolini et al., 1997\] it is shown that it essentially averages out over an orbital revolution.

In Tab.2 we summarize the results obtained with the proposed combination using LAGEOS and LAGEOS II for a 8-year time span with orbital fits, arcs, each of 7 days. In assessing the total systematic error we have accounted for the fact that the gravitational errors are not independent, indeed we have simply summed them up; then we have added these uncertainties and the other independent systematic errors in a root-sum-square fashion. In Tab.2 we can observe that the error budget is mainly dominated by the even zonal harmonics and by the LAGEOS II perigee measurement error. In regard to the non–gravitational errors, according to the results of \[Lucchesi, 2001\] for the direct solar radiation pressure and the Earth’s albedo and to the results of \[Lucchesi, 2002\] for the thermal thrust perturbations and the asymmetric reflectivity, the corresponding uncertainty would amount to almost $`1\times 10^2`$ over 7 years. In obtaining this result a very pessimistic approach has been adopted by assuming a mismodeling of 20$`\%`$ for all the perturbing effects except for the direct solar radiation pressure which has been assumed to have an uncertainty at the 0.5$`\%`$ level. However, we stress that only the Earth’s thermal thrust, or Yarkovski-Rubincam effect, induces a mismodeled linear trend whose impact would amount to about $`1\times 10^4`$: the other non–gravitational forces are time–varying with known periodicities and can be fitted and removed from the signal as done for the tidal perturbations.

Forthcoming improvements in the Earth gravity field modeling achievable with the CHAMP and GRACE data have been recently extensively studied and reported in \[Pavlis, 2002a\]. They will yield important improvements in the accuracy of our experiment. We finally observe that we already have the required laser-ranging data since, at present, we have almost 10 years of LAGEOS and LAGEOS II data.

### 3.3 The parameters $`\beta `$ and $`\gamma `$

The results obtained for the combination 1 examined here, together with $`\delta \eta =8\times 10^4`$ reported in \[Anderson and Williams, 2001\] for the combination $`\eta =4\beta \gamma 3`$ \[Nordtvedt, 1968; 1991\], would allow to measure $`\beta =\frac{2}{7}\eta +\frac{3}{7}\nu +\frac{4}{7}`$ at the level $`3\times 10^3`$ independently of other measurements of $`\gamma `$; this result, which is of the same order of magnitude of that obtained with radar ranging \[Shapiro, 1990\], should be compared with the most recent $`\delta \beta =4.7\times 10^4`$, obtained with LLR data \[Anderson and Williams, 2001\].

The parameter $`\gamma `$, that may be written as $`\gamma =\frac{1}{7}\eta +\frac{12}{7}\nu \frac{5}{7}`$, would be measured less precisely: $`\delta \gamma =1.2\times 10^2`$. However, $`\gamma `$ can be directly measured via light deflection or radar time delay. For example, in \[Frœschlé et al., 1997\] the result $`\gamma =0.997\pm 3\times 10^3`$ is based on the astrometric observations of electromagnetic waves deflection in the visible frequency. By using the Viking time delay \[Reasenberg et al., 1979\], $`\gamma =1.000\pm 2\times 10^3`$ was obtained; the quoted uncertainty allows for possible systematic errors. A more recent measurement based on the time delay with the NEAR spacecraft \[Elliott et al., 1998\] claims an uncertainty $`10^3`$.

For an updated overview of the present status of the measurements of the PPN parameters, see \[Will, 2001\]

## 4 The role of LARES and Starlette

The proposed LARES satellite \[Ciufolini, 1998\] will enable us to perform several general relativistic tests and geophysical measurements. In this section we investigate its role for the measurement of the gravitoelectric perigee advance that would be well detectable on LARES thanks to its proposed relatively large eccentricity: $`e_{\mathrm{LR}}0.04`$.

A possible observable is the following combination of orbital residuals

$$\delta \dot{\omega }^{\mathrm{II}}+c_1\delta \dot{\omega }^{\mathrm{LR}}+c_2\delta \dot{\mathrm{\Omega }}^{\mathrm{II}}=x_{\mathrm{GR}}\nu ,$$
($`4`$)
with

$`c_1`$ $`=`$ $`4.71,`$ ($`5`$)
$`c_2`$ $`=`$ $`2.26,`$ ($`6`$)
$`x_{\mathrm{GR}}`$ $`=`$ $`12,117.4\mathrm{mas}/\mathrm{y}.`$ ($`7`$)
The coefficients $`c_i`$ depend on the orbital parameters of the selected satellites and have been calculated corresponding to the proposed orbital parameters of LARES: $`a_{\mathrm{LR}}=12,270`$ km, $`i_{\mathrm{LR}}=70`$ deg and $`e_{\mathrm{LR}}=0.04`$. Let us now calculate the systematic errors induced by the even zonal harmonics of the geopotential, according to the EGM96 gravity model. Eq. ($`4`$) would allow to cancel out the errors due to the first four even zonal harmonics; the impact of the remaining ones, with degree higher than four, amounts to an uncertainty of about $`6\times 10^3`$. This error will further reduce when the new gravity models from the CHAMP and GRACE missions will be available. An interesting feature of eq. ($`4`$) is that the impact of the non–gravitational perturbations would produce an uncertainty of about $`3.9\times 10^3`$ over a time span of 7 years, according to the very conservative estimates discussed above; in particular the Earth thermal thrust, i.e. the so-called Yarkovski–Rubincam effect, would induce a mismodeled secular trend with an uncertainty of only $`6.4\times 10^5`$. The inclusion of LARES would thus represent an improvement with respect to the LAGEOS and LAGEOS II scenario previously outlined, especially when, in the near future, the impact of the mismodeled non–gravitational perturbations will increase relatively to the gravitational perturbations thanks to the more accurate Earth’s gravity field models.

Among the existing laser-ranged satellites, Starlette, with its relatively large eccentricity of 0.0204, seems to be a natural candidate for the measurement of the perigee advance. Unfortunately, it orbits at an altitude much lower than the LAGEOS satellites ($`a_{\mathrm{Str}}=7,331`$ km) and thus the mismodeling of the even zonal harmonics of the geopotential would induce large systematic errors in the measurement of the perigee rate. It turns out that the inclusion of the perigee of Starlette in the combined residuals if, on the one hand, provides a further observable and thus would allow to cancel out the error of an additional mismodeled even zonal harmonic, on the other hand it raises the total error due to the geopotential, even with the LARES data. Indeed, in a combination including the nodes of LAGEOS and LAGEOS II and the perigees of LAGEOS II and Starlette, the relative error due to the geopotential amounts to only $`6.2\times 10^2`$. This rules out the possibility of using the Starlette data, at present level of knowledge of the Earth’s gravitational field, to improve the measurement of the relativistic perigee advance.

## 5 Conclusions

The main features of our proposed measurement of perigee advance using laser-ranged satellites are the following:

* It will provide a test of a general relativistic effect, i.e. the pericenter advance, in the gravitational field of Earth, complementary to similar tests in the field of Sun. Indeed, such a measurement in the gravitational field of Earth would place stronger limits on possible, very weak, Yukawa-type components of the gravitational interaction with a finite range of the order of $`10^4`$ km \[Iorio, 2002b\].
* For a suitable combination of the orbital residuals of the nodes of LAGEOS and LAGEOS II and perigee of LAGEOS II, the relative error due to the mismodeled classical precessions by the static even zonal harmonics of the geopotential, which is the most insidious source of systematic error, amounts, according to the EGM96 Earth gravity model, to about $`6\times 10^3`$. The total error in our proposed combination of orbital elements, also due to the other time–dependent systematic errors, amounts, over an observational time span of eight years, to almost $`7\times 10^3`$.
* In the near future the knowledge of the Earth gravitational field will increase thanks to the data from the CHAMP and GRACE missions. Moreover, the rms accuracy in measuring the position of the LAGEOS satellites should reach the few mm level. Consequently, there will be an increase in the accuracy of our measurement.
* In regard to the terrestrial space environment, which induces a large number of orbital perturbations on the satellites, it is very well studied and the action of these perturbations, especially on the LAGEOS satellites, is very well modeled. The data processing would be reliably based on the large experience collected in several years of analysis of the LAGEOS orbits \[Pavlis, 2002b\]. Moreover, it would be relatively fast and easy to include the future improvements in the force modeling in our analysis, e.g. the improvements due to a better experimental and theoretical knowledge of the rotation rate and spin axis orientation of LAGEOS type satellites.
* The proposed measurement, as an additional outcome, from $`\frac{2+2\gamma \beta }{3}`$ in conjunction with $`\eta =4\beta \gamma 3`$ from Lunar Laser Ranging, would allow to obtain the PPN parameters $`\gamma `$ and $`\beta `$ in the gravitational field of Earth (see above) with an accuracy at the level of $`10^210^3`$. This accuracy would improve with the new Earth’s gravitational models from CHAMP and GRACE.
* Starlette data, although this satellite has a relatively eccentric orbit, should not be included in the combined residuals. Indeed, Starlette would raise the systematic error in the geopotential due to its lower altitude with respect to the LAGEOS satellites and thus due to its higher sensitivity to the even zonal harmonics of the terrestrial gravitational field.
* The use of the perigee of the proposed LARES satellite, for a suitable combination of observables including also the node and the perigee of LAGEOS II, would allow to reduce the systematic relative error due to the mismodeled even zonal harmonics of the geopotential to about 6$`\times 10^3`$, according to EGM96. This proposed measurement would also take advantage of the new Earth gravity models from CHAMP and GRACE. According to conservative estimates, the impact of the non–gravitational perturbations would amount to about $`4\times 10^3`$ over 7 years; in particular the impact of the mismodeled Earth thermal thrust, or Yarkovski-Rubincam effect, would amount to only about $`6\times 10^5`$.

## Acknowledgements

L. I. would like to thank L. Guerriero for his encouragement and support, C.M. Will and I.I. Shapiro for the useful information kindly supplied to him. E. C. Pavlis gratefully acknowledges NIMA support through the NURI grant NMA201-01-1-2008.


text:
# 1 Introduction

## 1 Introduction

As the hadrons with the largest gluon content (despite probably sizeable quarkonium admixtures), glueballs provide a privileged source of information on how glue manifests itself in hadronic bound states (in their rest frame). Progress in understanding glueball structure may therefore even help to unravel the elusive gluonic component of the classical hadrons.

The intricacies of hadronic glue are exemplified by the problems encountered when implementing various effective, gluonic degrees of freedom, such as constituent gluons, strings, or flux tubes, into hadron models. The resulting effects are generally less prominent and more ambiguous than those of the quark substructure, for which phenomenologically rather successful effective degrees of freedom are provided, e.g., by constituent quarks. The problems in modeling the gluon sector are partially related to the fact that gluons, unlike quarks, do not carry internal, global quantum numbers (like flavor or charge) to which external probes can couple. Unambiguous and transparent gluonic effects in hadrons are therefore much harder to identify.

Model-independent and more directly QCD-based approaches are thus called upon to clarify the structure of hadronic glue and to identify potentially dominant gluon field configurations. Promising candidates for the latter, especially in the $`0^{++}`$ glueball channel, are instantons  (i.e. the strong, coherent gluon fields which, by mediating tunneling processes, give rise to the $`\theta `$-angle of the QCD vacuum), as qualitative arguments  and instanton-liquid model results  suggest. In the following, we will give a brief summary of work in an analytical approach to the scalar glueball correlator at short distances , based on an instanton-improved operator product expansion (IOPE) and QCD sum-rule techniques, in which these issues can be addressed, and glueball properties calculated, in a largely model-independent fashion.

## 2 IOPE sum rules

Consider the correlation function

$$\mathrm{\Pi }\left(q^2\right)=id^4xe^{iqx}0|TO_S\left(x\right)O_S\left(0\right)|0$$
(1)
where the interpolating field

$$O_S=\alpha _sG_{\mu \nu }^aG^{a,\mu \nu }$$
(2)
carries the quantum numbers of the scalar ($`0^{++}`$) glueball. The standard OPE separates this correlator into contributions from hard field modes (with typical momenta $`k\mu 1/2`$ GeV), contained in perturbative Wilson coefficients, and from soft ($`k<\mu `$) modes contained in the vacuum expectation values of composite QCD operators, the so-called condensates. The perturbative Wilson coefficients of the low-dimensional operators can be found to $`O\left(\alpha _s\right)`$ in .

Remarkably, the nonperturbative power (i.e. condensates) corrections to this conventional OPE turn out to be small and are, except for the lowest-dimensional gluon condensate contribution, negligible. The bulk of the nonperturbative physics reponsible for the strong binding in the scalar glueball channel should therefore manifest itself in the Wilson coefficients, i.e. predominantly via direct instantons  which are neglected in the standard analyses . This expectation is strengthened by the fact that instantons couple particularly strongly to the gluonic $`0^{++}`$ interpolators.

In Ref. , we have evaluated the direct-instanton contributions and implemented them into the corresponding IOPE sum rules, which are based on the Borel-transformed moments

$$_k\left(\tau \right)=\widehat{B}\left[\left(Q^2\right)^k\mathrm{\Pi }\left(Q^2\right)\right]$$
(3)
(for the explicit form of the Borel operator $`\widehat{B}`$ see, e.g. ) with $`k\{1,0,1,2\}`$. Besides the standard OPE contributions $`_k^{\left(OPE\right)}`$ of Ref. , the IOPE includes the instanton contributions

$`_k^{\left(I+\overline{I}\right)}\left(\tau \right)`$ $`=`$ $`2^6\pi ^2\overline{n}\left({\displaystyle \frac{}{\tau }}\right)^{k+1}`$ (4)
$`\times \left\{\xi ^2e^\xi \left[\left(1+\xi \right)K_0\left(\xi \right)+\left(2+\xi +{\displaystyle \frac{2}{\xi }}\right)K_1\left(\xi \right)\right]2\right\},`$
($`k1`$, $`\xi \overline{\rho }^2/\left(2\tau \right)`$). The semiclassical result (4) neglects $`O\left(\mathrm{}\right)`$ corrections (which are suppressed by the large instanton action $`S_I\left(\overline{\rho }\right)10\mathrm{}`$) and multi-instanton correlations (since $`\left|x\right|\left|Q^1\right|0.2`$ fm is small compared with the average instanton separation $`\overline{R}1`$ fm). We have also removed the soft subtraction term $`\mathrm{\Pi }^{\left(I+\overline{I}\right)}\left(0\right)=2^7\pi ^2\overline{n}`$ to avoid double-counting with the condensates. Below, the average instanton size $`\overline{\rho }(1/3)`$ fm and density $`\overline{n}(1/2)`$ fm<sup>-4</sup> will be fixed at the values obtained (approximately) from instanton vacuum model  and lattice  simulations.

To obtain sum rules, the IOPE expressions are matched to their “phenomenological” counterparts, which are derived from the borelized dispersive representation

$$_k^{\left(ph\right)}\left(\tau \right)=\frac{1}{\pi }_0^{\mathrm{}}𝑑ss^k\mathrm{Im}\mathrm{\Pi }^{\left(ph\right)}\left(s\right)e^{s\tau }+\delta _{k,1}\mathrm{\Pi }^{\left(ph\right)}\left(0\right)$$
(5)
where the spectral function $`\mathrm{Im}\mathrm{\Pi }^{\left(ph\right)}\left(s\right)`$ contains a glueball pole contribution and an effective continuum from the dispersive cut of the IOPE, starting at an effective threshold $`s_0`$. This corresponds to

$$\mathrm{Im}\mathrm{\Pi }^{\left(ph\right)}\left(s\right)=\pi f_G^2m_G^4\delta \left(sm_G^2\right)+\mathrm{Im}\left[\mathrm{\Pi }^{\left(OPE\right)}+\mathrm{\Pi }^{\left(I+\overline{I}\right)}\right]\left(s\right)\theta \left(ss_0\right).$$
(6)
The instanton continuum contributions

$$\mathrm{Im}\mathrm{\Pi }^{\left(I+\overline{I}\right)}\left(s\right)=2^4\pi ^4\overline{n}\overline{\rho }^4s^2J_2\left(\sqrt{s}\overline{\rho }\right)Y_2\left(\sqrt{s}\overline{\rho }\right)$$
(7)
($`J_2`$ ($`Y_2`$) are Bessel (Neumann) functions), introduced in , will play an essential role in the subsequent analysis. The IOPE sum rules can then be written as

$$_k^{\left(IOPE\right)}\left(\tau \right)_k^{\left(I+\overline{I}\right)}\left(\tau \right)+_k^{\left(OPE\right)}\left(\tau \right)=_k^{\left(ph\right)}(\tau ,s_0).$$
(8)
The subtraction constant $`\mathrm{\Pi }^{\left(ph\right)}\left(0\right)`$ in the $`_1`$ sum rule (regularized by removing the high-momentum contributions) is related to the gluon condensate by the low-energy theorem

$$\mathrm{\Pi }\left(0\right)=\frac{32\pi }{b_0}\alpha G^2.$$
(9)
This relation provides an important consistency check for the sum-rule analysis, as discussed below.

## 3 Results and conclusions

The quantitative analysis of the sum rules (8) amounts to determining those values of the glueball parameters and $`s_0`$ in Eq. (6) for which both sides optimally match in the fiducial $`\tau `$ domain (determined such that the approximations on both sides of the sum rules are expected to be reliable ). The previously neglected instanton contributions turn out to be dominant and render (8) the first overall consistent set of QCD sum rules in the scalar glueball channel.

In particular, the IOPE resolves two long-standing flaws of the earlier sum rules : the mutual inconsistency between different Borel moments and the inconsistency with the low-energy theorem (9). Even the previously deficient and usually discarded lowest-moment ($`_1`$) sum rule is rendered consistent both with the higher-moment sum rules and the low-energy theorem. (A subsequent analysis of the related gaussian variant of these sum rules can be found in Ref. .) No evidence for a low-lying ($`m1`$ GeV) gluonium state (or any state strongly coupled to gluonic interpolators), which had been argued for on the basis of this sum rule , remains.

All four IOPE sum rules show an unprecedented degree of consistency and allow for a simultaneous 3-parameter fit to the glueball mass, its coupling, and the continuum threshold $`s_0`$. In Fig. 1 the glueball pole contribution (solid line) to the optimized $`_0`$ sum rule is compared with the remaining components (OPE with subtracted OPE-continuum (dash-double-dotted), instanton contribution (dashed), instanton continuum (dash-dotted), and their sum (dotted)). Both parts match almost perfectly over the whole fiducial region, mostly due to the dominant instanton contributions. The remaining sum rules (including the $`_1`$ sum rule) show a similarly high degree of stability. Together with the mutual agreement of all four IOPE sum rules (in the typical range of uncertainty) and their consistency with the low-energy theorem, this indicates that the IOPE provides a sufficiently complete description of the short-distance glueball correlator. The most dramatic quantitative effect of the direct instantons is to increase the residuum of the glueball pole, $`f_G^2`$, by about a factor of five. From the $`_2`$ sum rule (likely to be the most reliable one since it receives the strongest relative pole contribution) we obtain $`m_G=1.53\pm 0.2`$ GeV and $`f_G=1.01\pm 0.25`$ GeV.

The instanton continuum contributions, Eq. (7), are indispensable for the overall consistency and stability of the sum rules and shed new light on the spectral content of the glueball correlator. For once, they compensate, together with the perturbative terms, the pole contribution and thereby lead to an improved description of the correlator towards low momenta. This reconciles the sum rules with the low-energy theorem in the $`Q0`$ limit. In the opposite limit, i.e. for $`\tau 0`$, the instanton part of the continuum remains effective despite the suppression incurred by funnelling a sizeable momentum through a coherent field. This indicates that small-instanton physics accounts for part of the higher-lying strength in (1) and may play a rather prominent role in excited glueball states.

In contrast to previously studied IOPE sum rules , built on quark-based correlators, those for the scalar glueball are the first where i) the instanton contributions do not enter via topological quark zero-modes and where ii) the sum rules reach a satisfactory (though not excellent) level of consistency even without any perturbative and soft contributions, i.e. with the instanton terms alone. The latter result explains why instanton models find scalar glueball properties similar to those obtained above . It also illustrates that the quantitative impact of the instanton contributions can be judged only if all contributions are consistently taken into account. At present, the IOPE seems to be the only controlled and analytical framework in which this is possible.

The predominance and approximate self-sufficiency of the instanton contributions indicates that instantons may generate the bulk of the forces which bind the scalar glueball. A further, striking consequence of the exceptionally strong instanton contributions is that the scales of the predicted $`0^{++}`$ glueball properties are approximately set by the bulk features of the instanton size distribution. More specifically, neglecting the standard OPE contributions, one finds

$`m_G`$ $``$ $`\overline{\rho }^1,`$ (10)
$`f_G^2`$ $``$ $`\overline{n}\overline{\rho }^2.`$ (11)
Conceptually, the main virtue of these scaling relations lies in establishing an explicit link between fundamental vacuum and hadron properties. They could also be of practical use, e.g. for the test of instanton vacuum models, to provide constraints for glueball model building, or generalized to finite temperature and baryon density, where the instanton distribution changes. If the glueball size $`r_G`$ scales like its Compton wavelength, one furthermore has $`r_G\overline{\rho }`$, in agreement with the lattice calculations  which find $`r_G0.2`$ fm.

It is a pleasure to thank the organizers for this interesting workshop and the Deutsche Forschungsgemeinschaft for financial support under Habilitation Grant Fo 156/2-1.