text:
# Production of intermediate-mass dileptons in relativistic heavy ion collisions

## I Introduction

The study of relativistic heavy ion collisions derives its unique importance from the opportunity of producing hot and dense strongly interacting matter. More specifically, one of the hopes is the creation and study of a quark-gluon plasma (QGP), a prediction of QCD , which filled the nascent Universe roughly a microsecond after the Big Bang. This field of endeavour currently represents a cutting edge of research in subatomic physics and has generated tremendous activity, both in theory and in experiment . Good penetrating probes of in-medium phenomena in heavy ion collisions are electromagnetic observables, owing to their small rescattering cross sections and to their privileged coupling to vector mesons . For collisions at ultrarelativistic energies, the lepton pair invariant mass (M) spectrum can roughly be divided in three regions: low mass (M $`<m_\varphi `$), high mass (M $`>m_{J/\psi }`$), and intermediate mass. The high mass region is Drell-Yan-dominated at the CERN SPS, and this contribution is in principle calculable within the realm of perturbative QCD. In connection with the low mass sector in heavy ion collisions at similar energies, the measurements of the CERES collaboration  have been seminal in building the case for an observation of in-medium effects on the vector meson spectral densities . There, the experimental collaboration has identified an excess over the sources that were known to be sufficient for the understanding of the spectra in pA collisions. This modification of the vector meson properties in matter represents an exciting discovery and is a genuine consequence of many-body physics.

In the intermediate invariant mass region, an excess of dimuons over the sources expected from pA measurements has also been identified experimentally by the Helios/3  and NA50  collaborations. Some excitement and interest have also been associated with this observation as this mass region is the one originally predicted to contain a plasma signal . Theoretical models have been devised to understand this dimuon excess. Those include rescattering of the open charm mesons with the surrounding hadronic constituents . This approach can not however reproduce the dimuon mass and transverse momentum spectra for central nuclear collisions . An overall charm production enhancement has also been invoked as a prospective interpretation . Whereas this scenario raises an interesting possibility, it has to be reconciled quantitatively with the observation of $`J/\psi `$ suppression. This issue can however be settled by a direct measurement . An earlier theoretical study of the Helios/3 measurements has attributed the excess to secondary hadronic reaction processes . This result has some consequences that are important to put on a firmer basis, by comparing with other newer data. In this strain, more recent NA50 results have also been interpreted by thermal models, with some quark-gluon content .

Our goal in this paper is to consider the NA50 intermediate mass dimuon results in the light of a hydrodynamic modeling of the nuclear dynamics, of a detailed analysis of the hadronic dilepton rates, and of a precise simulation of the detector acceptance and resolution. The issue of centrality dependence of the signal is also addressed. The next section briefly describes the hydrodynamic approach utilized, along with its reproduction of hadron spectra. The dilepton rate calculation is then described and the numerical modeling of the detector cuts and resolution is introduced. After the dilepton spectra in invariant mass and transverse momentum have been discussed, the centrality dependence is considered. The paper closes with a summary and a conclusion.

## II Hydrodynamic evolution and hadron spectra

It is well known that the final yield of the thermal dileptons has to be obtained by an integration over the space-time history of evolution of the interacting system. We shall assume that a thermally and chemically equilibrated quark-gluon plasma is produced in such collisions at a time $`\tau _0`$, which may be estimated from the condition of isentropic expansion of the plasma :

$$\frac{2\pi ^4}{45\zeta (3)}\frac{1}{A_T}\frac{dN}{dy}=4aT_0^3\tau _0$$
(1)
where $`dN/dy`$ is the particle rapidity density for the collision, and $`a=42.25\pi ^2/90`$ for a plasma of massless u, d, s quarks, and gluons. The number of flavours is taken to be $`N_f2.5`$ to account for the mass of strange quarks. We see that once the transverse area of the colliding system $`A_T`$ is known, along with $`dN/dy`$, the above relation uniquely relates $`T_0`$ to $`\tau _0`$. While one may get arbitrary large value of $`T_0`$ by choosing a smaller formation time of the plasma, the uncertainty principle provides a lower physical limit on $`\tau _01/3T_0`$ .

While considering Pb + Pb collisions at SPS energies, we have taken the average particle rapidity density as 750 for the 10% most central Pb + Pb collisions at the CERN SPS energy, as measured. We estimate the average number of participants for the corresponding range of impact parameters ($`0b\mathrm{\hspace{0.17em}4.5}`$ fm) as about 380, compared to the maximum of 416 for a head-on collision. We thus use a mass number of 190 to get the radius of the transverse area of the colliding system and neglect its deviations from azimuthal symmetry, for simplicity. As the deviation measured in terms of the number of participants is marginal ($`<`$ 9%), we expect the error involved also to be small. We also recall that the azimuthal flow is minimal for near central collisions.

The plasma is then assumed to undergo a boost-invariant longitudinal expansion and an azimuthally symmetric radial expansion, with a transition to a hot hadronic gas consisting of all hadrons having M$`<`$ 2.5 GeV, in a thermal an chemical equilibrium at a transition temperature $`T_c`$. This makes for a rich equation of state. Once all parton matter is converted into hadronic matter, we assume the hot hadronic system to continue to expand until it undergoes a freeze-out at some temperature $`T_F`$. During this evolution, the speed of sound of the matter is consistently calculated at every temperature to be used in the equation of state needed for solving the hydrodynamic equations .

It is important to realize that the (pressure) gradients in the system cause and control the radial expansion of the plasma, which in turn leads to a more rapid cooling and to the development of a transverse velocity for the fluid. Support for these ideas comes from intensity interferometry and from measured particle spectra. It is thereby necessary to have a proper initial energy density profile as it affects the hydrodynamic developments by introducing additional gradients. We assume it to follow the so-called “wounded-nucleon” distribution, which for central collision of identical nuclei leads to:

