text:
# Lattice Gauge Fixing and the Violation of Spectral Positivity

## Abstract

Spectral positivity is known to be violated by some forms of lattice gauge fixing. The most notable example is lattice Landau gauge, where the effective gluon mass is observed to rise rather than fall with increasing distance. We trace this violation to the use of quenched auxiliary fields in the lattice gauge fixing process, and show that violation of spectral positivity is a general feature of quenching. We illustrate this with a simple quenched mass-mixing model in continuum field theory, and with a quenched form of the Ising model. For lattice gauge fixing associated with Abelian projection and lattice Landau gauge, we show that spectral positivity is violated by processes similar to those found in quenched QCD. For covariant gauges parametrized by a gauge-fixing parameter $`\alpha `$, the $`SU(2)`$ gluon propagator is well described by a simple quenched mass-mixing formula. The gluon mass parameter appears to be independent of $`\alpha `$ for sufficiently large $`\alpha `$.

Although many observables can be determined in lattice gauge theories without gauge fixing, there are several reasons why gauge fixing is desirable in lattice simulations. Gauge fixing is necessary to make the connection between continuum and lattice gauge fields. Continuum theories of the origin of confinement often make predictions about the gauge field propagator. Gauge fixing has also been a key technique in lattice studies of confinement as well.Greensite:2003bk  Important properties of the quark-gluon plasma phase of QCD, such as screening masses, are contained in the finite-temperature gluon propagator.

Techniques for lattice gauge fixing have been known for some time.Mandula  It has been clear from the beginning that non-Abelian lattice gauge field propagators show a violation of spectral positivity. This is readily seen from the effective mass: for a normal operator which connects only states of positive norm to the vacuum, the effective mass monotonically decreases with distance to the lightest mass state coupling to the operator. Covariant gauge gluon propagators have an effective mass increasing with distance. In one sense, this is not surprising. We know from perturbation theory that covariant gauges contain states of negative norm. However, that knowledge has neither explained the form of the lattice gluon propagator nor aided in the interpretation of the mass parameters measured from it. In fact, no similar violation of spectral positivity is observed in the $`U(1)`$ case Coddington:yz , which has negative-norm states in covariant gauges.

In lattice simulations, gauge fixing has typically involved choosing a particular configuration on each gauge orbit. A brief review of this approach is given in Ref. Mandula:nj . In the continuum, on the other hand, gauge fixing usually includes a parameter that causes the functional integral to peak around a particular configuration on the gauge orbit. As shown below, the extension of this idea to lattice gauge theories makes clear that lattice gauge fixing is a form of quenching, with the gauge transformations acting as quenched fields. As has been demonstrated in quenched QCD, quenching can violate spectral positivity, with significant effects on many observables.Golterman:1994jr ; Bardeen:2001yz

We begin with a review of lattice gauge fixing, including the generalization of lattice Landau gauge to covariant gauges with a gauge parameter.Zwanziger:tn ; Parrinello:1990pm ; Fachin:1991pu  This generalization will be directly interpreted as a quenched Higgs theory. We then explore the origin of violations of spectral positivity in some simple lattice and continuum models of quenching. Simulation results for the effective mass of an $`SU(2)`$ lattice gauge field will show behavior very similar to these models as the gauge fixing parameter is varied. We will argue that spectral positivity violations in both lattice covariant gauges and in studies of Abelian projection originate in the quenching process.

The standard approach to lattice gauge fixing is a two step process.Mandula:nj  An ensemble of lattice gauge field configurations is generated using standard Monte Carlo methods, corresponding to a functional integral

$$Z_U=\left[dU\right]e^{S_U\left[U\right]},$$
(1)
where $`S_U`$ is a gauge-invariant action for the gauge fields, e.g., the Wilson action. The gauge action $`S_U`$ is invariant under gauge transformations of the form $`U_\mu \left(x\right)g\left(x\right)U_\mu \left(x\right)g^+\left(x+\mu \right)`$.

In order to measure gauge-variant observables, each field configuration in the $`U`$-ensemble may be placed in a particular gauge, i.e., a gauge transformation is applied to each configuration in the $`U`$-ensemble which moves the configuration along the gauge orbit to a gauge-equivalent configuration satisfying a lattice gauge fixing condition. The simplest gauge choice is defined by maximizing $`_{x,\mu }Tr\left[U_\mu \left(x\right)+U_\mu ^+\left(x\right)\right]`$ for each configuration over the class of all gauge transformations. Any local extremum of this functional satisfies a lattice form of the Landau gauge condition:

$$\underset{\mu }{}\left[A_\mu \left(x+\mu \right)A_\mu \left(x\right)\right]=0$$
(2)
where $`A_\mu \left(x\right)`$ is a lattice approximation to the continuum gauge field, given by

$$A_\mu \left(x\right)=\frac{U_\mu \left(x\right)U_\mu ^+\left(x\right)}{2i}\frac{1}{N}Tr\left[\frac{U_\mu \left(x\right)U_\mu ^+\left(x\right)}{2i}\right]\text{.}$$
(3)
Other gauge-fixing conditions may also be used Giusti:2001xf , and lattice improvement techniques can be applied to the definition of $`A_\mu `$ to reduce discretization errors as well. The global maximization needed is often implemented as a local iterative maximization. The issue of Gribov copies arises in lattice gauge fixing because such a local algorithm tends to find local maxima of the gauge-fixing functional. There are variations on the basic algorithm that ensure a unique choice from among local maxima.Giusti:2001xf

For analytical purposes, it is necessary to generalize this procedure Fachin:1991pu , so that a given single configuration of gauge fields will be associated with an ensemble of configurations of $`g`$-fields. We will generate this ensemble using

$$S_{gf}[U,g]=\underset{l}{}\frac{\alpha }{2N}Tr\left[g\left(x\right)U_\mu \left(x\right)g^+\left(x+\mu \right)+g\left(x+\mu \right)U_\mu ^+\left(x\right)g^+\left(x\right)\right]$$
(4)
as a weight function to select an ensemble of $`g`$-fields. The sum over $`l`$ is a sum over all links of the lattice. The normal gauge-fixing procedure is formally regained in the limit $`\alpha \mathrm{}`$. Computationally, this can be implemented as a Monte Carlo simulation inside a Monte Carlo simulation.

Note that the $`g`$-fields must be thought of as quenched variables, since they do not affect the $`U`$-ensemble. The expectation value of an observable $`O`$, gauge-invariant or not, is given by

$$O=\frac{1}{Z_U}\left[dU\right]e^{S_U\left[U\right]}\frac{1}{Z_{gf}[U]}\left[dg\right]e^{S_{gf}[U,g]}O,$$
(5)
where

$$Z_{gf}[U]=\left[dg\right]e^{S_{gf}[U,g]}\text{.}$$
(6)
Formally, the field $`g`$ is a quenched scalar field with two independent symmetry groups, $`G_{global}G_{local}`$, so that it appears to be in the adjoint representation of the gauge group, but the left and right symmetries are distinct. The generating functional $`Z_{gf}[U]`$ is in some ways a lattice analog of the inverse of the Fadeev-Popov determinant.Bock:2000cd  However, there are important differences. Note immediately that $`Z_{gf}[U]`$ depends on the gauge-fixing parameter $`\alpha `$. More fundamentally, the lattice formalism resolves the Gribov ambiguity. By construction, gauge-invariant observables are evaluated by integrating over all configurations. Gauge-variant quantities receive contributions from Gribov copies, always with positive weight. Thus the connection between this formalism for lattice gauge fixing and gauge fixing in the continuum is not simple. Furthermore, alternative lattice gauge fixing procedures have been proposed, along with new gauge choices specific to the lattice. A comprehensive review of gauge fixing technology is available.Giusti:2001xf

We begin our analysis of spectral positivity violation with the simplest model of quenching possible: two free, real scalar fields with a non-diagonal mass matrix. The Lagrangrian is

$$L=\frac{1}{2}\left[\left(\varphi _1\right)^2+m_1^2\varphi _1^2\right]+\frac{1}{2}\left[\left(\varphi _2\right)^2+m_2^2\varphi _2^2\right]\mu ^2\varphi _1\varphi _2.$$
(7)
We treat the quenched approximation of this model in a manner completely parallel to our discussion of lattice gauge fixing above. We divide the action into three parts $`S=S_1+S_2+S_{12}`$, where $`S_1`$ and $`S_2`$ are functionals only of $`\varphi _1`$ and $`\varphi _2`$, respectively, and $`S_{12}`$ contains the mixing term. We quench the field $`\varphi _2`$. Although there are no loops in this simple theory, quenching implies that $`\varphi _2`$ cannot appear as an internal line in the complete propagators. The generating functional in the quenched approximation, including sources $`J_1`$ and $`J_2`$ is

$$Z=\left[d\varphi _1\right]e^{S_1+{\scriptscriptstyle J_1\varphi _1}}\frac{\left[d\varphi _2\right]e^{S_2S_{12}+{\scriptscriptstyle J_2\varphi _2}}}{\left[d\stackrel{~}{\varphi }_2\right]e^{S_2S_{12}}}\text{.}$$
(8)
where we have introduced a kind of ghost variable $`\stackrel{~}{\varphi }_2`$; space-time variables are implicit.

From the generating functional we can obtain the $`\varphi _1\varphi _1`$ and $`\varphi _2\varphi _2`$ propagators. In momentum space, the $`\varphi _1\varphi _1`$propagator is $`1/(p^2+m_1^2)`$, since $`\varphi _1`$ is unaffected by $`\varphi _2`$. On the other hand, the $`\varphi _2\varphi _2`$ propagator is

$$\frac{1}{p^2+m_2^2}+\frac{1}{p^2+m_2^2}\mu ^2\frac{1}{p^2+m_1^2}\mu ^2\frac{1}{p^2+m_2^2}\text{.}$$
(9)
An alternative diagrammatic procedure is to sum Dyson’s series, as shown in Fig. 1, noting that the $`\varphi _2\varphi _2`$ propagator is truncated at two terms. The propagator has a structure similar to the $`\eta ^{}`$ propagator in quenched QCD Golterman:1994jr ; Bardeen:2001yz ; the $`\eta ^{}`$ has a double pole form in quenched QCD when singlet self-energy graphs are approximated by a constant. The $`\varphi _2\varphi _2`$ propagator also may be written as