$$ϵ(\tau _0,r)_{\mathrm{}}^{\mathrm{}}\rho (\sqrt{r^2+z^2})𝑑z$$
(2)
where $`\rho `$ is the (Woods-Saxon) distribution of nucleons and $`r`$ is the transverse distance. This is prompted by the experimental observation that transverse energy deposited in these collisions scales with the number of participants. The normalization is then determined through numerical integration such that

$$A_Tϵ_0=\mathrm{\hspace{0.17em}2}\pi rϵ(r)𝑑r.$$
(3)
A proper profile is also needed to get a quantitative description of the dilepton yields and through this to infer the initial conditions. It is clear that the energy densities at larger $`r`$ would be smaller, implying lower temperatures and a shorter life-time till the fluid there undergoes freeze-out. Now recall that the four-volume $`d^4x=\tau d\tau rdrd\eta d\varphi `$ and thus the contribution of a given fluid element rises linearly with its transverse distance, which further “magnifies” the contribution from large $`r`$. One can immediately see that a uniform energy density profile often employed in the literature would lead to an erroneously large contribution from large $`r`$. This fact plays an important role in arriving at the required initial conditions.

We now have all the ingredients to solve the hydrodynamic equations. This is done through the procedures described in Refs. .

### A Initial conditions at the SPS

On the basis of the rapidity density etc. observed at SPS energies in central collisions involving lead nuclei, it was shown recently  that a good empirical description of the single photon data measured by the WA98 collaboration was obtained when the plasma formation time was taken to be $`\tau _0=`$ 0.2 fm/$`c`$. This corresponded to an average initial temperature of about 330 MeV. In the present work, we shall use these values but also explore the consequences of varying $`\tau _0`$ and $`T_0`$, keeping the resulting $`dN/dy`$ fixed.

The phase transition is assumed to take place at $`T=`$ 180 MeV and the freeze-out at 120 MeV. This value of the critical temperature is motivated by lattice QCD results which give values of about 170 – 190 MeV , and the thermal model analyses of hadronic ratios which suggest that the chemical freeze-out in such collisions takes place at about 170 MeV (recent analyses yield a value of $`158.1\pm \mathrm{\hspace{0.17em}3.2}`$ MeV  for the chemical freeze-out temperature). The phase transition should necessarily take place at a higher temperature.

### B Hadronic Spectra

As a first step we look at the spectra of several hadrons for the central collision of lead nuclei measured at the SPS.

In Fig. 1 we show our results when the formation time is take as $`\tau _0=`$ 0.2 fm/$`c`$ as in Ref. . The data are from Refs. . Even though it is known for quite some time that a reasonable variation of the initial temperature and the formation time, keeping the corresponding $`dN/dy`$ fixed, affects the flow only marginally, we study the consequences of increasing $`\tau _0`$ to 1 fm/$`c`$ in Fig. 2. We see that the hadronic spectra shapes are relatively insensitive to such variations in the initial conditions, as the flow takes some time to develop. In all cases, the agreement is quite good. This fact lends some credibility to the dynamics advocated here.

## III Dilepton Spectra at SPS energy

### A Dilepton emission rates

For the dilepton emission from the deconfined sector of QCD, we have assumed that the dominating process is $`q\overline{q}\mathrm{}^+\mathrm{}^{}`$. The emission rates that correspond to this process have been calculated many times. An example of their quantitative contribution can be found in Ref. . Even though those rates are the Born terms in the parton sector, a very recent lattice calculation has found that those were surprisingly accurate .

In the evaluation of hadronic rates most calculations have relied on effective Lagrangian techniques, where the available parameters are fitted to empirical measurements. An example of the application of this line of thought to the low mass dilepton sector can be found in Refs. . A problem which arises when one tries to extend those techniques to the intermediate mass region is the appearance of off-shell effects. Indeed, different approaches that agree in the soft sector can yield widely different results in higher invariant mass extrapolations . Fortunately, constraints on the hadronic processes can be obtained through the wealth of data of the type $`e^+e^{}`$ hadrons . In addition, those measurements cover exactly the same invariant mass range as the one which concerns this work. Those have been used, together with $`\tau `$-decay data, to construct the vector and axial-vector spectral densities . Similarly, the intermediate invariant mass initial state $`e^+e^{}`$ data has been analyzed specifically in the most important exclusive channels. This information can then be used to construct rates for hadrons $`e^+e^{}`$ . It is this procedure that is followed in the current work. The contributing channels producing lepton pairs in the invariant mass range 1 GeV $`<`$ M $`<`$ 3 GeV have been found to correspond to the initial states: $`\pi \pi `$, $`\pi \rho `$, $`\pi \omega `$, $`\eta \rho `$, $`\rho \rho `$, $`\pi a_1`$, $`K\overline{K}`$, $`K\overline{K^{}}+c.c.`$ . A detailed channel-by-channel discussion is too long to be had here, but a quantitative assessment of the net dilepton rate from a gas of mesons at T = 150 MeV is shown in Fig. 3. Also shown is the rate extracted from the spectral function evaluation . In principle, the latter contains all the hadronic sources. The correspondence between those two results is a confirmation that the important channels have indeed been identified in the hadronic scenario.

Using the rates of dilepton production from quark and hadronic matter, it is straightforward to calculate the spectra for dileptons, say for the central collision of lead nuclei at SPS energy. We add the contributions of the quark matter from the QGP phase and the quark matter part of the mixed phase and call it QM. Similarly we add the contributions of the hadronic matter part of the mixed phase and the hadronic phase and call it HM. The contribution of the Drell-Yan process is obtained by a scaling of nucleon-nucleon estimates to the case of nucleus-nucleus collisions. We get the results shown in Fig. 4, no acceptance or resolution corrections are implemented. It is clear that in principle, if the formation time is small (i.e., the initial temperature is large) then there could be a substantial contribution of the dileptons having their origin in the quark matter, in the intermediate mass window. Note that the HM contribution is only very marginally affected by variations in the $`\tau _0`$. We shall proceed to model the effect of detector resolution and to compare with experimental measurements.

### B Detector acceptance and dilepton spectra

It is vital to account for the finite acceptance of the detectors and for their resolution when comparing the results of theoretical calculations with measured experimental data. Those effects are important in the NA50 experiment . One approach to this problem is to model approximately and analytically the acceptance. While this can be readily implemented , a legitimate doubt can subsist about the accuracy of the experimental representation, especially in regions where edge effects might be important. In order to circumvent this problem we use here for the first time a numerical subroutine developed to reproduce the NA50 acceptance cuts and finite resolution effects in the measurements of dimuon pairs in Pb + Pb collisions at the CERN SPS .

We compute the invariant mass distribution of the dileptons in our hydrodynamic approach and then we run our pairs through the numerical detector simulation. The normalization is determined by a fit to the Drell-Yan data using the MRSA parton distribution functions as in the NA50 analysis. For getting the $`p_T`$ distribution we supplement our $`dN/dM^2`$ estimates for the Drell-Yan with a gaussian distribution in $`p_T`$ as in Ref.  which very closely reproduced the estimates obtained by the NA50 collaboration. The result is shown in Fig. 5. We also compute the $`p_T`$ distribution for the muon pairs, in the invariant mass window 1.5 $`<`$ M $`<`$ 2.5 GeV. This spectrum is broken down in the same sources as previously and is shown in Fig. 6.

At this point it is appropriate to consider the following question: which initial temperature is demanded by the intermediate mass dilepton data? A critical and quantitative assessment of this issue can be obtained by examining a linear plot of the lepton pair spectrum in the mass region under scrutiny. This is show in Fig. 7. From this figure, it appears that the best fit is provided by $`\tau _0`$ = 0.2 fm/c, and that the second best (less than two standard deviations away for most of the data points) belongs to $`\tau _0`$ = 0.4 fm/c. In terms of initial temperature those two values of the formation time correspond to $`T_0`$ 330 and 265 MeV respectively. A conservative and reasonable point of view is that it is probably not fair in such a challenging and complex environment as that of ultrarelativistic heavy ion collisions to ask for an agreement that is better than two standard deviations, considering the inherent uncertainties. This circumvents the issue of initial temperature determination. Also, there are questions that have to do with the specific dynamics used here to model the nuclear collisions. However in this respect, the hydrodynamic model used in this work is solved self-consistently and its parameters are constrained by the hadronic data as well. In connection with the important discussion on plasma signals the quark matter contribution is $``$ 23% for $`\tau _0`$ = 0.2 fm/c, and $``$ 19% for $`\tau _0`$ = 0.4 fm/c, around the lepton pair invariant mass of 1.5 GeV.

### C Centrality dependence

Our discussion so far was concerned only with the high multiplicity bin or, with mostly central collisions. To extend the hydrodynamic model to non-central events and to properly treat the azimuthal anisotropy is not a simple task. However, one can get a fair estimate of the centrality dependence by ignoring the broken azimuthal symmetry and by approximating the region of nuclear overlap by a circle of radius $`R1.2(N_{\mathrm{part}}/2)^{1/3}`$, where $`N_{\mathrm{part}}`$ is the number of participants . Assuming the uncertainty relationship, $`\tau _0=1/3T_0`$, one can track the changing multiplicities without introducing additional parameters. The results of this exercise are shown in Fig. 8, for S + U and Pb + Pb collisions. Also shown are the contributions from the “enhanced production of charm” estimated by the NA50 collaboration, normalized to the dileptons estimated by us. We estimate the contribution of the “excess charm” $`\mathrm{\Delta }N_{ch}`$ from

$$\mathrm{\Delta }N_{ch}(E1)N_{coll}$$
(4)
where $`N_{coll}`$ is the average number of collisions for the given centrality and $`E`$ is the ‘enhancement’ factor given by the NA50 measurement . The constant of proportionality is identical for all the bins for a given system. It is seen that the agreement with the measured data is quite good, and our approach gives a fair description of the centrality dependence of the excess dilepton measurement.

## IV Summary and conclusion

In the light of the findings reported on in this work, it does appear that the case of attributing the intermediate mass dilepton excess to an enhanced charm production rests on evidence that is no longer very compelling. The lepton pairs coming from a thermal hadronic medium comprise an important part of the signal, in agreement with earlier analyses of the Helios/3 measurements , and with those of more recent studies . It is satisfying to notice that the behaviour of the centrality dependence is also reproduced by the modeling and the dynamics used here. We have also considered the constraints on the initial temperature set by the data in the framework of our approach. There is a QGP component in the treatment used in this work. Unfortunately, its numerical value is not important enough to provide irrefutable evidence. The future of this field is bright: RHIC electromagnetic measurements should soon signal bold incursions into a new territory.

###### Acknowledgements.

It is a pleasure to acknowledge useful conversations with L. Capelli, O. Drapier, and L. Kluberg. This work was supported in part by the Natural Sciences and Engineering Research Council of Canada, in part by FCAR of the Quebec government, and in part by the US Department of Energy under grant DE-FG02-96ER40945..


format: text

text:
# 𝐵 Mesons on the Transverse Lattice

## 1 Introduction

Parton distributions measured in deep inelastic scattering, as well as many other high-energy observables (e.g. the deeply virtual Compton scattering amplitudes), are dominated by correlations along the light-cone ($`x^2=0`$). This simple fact poses a big obstacle for non-perturbative calculations of these important observables. For example, this makes direct evaluations of parton distributions on a Euclidean lattice, where all distances are space-like, impossible and calculations performed in a Euclidean framework usually try to reconstruct parton distribution functions from their moments. Furthermore, in an equal time quantization scheme, deep inelastic structure functions are described by real time response functions which are not only very difficult to interpret but also to calculate.