$$\left(1\frac{\mu ^4}{\left(m_2^2m_1^2\right)^2}\right)\frac{1}{p^2+m_2^2}+\frac{\mu ^4}{\left(m_2^2m_1^2\right)^2}\frac{1}{p^2+m_1^2}+\frac{\mu ^4}{m_1^2m_2^2}\frac{1}{\left(p^2+m_2^2\right)^2}\text{.}$$
(10)
This propagator always violates spectral positivity because of the double pole term, $`1/\left(p^2+m_2^2\right)^2`$, which has a coefficient whose sign depends on $`m_1^2m_2^2`$. Another possible violation of spectral positivity occurs for sufficiently strong mixing: if $`\mu ^4>\left(m_2^2m_1^2\right)^2`$, there is a simple pole at $`p^2=m_2^2`$ with negative residue.

The form of the $`\varphi _2\varphi _2`$ propagator in coordinate space is very interesting, and forms the basis for our study of other quenched theories. In any number of dimensions, we can consider propagators using wall sources, i.e., of co-dimension $`1`$. This has the effect of setting the momentum equal to zero in all the directions of the wall. For wall sources, we have the propagator

$`G(x)`$ $`=`$ $`\left(1{\displaystyle \frac{\mu ^4}{\left(m_2^2m_1^2\right)^2}}\right){\displaystyle \frac{1}{2m_2}}e^{m_2\left|x\right|}+{\displaystyle \frac{\mu ^4}{\left(m_2^2m_1^2\right)^2}}{\displaystyle \frac{1}{2m_1}}e^{m_1\left|x\right|}`$ (11)
$`+{\displaystyle \frac{\mu ^4}{m_1^2m_2^2}}{\displaystyle \frac{1}{4m_2^3}}e^{m_2\left|x\right|}\left(1+m_2\left|x\right|\right)\text{.}`$
The factor $`m_2\left|x\right|e^{m_2\left|x\right|}`$ shows an initial rise rather than a decay with increasing $`\left|x\right|`$, violating spectral positivity.

We define an effective mass associated with the $`\varphi _2`$ field as

$$m_{eff}(x)=\underset{a0}{lim}\frac{1}{a}\mathrm{ln}\left(\frac{G(x)}{G(x+a)}\right)=\frac{d}{dx}\mathrm{ln}(G(x))\text{.}$$
(12)
One can easily check explicitly that $`m_{eff}(x)\mathrm{min}(m_1,m_2)`$ as $`x\mathrm{}`$. For any field theory which obeys spectral positivity, $`m_{eff}(x)`$ monotonically decreases to its limiting value. Theories violating spectral positivity may display a complicated behavior for $`m_{eff}(x)`$ before the eventual onset of asymptotic behavior.

We have identified three different possible behaviors for $`m_{eff}(x)`$ in this simple quenched model. If the mixing parameter $`\mu `$ is sufficiently small and $`m_1<m_2`$, $`m_{eff}(x)`$ monotonically decreases to its value at infinity, as in a normal field theory which obeys spectral positivity, as shown in Fig. 2(a). As $`\mu `$ is increased relative to $`m_1`$ and $`m_2`$, $`m_{eff}(x)`$ may develop a minimum, as displayed in Fig. 2(b). On the other hand, if $`m_2<m_1`$, the behavior seen in Fig. 2(a) is not possible, and only the behaviors seen in Fig. 2(c) and Fig. 2(d) are possible. In Fig. 2(d), the minimum has moved to $`x=0`$. For sufficiently small $`\mu `$, these effects are difficult to observe, and $`m_{eff}(x)`$ is essentially equal to $`m_2`$ for all $`x`$. Regardless of the relative size of $`m_1`$ and $`m_2`$, an observable violation of spectral positivity associated with $`m_{eff}(x)`$ not monotonically decreasing indicates a significant mixing parameter $`\mu `$.

Similar behavior can be observed in a very simple lattice model based on the Ising model, where real-space arguments can be used to find an approximate propagator. We consider two coupled one-dimensional Ising models, with spins $`\mu _i`$, $`\sigma _i\{1,+1\}`$ and respective nearest-neighbor couplings $`J`$ and $`K`$. The $`\sigma `$ spins are coupled to the $`\mu `$ spins via an interaction of the form $`_iL\sigma _i\mu _i`$, and the $`\sigma `$’s are quenched. This simple model is a form of spin glass, with the averaging over the ensemble of $`\mu `$ spins representing the “quenching” process.

The $`\sigma `$propagator is given by

$$\sigma _0\sigma _n=\frac{1}{Z_\mu }\underset{\left\{\mu \right\}}{}\mathrm{exp}\left[\underset{i}{}J\mu _i\mu _{i+1}\right]\frac{1}{Z_\sigma \left[\mu \right]}\underset{\left\{\sigma \right\}}{}\sigma _0\sigma _n\mathrm{exp}\left[\underset{i}{}\left(K\sigma _i\sigma _{i+1}+L\sigma _i\mu _i\right)\right],$$
(13)
where $`Z_\mu `$ is the partition function for $`\mu `$ and $`Z_\sigma [\mu ]`$ is the partition function for $`\sigma `$ in the presence of a particular $`\mu `$ background. The parameter $`L`$is a mixing parameter. We can approximately evaluate the $`\sigma `$propagator for $`J`$, $`K`$, and $`L`$ sufficiently small by considering the direct contribution $`\left(\mathrm{tanh}K\right)^n`$ combined with mixing of $`\sigma `$with $`\mu `$. This indirect term can be written as (compare Fig. 1)

$$\underset{p=1}{\overset{n}{}}\underset{m=0}{\overset{np}{}}\left(\mathrm{tanh}K\right)^{np}\left(\mathrm{tanh}L\right)^2\left(\mathrm{tanh}J\right)^p.$$
(14)
After performing the summations, the propagator is given approximately as

$$\sigma _0\sigma _n\left(\mathrm{tanh}K\right)^n+\left(\mathrm{tanh}L\right)^2\left(\mathrm{tanh}K\right)^n\left[\frac{nx}{\left(1x\right)}\frac{x^2(1x^n)}{\left(1x\right)^2}\right].$$
(15)
where $`x=`$ $`\mathrm{tanh}J/\mathrm{tanh}K`$. The $`n(\mathrm{tanh}K)^n`$ factor signals a violation of spectral positivity, just as the $`m\left|x\right|\mathrm{exp}\left(m\left|x\right|\right)`$ term did in the mixing model. Of course, the arguments which led to Eq. (9) and Eq. (15) are essentially the same, but carried out in momentum space and real space, respectively. For small $`J`$, $`K`$, and $`L`$, Eq. (15) fits lattice simulations of the $`\sigma _0\sigma _n`$ propagator well.

In Fig. 3, we show the effective mass determined from the $`\sigma _0\sigma _n`$ and $`\mu _0\mu _n`$ propagators for the parameter set $`J=0.7`$, $`K=0.9`$, and $`L=0.3`$ for a one-dimensional lattice of size $`26`$. The propagators were obtained from $`40000`$ heat bath sweeps of the $`\mu `$ variables; after each such sweep, $`100`$ heat bath sweeps of the $`\sigma `$ variables were carried out. The parameters $`K`$ and $`L`$ were chosen empirically so as to display a clear violation of spectral positivity. The $`\mu `$ mass fits very well with the analytical solution $`m=\mathrm{ln}\mathrm{tanh}(J)`$ for the $`d=1`$ Ising model out to a distance of $`8`$. Note that the $`\sigma `$ reaches its asymptotic value of $`\mathrm{ln}(\mathrm{tanh}K)`$ from below, and only at $`n8`$. The similarity to the simple field theoretic model of quenching is clear.

We will now show that the $`SU(2)`$ lattice gluon propagator regarded as a function of $`\alpha `$ shows behavior similar to that of the other, simpler quenched models studied above. Simulations of this type of lattice field theory, with stochastic quenched gauge fixing fields, were first performed by Henty et al. Henty:1996kv , who studied the case of $`SU(3)`$ as a function of $`\alpha `$ at $`\beta =5.7`$ on $`8^4`$ lattices. They found evidence for a first-order phase transition as $`\alpha `$ was varied, but did not determine the full phase diagram in the $`\alpha `$-$`\beta `$ plane. They also found that the gluon propagator was dependent on $`\alpha `$, a result which could be anticipated from the strong-coupling expansion.Fachin:1991pu

Let us consider for the moment the unquenched version of the gauge fixing model. This is a model with scalar fields in the fundamental representation of the gauge group in addition to the gauge fields. The scalar fields explicitly break the $`Z(N)`$ global symmetry associated with confinement in the pure gauge case, and external color charges are screened. As first shown by Fradkin and Shenker Fradkin:1978dv , this leads to a connection between the strong-coupling, confining phase and the Higgs phase, so the two phases are not actually distinct. We have verified that this phase structure is preserved in the quenched form of the model. For $`\beta `$ sufficiently large, there is a line of first-order phase transitions in the $`\beta `$-$`\alpha `$ plane. It is very reasonable that such a line exists in the quenched model, since it can be thought of as the continuation of the critical point of a pure spin model at $`\beta =\mathrm{}`$. However, this line terminates at a critical end point; for sufficiently small $`\alpha `$, the nominal confining phase ($`\beta `$ small) and Higgs phase ($`\beta `$ large) are directly connected. This observation forms the starting point for a detailed analysis of the model.us\_longer\_ms

We have performed simulations of $`SU(2)`$gauge theory at $`\beta =2.6`$ with $`\alpha `$ ranging from $`1.0`$ to $`3.0`$ on a $`12^3\times 16`$ lattice. At this value of $`\beta `$, there is a first-order phase transition at $`\alpha 0.83`$. We have fit the data using a simple generalization of the quenched mixing model. The coordinate space propagator has the form

$$\left(A+Bm_2\left|x\right|\right)e^{m_2\left|x\right|}+Ce^{m_1\left|x\right|}\text{.}$$
(16)
This form for the propagator follows from the replacement of the mixing parameter $`\mu ^4`$ in Eq. (9) by the more general form $`\mu ^2(p^2+m_3^2)`$.

In Fig. 4, we show the effective mass as a function of $`x`$for $`\alpha =1`$, $`2`$, and $`3`$. The solid lines are obtained from fits to the propagator using Eq. (16). The similarity to the other quenched models is quite clear. For small $`\alpha `$ there is an initial decrease and then rise of the effective mass, much like Fig. 2(c); as $`\alpha `$ increases, this minimum vanishes. We plot in Fig. 5 the best-fit values of $`m_1`$ and $`m_2`$ as a function of $`\alpha `$. Note that the behavior of $`m_2`$ is consistent with it being constant in this region, while $`m_1`$ appears to be decreasing to a constant limit as $`\alpha `$ increases.