Light-Front (LF) quantization seems is a promising tool to describe the immense wealth of experimental information about structure functions for a variety of reasons:

* correlations along the light-cone become “static” observables in this approach \[i.e. equal $`x^+(x^0+x^3)/\sqrt{2}`$ observables\]
* structure functions are easy to evaluate from the LF wavefunctions
* structure functions are easily interpreted as LF momentum densities

Further advantages of the LF formalism derive from the simplified vacuum structure (nontrivial vacuum effects can only appear in zero-mode degrees of freedom) which provides a physical basis for the description of hadrons that stays close to intuition .

## 2 The Transverse Lattice

Before one can apply the LF formalism to QCD one has to remove the divergences first (i.e. regularize and renormalize). Then one has to cast bound state problems into a form that can be solved numerically with a reasonable effort.

The basic idea of the transverse lattice is very simple: one keeps two directions (the time and the z-direction) continuous but discretizes the transverse space coordinates (Fig. 1). The metric is Minkowskian. Two immediate advantages of this construction are

* manifest boost and translational invariance in the longitudinal direction — thus keeping parton distributions easily accessible
* a gauge invariant cutoff for divergences associated with large transverse momenta

Because of these features, the transverse lattice seems to be ideally suited for a light-cone formulation of QCD.

The gauge degrees of freedom on the transverse lattice are described by non-compact $`A`$-fields in the continuous longitudinal directions and by compact link fields $`U`$ in the transverse directions. The $`A`$ fields are defined on the sites of the lattice and the $`U`$s on the links. All degrees of freedom depend on two continuous and two discrete space-time variables.

For the canonical lattice approximation of $`G_{\mu \nu }G^{\mu \nu }`$ on the transverse lattice one needs to distinguish the following cases:

* Both $`\mu `$ and $`\nu `$ are longitudinal. Here the lattice representation for $`G_{\mu \nu }G^{\mu \nu }`$ is formally identical to the continuum representation.
* Both $`\mu `$ and $`\nu `$ are transverse. Here the lattice representation is just the plaquette interaction, which is familiar from Euclidean lattice gauge theory
* In the mixed case, i.e. for example when $`\mu `$ is longitudinal and $`\nu `$ is transverse, the lattice representation for $`G_{\mu \nu }G^{\mu \nu }`$ resembles the kinetic term of a gauged nonlinear sigma model

$`G_{\mu \nu }G^{\mu \nu }D_\mu U^{}D^\mu U.`$ (1)

The transverse lattice action was first introduced in Ref.  where it was already realized that a formulation with compact link-fields is not suitable for light-cone quantization. 111For a discussion of these technical difficulties, see for example Refs. . This has led to the color dielectric formulation of the transverse lattice where ‘macroscopic’ or smeared degrees of freedom are introduced to represent linearized link degrees of freedom on a coarse lattice. The effective action in the color dielectric formulation is obtained by making an ansatz which does not break the unbroken symmetries of the transverse lattice (e.g. gauge invariance, longitudinal boost invariance) and the coefficients in this ansatz are then obtained by seeking regions of enhanced Lorentz symmetry (e.g. where the static $`Q\overline{Q}`$ potential is ‘round’ and and where light glueballs have a dispersion relation with the same transverse speed of light). In Ref.  such an ansatz which included terms up to $`4^{th}`$ order in the link fields yielded glueball masses that are consistent with Euclidean results.

## 3 Spectrum and structure of light mesons

We used the effective link-field interactions, as determined in previous pure glue calculations within the color-dielectric formulation of $``$ lattice QCD . The fermion action was based on a $``$ lattice generalization of Wilson fermions . In this formulation, there are two hopping terms for fermions (hopping with and without spin flip for the quarks), which represent also the most general hopping terms that are possible to the same order in the fields which are consistent with the residual symmetries of the $``$ lattice. 222An alternative formulation, which represents a generalization of staggered fermions to the $``$ lattice, has been described in Ref. .

In addition to these hopping terms, one needs to introduce kinetic energy terms for the quarks as well as a coupling of the fermions to the longitudinal gauge degrees of freedom $`A^{}`$.

Because of gauge invariance, the quarks and $``$ link fields couple to the longitudinal gauge field with the same strength. Since gauge coupling for the $``$ link fields have already been determined in studies of glueball spectra as well the static $`Q\overline{Q}`$ potential, this coupling is no longer a free parameter and the only new parameters are the coefficients of the hopping terms as well as the (kinetic) quark masses. In the spirit of the color dielectric formulation, we determined these parameters by looking for regions in parameter space with enhanced Lorentz symmetry. The criteria that we used to test the violation of Lorentz symmetry were

* $`𝐤_{}`$ dependence in the dispersion relations of $`\pi `$ and $`\rho `$ mesons: using the relation $`𝐤_{}=a𝐩_{}`$ between the lattice and physical $``$ momentum, one can extract the effective $``$ lattice spacing $`a`$ for each hadron individually by performing a Taylor series expansion of its numerically determined dispersion relation

$`2p^+p_n^{}`$ $`=`$ $`m_n^2+c_n^2𝐤_{}^2+𝒪(𝐤_{}^4)`$ (2)
$`=`$ $`m_n^2+a_n^2c_n^2𝐩_{}^2+𝒪(𝐩_{}^4).`$
The covariant dispersion relation $`2p^+p^{}=m^2+𝐩_{}^2`$ is satisfied if the $``$ lattice spacing satisfies $`a_nc_n=1`$. The non-perturbative renormalization condition that we used in this work was to demand that the respective lattice spacings for $`\pi `$ and $`\rho `$ mesons are the same and also agree with $``$ lattice spacings determined within the pure glue sector.
* mass splitting within the $`\rho `$ meson spin multiplet

In a first principle calculation, it would be sufficient to input the $`m_\pi `$ to set all scales in the light quark sector. However, since the role of chiral symmetry and chiral symmetry breaking is still not very well understood within the LF framework, we needed to input the chiral symmetry breaking scale as one phenomenological parameter (on top of the string tension and $`m_\pi `$). We did this by using $`m_\rho `$ as an additional input parameter.

In the numerical calculations  we restricted ourselves to the $`N_C\mathrm{}`$ limit, which limited not only the number of possible terms in $`P_{eff}^{}`$ but also simplified the classification of states. For the numerical calculations we restricted the number of link fields to the minimal number which allows $``$ propagation of the entire hadron: $`q\overline{q}`$ on the same $``$ site as well as $`q\overline{q}`$ separated by one link with a link field $`U`$ connecting the two.

Typical results for the meson dispersion relations are displayed in Fig. 2. As one can read off from Fig. 2, it was not possible to restore full rotational invariance for the $`\rho `$ multiplet, but we hope that future calculations including additional terms in $`P_{eff}^{}`$ as well as higher Fock components can improve this situation.

In the continuum limit the dispersion relations should all be parabolas with the same curvature at the bottom. Near $`𝐤_{}=0`$ this is reasonably well satisfied, but of course there are larger violations of Lorentz invariance as the inverse momentum becomes comparable to the lattice spacing near the boundaries of the Brillouin zone. The level crossing is a remnant of species doubling since without $`r`$ term species doubling manifests itself on the transverse lattice at the hadronic level by giving rise to $`\rho `$ mesons at the boundary that are degenerate with $`\pi `$ mesons in the center of the Brillouin zone (and vice versa).

Results for the $`\pi `$ distribution amplitude are shown in Fig. 3. Although $`\varphi ^\pi (x)`$ resembles the asymptotic distribution $`\varphi _{asy}^\pi (x)=6x(1x)`$, our numerically determined result is somewhat broader, but clearly does not exhibit any ‘double hump’ feature. The $`\rho `$ meson distribution amplitude looks similar although it is slightly more peaked, which reflects the weaker binding of the quarks in the $`\rho `$.

For its normalization, i.e. for $`f_\pi `$ we found a result that is about a factor 2 larger than the experimental value. This discrepancy is most likely caused by the Fock space truncation, since including more higher Fock components tends to decrease the probability to find a hadron in its lowest Fock component and therefore also the normalization of the distribution amplitude.

## 4 B mesons on the $``$ lattice

In the limit where the $`b`$ quark is infinitely heavy, it acts as a static color source to which the light quark is bound. The extension of our light meson calculations to such a heavy-light system is straightforward and since the static source does not propagate, no new parameters appear in the Hamiltonian for such a system .

For the decay constant, which also plays an important role in mixing phenomenology, we find $`f_B240MeV\pm 20MeV`$, i.e. a value that is somewhat larger than those obtained using Euclidean lattice gauge theory or QCD sum rules, but here the discrepancy is much smaller than for $`f_\pi `$. This result is consistent with our observation that the Fock expansion also seems to converge much more rapidly for $`B`$ mesons.

For phenomenological applications  it is useful to have numerical estimates for the moments (normalization: $`_0^{\mathrm{}}𝑑z\varphi _{\mathrm{}}(z)𝑑z=1`$)333The error bars include only the estimated error from the numerical extrapolation to the heavy quark limit and the truncation of the Hilbert space, but not the systematic errors from the extrapolation in the Fock space and the transverse lattice spacing.

$`{\displaystyle _0^{\mathrm{}}}𝑑zz\varphi _{\mathrm{}}(z)`$ $`=`$ $`1.51\pm 0.1GeV`$
$`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{z}}\varphi _{\mathrm{}}(z)`$ $`=`$ $`1.22\pm 0.1(GeV)^1.`$ (3)
The numerical values for the moments indicate a larger momentum scale compared to other calculations. This result is consistent with a rather large value of the $`B`$-meson ‘binding energy’ $`\overline{\mathrm{\Lambda }}0.91.0GeV`$ obtained from the same transverse lattice eigenstates by calculating the expectation value of the $`p^+`$ momentum of all light degrees of freedom in the $`B`$ meson. We expect that including more Fock components will lead to a lowering of $`\overline{\mathrm{\Lambda }}`$.

An important observable in $`B`$-physics is the Isgur-Wise (IW) form factor, because of its use in the extraction of the CKM matrix element $`V_{bc}`$ from decays like $`B\overline{D}^{}l\nu `$. We work in the limit $`m_c`$, $`m_b\mathrm{}`$, where $`B^{}|\overline{b}\gamma ^\mu b|B`$, $`D^{}|\overline{c}\gamma ^\mu b|B`$ are all described by the same universal form factor

$`B^{}|\overline{b}\gamma ^\mu b|B=m_B\left(v_{}^{}{}_{}{}^{\mu }+v^\mu \right)F(vv^{})`$ (4)
For $`m_b,m_c\mathrm{\Lambda }_{QCD}`$, heavy quark pair creation is suppressed, i.e. relevant matrix element is diagonal in Fock space and an overlap representation exists for $`F(vv^{})`$

$`F(vv^{})=F^{(2)}(vv^{})+F^{(3)}(vv^{}),`$ (5)
where

$`F^{(2)}(x)`$ $`=`$ $`{\displaystyle \frac{2}{2x}}{\displaystyle \underset{s}{}}{\displaystyle _x^1}𝑑z\psi _s(z)\psi _s^{}\left({\displaystyle \frac{zx}{1x}}\right)`$ (6)
$`F^{(3)}(x)`$ $`=`$ $`{\displaystyle \frac{2}{2x}}{\displaystyle \frac{1}{\sqrt{1x}}}{\displaystyle \underset{s}{}}{\displaystyle _x^1}𝑑z`$
$`{\displaystyle _0^{1z}}𝑑w\psi _s(z,w)\psi _s^{}({\displaystyle \frac{zx}{1x}},{\displaystyle \frac{w}{1x}})`$
Here, $`\psi _s(z)`$ and $`\psi _s(z,w)`$ are the wave functions in the 2 and 3 particle Fock component and $`s`$ represents the spin/orientation labels.