Our results are roughly consistent with the work of Leinweber et al., who performed high-precision studies of the $`SU(3)`$ gluon propagator at $`\alpha =\mathrm{}`$.Leinweber:1998uu  Among a large variety of possible functional forms for the gluon propagator, they found that their data was best fit by the functional form

$$G(k)=Z\left[\frac{AM^{2\delta }}{(k^2+M^2)^{1+\delta }}+\frac{1}{k^2+M^2}L(k^2,M)\right],$$
(17)
with $`L(k^2,M)`$ an infrared-regulated version of the asymptotic behavior of the renormalized gluon propagator in the continuum. Their best fit was achieved with the parameters $`\delta =2.2_{0.20.3}^{+0.1+0.2}`$, $`M=(1020\pm 100\pm 25)\text{MeV}`$, and $`A=9.8_{0.9}^{+0.1}`$. Many other functional forms were ruled out.

Our results suggest that the lighter mass parameter $`m_2`$ is independent of $`\alpha `$, at least for large $`\alpha `$ (in the Higgs phase). If $`m_2`$ is indeed independent of gauge choice, as least within the class of covariant gauges considered, it seems natural to identify it as the gluon mass. As a consequence of the quenched character of lattice gauge fixing, this state partially mixes with another, heavier state, with a mass on the order of the scalar or vector glueball.Teper:1998kw

Note that the value of the lightest mass in the propagator may be difficult to extract from the effective mass. While it is true that the effective mass tends asymptotically to the lightest mass, the approach to the limit can be much slower than in a conventional field theory obeying spectral positivity. For example, at $`\alpha =3.0`$, $`m_{eff}`$ at $`x=7`$ is substantially lighter than $`m_2`$. Having a theoretical basis for the form of the propagator is crucial in estimating the mass.

Another application of lattice gauge fixing is Abelian projection, a method for investigating the confining properties of gauge theories. In lattice gauge theories, Abelian projection is implemented as an algorithm for extracting an ensemble of Abelian gauge field configurations from an ensemble of non-Abelian configurations. A notable success of lattice studies of Abelian projection Suzuki:1989gp ; Hioki:1991ai  has been the correlation of the string tension of the projected theory with the string tension of the underlying non-Abelian theory.

These lattice studies of Abelian projection typically use gauge fixing in an integral way. Taking for clarity the case of $`SU(2)`$, the gauge fixing functional is

$$S_{gf}=\underset{x,\mu }{}\frac{\alpha }{2}Tr\left[g\left(x\right)U_\mu \left(x\right)g^+\left(x+\mu \right)\sigma _3g\left(x+\mu \right)U_\mu ^+\left(x\right)g^+\left(x\right)\sigma _3\right],$$
(18)
which is conventionally maximized over the gauge orbit, corresponding to the limit $`\alpha \mathrm{}`$ in the formalism used here. An initial study of the phase structure in the $`SU(2)`$ case finds evidence for a first order phase transition as $`\alpha `$ is varied at $`\beta =2.4`$.Mitrjushkin:2001hr  The aim of this procedure is to transform an $`SU(2)`$ configuration into a gauge-equivalent one which lies mostly in a given $`U(1)`$ subgroup. After this gauge-fixing, the actual projection to $`U(1)`$ is performed.

In the case where no gauge fixing is done ($`\alpha =0`$), and only projection occurs, Faber et al. Faber:1998en  and Ogilvie Ogilvie:1998wu  have proved that the asymptotic string tension measured in the projected and underlying theories are the same. Furthermore, Ogilvie Ogilvie:1998wu  has proven that this result should continue to hold for small $`\alpha `$, under the assumption that the gauge fixing does not violate spectral positivity. However, the fact that the string tension evaluated using various forms of Abelian projection with gauge fixing is consistently slightly different from the actual non-Abelian string tension Stack:2000kf ; Bornyakov:2000cd  suggests that a violation of spectral positivity may indeed be occurring.

We identify the origin of this violation as the presence of a quenched scalar field. The case of $`SU(2)`$ is particularly clear. Note that the combination $`g^+\left(x\right)\sigma _3g\left(x\right)`$ occurring in $`S_{gf}`$ can be written as a Hermitian scalar field $`\varphi (x)`$, where $`\varphi `$ transforms as the adjoint representation of the gauge group. The field $`\varphi `$ is traceless, $`Tr(\varphi )=0`$, and satisfies $`Tr(\varphi ^2)=2`$. The gauge fixing action is thus equivalent to an adjoint scalar action of the form

$$S_{gf}=\underset{x,\mu }{}\frac{\alpha }{2}Tr\left[\varphi \left(x\right)U_\mu \left(x\right)\varphi \left(x+\mu \right)U_\mu ^+\left(x\right)\right]\text{.}$$
(19)
As we have seen, such quenched fields naturally lead to violations of spectral positivity. Suppose we wish to measure a $`U(1)`$ projected Wilson loop. This may be obtained from the expectation value of

$$Tr\underset{j}{}\frac{1}{2}(1+\sigma _3)g_jU_jg_{j+1}^+=Tr\underset{j}{}\frac{1}{2}(1+\varphi _j)U_j,$$
(20)
where the product is ordered along a closed path labeled by the index $`j`$. The $`U(1)`$ projected loop is represented in the underlying quenched Higgs theory as a sum of Wilson loops with all possible insertions of $`\varphi `$ at lattice sites on the path. When four or more $`\varphi `$ fields are inserted, problematic subdiagrams appear of the type shown in Fig. 6. Such terms lead to a violation of spectral positivity: there are no internal $`\varphi `$ loops in the quenched approximation, and an infinite set of diagrams occurring in the full, unquenched theory is omitted. This exactly parallels quenched QCD.

###### Acknowledgements.

This work was partially supported by the U.S. Department of Energy under grant number DE-FG02-91ER40628.


format: text

text:
# The first deployment of workload management services on the EU DataGrid Testbed: feedback on design and implementation.

## I Introduction

The workload management task (Work Package 1, or WP1) wp1  of the EU DataGrid project DataGrid  (also known, and referred to in the following text, as EDG) is mandated to define and implement a suitable architecture for distributed scheduling and resource management in the Grid environment. During the first year and a half of the project (2001-2002), and following a technology evaluation process, EDG WP1 defined, implemented and deployed a set of services that integrate existing components, mostly from the Condor Condor  and Globus Globus  projects. This was described in more detail at CHEP 2001 WP1Chep01 . In a nutshell, the core job submission component of CondorG (CondorG ), talking to computing resources (known in DataGrid as Computing Elements, or CEs) via the Globus GRAM protocol, is fundamentally complemented by:

* A job requirement matchmaking engine (called the Resource Broker, or RB), matching job requests to computing resource status coming from the Information System and resolving data requirements against the replicated file management services provided by EDG WP2.
* A job Logging and Book-keeping service (LB), where a job state machine is kept current based on events generated during the job lifetime, and the job status is made available to the submitting user. The LB events are generated with some redundancy to cover various cases of loss.
* A stable user API (command line, C++ and JAVA) for access to the system.

Job descriptions are expressed throughout the system using the Condor Classified Ad language, where appropriate conventions were established to express requirement and ranking conditions on Computing and Storage Element info, and to express data requirements. More details on the structure and evolution of these services and the necessary integration scaffolding can be found in various EDG public deliverable documents.

This paper focuses on how the experience of the first year of operation of the WP1 services on the EDG testbed was interpreted, digested, and how a few design principles were learned (possibly the hard way) from the design and implementation shortcomings of the first release of WP1 software.

These principles were applied to design and implement the second major release of WP1 software, that is described in another CHEP 2003 paper (WP1Chep03Sga ).

To illustrate the logical path that leads to at least some of these principles, we start by exploring the available techniques to model the behaviour and throughput of the integrated workload management system, and identify two factors that significantly complicate the system analysis.

## II The Workload Management System as a network of queues

The workload management system provided by EDG-WP1 is designed to rely as much as possible on existing technology. While this has the obvious advantages of limiting effort duplication and facilitating the compatibility among different projects, it also significantly complicates troubleshooting across the various layers of software supplied by different providers, and in general the understanding of the integrated system. Also, where negotiations with external software providers couldn’t reach an agreement within the EDG deadlines, some of the interfaces and communication paths in the system had to be adapted to fit the existing external software incarnations.

To get a useful high-level picture of the integrated Workload Management system, beyond all these practical constraints, we can model it as a queuing system, where job requests traverse a network of queues, and the “service stations” connected to each queue represent one of the various processing steps in the job life-cycle. A few of these steps are exemplified in Figure 1.

Establishing the scale factors for each service in the WP1 system (e.g.: how many users can a single matchmaking/job submission station serve, how many requests per unit time can a top-level access point to the information system serve, what is the sustained job throughput that can be achieved through the workload management chain, etc.) is one of the fundamental premises for the correct design of the system. One could expect to obtain this knowledge either by applying queuing theory to this network model (this requires obtaining a formal representation of all the components, their service time profiles and their interconnections) or by measuring the service times and by identifying where long queues are likely to build up when a “realistic” request load is injected in the system. This information could in principle also be used to identify the areas of the system where improvement is needed (sometimes collectively called bottlenecks).

Experience with the WP1 software integration showed that both of these approaches are impractical for either dimensioning the system or (possibly even more important) for identifying the trouble areas that affect the system throughput. We identified two non-linear factors that definitely work against the predictive power of queuing theory in this case, and require extra care even to apply straightforward reasoning when bottlenecks are to be identified to improve system throughput. These are the consequence of common programming practice (and are therefore easy to be found in the software components that we build or are integrating) and are described in the following Section.

## III Troubleshooting the WMS

One of the most common (and most frustrating, both to developers and to end users) experiences in troubleshooting the WP1 Workload Management system on the EDG testbed has been the fact that often, perceived improvements to the system (sometimes even simple bug fixes) result in a decrease in the system stability, or reliability (fraction of requests that complete successfully). The cause is often closely related to the known fact that removing a bottleneck, in any flow system, can cause an overflow downstream, possibly close to the next bottleneck. The complicating factor is that there are at least two characteristics that could (and possibly still can) be found in many elements of our integrated workload management queuing network, that can cause problems to appear even very far from the area of the network where an improvement is being attempted:

* Queues of job requested can form where they can impact on the system load.
Different techniques can be chosen or needed to pass job requests around. Sometimes a socket connection is needed, sometimes sequential request processing (one request at a time in the system) is required for some reason, and multiple processes/threads may be used to handle individual requests. Having a number of tasks (processes/threads) wait for a socket queue or a sequential processing slot is one way to “queue” requests that definitely generates much extra work for the process scheduler, and can cause any other process served by the same scheduler to be allocated less and less time. Queues that are unnecessarily scanned while waiting for some other condition to allow the processing of their element can also impact on the system load, especially if the queue elements are associated to significant amounts of allocated dynamic memory.
* Some system components can enforce hard timeouts and cause anomalies in the job flow.
When handling the access (typically via socket connections) to various distributed services, provisions typically need to be made to handle all possible failure modes. “Reasonably” long timeouts are sometimes chosen to handle failures that are perceived to be very unlikely by developers (failure to establish communication to a local service, for instance). This kind of failures, however, can easily materialise when the system resources are exhausted under a stress test or load peak.

Figure 2 illustrates how these two effects can conspire to frustrate a genuine effort to remove what seems the limiting bottleneck in the system (the example in the Figure does nor refer to any real case or component): removing the bottleneck (1) causes a request queue to build up at the next station (2), and this interferes via the system load to cause hard timeouts and job failures elsewhere (3). This example is used to rationalise some of the unexpected reactions that, in many cases, were found while working on the WP1 integrated system. The experience on practical troubleshooting cases similar to this one, while bringing an understanding of the difficulties inherent in building distributed systems, also drove us to formulate some of the principles that are presented in the next section.

## IV Principles that were learnt (and applied to improve the design)

The attempts at getting a deeper understanding of the EDG-WP1 Workload Management System and their failures led us to formulate a few design principles and to apply them to the second major software release. Here are the principles that descend from the paradigm example described in Section III:

1. Queues of various kinds of requests for processing should be allowed to form where they have a minimal and understood impact on system resources.
Queues that get ‘filled’ in the form of multiple threads or processes, or that allocate significant amounts of system memory should be avoided, as they not only adversely impact system performance, but also generate inter-dependencies and complicate troubleshooting.
2. Limits should always be placed on dynamically allocated objects, threads and/or subprocesses.
This is a consequence of the previous point: every dynamic resource that gets allocated should have a tunable system-wide limit that gets enforced.
3. Special care needs to be taken around the pipeline areas where serial handling of requests is needed.
The impact of any contention for system resources becomes more evident near areas of the queuing system that require the acquisition of system-wide locks.

So far we concentrated on a specific attempt at modeling and understanding the workload management system that led to an increased attention to the usage of shared resources. There were other specific practical issues that emerged during the deployment and troubleshooting of the system and that led to the awareness of some fundamental design or implementation mistake that was made. Here is a short list, where the fundamental principle that should correct the fundamental mistake that was made is listed:

1. Communication among services should always be reliable:

* Always applying double-commit and rollback for network communications.
    * Going through the filesystem for local communications.

In general, forms of communication that don’t allow for data or messages to be lost in a broken pipe lead to easier recovery from system or process crashes. Where network communication is necessary, database-like techniques have to be used.
2. Every process, object or entity related to the job lifecycle should have another process, object or entity in charge of its well-being.
Automatic fault recovery can only happen if every entity is held accountable and accounted for.
3. Information repositories should be minimized (with a clear identification of authoritative information).
Many of the software components that were integrated in the EDG-WP1 solution are stateful and include local repositories for request information, in the form of local queues, state files, database back-ends. Only one site with authoritative information about requests has to be identified and kept.
4. Monolithic, long-lived processes should be avoided.
Dynamic memory programming, using languages and techniques that require explicit release of dynamically allocated objects, can lead to leaks of memory, descriptors and other resources. Experimental, R&D code can take time to leak-proof, so it should possibly not be linked to system components that are long-lived, as it can accelerate system resource starvation. Short-lived, easy-to-recover components are a clean and very practical workaround in this case.
5. More thought should be devoted to efficiently and correctly recovering a service rather than to starting and running it.
This is again a consequence of the previous point: the capability to quickly recover from failures or interruption helps in assuring that system components ‘can’ be short-lived, either by design or by accident.

## V Conclusions

EDG-WP1 has been distributing jobs over the EDG testbed in a continuous fashion for one and a half years now, with a software solution where existing grid technology was integrated wherever possible.

The experience of understanding the direct and indirect interplay of the service components could not be reduced to a simple scalability evaluation. This because understanding and removing bottlenecks is significantly complicated by non-linear and non-continuous effects in the system. In this process, few principles that apply to the very complex practice of distributed systems operations were learned the hard way (i.e. not by just reading some good book on the subject). EDG-WP1 tried to incorporate these principles in its second major software release that will shortly face deployment in the EDG testbed.

###### Acknowledgements.

DataGrid is a project funded by the European Commission under contract IST-2000-25182. We also acknowledge the national funding agencies participating to DataGrid for their support of this work. We wish to thank the Condor development team for many valuable discussions.


text:
# A Monte Carlo Calculation of the Pionium Break-up Probability with Different Sets of Pionium Target Cross Sections

## 1 Introduction

Pionium is the hydrogen-like electromagnetic bound state of a $`\pi ^+\pi ^{}`$ pair. Its lifetime is determined by the strong $`\pi ^+\pi ^{}\pi ^0\pi ^0`$ annihilation in the low relative momentum regime where Chiral Perturbation Theory applies. The theory predicts the pionium lifetime value of $`(2.9\pm 0.1)\times 10^{15}\text{ s}`$ . Experimentally, this value is being measured in a model-independent way by the DIRAC experiment at CERN .

Under the experimental conditions of DIRAC the pionic atoms are created in the inelastic scattering between $`24\text{GeV}/c`$ protons, supplied by the PS at CERN, and the nuclei of the target , a chemically pure material of well-determined thickness.

A pionic atom propagating at relativistic speed inside the target interacts mainly electromagnetically with the target atoms. The electromagnetic cross section is on the order of a megabarn. These interactions can lead either to a transition between two bound states of pionium or to a dissociation (break-up). The scattering with the target atoms competes with the process of annihilation since a shorter pionium lifetime leads to a smaller break-up probability.

The main thrust of DIRAC’s experimental technique is the detection of $`\pi ^+\pi ^{}`$ pairs with very low relative momentum. From the low-momentum part of the spectrum, the break-up probability of pionium is then determined.

The goal of this study is to establish the theoretical dependence of the break-up probability on the lifetime with a very high accuracy. We modeled the dynamics of pionium in different targets taking into account a large number of atomic shells that become populated as the pionic atom evolves in the target.

We have performed the calculations by using three sets of pionium target interaction cross sections. The original studies by Afanasyev and Tarasov  made use of the Born approximation and pure electrostatic interaction. In our work we have applied new corrected cross sections that take into account relativistic effects  and the multiphoton exchanges . In addition, these cross section sets include a more accurate description of the target atomic form factor. As an error in the break-up probability translates directly to an error in lifetime, we have checked whether the results obtained in the new calculations significantly deviate from the previous ones. We conclude that the magnetic and relativistic corrections together with the target form factor choice are non-negligible mainly for the small $`Z`$ targets, whereas the multiphoton exchange should be considered for the large $`Z`$ ones.

## 2 Monte Carlo Simulation of Pionium in the Target

As we have noted in the introduction, to determine the pionium lifetime, the experiment DIRAC needs a precise theoretical calculation of the break-up probability $`P_{\mathrm{br}}`$ of pionium due to the electromagnetic interaction with target atoms. This calculation can be done by means of a Monte Carlo transport code that simulates the evolution of pionium from its creation to either its annihilation or its break-up under the given experimental conditions.

### 2.1 Pionium Production

The pionic atoms are formed as a consequence of the Coulomb final state interaction of two oppositely charged pions. These pions are created in an inelastic collision of a $`24\text{GeV}/c`$ proton and one of the target nuclei. The pionic atom is described by six quantities, accounting for the six degrees of freedom of a two body system. A particularly convenient choice of these quantities is given by the laboratory momentum of the center of mass of the atom $`\stackrel{}{P}`$ and the quantum numbers of the created bound state $`n`$, $`l`$, and $`m`$ in the spherical coordinate representation 111Also parabolic quantum numbers have been used elsewhere ..

The probability of pionium being created is given by

$$\frac{\sigma _{nlm}^A}{\stackrel{}{P}}=(2\pi )^3\left|\psi _{nlm}(0)\right|^2\frac{E}{M}\frac{\sigma _s^0}{\stackrel{}{p}\stackrel{}{q}}|_{\stackrel{}{p}=\stackrel{}{q}=\stackrel{}{P}/2}.$$
(1)
The two terms on the right-hand side of the equation illustrate the final state interaction mechanism. The rightmost factor is the doubly inclusive cross section of $`\pi ^+`$ and $`\pi ^{}`$ pairs at equal momenta ($`\stackrel{}{p}=\stackrel{}{q}`$) without considering the final state interaction, as indicated by the superscript $`0`$. The subscript $`s`$ means that only pions created from direct hadronic processes and decays of resonances with a very short lifetime are considered, because the Coulomb interaction of pions from long-lived sources (e.g. $`\eta `$, $`K_S^0`$ and $`\mathrm{\Lambda }`$) is negligible and hence they do not contribute to pionic atoms production. The effect of the final state Coulomb interaction is to create a bound state with quantum numbers $`n`$, $`l`$, and $`m`$; it is given by the squared wave function at the origin.

The doubly inclusive cross section can be obtained from the direct measurements of time correlated $`\pi ^+\pi ^{}`$ pairs in DIRAC, according to the following reasoning:

* The final state Coulomb interaction for short-lived sources is given, as in the case of the creation of a bound state, by a multiplicative factor depending only on $`Q`$, the magnitude of the relative momentum between the two pions. This is the so-called Coulomb or Gamow factor

$$\frac{\sigma _s}{\stackrel{}{p}\stackrel{}{q}}=A_C(Q)\frac{\sigma _s^0}{\stackrel{}{p}\stackrel{}{q}};A_C(Q)=\frac{2\pi M_\pi \alpha /Q}{1\mathrm{}^{2\pi M_\pi \alpha /Q}},$$
(2)
where $`\alpha `$ is the fine structure constant.
* The contribution to the doubly inclusive cross section of pairs containing at least one pion from a long-lived source, $`\omega _l(\stackrel{}{P})`$, can be calculated with a hadron physics Monte Carlo simulation. In our case we have used FRITIOF6 . This function has been shown to depend only on $`P`$ , the magnitude of the total momentum of the pion pair. Taking this into account together with (2) we find