Numerical results for the shape of the IW form factor, obtained from our numerically determined eigenstates on the $``$ lattice are consistent with experimental results.


text:
# Report on A5. Computer Methods

## 1 Introduction

In this Chapter, I shall describe the current status of computer simulation and algebraic computing for general relativity and gravitation. Most of the effort, both as represented in this Workshop and world-wide, is focused on “numerical relativity,” the name usually given to large scale simulations of isolated gravitational systems. Elsewhere in these proceedings is the contribution by Lehner based on his plenary talk on the status of the colliding black hole (CBH) problem. This problem is extremely difficult but urgently needs to be solved. The urgency arises because one desired output of CBH simulations is a set of gravitational wave forms that can be used by the ground based gravitational wave detectors (LIGO, VIRGO, GEO600, and TAMA) now coming into operation. The difficulties arise on many fronts. These include (1) the computational spacetime grid is not the physical spacetime grid; (2) CBH spacetimes are singular within the horizons of the black holes. Note that simulations cannot be followed (much) past the time at which NaN (not a number) or INF (infinity) is recorded as a value anywhere in the computational grid; (3) Setting up initial data is highly nontrivial. In the process of solving the constraints, the connection between the actual initial system and the desired physical initial state is lost. For example, the “no incoming radiation” boundary condition is problematic; (4) Existing 3D CBH codes which have made considerable progress in dealing with (1) and (2) are plagued by instabilities which cause them to crash long before any physically relevant results are obtained. Note that a code that runs long enough can overcome the problem of radiation in the initial data by evolving until after it dissipates; (5) Identification of outgoing gravitational radiation is also nontrivial because radiation is a small effect thus requiring accurate simulations for a credible extraction and because many computational grids must terminate before the radiation zone is reached. All of these issues are under intensive investigation and were discussed in this Workshop. Similar techniques have been developed to study the coalescence of neutron star binaries. Such simulations do not have to face the issues of singularity and horizon but do have the complication of matter.

The travails of the CBH simulations should not be taken to mean that all computer simulations in general relativity and gravitation have not yet been completely successful. In fact, rather the opposite is true. Most spectacular, of course, was the discovery by Choptuik of critical phenomena in the collapse of a spherically symmetric self-gravitating scalar field. Careful numerical techniques (including adaptive mesh refinement) allowed study of the critical system dividing scalar field configurations which eventually disperse from those which eventually form BHs. Scaling typical of that found in other phase transitions was observed. The critical solution itself (a zero mass naked singularity) displays periodic self-similarity — a phenomenon unknown outside gravitating systems. Choptuik’s original discovery sparked a large number of further and ongoing numerical and analytic studies of critical phenomena in gravitational collapse. Some of these were discussed in this Workshop.

Another area where there has been some success has been the numerical investigation of cosmological singularities. While an actual singularity cannot be handled on the computer, the regular behavior in the approach to it may be studied. If, asymptotically as the singularity is approached, some types of terms in Einstein’s equations dominate over others, it might be possible to characterize the approach to the singularity in a simple way. This was argued heuristically long ago by Belinskii, Lifshitz, and Khalatnikov (BKL) with their predictions of local Mixmaster dynamics (LMD) recently supported in the cosmological context with numerical experiment and, in some cases, mathematical proof. In fact, the synergy between the mathematical and numerical approaches in this area has been useful in understanding the relevant phenomenology.

Originally, this Workshop was designed to feature algebraic computing as well as computer simulation. However, the co-organizer, Jim Skea, was unable to attend. In the end, there was only a single talk on this subject where considerable progress was described toward the goal of automated classification of exact or perturbative solutions of Einstein’s equations.

In the rest of this Chapter, the highlights of the talks in each area will be given. For more details please see the indicated references.

## 2 Colliding Black Holes and Neutron Stars

Lousto presented an overview of the Lazarus Project to combine 3D simulations and perturbation approaches in an ambitious attempt to follow a BH collision and extract waveforms. This required the development of techniques to convert the coordinates suitable to the 3D simulation to those appropriate for the single black hole close-limit approximation. These simulations were checked by interchanging the use of 3D simulation and perturbation where both were valid. In addition, the method was applied to a single Kerr BH where computed spurious radiation indicates the degree of error. This approach has made it possible to compute the first complete waveforms covering the post-orbital dynamics of a binary black hole system.

A number of talks were given on various techniques needed for standard (i.e. $`3+1`$ formulation) 3D codes. It is currently believed that an important source of instability in 3D codes is constraint violating modes. Einstein’s equations may be written in an infinite number of ways by adding arbitrary multiples of the constraints to the dynamical equations of motion. The goal is to identify a formulation which will drive the system to the constraint hypersurface and thus eliminate these instabilities. One such approach was discussed by Shinkai.

Pollney discussed work (with AEI collaborators) on BH excision and gauge techniques. They describe recent progress in an implementation of excision for 3D black holes along with a set of gauge conditions which respond naturally to spacetime dynamics. Through the combined use of these techniques, they are able to produce accurate and long-lived evolutions of highly distorted, rotating black hole spacetimes. Accurate waveforms can be extracted.

Diener presented the results of the head-on CBH code test for 3D simulations. Koppitz reported on preliminary work in the construction of initial data sets. The objective is to consider thin sandwich, Kerr-Schild, and post-Newtonian approximation initial data using various implementations.

Bondarescu described how one might visualize BH horizons as they evolve in CBH simulations. The 3D codes allow one to locate and track the evolution of apparent and event horizons in the coordinates of the simulation. One typically has little information about their real geometry. Previous studies in axisymmetric spacetimes have visualized horizons via an embedding in flat space. A new method which goes beyond axisymmetry was discussed. The method correctly reproduces known results and, for the first time, has allowed construction of embeddings of non-axisymmetric, distorted black holes.

Ashtekar described a method to extract physics from strong field regimes in numerical simulations using the isolated horizon (IH) framework. If the IH (but not the spacetime) is axisymmetric, e.g., then, knowing just the intrinsic metric and one Newman-Penrose spin coefficient, the IH framework enables one to calculate the angular momentum and mass of the final black hole. These quantities refer only to the black hole in equilibrium and do not include the angular momentum and mass in matter fields or radiation outside the IH and thus could be found numerically.

Grandclément described a new approach to binary black hole evolution prior to the plunge phase. One important feature of this phase of the evolution is the location of the innermost stable circular orbit (ISCO). This has been calculated both using post-Newtonian approximations and numerically from quasi-static equilibria with different results. If two orbiting bodies are far apart, the time-scale associated with the gravitational radiation is much longer than the orbital period allowing one to assume that the two black holes are on exact circular orbits. Using multi-domain spectral methods, a sequence of two identical corotating black holes is computed. The ISCO appears as a turning point in the total energy and angular momentum curves. Its position agrees well (for the first time) with the post-Newtonian values.

Several talks were given on characteristic codes. For a review of this approach, see \[\]. Bishop described the initial value problem for a 3D characteristic code. The code is robust and stable and has been extended to include matter, although in a simple way that would not handle shocks. The code cannot evolve a whole spacetime, but requires data on an inner world tube. For this reason, it is important to identify initial data of the required type to represent a physical situation with strong gravitational fields.

D’Inverno and Pollney discussed lower dimensional code tests for a Cauchy characteristic matching (CCM) code. See (e.g.) \[\]. The advantages of such a code are that it dispenses with an outer boundary condition and (since it uses a compactified coordinate) can yield global solutions in which gravitational waves can be identified unambiguously at future null infinity. Recently, a master vacuum axisymmetric CCM code has been completed. The main motivation for this work is to construct a three dimensional code possessing the characteristic, injection and extraction modules (now present in the 2D code) which can be attached to existing interior codes based on a finite grid.

A completely different approach involves a 3D code based on the conformal field equations developed by Friedrich. This work was presented by Husa. He discussed the evolution of linear and (mildly) nonlinear gravitational waves. This approach aims at the numerical computation of the global structure of generic isolated systems in GR using well defined and unambiguous algorithms even out to and beyond future null infinity. The basic ideas and algorithms were developed by Hübner and Weaver in recent years and implemented by Hübner. Hübner obtained the complete future of (the physical part of) the initial slice, illustrating a theorem by Friedrich, which states that for sufficiently weak initial data a regular point $`i^+`$ exists. However Husa found that for higher amplitudes the gauge chosen in by Hübner results in code crashes which can be cured by some ad-hoc modification of the lapse. For still stronger (but still only mildly nonlinear) data, evolution up to $`i^+`$ would require the tuning of too many parameters to compensate for the lack of symmetry. Evolutions of axially symmetric initial data modeled in the spirit of “Brill waves”, and evolutions of Minkowski space in a static gauge were discussed. While the former evolution is numerically stable and quite robust, the latter exhibit instabilities which seem to be rooted in exponentially growing constraint violating modes inherent in the analytical formulation.

Siebel presented the results of a study using the characteristic approach to identify gravitational radiation. The implementation follows the characteristic initial value formulation of general relativity based on Bondi’s radiative metric. The study focuses on the evolution of neutron stars modeled by a $`N=1`$ polytrope with fully relativistic hydrodynamics. The code was used to study spherical relativistic stars with a scalar field to model the gravitational degrees of freedom. Depending on the specific neutron star model, the scalar field either induces oscillations of the star or makes it collapse to a black hole. The extracted frequencies of the oscillations agree very well with linear studies, except, as expected, at the threshold of black hole formation. The code was also applied to evolutions of neutron stars in axisymmetry and to their interaction with gravitational waves.

Stark (with Lun) considered an alternative $`3+1`$ scheme as applied to spherical perfect fluid collapse. The scheme utilizes the constraint equations contained in the Einstein field equations with those in the Bianchi identities to determine the 3-curvature tensor and the 3-covariant derivatives of the extrinsic curvature of the hypersurface. This yields a hyperbolic system of equations for the evolution of the gravitational potentials and hydrodynamical quantities. A numerical implementation of this system was used to model the Oppenheimer-Volkoff solution and the gravitational collapse of an initially homogeneous dust sphere to a BH.

Moreschi discussed modeling the collision of black holes with angular momentum with Robinson-Trautman (RT) spacetimes. These spacetimes contain purely outgoing gravitational radiation and decay to leave a Schwarzschild-like horizon. These RT spacetimes are natural candidates to study systems settling to a single black hole and may be used as background spacetimes for a perturbation treatment. Previous studies of this type give excellent agreement with numerical relativity results for the head-on collision, zero spin case. They are then extended to include angular momentum.

Reula presented results using the collapse of spherically symmetric self-gravitating scalar fields as a model problem to study numerical instabilities. There are instabilities present in this simple problem, and furthermore there are strong indications that, more generally, the instabilities are predominantly due to a longitudinal or Newtonian mode, namely the only one present in the spherically symmetric case. Free and constrained evolutions are compared. In \[\] a symmetric hyperbolic formulation is used, and certain freedom, still available in that setting, is used to suppress the main instability found there, to allow a stable propagation. However, their discretization scheme has first order numerical dissipation and the boundary-value problem for the system is not known to be well posed. Reula discussed an improved system that overcomes these problems, preserves the constraints, and reproduces several known results.

Gentle applied Regge calculus methods to the evolution of Brill waves. This is the first significant test of Regge calculus in a highly dynamic setting. An axisymmetric lattice is obtain by triangulating and collapsing a hypercubic grid which is aligned with a cylindrical polar coordinate system (see \[\] and references therein).