$$\frac{\sigma }{\stackrel{}{p}\stackrel{}{q}}=\frac{\sigma _s}{\stackrel{}{p}\stackrel{}{q}}+\frac{\sigma _l}{\stackrel{}{p}\stackrel{}{q}}=A_C(Q)\frac{\sigma _s^0}{\stackrel{}{p}\stackrel{}{q}}+\omega _l(P)\frac{\sigma }{\stackrel{}{p}\stackrel{}{q}},$$
(3)
thus relating $`\sigma `$ and $`\sigma _s^0`$.
* Finally, we have found that the $`\stackrel{}{P}`$-dependence of the doubly inclusive cross section is not correlated to $`\stackrel{}{Q}`$, given that $`\stackrel{}{Q}30\text{MeV}/c`$.

These findings allow us to relate the $`\stackrel{}{P}`$-dependence of $`\sigma `$ and $`\sigma _s^0`$ by

$$\frac{\sigma _s^0}{\stackrel{}{p}\stackrel{}{q}}|_{\stackrel{}{p}=\stackrel{}{q}=\stackrel{}{P}/2}_0^{Q2\mathrm{MeV}/c}(1\omega _l(P))\frac{\sigma }{\stackrel{}{p}\stackrel{}{q}}\stackrel{}{Q},$$
(4)
where the $`\stackrel{}{P}`$ distribution is obtained from the direct measurement of the laboratory momentum of low relative momentum $`\pi ^+\pi ^{}`$ pairs in DIRAC. In figure 1 we show the distribution of the magnitude of the momentum $`P`$ and the angular distribution relative to the proton beam axis for low relative momentum $`\pi ^+\pi ^{}`$ pairs.

The initial quantum number distribution depends on the value of the wave function at the origin. It has been shown  that the effect of the strong interaction between the two pions of the atom significantly modifies $`|\psi _{nlm}(0)|`$ in comparison to the pure Coulomb wave function. However, the ratio between the production rate in different states has been demonstrated to be kept as for the Coulomb wave functions . Then, considering that the Coulomb functions obey