Neutron stars might act as a source of gravitational waves. Sperhake presented a new numerical approach to nonlinear oscillations in neutron stars. This approach rewrites the equations to cancel out the background terms in the numerical evolution of the perturbations. A comparison of the resulting perturbative code with a “standard” non-perturbative method was made. The perturbative scheme reproduces the expected harmonic oscillations with high accuracy, while the non-perturbative scheme produced spurious behavior. The new method is then used to accurately evolve initial data of arbitrary amplitude to investigate non-linear coupling of eigenmodes for a simplified neutron star model.

Hawley (with Choptuik) explore the genericity of initial data for multi-scalar stars. They consider a class of general relativistic soliton-like solutions composed of multiple minimally coupled, massive, real scalar fields without self-interaction.

## 3 Numerical Investigation of Singularities

For reviews of critical phenomena see \[\]. Gundlach described what happens to the critical behavior of the gravitational collapse of a perfect fluid when there are departures from spherical symmetry. Does the type II critical behavior in the gravitational collapse of a $`p=k\rho `$ perfect fluid persist beyond the restriction to spherical symmetry? It was found that for $`1/9<k<0.49`$, all non-spherical perturbations decay. For $`0.49<k<1`$, there are one or several growing modes in the polar perturbations. For $`0<k<1/9`$, there is exactly one non-spherical growing mode. It is an $`l=1`$ axial mode, linked to infinitesimal rotation. BH mass and BH angular momentum become a function of 2-parameter families of initial data, chosen so that one of the parameters is linked to rotation (and thus can be tuned to zero by axisymmetry). At the BH threshold they depend on the initial data through one universal function of one variable,a combination of these two parameters. The formula obtained is completely analogous to the dependence of the magnetization of a ferromagnet as a function of temperature and external magnetic field near its critical point.

Aichelburg presented results (including those of Lechner) on a new transition between discrete and continuous self-similarity in collapsing $`SU(2)`$ $`\sigma `$-models parametrized by a dimensionless coupling constant. In the intermediate range, a competition between CSS and DSS solutions gives rise to new phenomena in the dynamical evolution of a critical search: the appearance of episodes of approximate CSS, where repeatedly the evolution approaches and departs from CSS before leading to black hole formation or dispersion. This picture is supported by the explicit numerical construction of the CSS and DSS solutions and its comparison with critical time evolutions.

Szpak studied critical behavior in a model nonlinear wave problem. The motivation was to study the existence of an universal intermediate attractor in the dynamics at the threshold of forming singularity.

Hobill (with Webster) presented studies of trapped surface formation in Brill wave evolution to see if a naked singularity could form. The oblateness or prolateness of the system is determined by measuring the ratio of polar circumference to equatorial circumference for a surface with constant logarithmic radial coordinate. It is found that for initial data with positive amplitude the systems are prolate and for negative amplitude the systems are oblate. It is seen that, as the Brill wave amplitude is increased (thereby decreasing the radius of the outermost minimal surface), the location of the outermost trapped surface eventually coincides with the minimal surface thereby preventing the formation of a naked singularity. This behavior is in agreement with recent results of Garfinkle and Duncan.

I presented a summary of our (with D. Garfinkle, J. Isenberg, V. Moncrief, M. Weaver) program of investigation of collapsing spatially inhomogeneous cosmologies. Significant recent progress has been made in rigorous studies of generic collapse in cases where one expects the approach to the singularity to be asymptotically velocity term dominated. However, one expects that an even larger class of spacetimes exhibits local Mixmaster dynamics (LMD) in the vicinity of the singularity. While numerical evidence for LMD in spatially inhomogeneous spacetimes with $`T^2`$ and $`U(1)`$ symmetries is compelling, rigorous results for spatially inhomogeneous cosmological spacetimes with LMD do not exist. New mathematical methods to handle LMD are needed.

## 4 Algebraic Computing: Invariant Classification in Maple

D’Inverno provided a summary of the program to classify solutions to Einstein’s equations automatically. Recent developments are in collaboration with Pollney, Skea, Araújo, de Albuquerque, and Roveda. A short review of the equivalence problem in general relativity was presented which included a description of the computer database of exact solutions whose central site is located at http://www.jim.dft.uerj.br. The database contains over 200 solutions of Einstein’s equations which have been classified with CLASSI, an extension of the general relativity system called SHEEP. In the next stage of this collaborative project attention has been focused on setting up a database which can easily be accessed and updated by the user community. The choice of platform is the general relativity package GRTensor which is an applications package of the computer algebra system Maple. The underlying algorithm is based on the Cartan-Karlhede invariant classification of geometries. In particular, new algorithms have been obtained for putting the Weyl spinor, Ricci spinor and general spinors into standard form and the derivative operators needed in the classification algorithm and the behavior of symmetric spinors under frame rotations have been investigated. More recent work was reported on a refinement of the JMS (Joly-MacCallum-Seixas) algorithm for completely determining the Segre type of a Ricci spinor. To date 30 of the 45 possible sub-cases have been completed and the algorithms have been implemented in GRTensor. In addition, the status of the work on obtaining the isometry group of a space-time (rather than simply its dimension as in the current database) was reported. It is currently possible to obtain the isometry group for spacetimes admitting a G2, G3 and some G4 symmetries. This may require the user to input limits on parameters and coordinates appearing in the line element, to overcome an ambiguity in determining the signature of a matrix occurring in the Bianchi classification. The goal is to complete the G4 case and extend the work to the G5, G6 and G7 cases. Work was also reported on the Maple package SPIDOR, which has been developed for obtaining spinor contractions. An example of its application is its ability to automatically generate the spinor curvature invariants built out of the Weyl spinor, its complex conjugate and the Ricci spinor up to the 20th order on a laptop. This package may be useful in relation to formalisms which have been developed to address particular issues in exact solutions.