$$\left|\psi _{nlm}^{(\mathrm{C})}(0)\right|^2=\{\begin{array}{cc}0& \text{if }l0\text{,}\hfill \\ \frac{(\alpha M_\pi /2)^3}{\pi n^3}& \text{if }l=0\text{,}\hfill \end{array}$$
(5)
we note that only $`S`$ states are created, according to a $`1/n^3`$ distribution.

Another quantity to be specified is the position where the proton–target interaction took place ($`\stackrel{}{R}`$). This is also the position where the atom was created. Since the target thickness is chosen much smaller than the nuclear interaction length of the target material, the $`\pi ^+\pi ^{}`$ atoms are supposed to be uniformly generated all across the target thickness. The position in the transverse coordinates is unimportant, but it can be generated according to the beam profile characteristics, too.

We have now gathered all the information needed to simulate the creation of pionium atoms with a center of mass momentum $`\stackrel{}{P}`$ according to the experimentally measured momentum and angular distributions, at a position $`\stackrel{}{R}`$ uniformly distributed through the target, and in an initial $`S`$-wave state with the principal quantum number $`n`$ distributed according to a $`1/n^3`$ distribution.

### 2.2 Pionium Annihilation

Once an atom has been created in its initial state, specified by $`\stackrel{}{P}`$, $`n`$, $`l`$, and $`m`$, its dynamics are those of a free system that can either be annihilated, mainly via the $`\pi ^+\pi ^{}\pi ^0\pi ^0`$ channel, or be electromagnetically scattered by one of the target atoms.

The strong interaction decay to two neutral pions determines the lifetime and is related to the $`a_0^0a_0^2`$ scattering lengths difference and to the wave function at the origin by

$$\frac{1}{\tau _{nlm}}=\frac{16\pi }{9}\frac{\sqrt{M_\pi ^2M_{\pi ^0}^2\frac{1}{4}M_\pi ^2\alpha ^2}}{M_\pi }(a_0^0a_0^2)^2(1+\delta _\mathrm{\Gamma })\left|\psi _{nlm}^{(\mathrm{C})}(0)\right|^2,$$
(6)
where $`M_\pi `$ and $`M_{\pi ^0}`$ are the masses of the charged and the neutral pion, respectively, and $`\delta _\mathrm{\Gamma }`$ is the correction to next-to-leading order ($`\delta _\mathrm{\Gamma }=0.058`$) that includes the effect of the strong interaction between the two pions. Using Chiral Perturbation Theory, Colangelo et al  have been able to calculate the most precise value of the scattering lengths difference to date ($`a_0^0a_0^2=0.265\pm 0.004`$). Employing this value in (6) yields

$$\tau _{100}=(2.9\pm 0.1)10^{15}\text{s}.$$
(7)
Note, however, that due to (5) pionium may only decay from $`S`$ states. Moreover, the lifetime of any $`S`$ state is related to the lifetime of the ground state, by

$$\tau _{n00}=n^3\tau .$$
(8)
For the purpose of simulating pionium in the target, we shall from now on refer to $`\tau `$ as the *pionium lifetime*.

Hence, the probability for a $`\pi ^+\pi ^{}`$ atom to annihilate per unit length, after the Lorentz boost transformation to the laboratory system, is given by

$$p_{nlm}^{\mathrm{anh}}=\frac{1}{\lambda _{nlm}^{\mathrm{anh}}}=\{\begin{array}{cc}\frac{2M_\pi }{Pn^3\tau }& \text{if }l=0\text{,}\hfill \\ 0& \text{other cases,}\hfill \end{array}$$
(9)
where $`\lambda _{nlm}^{\mathrm{anh}}`$ is the annihilation mean free path.

### 2.3 Electromagnetic Interaction of Pionium with the Target

The electromagnetic pionium–target interaction of a pionic atom in an initial $`nlm`$ bound state can induce a transition to another $`n^{}l^{}m^{}`$ bound state. The probability of such an interaction per unit length is given by

$$p_{nlm}^{n^{}l^{}m^{}}=\frac{\rho N_0}{A}\sigma _{nlm}^{n^{}l^{}m^{}}$$
(10)
where $`\rho `$ is the target density, $`A`$ its atomic weight, $`N_0`$ is the Avogadro number, and $`\sigma _{nlm}^{n^{}l^{}m^{}}`$ are the discrete (bound–bound) transition cross sections.

The break-up mechanism is analogous to the discrete one; the break-up probability per unit length of an atomic bound state $`nlm`$ is given by

$$p_{nlm}^{\mathrm{br}}=\frac{1}{\lambda _{nlm}^{\mathrm{br}}}=\frac{\rho N_0}{A}\sigma _{nlm}^{\mathrm{br}}$$
(11)
where $`\sigma _{nlm}^{\mathrm{br}}`$ is the break-up (ionization) cross section.

Finally, the total cross section gives the probability of an atom to undergo an electromagnetic interaction and of course fulfills

$$\sigma _{nlm}^{\mathrm{em}}=\underset{n^{}l^{}m^{}}{}\sigma _{nlm}^{n^{}l^{}m^{}}+\sigma _{nlm}^{\mathrm{br}}.$$
(12)
The total probability for a pionic atom to suffer an electromagnetic collision per unit length is then given by

$$p_{nlm}^{\mathrm{em}}=\frac{1}{\lambda _{nlm}^{\mathrm{em}}}=\frac{\rho N_0}{A}\sigma _{nlm}^{\mathrm{em}},$$
(13)
where $`\lambda _{nlm}^{\mathrm{em}}`$ is the mean free path before an electromagnetic interaction takes place. Exploiting the completeness of the eigenstates of the Coulomb Hamiltonian the total electromagnetic cross sections can be calculated directly  and not just via (12222Note, however, that this is strictly true only within the framework of the sudden approximation ..

The electromagnetic cross sections have been obtained with different approaches in . We will devote section 5 to discussing the different break-up probabilities they lead to.

To get an insight into the magnitude of these interaction probabilities we show in figure 2 the average values of the annihilation, ionization, excitation and de-excitation probabilities per unit length. The average is taken over the even $`z`$-parity states (i.e., $`lm`$ even) for fixed $`n`$. The atoms are created in even $`z`$-parity states ($`l=m=0`$) and the transitions to odd $`z`$-parity ones are strongly suppressed. The figure shows the probabilities using the coherent (interaction with the atom as a whole) contribution of the Born2 set of cross sections. This cross section set will be described in section 4. Any other choice among the cross section sets described in section 4 would lead to very similar results. The averages are defined as

$`\overline{p}_n^{\mathrm{anh}}={\displaystyle \frac{1}{n(n+1)/2}}{\displaystyle \underset{lm}{}}p_{nlm}^{\mathrm{anh}},`$ (14)
$`\overline{p}_n^{\mathrm{br}}={\displaystyle \frac{1}{n(n+1)/2}}{\displaystyle \underset{lm}{}}p_{nlm}^{\mathrm{br}},`$ (15)
$`\overline{p}_n^{n^{}<n}={\displaystyle \frac{1}{n(n+1)/2}}{\displaystyle \underset{lm}{}}{\displaystyle \underset{n^{}<n,l^{}m^{}}{}}p_{nlm}^{n^{}l^{}m^{}},`$ (16)
$`\overline{p}_n^{n^{}>n}={\displaystyle \frac{1}{n(n+1)/2}}{\displaystyle \underset{lm}{}}\left(p_{nlm}^{\mathrm{em}}p_{nlm}^{\mathrm{br}}{\displaystyle \underset{n^{}n,l^{}m^{}}{}}p_{nlm}^{n^{}l^{}m^{}}\right),`$ (17)
where $`n(n+1)/2`$ is the number of even $`z`$-parity states for a given $`n`$.

### 2.4 Pionium Evolution in the Target

To simulate the evolution of a pionic atom we use the following algorithm:

1. We generate a laboratory momentum $`\stackrel{}{P}`$, an initial set of quantum numbers and an initial position $`\stackrel{}{R}`$ for the atom as described in subsection 2.1.
2. We generate a free path according to:

$$p(x)x=\frac{1}{\lambda _{nlm}}\mathrm{}^{x/\lambda _{nlm}}x$$
(18)
where $`\lambda _{nlm}=(1/\lambda _{nlm}^{\mathrm{em}}+1/\lambda _{nlm}^{\mathrm{anh}})^1`$ is the mean free path before either an electromagnetic interaction or the annihilation takes place.
3. We displace the atom by the distance $`x`$:

$$\stackrel{}{R}^{}=x\frac{\stackrel{}{P}}{P}+\stackrel{}{R}.$$
(19)
4. We determine whether the atom has been annihilated, excited (or de-excited) in a discrete collision, or broken up. The relative weights of the respective branches of the evolution are given by the probabilities of equations (9), (10), and (11).
5. If the atom has been scattered and suffered a discrete transition we return to step (ii) using the new quantum numbers $`n^{}`$, $`l^{}`$ and $`m^{}`$ and the new position $`\stackrel{}{R}^{}`$ as the initial values.

More details on this model may be found in .

## 3 Break-up Probability Calculation

In principle, the break-up probability calculation of pionium should be straightforward once we have established the Monte Carlo model. The rest would be a matter of generating an atom sample and computing how many of them break up in the target. However, two main difficulties arise when trying to implement the algorithm.

The first difficulty is due to the presence of an infinite number of atomic bound states in the calculations. Clearly, only a finite number of states can be taken into account in the simulation of the evolution of pionium. In our calculations we have imposed a cut on the states with $`nn_{\mathrm{max}}`$. This would not pose a serious problem if the atoms, being created mainly in very low $`n`$ states, could not get highly excited. Unfortunately, excitation to ever higher lying bound states constitutes a major branch in the evolution of pionium. As a consequence we cannot directly calculate the break-up probability as outlined in the previous paragraph.

The other difficulty lies in the fact that for some of the cross section sets to be studied in section 4, the break-up cross sections have not been calculated. In this case it is imperative to find an indirect way to compute the break-up probability.

We have discussed in the previous section that pionium terminates its evolution in the target by being either annihilated or broken up. However, the atom can also leave the target in a bound state. This would happen if one of the generated free paths in the Monte Carlo procedure carries it to a position outside the target. The break-up probability ($`P_{\mathrm{br}}`$), the annihilation probability ($`P_{\mathrm{anh}}`$), and the probability to leave the target in a discrete state ($`P_{\mathrm{dsc}}`$) are related by:

$$1=P_{\mathrm{br}}+P_{\mathrm{anh}}+P_{\mathrm{dsc}}.$$
(20)
This equation allows us to compute the break-up probability indirectly.

### 3.1 Computation Difficulties Due to Physical Characteristics of the Problem

The probability to generate an atom in a specific shell decreases as $`1/n^3`$. This means that the number of atoms created with $`n4`$ is very small. If the atoms could not get excited to states with large $`n`$, we could safely solve the evolution system by setting $`n_{\mathrm{max}}>4`$. However, as we saw in figure 2, the atoms have a tendency to be excited, as $`n`$ increases, rather than being annihilated or ionized.

Hence we expect a significant fraction of atoms excited into $`n>n_{\mathrm{max}}`$ shells even for large values of $`n_{\mathrm{max}}`$. The probability of an atom in a $`nlm`$ state to be excited into a state beyond the cut, i.e. with $`n>n_{\mathrm{max}}`$, is given by

$$p_{nlm}^{n^{}>n_{\mathrm{max}}}=p_{nlm}^{\mathrm{em}}p_{nlm}^{\mathrm{br}}\underset{n^{}n_{\mathrm{max}},l^{}m^{}}{}p_{nlm}^{n^{}l^{}m^{}}$$
(21)
where we have used (10), (11) and (13). However, once the atom jumps into one of these states we loose control over it and we have to stop its evolution.

To analyze the change of the Monte Carlo results with $`n_{\mathrm{max}}`$ we have modeled the evolution of a sample of atoms by changing $`n_{\mathrm{max}}`$ from 7 to 9. We observed three main effects:

* The fraction of annihilated atoms ($`P_{\mathrm{anh}}(nn_{\mathrm{max}})`$) does not change significantly.
* The portion of atoms leaving the target in discrete states ($`P_{\mathrm{dsc}}(nn_{\mathrm{max}})`$) changes only slightly.
* The fraction of dissociated atoms ($`P_{\mathrm{br}}(nn_{\mathrm{max}})`$) changes significantly.

This effect can be understood by checking the dependence of the annihilation, the discrete and the break-up probabilities on $`n`$, the principal quantum number of the state from which the atom was annihilated, broken-up, or in which it left the target. In figure 3 we show the result of the Monte-Carlo simulation with $`n_{\mathrm{max}}=8`$ for a sample of one million atoms using the Born2 cross section set that also includes cross sections for the ionization (refer to section 4). For the annihilated atoms we can see that $`P_{\mathrm{anh}}(n)`$ is negligible for values of $`n4`$. The $`P_{\mathrm{dsc}}(n)`$ dependence also shows a fast, but less drastic, decrease with $`n`$. Only the solution for the states with $`n=n_{\mathrm{max}}1`$ or $`n=n_{\mathrm{max}}2`$ is unstable under variation of $`n_{\mathrm{max}}`$. For $`n_{\mathrm{max}}=8`$ this is a small contribution to the total $`P_{\mathrm{dsc}}`$ value. Finally, $`P_{\mathrm{br}}(n)`$ decreases very slowly with $`n`$, showing that there is a significant fraction of atoms broken up from states with $`n>n_{\mathrm{max}}`$. The probability of an atom to be excited into such a state with $`n>n_{\mathrm{max}}`$ is also shown as a function of the principal quantum number of the last state before the excitation. Obviously, this effect is non-negligible.

For the cross section sets without break-up cross sections, we can calculate directly only the total probability for all electromagnetic processes and the probabilities for discrete transitions to states with $`n^{}<n_{\mathrm{max}}`$. In these cases we can therefore not distinguish whether an atom has been broken up or excited into a state with $`n^{}>n_{\mathrm{max}}`$, that is, we can only determine the *combination* of probabilities

$$p_{nlm}^{n^{}>n_{\mathrm{max}}}+p_{nlm}^{\mathrm{br}}=p_{nlm}^{\mathrm{em}}\underset{n^{}n_{\mathrm{max}},l^{}m^{}}{}p_{nlm}^{n^{}l^{}m^{}}.$$
(22)
This is, of course, equivalent to (21), but in this case $`p_{nlm}^{\mathrm{br}}`$ is unknown. Thus, for these cross section sets not even the break-up probability for low $`n`$ states could be directly calculated and we are forced to use the procedure described below.

### 3.2 Calculation Procedure

Based on the fast decrease of $`P_{\mathrm{anh}}`$ and $`P_{\mathrm{dsc}}`$ as a function of $`n`$ we can assume that almost every atom excited to a state $`n>n_{\mathrm{max}}`$ will be eventually broken up. This will be true even though the excitation probability per unit length of a given bound state is significantly larger than the break-up probability per unit length. We can explain it as follows. The mean free path of the excited atoms strongly decreases with increasing $`n`$. For $`n8`$ the mean free path is $`0.1\mu `$m. An excited atom will thus interact many times within a very short distance. In every scattering the atom will have some small probability to break up, thereby terminating its evolution. In summary, the most probable evolution of an atom that has been excited to any state with $`n4`$ is a sequence of excitations (and less frequent de-excitations) terminated by break-up.

However, while we are neglecting the atomic annihilation from states with $`n>8`$ and thus setting $`P_{\mathrm{anh}}=P_{\mathrm{anh}}(n8)`$, we can estimate $`P_{\mathrm{dsc}}(n>8)`$ by means of a fit to the $`P_{\mathrm{dsc}}(n)`$ histogram as recommended in

$$P_{\mathrm{dsc}}(n>8)=\frac{a}{n^3}+\frac{b}{n^5}.$$
(23)
Hence, taking into account (20) we obtain

$$P_{\mathrm{br}}=1P_{\mathrm{dsc}}P_{\mathrm{anh}},$$
(24)
where $`P_{\mathrm{dsc}}`$ consists of two parts,

$$P_{\mathrm{dsc}}=P_{\mathrm{dsc}}(n8)+P_{\mathrm{dsc}}(n>8),$$
(25)
of which $`P_{\mathrm{dsc}}(n8)`$ is computed directly and $`P_{\mathrm{dsc}}(n>8)`$ is calculated from (23). In this manner we can calculate the break-up probability even without ionization cross sections as input.

In table 1 and in figure 4 (top left) we show a few results for the probability for different lifetime values in a $`95\mu `$m Ni target. The target choice coincides with that of the DIRAC experiment. We observe that the result of $`P_{\mathrm{dsc}}(n>8)`$ adds only a small correction. In figure 4 we also show the ionization and annihilation distributions as a function of the target position, and finally the creation position for those atoms that managed to emerge from the target in a bound state. As emphasized in subsection 3.1, with increasing $`n`$ only the atoms very near the target end will be able to leave the target in a discrete state.

## 4 Cross Sections Sets

In our calculations of the break-up probability we employed three different sets of cross sections. The first two have been calculated in the framework of the Born approximation. We assign the labels Born1 to the calculations made in reference  and Born2 to those of . The two sets differ in four main points:

* The Born1 set neglects the contribution of incoherent scattering (collisions leading to an excitation of the target atom), thus considering the coherent interaction only (collisions with the target atom as a whole), i.e. the leading term. By contrast the Born2 set accounts for target excitations.
* The Born1 set uses Molière’s parameterization  for the Thomas-Fermi equation solution as the target atom form factor of the pure electric interaction, whereas the Born2 set takes electron orbitals determined numerically within the Hartree-Fock framework for the same purpose. The Thomas-Fermi-Molière parameterization of the atomic form factor is accurate for low momentum exchange, but gives a small excess for harder scattering.
* The Born1 set considers the sudden approximation (no recoil energy for the target and the pionic atom) and neglects the energy difference between the initial and the final state, while the Born2 set accounts for these two effects.
* Finally, Born2 set also considers the effect of magnetic and relativistic terms.

In principle it has been concluded  that accounting for second order effects like the magnetic terms of the Hamiltonian, the recoil energy of the atoms, or the relativistic terms generally leads to an overall decrease of the sudden approximation pure electrostatic coherent cross section value due to destructive interference with the leading orders. Moreover, employing atomic orbitals obtained in the Hartree-Fock approach for the form factor used to compute these cross sections leads to lower values than those of the Born1 set since the Molière parameterization of the solution to the Thomas-Fermi equation is excessive for mean and large values of the photon momentum transfer. This last issue leads to discrepancies that increase for large $`n`$ states and decrease for large $`Z`$ target atoms. The difference appears to be balanced by neglecting the incoherent contribution to the cross section in the Born1. This results in a systematically smaller ground state cross section of Born1 set whereas for larger $`n`$ Born1 cross sections are larger (up to $`10\%`$ discrepancy) or compatible with Born2 results. We can observe the comparison for three target materials in figure 5 and we shall analyze the disagreement in the break-up probability resulting from this effect.

Finally we have also used a set of cross sections where the Glauber formalism has been applied to calculate the coherent contribution to the cross section value. The details are given in . This calculation technique accounts for multi-photon exchange in the pionium–target atom collision. Contrary to what one would expect the consideration of more than one photon being exchanged diminishes the values of the cross sections due to a destructive interference of the $`n`$-photon exchange contributions (this happened also when accounting for magnetic terms in the Born2 set). The leading order of the Glauber result matches the sudden approximation of the Born cross sections (since both neglect the difference between the initial and the final state energies). However, this cross section set uses a parameterization for the target atom form factors similar to the ones used in Born2. This explains the disagreement with respect to the Born1 and the agreement with Born2 set for low $`Z`$ targets, as can be seen in figure 5. The corrections due to multiphoton exchange are important for large $`Z`$ targets  and this explains the large discrepancies obtained for Platinum.

## 5 Results and Conclusions

After the discussions of the previous sections, we finally present the results of the break-up probability calculation. The DIRAC experiment has the possibility of choosing between several targets. The design of these different targets has been made to achieve the maximum break-up probability resolution in different lifetime ranges. Large $`Z`$ targets with larger interaction cross sections are better suited for small lifetime values whereas lower $`Z`$ materials are more sensitive to larger lifetime values. Three of these targets are the Pt $`28\mu `$m target, suitable for lifetime ranges $`\tau <110^{15}`$s, the Ni $`95\mu `$m target for $`\tau 310^{15}`$s and the Ti $`251\mu `$m target for $`\tau 410^{15}`$s. The target thickness was chosen so as to have the same radiation length and hence equivalent multiple scattering effects for all three targets. The Nickel target constitutes DIRAC’s main target with which 90% of the data have been collected, as it is optimal for the theoretically predicted lifetime value.

In figure 6 the break-up probability curves are shown for these three targets. The calculation has been carried out for samples of ten million events, with a statistical error less than $`0.08\%`$.

One can clearly see that for the Ti and Ni targets the Glauber and Born2 sets lead to similar results whereas the Born1 set shows a $`8\%`$ disagreement. For the large $`Z`$ target (Pt) both Born1 and Born2 are biased toward large values. In this case, the multi-photon contributions to the total cross sections are not negligible. In any case the discrepancies between the break-up probability results are at the level of the discrepancies between ground state cross sections and much smaller than the differences between cross sections sets for medium or highly excited states. We can understand this based on the fact that the probability for the atoms to leave the target in a discrete state other than the ground state and maybe the first excited shell $`P_{\mathrm{dsc}}(n2)`$ is of the order of, or smaller than, 5%. Hence, even large uncertainties in this magnitude (up to $`1015\%`$) lead to very small changes in the break-up probability result. Only discrepancies in the ground state population and maybe the first excited shell, where most of the atoms are created, would lead to significant differences between the break-up probability results of the different sets.

Graphically we can view the atom as a balloon being inflated in every collision with the target. The different sets will lead to similar results of the size increase rate as long as the atom remains in a low excited state. However, as the atom grows (inflates) it will no longer be able to advance as easily in the target due to its large size and will finally break-up (explode). Large discrepancies in the excitation and break-up rate of the excited atom will not be important given that the mean free paths for the excited states are very small compared to the target dimensions.

In summary, we recall that the high precision measurement attempted by the DIRAC collaboration also requires an accuracy to better than $`1\%`$ in our theoretical break-up probability calculations. We note that the seemingly large discrepancies among our different cross section sets particularly for pionium transitions starting from highly excited states do not lead to significant differences in the theoretical break-up probabilities. The discrepancies between break-up probabilities are coming almost entirely from differences in the cross sections for the lowest lying states, where both the atomic structure of the target and the multi-photon transitions need to be treated as accurately as possible. This challenge, however, has already been mastered in our previous works  where we showed that the required 1% accuracy can be achieved with our calculations, albeit only with the Born2 and the Glauber sets for low $`Z`$ and with the Glauber set for large $`Z`$ targets. The important conclusion of the present investigation is the finding that the (infinitely many!) highly excited states of pionium do *not* limit the validity of our approach even though we can include explicitly only a moderate number of these states in our simulations.

We would like to thank B. Adeva, L.L. Nemenov, L. Tauscher and D. Trautmann for their support. We are also indebted to M. Pló, J.J. Saborido and A.V. Tarasov for their invaluable help.

## References


text:
# FUSE observations of H2 around the Herbig AeBe stars HD 100546 and HD 163296

## 1 Introduction

HD 100546 and HD 163296 are two Herbig AeBe stars harboring young circumstellar disks. Herbig AeBe stars are intermediate-mass (2–10 M) young emission-line stars still in the pre-main sequence phase (Herbig Herbig1960 (1960)). Material surrounding Herbig stars is thus considered as the constituent of protoplanetary disks, similar to the disks seen around the less massive T-Tauri stars. These disks are the material from which (and the places where) planets are supposed to form according to the standard model of planetary formation (Lissauer Lissauer1993 (1993)).

HD 100546 (B9Vne, $`d`$=103$`{}_{6}{}^{}{}_{}{}^{+7}`$ pc, $`t`$$`>`$10 Myr; van den Ancker et al. Ancker1998 (1998)) is older than HD 163296 (A1Ve, $`d`$=122$`{}_{13}{}^{}{}_{}{}^{+17}`$ pc, $`t`$=$`4_{2.5}^{+6}`$ Myr; van den Ancker et al. Ancker1998 (1998)). However, both stars present very similar characteristics. In the two cases, ISO detected silicate features resembling cometary materials (Malfait et al. Malfait1998 (1998); van den Ancker et al. Ancker2000 (2000); Bouwman et al. Bouwman2000 (2000), Bouwman2001 (2001)). Disks have recently been imaged around both stars (Pantin et al. pantin (2000); Grady et al. grady2000 (2000), grady2001 (2001); Augereau et al. augereau (2001)). The images show dusty disks at intermediate inclination ($``$60°). Absorption spectroscopy also revealed the presence of circumstellar gas around HD 100546, including spectral variabilities which are interpreted as accreting material (Grady et al. grady1996 (1996), Grady1997 (1997)).

In the case of HD 163296, a CO disk of about 310 AU in semi-major axis has been detected (Mannings & Sargent Mannings (1997)). H2 infrared emission lines have also been reported by Thi et al. (thi (2001)) at 17 and 28$`\mu `$m. From their ISO observations and assuming thermal equilibrium up to $`J=3`$, they derive a total H2 mass of $`0.4(\pm 0.2)\times 10^3`$ M, and the ratio of the two emission lines ($`v=0`$, $`J`$=2$``$0 and $`J`$=3$``$1) gives a temperature of $`220(\pm 22)`$ K. However, this detection has been challenged by more recent ground based observation at 17$`\mu `$m which shows no detection with three times better sensitivity (Richter et al. richter (2002)).

These controversial observations raise the question of the molecular content of such disks, where most of the mass is in the form of molecules, and particularly H2. The formation of giant planets requires a large reservoir of molecular gas. Moreover some of these planets migrate close to their parents where they are observed. This migration also needs a massive disk to allow the angular momentum exchange between the migrating planet and the disk. The H2 content is thus a key ingredient in the recipe for the formation of giant planets. But H2 is a symetrical molecule, and infrared emission by quadrupole rotational transitions is very inefficient. In contrast, when seen against a UV bright source, far-UV absorption lines due to electronic transitions allow sensitive observation of H2 even at low temperatures. FUSE offers a unique opportunity to scrutinize in detail the H2 content of protoplanetary disks around Herbig stars.

## 2 Observation and data analysis

HD 100546 was observed with FUSE through the LWRS aperture ($`30\mathrm{}`$$`\times `$$`30\mathrm{}`$) for a total time of 5.8 hours on March 26, 2000 (Program P1190303) and March 3, 2002 (Program P2190401) (for an overview of FUSE, see Moos et al. moos (2000) and Sahnow et al. sahnow (2000)). HD 163296 was observed twice, on April 27, 2001 and April 29, 2001 for a total time of 8.9 hours (Programs P2190601 and Q2190101, respectively). The data of both targets were reprocessed with the version 2.0.5 of the CALFUSE pipeline. The output of the pipeline is a total of 5 and 10 sub-exposures for HD 100546 and HD 163296, respectively. The sub-exposures have been aligned and coadded resulting in a set of four independent spectra, one for each FUSE channel (2 LiF spectra and 2 SiC spectra).

The version 2.0.5 of the CALFUSE pipeline is known to slightly over-estimate the tabulated errors on each pixel. We compared the data used in the present work with data of the same observations but obtained with the version 2.2.1 of the CALFUSE pipeline; we conclude that the error bars given below are not significantly affected by the improvement in the error propagation of the different pipelines.

Apart from the observation of the molecular hydrogen, the spectra are rich in emission and absorption lines from atomic and ionic species. For instance, a large number of Fe ii lines from the ground level as well as from excited levels are clearly detected. As already observed in AB Aur and $`\beta `$ Pic (Roberge et al. Roberge2001 (2001); Deleuil et al. deleuil2001 (2001)), bright emission from C iii and O vi are also detected in the two spectra. All these features indicate the presence of circumstellar material, chromospheric activity and/or accretion (Bouret et al. bouret2002 (2002)). A detailed analysis of the whole FUSE spectrum will be made in a forthcoming paper (Deleuil et al. deleuil2003 (2003)). Here we focus on the analysis of the H2 lines, probing the molecular portion of the disk.

In the far-UV, HD 100546 is brighter than HD 163296 and the observed spectrum has a better S/N ratio allowing detection of H2 in its pure rotational level ($`v=0`$) up to $`J=9`$ (a line of $`v=0`$, $`J=10`$ is also marginally detected at 1058.6 Å) and in its first vibrational level ($`v=1`$) up to $`J=5`$ (Fig. 1). In the HD 163296 spectrum, H2 lines are detected up to $`v=0`$, $`J=4`$. Higher $`J`$-levels are beyond the detection limit. The data analysis and profile fitting has been done using the Owens code kindly made available to us by Dr. M. Lemoine (see for example Lemoine et al. lemoine (2002) and Hébrard et al. hebrard (2002)). For the electronic transitions, we used the wavelengths and oscillator strengths tabulated by Abgrall et al. (1993a  and 1993b  for the Lyman and the Werner system, respectively), and the inverses of the total radiative lifetimes tabulated by Abgrall et al. (abgrall2000 (2000)).

With low H2 column densities in HD 100546 ($`N10^{16}`$cm-2 for each $`J`$-level), we can choose to fit only the unsaturated lines having low oscilator strength ($`f10^2`$). This allows us to avoid systematic errors in the estimates of column densities which could be included by the fit of saturated lines with uncertain instrumental line spread function. Note, however, that additional fits have been performed with saturated lines ($`f10^2`$). They give very similar column densities and allow the determination of the intrinsic line width ($`b`$). In HD 163296 the column densities are larger, and all the H2 lines are saturated. Therefore the error bars on the column densities are also larger than in HD 100546. An example of the fit of two H2 lines ($`J=3`$ and $`J=4`$) which are superimposed on the O vi stellar emission line is given in Fig. 2.

The estimated column densities are tabulated with $`2\sigma `$ error bars (Table 1). We obtain the total column densities

$$N_{\mathrm{H}_2}(\mathrm{HD}\mathrm{\hspace{0.17em}100546})2.8\times 10^{16}\mathrm{cm}^2$$

and

$$N_{\mathrm{H}_2}(\mathrm{HD}\mathrm{\hspace{0.17em}163296})1.4\times 10^{18}\mathrm{cm}^2.$$

The error bars are estimated by the classical method of the $`\mathrm{\Delta }\chi ^2`$ increase of the $`\chi ^2`$ of the fit (see Hébrard et al. (hebrard (2002)) for a full discussion on the fit method and error estimates with owens profile fitting). These error bars include the uncertainties in the continuum fits, intrinsic line widths ($`b`$) and the instrumental line spread function. We note that the observation of different lines with different oscillator strengths for the same species allows us to constrain all these quantities. In particular the fit to saturated lines constrains the line widths and the instrumental line spread function. This last one is allowed to slowly vary with the wavelength. The results are obtained from a final self-consistent fit to all the data in which the continuum, the instrumental line spread function, and physical parameters of the absorption lines are free parameters and estimated by the determination of the best $`\chi ^2`$ in the parameters’ space. For the intrinsic line widths, we obtain $`b_{\mathrm{HD}\mathrm{\hspace{0.17em}100546}}=3.4_{0.5}^{+0.8}`$km s-1 and $`b_{\mathrm{HD}\mathrm{\hspace{0.17em}163296}}=2.2_{0.7}^{+0.5}`$km s-1.

## 3 Results

### 3.1 Excitation diagrams

The excitation diagrams corresponding to the derived H2 column densities are presented in Fig. 3 and 4, for HD 100546 and HD 163296, respectively. For the $`v=1`$ levels observed toward HD 100546, we see that the column densities are above a simple extrapolation of the $`v=0`$ levels as plotted on the excitation plot (Fig. 3). This effect is also observed in the vibrationally excited interstellar H2 detected toward HD 37903 (Meyer et al. meyer (2001)), where the high energy levels are not thermally populated. Note also that the level degeneracies ($`g`$) of the $`v=1`$ levels are smaller than the degeneracies of the $`v=0`$ levels with similar energy, hence in the same region of the plot. The plotted values of $`N/g`$ are thus larger for these $`v=1`$ levels than for the nearby $`v=0`$ levels (Fig. 3).

### 3.2 Turbulence

The simultaneous determination of temperature and line broadening can be used to constrain the turbulence of the gas. The observed line width $`b`$ is a combination of the thermal broadening, the intrinsic turbulence $`v_{\mathrm{turb}}`$, and the projected radial component of the gas motion seen along the star diameter. With the knowledge of the temperature, and considering that the motion can only increase the line width, we can estimate a lower limit on the turbulence velocity of the gas. We have

$$b\sqrt{\frac{2kT}{\mu _{\mathrm{H}_2}}+v_{\mathrm{turb}}^2}$$

where $`\mu _{\mathrm{H}_2}`$ is mass of molecular hydrogen. This can be written as

$$\left(\frac{v_{\mathrm{turb}}}{\mathrm{km}\mathrm{s}^1}\right)^2\left(\frac{b}{\mathrm{km}\mathrm{s}^1}\right)^28.31\times 10^3\frac{T}{\mathrm{K}}.$$

Using the temperature given above and the 2$`\sigma `$ upper limits on the measured line widths ($`b_{\mathrm{HD}\mathrm{\hspace{0.17em}100546}}^{2\sigma }4.2`$ km s-1 and $`b_{\mathrm{HD}\mathrm{\hspace{0.17em}163296}}^{2\sigma }2.7`$ km s-1), we find that $`v_{\mathrm{turb},\mathrm{HD}\mathrm{\hspace{0.17em}100546}}`$$``$$`3.4`$ km s-1 and $`v_{\mathrm{turb},\mathrm{HD}\mathrm{\hspace{0.17em}163296}}`$$``$$`1.9`$ km s-1.

### 3.3 Nature and origin of the H2 gas

Many characteristics show that the detected H2 is not interstellar but circumstellar; this favours the interpretation that the detected lines arise from protoplanetary material. First, in both cases the radial velocity of H2 is similar to the radial velocity of the atomic lines due to the circumstellar gas. Lines from exited levels of Fe ii (like Fe ii* and Fe ii\**, Deleuil et al. deleuil2003 (2003)) cannot be due to the interstellar medium. They must be linked to dense circumstellar material. The similarity of their radial velocity with the radial velocity of H2 is a first clue that the H2 is linked with the circumstellar matter.

Importantly, the excitation diagrams show that the H2 is thermalized up to about $`J=4`$, with high temperatures. Hot H2 has also been observed around the T-Tauri star TW Hya ($`T`$=3000 K, Herczeg et al. herczeg2002 (2002)). Here we obtain

$$T740\pm 30\mathrm{K}\mathrm{for}\mathrm{HD}\mathrm{\hspace{0.17em}100546}$$

$$\mathrm{and}T430\pm 20\mathrm{K}\mathrm{for}\mathrm{HD}\mathrm{\hspace{0.17em}163296}.$$

Seen in absorption, it is very unusual to observe such a thermalization at high temperature up to $`J=4`$ (to our knowledge, this is the first case). This shows that the absorbing medium has very particular physical conditions, different from what is commonly seen in the diffuse interstellar medium. Indeed this requires a line of sight which presents simultaneously high density for the thermalization and low thickness to have column densities low enough for a detection of the far-UV stellar continuum. In the diffuse and transluscent interstellar media, levels higher than $`J=1`$ or $`J=2`$ are populated by mechanisms other than the collisional pumping. In these low density media, the collisional pumping is less efficient than the UV-pumping and mecanical processes such as C shocs or turbulence. Here levels up to $`J=4`$ show a dominant collisional pumping. This is a clue of a high density medium.

Note here that for the levels higher than $`J=4`$ as detected in HD 100546, it is not possible to discriminate between a collisional pumping within a higher temperature component of the absorber and other mecanical or radiative processes.

The observed high temperatures give additional clues that the detected H2 is circumstellar and close to the exciting stars. Using the numerical model developed by Le Bourlot et al. (lebourlot1993 (1993)), we find that the incident radiation on the H2 must be $`10^4`$ times larger than the mean galactic UV radiation field. This shows that the observed H2 is circumstellar and really close to its central star. These temperatures correspond to $``$1.5 AU from HD 100546 and $``$4 AU from HD 163296.

Most importantly, the thermalization allows us to estimate the H2 volume density ($`n_{\mathrm{H}_2}`$). Using the critical density given by Le Bourlot et al. (lebourlot1999 (1999)), we find :

$$n_{\mathrm{H}_2}(\mathrm{HD}\mathrm{\hspace{0.17em}100546})10^4\mathrm{cm}^3$$

$$\mathrm{and}n_{\mathrm{H}_2}(\mathrm{HD}\mathrm{\hspace{0.17em}163296})8\times 10^4\mathrm{cm}^3.$$

Assuming these volume densities, we can estimate the typical size of the absorbing layer responsible for the observed absorptions. Note that with low column densities and large volume densities, the absorbing cloud must be thin. The ratios of the H2 column densities to the volume densities give a typical thickness of $``$0.2 AU for HD 100546 and $``$1.2 AU for HD 163296.

All the observed and estimated quantities are consistent with absorption of the star light by a thin outer layer of a disk seen at an intermediate inclination. Because the disks of HD 100546 and HD 163296 are both inclined by about 60°, it is likely that we observed a line of sight grazing the disk surface. In any case, this shows that the H2 gas is not confined to a very flat disk. Alternatively, there is still a possibility that there might be shocked molecular gas in the vicinity of these stars, due to the combination of a bipolar jet (as observed around HD 163296) and a circumstellar envelope. However it is not possible to discriminate between these two possibilities with solely spectroscopic observation with no spatial information.

## 4 Conclusions

We observed the molecular hydrogen toward two Herbig stars surrounded by young protoplanetary disks. The observation of an extremely large number of rotational and vibrational levels allows for the first time the determination of the H2 physical conditions within circumstellar gas around two HAeBe stars. With this detailed view of the population of the H2 excited levels, we now need a full modeling in which all the cooling and heating terms must be taken into account. Such a modeling is under development. This opens a new window into the physical conditions operating in the protoplanetary disks surrounding massive stars.

The present observations of material surrounding stars with typical age between 106 and 107 years give new constrains on the relation between the radiation and the molecular protoplanetary disk in which planets are supposed to form and migrate (Terquem et al. Terquem2000 (2000)). This gives information on the interface between the dense part of the disks and the interstellar medium, and can help to better understand the late stages of the planetary formation. The radiation is indeed a key element in the evaporation of the protoplanetary disks whose lifetime seems to be limited (Zuckerman et al. Zuckerman1995 (1995)) although the process responsible for the clearing of the disks remains unclear.

###### Acknowledgements.

This work has been done using the profile fitting procedure developed by M. Lemoine and the FUSE French Team. The data were obtained for the Guaranteed Time Team by the NASA-CNES-CSA FUSE mission operated by the Johns Hopkins University. Financial support to French participants has been provided by CNES. Financial support to U. S. participants has been provided by NASA contract NAS5-32985. We warmly thank E. Roueff for providing H2 transition data in electronic format and for fruitful discussions. We thank A. Dutrey, G. Herczeg and J. Linsky for fruitful comments.