text:
# Computational Aspects of Lattice QCD

## 1 Introduction

Quantum chromodynamics (QCD) is the theory describing the strong interactions, which occur between hadrons (e.g. protons and neutrons) . The description is based on a model of elementary particles — the quarks — possessing “color charge” and interacting through the exchange of gauge fields — the gluons. QCD is a quantum field theory, with local $`SU(3)`$ gauge symmetry, corresponding to three possible colors. The theory is written in simple and elegant form. Its only parameters are the masses of the various types (called “flavors”) of quarks considered and the value of the strong coupling constant. Except for the symmetry under the $`SU(3)`$ gauge group \[instead of the $`U(1)`$ group\], the form of the QCD Lagrangian is the same as the one of quantum electrodynamics (QED), with the quarks corresponding to the electrons and the gluons to the photons. (The former two are spin-$`1/2`$ fermions and the latter two are massless vector bosons.) Analogously, the strong coupling constant $`\alpha _s`$ corresponds to the fine structure constant $`\alpha 1/137`$. The fact that the gauge group of QCD is non-Abelian introduces, however, qualitative differences between the two theories, reflecting the differences between the strong interactions and the electromagnetic interactions. In particular, one obtains that the gluons possess color charge and therefore interact with each other, as opposed to the photons.

An important characteristic of the strong interaction is that the coupling constant $`\alpha _s`$ becomes negligible only in the limit of small distances, or equivalently in the limit of high energy or momentum. This property is called asymptotic freedom. At larger distances (i.e. smaller energies) there is an increase in the intensity of the interaction and it is believed that at large distances the force of attraction between quarks is constant, determining the confinement of quarks and gluons inside the hadrons. The fact that $`\alpha _s`$ is not negligible at low energies makes the study of important phenomena such as the mechanism of quark confinement, the hadron mass spectrum and the deconfining transition at finite temperature inaccessible to calculations using perturbation theory, which is based on a weak-coupling expansion. These phenomena must therefore be studied in a nonperturbative way.

The nonperturbative study of QCD is possible in the lattice formulation of the theory . In this formulation — which consists in quantization by means of path integrals, in the continuation to imaginary or Euclidean time and in the lattice regularization (given by the discretization of space-time) — the theory becomes equivalent to a model in classical statistical mechanics. The continuum limit, in which physical results are obtained, is given by the critical point of this model, which may be studied through the usual methods of statistical mechanics. In particular, one may perform numerical simulations by Monte Carlo methods, which are based on a stochastic description of the systems considered . Due to the greater complexity of the interaction and to the larger number of degrees of freedom, these simulations are much more elaborate for QCD than for the usual models in statistical mechanics, requiring considerable computational resources. In fact, one usually needs to simulate on powerful parallel supercomputers, some of which were designed and built specifically for the study of lattice QCD, such as the QCDSP in the USA, the Hitachi/CP-PACS in Japan and the APE-Mille in Europe, all with performance in the teraflops range. Only recently has simulation of QCD on systems of small computers, the so-called PC clusters, become possible. These systems do not yet provide the same efficiency in parallelization as the machines with parallel architecture, but their cost is much lower. In addition to the computational power, the numerical and analytical techniques used in the simulations and in the interpretation of the produced data are of great importance in the field. Significant progress has been achieved through the development of more efficient simulation algorithms, new methods for interpolation and extrapolation of the numerical data and a better understanding of the systematic effects to which the simulation may be subject, such as finite-volume effects and discretization errors.

There is presently great interest in the results of the simulations described above and one hopes to be able to solve many theoretical questions about QCD and the standard model . In fact, despite the great computational difficulty, numerical studies of QCD have provided important contributions recently, such as accurate calculations of the strong coupling constant  and of the hadronic mass spectrum . In particular, lattice simulations constitute the only known evidence for the quark-deconfining transition at finite temperature  and its predictions are of direct interest for the current experiments in search of new states of matter in the laboratories Brookhaven and CERN.

Simulations of the so-called full QCD — i.e. including effects of dynamical fermions — for quark masses in the region of physical values are still extremely slow. They must in general be carried out on supercomputers as the ones mentioned above and involving the effort of large collaborations, such as the UKQCD in the United Kingdom and the JLQCD in Japan. The methods used in these simulations, which take on average several months or even a few years, are often developed in studies of simplified versions of the theory, such as pure QCD — the so-called quenched approximation, in which effects of dynamical fermions are neglected — and the pure $`SU(2)`$ theory, or models in lower dimensions. The consideration of this type of problems and the use of PC clusters is the goal of our research group in the IFSC–USP.

## 2 The lattice formulation

A difficulty in the study of QCD, common to virtually all quantum field theories, is the appearance of ultraviolet divergences (i.e. divergences for high energies or short distances) in the calculation of physical quantities . Only after these “infinities” are removed by means of some renormalization procedure can one obtain finite results, which may be compared to experiment. It is therefore necessary to regularize the theory first, writing it in such a way that the singularities are isolated, and then to remove these singularities through a redefinition of the parameters in the Lagrangian. The lattice QCD formulation, introduced in 1974 by Wilson , offers a convenient nonperturbative regularization, preserving the theory’s gauge invariance. Quarks are represented at lattice sites, while gluons are represented on the links between neighboring sites. The gluonic fields are given by $`SU(3)`$ matrices. The lattice action is written in terms of products of link variables along closed loops, so that the gauge symmetry of the original action is preserved. An excellent introduction to lattice QCD is Ref. . The essential ingredients for the lattice formulation are:

1. Feynman’s path integral formalism, in which expectation values are written for the observables of interest as integrals over all the degrees of freedom of the problem, with a statistical weight given by the exponential of the theory’s classical action.
2. The Euclidean formulation, obtained by analytic continuation of the time variable to imaginary times. In this way the (complex) oscillatory exponential present in the integrals described above becomes real and may be interpreted as a probability distribution.
3. The introduction of a discrete lattice for the space-time. Correspondingly, differential operators are rewritten as finite differences of the discretized fields.

The combination of the first two ingredients highlights the equivalence between quantum field theories and classical statistical mechanics: in Euclidean space a path integral for the quantum theory is equivalent to a thermodynamic average for the corresponding statistical mechanical system. For QCD, the square of the bare coupling constant $`g_0`$ of the field theory corresponds directly to the temperature $`1/\beta `$ of the statistical mechanical model.

The third ingredient — the lattice discretization — represents an ultraviolet regularization. In fact, the lattice spacing $`a`$ corresponds to a high-momentum cutoff, since momenta higher than $`1/a`$ cannot be represented on the lattice. In this way the modes causing divergences are suppressed and the theory is well defined. In order to recover the continuum-space theory we must take the limit $`a0`$. In this process it is necessary to “tune” the bare parameters of the theory — e.g. the bare coupling constant $`g_0`$ — in such a way that physical quantities (correlation functions, masses, etc.) converge to finite values, which can then be compared to experiment. In particular, a correlation length $`\xi `$ (corresponding to an inverse mass) measured in physical units — e.g. fermi — must approach a finite limit when the lattice spacing $`a`$ (measured in fermi) goes to zero. This means that the correlation length measured in units of the lattice spacing $`\xi /a`$ must go to infinity. In other words, the lattice theory considered must approach a critical point, i.e. a second order phase transition. Thus, the study of the continuum limit in quantum field theories on the lattice is analogous to the study of critical phenomena in statistical mechanics. The correspondence between Euclidean field theories and classical statistical mechanics allows the application of usual statistical-physics methods to the study of QCD. One may use, for example, high- and low-temperature expansions, corresponding respectively to strong- and weak-coupling expansions for the field theory. Another example of interplay between field theories and statistical mechanics is the renormalization-group method, developed in parallel for both fields . A particularly important statistical technique, especially for QCD, is Monte Carlo simulation, which allows a nonperturbative study of the models considered.

## 3 Monte Carlo simulations

Monte Carlo methods are generally used to sample the Boltzmann distribution for a statistical system in a stochastic way . One generates on the computer $`N`$ configurations for the system, in such a way that each configuration is chosen with a probability given by its Boltzmann weight. The average of an observable over the $`N`$ configurations produced converges in the $`N\mathrm{}`$ limit to the expectation value (or thermodynamic average) of this observable. For large (finite) values of $`N`$ we get central values for averages of the observables of interest, with a statistical error proportional to $`1/\sqrt{N}`$. It is therefore possible to obtain an arbitrarily small error by increasing the number of configurations produced, i.e. by increasing the computational effort. The Boltzmann distribution may be defined through a model for a physical system, especially if this system can be discretized, either naturally or by an approximation. In many examples the art of simulation is in the appropriate choice (or invention) of a model, for example in the modeling of polymers by random walks. In the case of lattice gauge theories, as explained in the previous section, the Boltzmann distribution is obtained directly from the Lagrangian (or equivalently the action) of the theory, without any approximations other than the discretization of space-time, allowing therefore a nonperturbative study from first principles.

In order to generate configurations with the desired statistical weight, one usually introduces a Markov dynamics for the system considered, in such a way that the resulting Markov chain has the appropriate Boltzmann distribution as its equilibrium probability distribution. In this approach the thermodynamic averages described above are calculated as time averages in the chosen dynamics. At each “instant of time” we generate a new configuration of the system, in a manner that respects the appropriate (equilibrium) distribution. Except for this restriction, the updates that determine the dynamics may be chosen in the most convenient way, without necessarily coinciding with the physical dynamic behavior of the system out of equilibrium. We start from a general initial condition and follow the time sequence by successively applying the updating procedure, which generates a new configuration starting from the present one. It is common therefore to think of the simulated system as evolving by itself and one frequently says that the observables are being “measured” instead of calculated. There is also a statistical error associated with the stochastic method, as mentioned above, and one must apply error analysis to determine the final precision of the results. These are characteristic features of experimental studies and the methods used for data analysis are often the same. We must remember however that we are dealing with numerical “experiments”, obtained from a theory (in the case of QCD) or from a model for a physical system (in the case of statistical mechanics).

The simulation of lattice gauge theories, QCD in particular, constitutes one of the most intensive fields of application of Monte Carlo simulations \[3, Cap. 11\]. As said in the Introduction, this type of study is crucial for QCD, since a perturbative treatment of the theory is not possible for many relevant energy regions. Despite the similarity of the methods, the Monte Carlo simulation of gauge theories is much more complex than in the case of the usual statistical mechanical models, requiring great computational effort and specific numerical techniques for the production of the data. Moreover, the physical interpretation of the generated data depends on a correct extrapolation to the continuum limit, i.e. it is necessary to “go back” to the continuum space after the lattice simulation. More specifically, we must consider three limits in order to obtain the desired physical results from the simulation data:

* The infinite-volume limit (or thermodynamic limit): Just as in statistical mechanics, simulations of QCD are carried out for finite lattice volumes, since the computer’s memory is finite. The lattice volume must therefore be sufficiently large with respect to the physical distance that is relevant to the problem at hand, so that finite-volume effects are not significant. (Correspondingly, one may not consider energies or momenta that are too low, since the finite lattice is equivalent to an infrared cutoff.) Finite-volume effects can in general be estimated through a finite-size-scaling analysis of the data. For QCD simulations it generally would suffice to consider lattice sizes $`L7fm`$.
* The continuum limit: In order to recover the original continuum physics it is necessary to take the limit $`a0`$. At the same time, the calculated quantities must be renormalized, i.e. redefined so as to generate finite results in the continuum limit. This can be done nonperturbatively by using the physical values of some observables, which are known experimentally. For example, writing the pion mass calculated on the lattice as $`m_\pi a`$ (where $`m_\pi `$ is the physical mass), we obtain the value of the ultraviolet cutoff $`a`$ in physical units. The other calculated quantities may then be written in the same form in terms of $`a`$ (which tends to zero) and “translated” into physical units, generating (physical) finite values. In practice, the value of $`a`$ must be sufficiently small when compared to the relevant distance for the problem, for example $`a0.05fm`$. It is important to notice that several discretizations of the action are possible and that one may consider the so-called improved actions, which converge to the continuum limit faster, i.e. for larger values of $`a`$. (One may also consider several different discretizations for the fermionic fields.)
* The chiral limit: It is very hard to consider physical values for the quark masses in the numerical simulations. That happens especially for the light quarks (up and down), whose masses are close to zero, the so-called chiral limit. The simulations are usually done for larger masses and the results are then extrapolated using chiral perturbation theory.

The above limits are not independent, since to get to the continuum limit and to be able to consider small masses for the quarks one needs a sufficiently large number of lattice points (corresponding to a small enough lattice spacing and to a large enough physical size of the lattice), which increases considerably the computational effort. For example, with the present algorithms, a numerical simulation at the ideal values of $`L`$ and $`a`$ given above and for physical values of the quark masses would take approximately  300 years on a super-computer with $`1Tflops`$ power, corresponding to the fastest computers available today.111One teraflops ($`Tflops`$) is equivalent to $`10^{12}`$ floating-point operations per second. Employing improved actions and making use of the chiral-limit extrapolation, the same simulation can be done in 2 years. These simulations are much faster for the so-called quenched case, in which the configurations are produced considering quarks of infinite mass, i.e. without considering effects of dynamical quarks. (Note that the observables calculated for each configuration may still include quarks with the desired masses, which are then called valence quarks.) Despite being a rough (and uncontrolled) approximation, one verifies that in many cases the quenched approximation shows small corrections with respect to the full QCD, indicating that in these cases the effect of dynamical fermions is small. Today, quenched simulations can be done with good precision and there are also several ongoing simulations of QCD with dynamical quarks, for several discretizations of the Dirac operator.

## 4 Progresses in lattice QCD

Progresses in the field are reported annually at the Lattice conference . We describe briefly below three important contributions of lattice QCD to the confirmation/prediction of experimental results.

* Strong coupling constant: The strong coupling constant $`\alpha _s(\mu _0)`$, taken at a fixed reference scale $`\mu _0`$, is the only free parameter of QCD and must therefore be known to the highest precision possible. Numerical simulations of QCD are now able to produce calculations of $`\alpha _s`$ with precision comparable to the experimental one or better. These results are presently included in the world average for this quantity . The value of $`\alpha _s`$ should also be determined with good precision over as large a range of values as possible, describing the behavior of the interaction from the nonperturbative (strong-coupling) to the perturbative (weak-coupling) regime. Several methods for the calculation of the running strong coupling are being pursued, both in the quenched approximation and in the full-QCD case. A recent review of these methods can be found in .
* Hadron spectrum: It is possible to obtain physical values of hadron masses from lattice QCD simulations, as described in the previous section. To this end one must tune the $`n_f+1`$ parameters of the theory, where $`n_f`$ is the number of quark flavors considered. Thus, one must use $`n_f+1`$ known experimental results as inputs and the subsequent calculations are physical predictions of the numerical simulation. The mass spectrum of the light hadrons (including the two light quarks and the strange quark) has been determined for the quenched case with great precision in . One does not obtain complete agreement with the experimental spectrum, but the observed discrepancies are of at most 10%. We conclude then that the difference introduced by the quenched approximation is only quantitative. Similar calculations are now being performed for the full-QCD case.
* QCD phase transition: The predicted phase transition for QCD at high temperatures is clearly observed in lattice QCD simulations . For the quenched case one studies the deconfining transition itself, while for the full-QCD case one must consider the transition associated with the restoration of the chiral symmetry. This is an exact symmetry of the QCD Lagrangian in the limit of zero masses for the quarks and is spontaneously broken at low temperatures. Studies of pure QCD (i.e. the quenched approximation) are done with high precision for the determination of the critical temperature and of the thermodynamic equation of state for the system. In this case one encounters qualitative differences between pure QCD and the full-QCD case. (In particular, the order of the phase transition is different for the two cases.)

## 5 Lattice QCD at the IFSC–USP

Since the beginning of 2001 we have been carrying out a project on numerical simulations of lattice gauge theories at the Physics Department of the University of São Paulo in São Carlos (IFSC–USP), funded by FAPESP . The project included the installation of a PC cluster with 20 processing nodes (16 nodes with 866 MHz Pentium III CPU, 4 nodes with 1.7 GHz Pentium 4 CPU, all with 256 MB of memory).222The resulting computer power is of approximately $`20Gflops`$ for peak performance. We have performed production runs since July of 2001 and have started intensive parallel simulations in 2003. We consider applications that require moderate computer power, such as several problems for the $`SU(2)`$ gauge theory. In particular, we propose a new method for the study of the running coupling constant $`\alpha _s`$, based on the calculation of gluon and ghost propagators . We also carry out numerical studies of the infrared properties of QCD (through the study of the behavior of gluon and ghost propagators in the infrared limit) , of gauge-fixing techniques , of the chiral phase transition of QCD with two dynamical fermions  \[for which one expects analogies with the $`O(4)`$ spin model\], of the electroweak phase transition and aspects of spin models and percolation.

## Acknowledgements

The author thanks Attilio Cucchieri for the help in the preparation of this manuscript. Research supported by FAPESP, Proc. No. 00/05047-5.


format: text

text:
# Talk presented at “Lattice 2003”, July 15-19, 2003, Tsukuba, JapanKEK-TH-915RBRC-334 Nucleon axial charge and structure functions with domain wall fermions

## Abstract

We report the current status of RBCK calculations on nucleon structure with both quenched and unquenched lattice QCD. The combination of domain wall fermions and DBW2 gauge action works well for isovector vector and axial charges, and moments of structure functions $`x_q`$, $`d_1`$, and $`1_{\delta q}`$.

Domain wall fermions (DWF)  have been successfully used by the RBC Collaboration for several applications . In the baryon sector we had reported that DWF with Wilson gauge action in the quenched approximation can reproduce the first negative-parity excited state of nucleon  which defied theoretical attempts for a long time. Here we report the current status of our calculations on the nucleon electroweak matrix elements, namely the isovector vector and axial charges and some moments of the structure functions, in both quenched and full-QCD, using DWF and the DBW2 gauge action .

Four form factors appear in neutron $`\beta `$ decay: the vector and induced tensor form factors from the vector current, and the axial-vector and induced pseudo-scalar form factors from the axial-vector current. In the forward limit, the 4-momentum transfer should be small because the mass difference of the neutron and proton is only about 1.3 MeV. This makes the limit $`q^20`$, where the vector and axial-vector form factors dominate, a good approximation. Their values in this limit are called the vector and axial charges of the nucleon: $`g__V`$ and $`g__A`$.

Lattice numerical calculation of these charges is relatively simple: prepare appropriate nucleon source and sink at rest, insert an appropriate current operator in between, and identify a plateau which gives a lattice estimate. There already had been several preceding works using Wilson or improved Wilson fermions . However since these fermions schemes do not respect the chiral symmetry, the lattice renormalizations of the relevant vector and axial currents are not the same. Worse, unwanted lattice artifact may result in unphysical mixing of chirally distinct operators. The combination of DWF and DBW2 preserves chiral symmetry well enough so that we are free of these difficulties: especially the relation $`Z__A=Z__V`$ is easily and well maintained, up to $`O(a^2)`$, or a few % at $`a0.15`$ fm. And in contrast to those preceding studies, we found a significant finite volume effect between the results calculated on lattices with $`(1.2\mathrm{fm})^3`$ and $`(2.4\mathrm{fm})^3`$ volumes : about 20% increase was observed. On the large volume we find $`g_A=1.212\pm 0.027(\mathrm{stat})\pm 0.024(\mathrm{norm})`$. The quoted systematic error is the dominant known one, corresponding to current renormalization. This theoretical first principles calculation, which does not yet include isospin breaking effects, yields a value of $`g_A`$ only a little bit below the experimental one, $`1.2670\pm 0.0030`$ . The dependence on the quark mass on the larger volume is very mild, while it is steep on the smaller.

First few moments of nucleon structure functions are as easy as these charges, in that they do not require finite momentum states. The structure functions such as $`F_1(x,Q^2)`$, $`F_2(x,Q^2)`$, $`g_1(x,Q^2)`$, and $`g_2(x,Q^2)`$ are measured in lepton deep inelastic scattering, and $`h_1(x,Q^2)`$ in the RHIC Spin experiment. Their moments are factorized in the operator product expansion into Wilson coefficients and non-perturbative quantities such as $`x^n_q`$, $`x^n_{\mathrm{\Delta }q}`$, $`d_n`$ or $`1_{\delta q}`$. The latter are nucleon forward matrix elements of certain local operators, such as $`\gamma D`$, $`\gamma _5\gamma D`$, and $`\gamma _5\sigma D`$, and hence are accessible by lattice QCD numerical methods similar to form factor calculations. There are calculations by QCDSF, QCDSF/UKQCD, and LHPC-SESAM collaborations with (improved) Wilson fermions . In Lattice 2002 we presented some preliminary quenched results from the present DWF/DBW2 calculations . The DWF/DBW2 combination indeed works better: good chiral behavior, ease of non-perturbative renormalization, and better approach to the continuum were seen. We present higher statistics quenched results here, and some preliminary full-QCD ones. We still lack non-perturbative renormalization for some observables.

In the on-going full-QCD calculations, a $`16^3\times 32\times 12`$ lattice is used with the DBW2 gauge coupling set at $`\beta `$=0.8 and the domain-wall height at $`M_5`$=1.8. Dynamical DWF quarks with mass 0.02, 0.03, and 0.04 (about strange mass)  in the lattice units are used. The nucleon structure calculations reported here are done only at the dynamical quark mass, and no partially quenched study has been made. About 900–200 conjugate-gradient iterations are necessary to achieve $`10^8`$ accuracy for the Dirac equation solver. Hybrid Monte Carlo acceptance of about 70-80 % is observed. About 50 configurations are used at each dynamical quark mass value. The residual mass of $`m_{\mathrm{res}}=0.00137`$, $`\rho `$ meson mass $`m_\rho =0.43(2)`$ which translates into the lattice cutoff estimate of $`a^1`$ $``$ 1.8 GeV, and the nucleon mass $`m_N=0.56(4)`$ or $``$1.01(7) GeV are obtained. Quality of plateaus in the three-point functions is similar to the quenched case at similar statistics.

In Fig. 1 we show the ratio of the isovector axial and vector charges as a function of the pion mass squared. The full-QCD results (diamonds) are consistent with the old quenched results (circles) which are also shown: no mass dependence is seen. The weighted average of the three full-QCD points is $`g__A/g__V`$=1.21(8) and is consistent with the quenched results.

In Fig. 2 we present the quark momentum fraction in quenched (circle) and full QCD (diamond). Though not yet renormalized, in the quenched approximation it does not show any curvature toward the chiral limit. The first moment of the transversity is shown in Fig. 3. In both the quenched and full QCD we obtain the chiral limit of the bare matrix element by linear extrapolation. We use the RI-MOM scheme with one loop running to non-perturbatively renormalize the matrix element. The quenched result is 1.194(25) ($`\overline{\mathrm{MS}}`$, 2 GeV) and the full QCD result will be reported soon.

We have also calculated the $`d_1`$ matrix element for which lack of chiral symmetry results in a mixing with a lower dimensional operator. We reported last year that this mixing is suppressed with domain wall fermions in the quenched approximation . Our full-QCD $`d_1`$ results, shown in Fig. 4, are again in sharp contrast with the unsubtracted Wilson fermion results indicating that DWF suppress the power divergent mixing.

To summarize, the quenched calculations with DBW2 and DWF works well for nucleon matrix elements. The axial to vector charge ratio, $`g__A/g__V`$, is well reproduced to within 5 % of the experimental value even on a modest volume and relatively heavy quark mass. The dependence on quark mass is very mild. Likewise the moments of structure functions $`x_q`$, $`d_1`$, and $`1_{\delta q}`$ all show mild linear dependence on quark mass, though some of them are yet to be renormalized. The full-QCD calculations of these quantities seem to work as well. Similar mild linear dependence on the quark mass was observed, extrapolating to estimates not contradicting with the quenched ones (where renormalization is done), albeit with limited statistics. With QCDOC computer , we will be adding more observables such as flavor-singlet quantities and form factors and higher moments of structure functions that require finite momentum, and moving to larger lattice volumes and smaller quark masses.


text:
# Nuclear profile dependence of elliptic flow from a parton cascade

## Abstract

The transverse profile dependence of elliptic flow is studied in a parton cascade model. We compare results from the binary scaling profile to results from the wounded nucleon scaling profile. The impact parameter dependence of elliptic flow is shown to depend sensitively on the transverse profile of initial particles, however, if elliptic flow is plotted as a function of the relative multiplicity, the nuclear profile dependence disappears. The insensitivity was found previously in a hydrodynamical calculation. Our calculations indicate that the insensitivity is also valid with additional viscous corrections. In addition, the minimum bias differential elliptic flow is demonstrated to be insensitive to the nuclear profile of the system.

Recently, the Relativistic Heavy Ion Collider has produced large amount of exciting new data. These new data give us valuable insight into the hot and dense nuclear matter. One of the important observables is elliptic flow which reflects the transverse anisotropy of particle momentum distribution. Elliptic flow has been studied by many theoretical models, including non-Abelian energy loss models , saturation models , parton recombination models , hydrodynamical models , and parton cascade models . In this paper, we will study the elliptic flow using a parton cascade model. We will first introduce the elliptic flow, and the parton cascade model used for this study. Then, we will use the parton cascade model to study the elliptic flow produced from two different initial transverse distributions, one proportional to the number of binary collisions, one proportional to the number of wounded nucleons. We demonstrate that even though, elliptic flow as a function of the impact parameter is very sensitive to the initial transverse distribution of particles, elliptic flow as a function of the relative multiplicity is almost independent of the transverse distribution. Furthermore, we show that the minimum bias differential elliptic flow is also insensitive to the initial transverse particle distribution.

Elliptic flow is the elliptic deformation in the particle transverse momentum distribution . It is usually characterized by the second Fourier coefficient of the particle azimuthal distribution . If we use $`f(\varphi )`$ for the azimuthal distribution, and choose the azimuthal angle of the reaction plane to be zero, then

$$f(\varphi )=v_0+2v_1\mathrm{cos}(\varphi )+2v_2\mathrm{cos}(2\varphi )+\mathrm{}.$$
(1)
The coefficient, $`v_2`$, is the elliptic flow observable. It is the average of $`\mathrm{cos}(2\varphi )`$ of produced particles. If the transverse components of the momenta are know, it can be calculated by:

$$v_2=\frac{p_x^2p_y^2}{p_x^2+p_y^2}.$$
(2)
In the above formula, $`\mathrm{}`$ denotes the average over particles. Since initial momentum distribution is isotropic, the elliptic flow is generated by final state interactions. Final state interactions (or pressure gradient) turn(s) the spatial anisotropy into momentum space anisotropy. It has been shown that the elliptic flow is very sensitive to the initial stage evolution, and can be used as a sensitive probe of early dynamics .

In the following, we are going to use Zhang’s Parton Cascade (ZPC)  to study elliptic flow at relativistic energies. The initial conditions are set up similar to those in a recent study by Molnar’s Parton Cascade (MPC) . In the local rest frame, the initial momentum distribution is thermal, with a temperature of 700 MeV. Particles are uniformly distributed between a space time rapidity of -5 and +5. The particle formation proper time is 0.1 fm/c. There are totally 2100 gluons per central event. As the momentum transport is determined by the momentum opacity, the following results are also correct if the total number of particles increases and the transport cross section decreases by the same factor. To efficiently simulate momentum transport, we use isotropic differential cross sections that preserve the reaction plane of a collision. The total parton-parton elastic cross section $`\sigma _{gg}`$ will be varied to have values of 40 mb, 20 mb, 10 mb to study the response of the system. These cross sections are effective cross sections as no radiative energy loss or parton recombinations are included in the calculations.

We will study the elliptic flow produced from two different initial transverse spatial distributions. One is proportional to the number of binary collisions per unit area. In this case, the particle number as a function of the impact parameter is also proportional to the number of binary collisions. The other distribution is proportional to the number of wounded nucleons per unit area and the particle number as a function of the impact parameter is also proportional to the number of wounded nucleons. In generating the above distributions, the nucleon-nucleon inelastic cross section $`\sigma _{NN}`$ is taken to be 40 mb, and the three parameter Woods-Saxon distribution is used for the nucleons inside one nucleus. These two spatial distributions are related to hard and soft particle production mechanisms, respectively .

The above initial spatial distribution and initial momentum distribution factorize. As local densities are sampled for the evolution of the expanding parton system according to the Boltzmann equation, the factorization is not automatically conserved. We also note that a geometry with sharp cylindrical nuclei always leads to larger elliptic flow values compared to the binary collision scaling case .

Fig. 1 gives the number of particles as a function of the impact parameter. While the two distributions have the same number of particles when $`b=0`$, the wounded nucleon scaling has more particles than the binary collision scaling case. In particular, when $`b=10`$ fm, the wounded nucleon scaling has about twice as many particles as the binary collision scaling case. Fig. 2 has the initial spatial ellipticity as a function of the impact parameter. The initial spatial ellipticity in the figure is defined through

$$ϵ=\frac{y^2x^2}{y^2+x^2}.$$
(3)
Note that this definition calculates the ratio first and then the average and the magnitude is smaller than the ratio of the averages.

We first study the impact parameter dependence of the elliptic flow. The set up is similar to the recent MPC model study . We approach the Boltzmann limit by increasing the number of particles and at the same time decreasing the cross section by the same factor . In the binary collision scaling case, the rescaling factors for $`b=0,2,4,6,8,10,12,14`$ fm are $`\lambda =100,100,100,220,450,1100,5000,50000`$. In the wounded nucleon scaling case, the rescaling factors are $`\lambda =100,100,100,200,300,600,2000,11000`$. The convergence is checked by comparing to calculations done with $`\lambda /2`$. The $`v_2`$ is calculated for particles with a rapidity range of $`|y|<2`$.

From Fig. 3, we observe that as the total cross section increases, or more precisely, as the transport cross section increases, the elliptic flow increases. The binary collision scaling case is larger than the wounded nucleon scaling case for small impact parameters and smaller than the wounded nucleon scaling case for large impact parameters. This follows the trend of the initial spatial ellipticity. However, the $`v_2`$ curves peak at smaller impact parameters than the $`ϵ`$ curves. This indicates that both initial ellipticity and initial particle density play roles in determining the elliptic flow. As the impact parameter increases, the elliptic flow increases with the initial ellipticity, however, after a point, the particle density is not high enough to generate enough response and the elliptic flow can not catch up with the initial ellipticity. It starts decreasing with increasing impact parameter.

An alternative way of characterizing centrality is to use the relative central rapidity density, which is the ratio of the central rapidity density to that in central collisions with $`b=0`$. If we plot the elliptic flow as a function of the relative central rapidity density as in Fig. 4, we see that for the same transport cross section, the two curves corresponding to the binary collision scaling and the wounded nucleon scaling overlap. In other words, the impact parameter dependence of elliptic flow is cancelled by the impact parameter dependence of the multiplicity. Similar observations have been made in a recent hydrodynamics study . This indicates that if we use the relative central rapidity density as a measure of centrality, the elliptic flow is not sensitive to whether the initial distribution is binary collision scaling or wounded nucleon scaling. It reflects the particle transport cross section, or in the case of hydrodynamics, the equation of state. Because of the viscous corrections, the cascade calculations have a bend over when the relative central rapidity density is small. In contrast, the hydrodynamic calculations have an almost straight line dependence and overshoot data when the relative rapidity density is small. We also note that the $`\sigma _{gg}=40`$ mb binary collision scaling case is consistent with set D of .

Now we turn to the study of $`p_t`$ differential elliptic flow which can give further information about the evolution . Hydrodynamic studies agree well with low $`p_t`$ data. At $`p_t>2`$ GeV, the data are consistent with a constant behavior while hydrodynamic results keep on increasing. In dynamic models, only when viscous effects are taken into account, is it possible to describe the deviation from the ideal hydrodynamical behavior. A recent hydrodynamic study demonstrates that the minimum bias $`p_t`$ differential flow is not sensitive to the initial nuclear profile. We want to know whether it is also true when viscous effects are taken into account. Fig. 5 shows the impact parameter averaged $`p_t`$ differential flow. The calculations with binary collision scaling agree well with those with wounded nucleon scaling. This is true not only for the low $`p_t`$ region, but also for the high $`p_t`$ region where viscous effects are important.

Another way of averaging over events is to calculate the multiplicity weighted average of cosine of the azimuthal angle. This gives the minimum bias elliptic flow. Results from the ZPC model are shown in Fig. 6. The binary collision scaling and the wounded nucleon scaling agree well with each other. This further demonstrates that the minimum bias $`p_t`$ differential elliptic flow is insensitive to the initial nuclear profile. A comparison of Fig. 5 and Fig. 6 shows that the minimum bias results are about 10% higher than the impact parameter averages. Hence the impact parameter averaged differential elliptic flow can be considered as a reasonably good approximation of the minimum bias differential elliptic flow. As pointed out in , the minimum bias elliptic flow weights in more central events. The central events in the wounded nucleon scaling case have lower elliptic flow than those in the binary collision scaling case. This can lead to a relatively smaller minimum bias elliptic flow in the wounded nucleon scaling case. However, the decrease in the relative amplitude can not be determined with the current statistics.

In summary, we demonstrate that the elliptic flow as a function of relative central rapidity density in the binary scaling case agrees well with that in the wounded nucleon case. In addition, the minimum bias $`p_t`$ differential elliptic flow in the binary scaling case also agrees well with that in the wounded nucleon case. This is true not only for the low $`p_t`$ region, but also for the high $`p_t`$ region where viscous effects are important and ideal hydrodynamics deviates from experimental data. As hadronization will not change the high $`p_t`$ elliptic flow , and recent research indicates the possibility of extracting the parton elliptic flow from hadron elliptic flow , the elliptic flow as a function of relative rapidity density and the minimum bias $`p_t`$ differential elliptic flow are promising observables for the extraction of information about final state interactions .

This work was supported by the U.S. National Science Foundation under Grant No. 0140046. We thank P. Kolb, Z.W. Lin, D. Molnar, and A. Poskanzer for helpful discussions. We also thank the Parallel Distributed System Facility at the National Energy Research Scientific Computer Center for providing computer resources.


text:
# 1 Introduction

## 1 Introduction

The Weak Equivalence Principle (WEP) and the Einstein Equivalence Prnciple (EEP) (Will 1993; Ciufolini and Wheeler 1995; Haugan and L$`\ddot{\mathrm{a}}`$mmerzahl 2001; Will 2001) are the cornerstones of Einstein General Theory of Relativity (GTR) and of all the other competing metric theories of gravity. The WEP states that non–rotating, uncharged bodies of different structure and compositions and with negligible amount of gravitational binding energy per unit mass fall with the same acceleration in a given gravitational field, provided that no other forces act on them. The EEP (called also Medium Strong Equivalence Principle) states that the outcome of any local, non–gravitational test experiment is independent of where and when in the gravitational field the experiment is performed. If the gravitational binding energy per unit mass of the freely falling bodies in a given external gravitational field is not negligible, as is the case for astronomical bodies, we have the Strong Equivalence Priciple (SEP. It is called also Very Strong Equivalence Principle). When the self–gravity of the falling bodies is accounted for it might happen, in principle, that their accelerations are different. This violation of the universality of free fall, called Nordvedt effect (Nordvedt 1968a), is predicted by all metric theories of gravity apart from GTR which, instead, satisfies the SEP. Lunar Laser Ranging (LLR) is able to test this effect for the motion of Earth and Moon in the gravitational field of Sun (Nordvedt 1968b). No different accelerations have been found at a 10-13 level (Anderson and Williams 2001). Of course, this is also a test of the WEP.

To go to space, where much larger distances, velocities, gravitational potential differences with respect to Earth, and, most important, free fall for an, in principle, infinitely long time are available, is important to perform very accurate tests of post–Newtonian gravity. In this paper we wish to investigate the level of accuracy which could be reached in testing the WEP, by analysing the orbital motion of a pair of Earth artificial satellites of different compositions for which the effects of self–gravity would be, of course, negligible: it is analogous to the analysis of the motion of Earth and Moon in the gravitational field of Sun in order to test the SEP. For the already performed experimental tests of the WEP on Earth see (Will 2001) and the references therein. The present–day level of accuracy is $`4\times 10^{13}`$ (B$`\ddot{\mathrm{a}}`$ssler $`etal`$ 1999). It has been obtained in the experiments of the so–called E$`\ddot{\mathrm{o}}`$t-Wash group of the Washington University by means of a sophisticated torsion balance in an Earth based laboratory set–up. Other proposed space–based experiments are the very complex and expensive STEP (Lockerbie $`etal`$ 2001), $`\mu `$SCOPE (Touboul 2001), GG (Nobili $`etal`$ 2000) and SEE (Sanders $`etal`$ 2000) missions111See also on the WEB http://einstein.stanford.edu/STEP/, http://www.onera.fr/microscope/ http://tycho.dm.unipi.it/$``$nobili/ggproject.html, and http://www.phys.utk.edu/see/. Notice that while the first three missions are in advanced stages of planning and hardware testing, and are expected to be launched in the next few years, SEE is still undergoing rigorous conceptual evaluation and is not yet a scheduled mission. whose goal is to reach the $`10^{15}`$$`10^{18}`$ level of accuracy.

The reason for searching for more and more accurate tests of the equivalence principle resides in the fact that all approaches to quantizing gravity and to unifying it with the other fundamental interactions currently under study are capable of predicting violations of the equivalence principle at some level. For example, departures from universal free fall accelerations of the order of $`10^{15}`$ have been calculated in (Damour and Polyakov 1994a; 1994b) in the string theory context. Violations of the WEP are predicted also by nonsymmetric theories of gravity (Will 1989). In (Moffat and Gillies 2002) a violation of the WEP of the order of 10-14 is predicted in the context of a non–local quantum gravity theory.

## 2 The orbital period

In this section we examine the possibility of testing the WEP by measuring the orbital periods of a pair of Earth artificial satellites of different compositions. It can be thought of as a comparison between two pendulums with enormously long threads swinging for an extremely long time222In the case of the SEP–LLR experiment the amplitude of the parallactic inequality long–period harmonic perturbation, proportional to $`\mathrm{cos}D`$, where $`D`$ is the synodic phase from New Moon, is sensitive to the possible different falling rates of Earth and Moon toward the Sun..

The orbital period of a satellite freely falling in Earth’s gravitational field can be written as

$$T=\frac{2\pi }{n(1+\mathrm{\Delta }n)}\frac{2\pi }{n}(1\mathrm{\Delta }n),$$
($`1`$)
in which the Keplerian unperturbed period is

$$T^{(0)}=\frac{2\pi }{n}=2\pi \sqrt{\frac{a^3}{GM}}\times \sqrt{\frac{m_\mathrm{i}}{m_\mathrm{g}}}=2\pi \sqrt{\frac{a^3}{GM}}\times \psi ,$$
($`2`$)
where $`G`$, $`M`$ and $`a`$ are the Newtonian gravitational constant, the mass of Earth and the satellite semimajor axis, respectively. We have explicitly written the square root of the ratio of the inertial to the (passive) gravitational mass of satellite $`\psi `$. The quantity $`\mathrm{\Delta }n`$ represents the various kind of perturbations, of gravitational and non–gravitational origin, which affect $`n`$. For example, the even zonal harmonic coefficients of the multipolar expansion of Earth gravitational potential, called geopotential, induce secular perturbations on $`n`$. The most important one is that due to the first even zonal harmonic $`J_2`$ and, for a circular orbit with eccentricity $`e=0`$, it is given by

$$\mathrm{\Delta }n_{\mathrm{obl}}^{(\mathrm{}=2)}=\frac{3}{4}J_2\left(\frac{R}{a}\right)^2(13\mathrm{cos}^2i),$$
($`3`$)
where $`R`$ is the Earth radius and $`i`$ is the inclination of the orbital plane to the equator. Also the time–varying part of Earth gravitational field should be considered, in principle, because the Earth solid and ocean tides (Iorio 2001a) induce long–period harmonic perturbations on $`n`$. For a given tidal line of frequency $`f`$ the perturbations induced by the solid tides, which are the most effective in affecting the orbits of a satellite, can be written as

$$\mathrm{\Delta }n_{\mathrm{tides}}^{\mathrm{}}=\underset{m=0}{\overset{\mathrm{}}{}}\left(\frac{H_{\mathrm{}}^m}{R}\right)\left(\frac{R}{a}\right)^{\mathrm{}}k_\mathrm{}m^{(0)}A_\mathrm{}m\underset{p=0}{\overset{\mathrm{}}{}}\underset{q=\mathrm{}}{\overset{+\mathrm{}}{}}F_{\mathrm{}mp}\left[2(\mathrm{}+1)G_{\mathrm{}pq}\frac{(1e^2)}{e}\frac{dG_{\mathrm{}pq}}{de}\right]\mathrm{cos}\gamma _{f\mathrm{}mpq},$$
($`4`$)
where $`H_{\mathrm{}}^m`$ are the tidal heights, $`k_\mathrm{}m^{(0)}`$ are the Love numbers, $`A_\mathrm{}m=\sqrt{\frac{(2\mathrm{}+1)(\mathrm{}m)!}{4\pi (\mathrm{}+m)!}}`$, $`F_{\mathrm{}mp}(i)`$ and $`G_{\mathrm{}pq}(e)`$ are the inclination and the eccentricity functions (Kaula 1966), respectively, and $`\gamma _{f\mathrm{}mpq}`$ is the frequency of the tidal perturbation built up with the frequencies of the lunisolar variables and of the satellite’s orbital elements.

In order to test the WEP, we propose to measure, after many revolutions, the difference of time spans which are multiple $`N`$ of the orbital periods $`\mathrm{\Delta }T_NT_N^{(2)}T_N^{(1)}`$ of a couple of satellites of different composition orbiting the Earth along circular orbits of almost same radius $`a`$

$$a^{(1)}a,a^{(2)}=a+d.$$
($`5`$)
The small difference $`d`$, which, in principle, should be equal to zero, is due to the unavoidable orbital injection errors. It can be made extremely small with a rocket launcher of good quality: for example, at the beginning of their mission the semimajor axes of the two GRACE spacecrafts were different by an amount of just333See on the WEB http://www.csr.utexas.edu/grace/newsletter/2002/august2002.html 0.5 km.

The observable quantity in which we are interested is $`\mathrm{\Delta }\psi `$ to be measured after $`N`$ orbital revolutions, where

$$\psi ^{(1)}\psi =1,\psi ^{(2)}=\psi +\mathrm{\Delta }\psi ,$$
($`6`$)
with $`\mathrm{\Delta }\psi /\psi 1`$. The quantity $`\mathrm{\Delta }\psi `$ can be expressed in terms of the standard E$`\ddot{\mathrm{o}}`$tv$`\ddot{\mathrm{o}}`$s parameter $`\eta `$. Indeed, the inertial mass of a body is composed by many types of mass–energy: rest energy, electromagnetic energy, weak–interaction energy, and so on. If one of these forms of energy contributes to $`m_\mathrm{g}`$ differently than it does to $`m_\mathrm{i}`$ we can put (Will 1993)

$$m_\mathrm{g}=m_\mathrm{i}+\underset{A}{}\eta ^A\frac{E^A}{c^2},$$
($`7`$)
where $`E^A`$ is the internal energy of the body generated by the interaction $`A`$ and $`\eta ^A`$ is a dimensionless parameter that measures the strength of the violation of WEP induced by that interaction. Then

$$\frac{m_\mathrm{i}^{(2)}}{m_\mathrm{g}^{(2)}}\frac{m_\mathrm{i}^{(1)}}{m_\mathrm{g}^{(1)}}\underset{A}{}\eta ^A\left[\frac{E_{(1)}^A}{m_\mathrm{i}^{(1)}c^2}\frac{E_{(2)}^A}{m_\mathrm{i}^{(2)}c^2}\right]\eta .$$
($`8`$)
From eq.($`8`$) it can be obtained

$$\psi ^{(2)}\sqrt{\frac{m_\mathrm{i}^{(2)}}{m_\mathrm{g}^{(2)}}}=\sqrt{\frac{m_\mathrm{i}^{(1)}}{m_\mathrm{g}^{(1)}}+\left[\frac{m_\mathrm{i}^{(2)}}{m_\mathrm{g}^{(2)}}\frac{m_\mathrm{i}^{(1)}}{m_\mathrm{g}^{(1)}}\right]}=\sqrt{\frac{m_\mathrm{i}^{(1)}}{m_\mathrm{g}^{(1)}}+\eta }\psi ^{(1)}+\frac{\eta }{2\psi ^{(1)}},$$
($`9`$)
so that, from eq.($`6`$),

$$\mathrm{\Delta }\psi =\frac{\eta }{2}.$$
($`10`$)
From eqs.($`1`$)-($`2`$) we can write

$$\mathrm{\Delta }T_N=N2\pi \left[\frac{1}{n^{(2)}}\frac{1}{n^{(1)}}+\frac{\mathrm{\Delta }n_{\mathrm{obl}}^{(1)}}{n^{(1)}}\frac{\mathrm{\Delta }n_{\mathrm{obl}}^{(2)}}{n^{(2)}}\right].$$
($`11`$)
In eq.($`11`$) we consider only the gravitational even zonal perturbations due to Earth oblateness, as usually done in orbital reduction programs in which a reference orbit including the $`J_2`$ effects is adopted444This approximation will be justified later.. They can be summarized as

$$\mathrm{\Delta }n_{\mathrm{obl}}=\underset{\mathrm{}=2}{}\left(\frac{R}{a}\right)^{\mathrm{}}𝒢_{\mathrm{}},$$
($`12`$)
where the $`𝒢_{\mathrm{}}=𝒢_{\mathrm{}}(i,e;J_{\mathrm{}})`$ functions include the even zonal harmonics coefficients $`J_{\mathrm{}}`$, the eccentricity $`e`$, the inclination angle $`i`$ and some numerical constants. For example, for $`\mathrm{}=2`$ and $`e=0`$, $`𝒢_2=\frac{3}{4}J_2(13\mathrm{cos}^2i)`$. By using the expansions of eqs.($`5`$)-($`6`$) it is possible to solve eq.($`11`$) for555Here we have considered $`i^{(1)}=i^{(2)}`$. Moreover, notice that if the two satellites would be inserted in counter–rotating orbits, in $`\mathrm{\Delta }T_N`$ it must be included also the time shift $`\mathrm{\Delta }T_N^{\mathrm{gm}}N\times 4\pi \frac{J}{c^2M}`$ due to the general relativistic gravitomagnetic clock effect (Mashhoon $`etal`$ 1999; 2001; Iorio $`etal`$ 2002) induced by the off-diagonal components of the metric proportional to the proper angular momentum $`J`$ of Earth. Indeed, it turns out (Iorio $`etal`$ 2002) that such time shift is independent of the ratio of the inertial to the gravitational masses of the satellite. It would act as a lower limit of the order of $`N\times 10^7`$ s to $`(\mathrm{\Delta }T_N)_{\mathrm{exp}}`$ if it was detectable. $`\mathrm{\Delta }\psi _N`$

$$\mathrm{\Delta }\psi _N=\frac{\frac{\sqrt{GM}}{N\pi }\mathrm{\Delta }T_N(A+B)}{(A+C)+(B+D)},$$
($`13`$)
with

$`A`$ $`=`$ $`3d\sqrt{a},`$ ($`14`$)
$`B`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}𝒢_{\mathrm{}}R^{\mathrm{}}\left[3\mathrm{}d^2a^{\left(\frac{1+2\mathrm{}}{2}\right)}+2\mathrm{}da^{\left(\frac{12\mathrm{}}{2}\right)}3da^{\left(\frac{12\mathrm{}}{2}\right)}\right],`$ ($`15`$)
$`C`$ $`=`$ $`2\sqrt{a^3},`$ ($`16`$)
$`D`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}𝒢_{\mathrm{}}R^{\mathrm{}}\left[2a^{\left(\frac{32\mathrm{}}{2}\right)}\right].`$ ($`17`$)
Eq.($`13`$), together with eqs.($`14`$)-($`17`$), allows to evaluate the accuracy obtainable in measuring the quantity $`\mathrm{\Delta }\psi _N`$. The error in measuring the difference of multiples of the orbital periods yields

$$[\delta (\mathrm{\Delta }\psi _N)]_{\mathrm{\Delta }T_N}=\frac{\sqrt{GM}}{N\pi }\frac{\delta (\mathrm{\Delta }T_N)}{(A+C)+(B+D)},$$
($`18`$)
while the error in Earth $`GM`$, which amounts to 8$`\times 10^{11}`$ cm3 s-2 (McCarthy 1996), yields

$$[\delta (\mathrm{\Delta }\psi _N)]_{GM}=\frac{(\mathrm{\Delta }T_N)_{\mathrm{exp}}}{N2\pi \sqrt{GM}}\frac{\delta (GM)}{(A+C)+(B+D)}.$$
($`19`$)
The uncertainty in Earth even zonal harmonics $`J_{\mathrm{}}`$, of which $`\delta J_2`$ amounts to 7.9626$`\times 10^{11}`$ (Lemoine $`etal`$ 1998), has an impact given by

$$[\delta (\mathrm{\Delta }\psi _N)]_J_{\mathrm{}}\underset{\mathrm{}=2}{}\left|\left\{\frac{\frac{B}{J_{\mathrm{}}}\times [(A+C)+(B+D)]+(A+B)\times \frac{(B+D)}{J_{\mathrm{}}}}{[(A+C)+(B+D)]^2}\right\}\right|\times \delta J_{\mathrm{}}.$$
($`20`$)
The errors in $`\mathrm{\Delta }\psi _N`$ due to the uncertainties in $`a`$ and in $`d`$ are

$`[\delta (\mathrm{\Delta }\psi _N)]_a`$ $`=`$ $`\left|\left\{{\displaystyle \frac{\frac{(A+B)}{a}\times [(A+C)+(B+D)]+(A+B)\times \frac{[(A+C)+(B+D)]}{a}}{[(A+C)+(B+D)]^2}}\right\}\right|\times \delta a,`$ ($`21`$)
$`[\delta (\mathrm{\Delta }\psi _N)]_d`$ $`=`$ $`\left|\left\{{\displaystyle \frac{\frac{(A+B)}{d}\times (C+D)}{[(A+C)+(B+D)]^2}}\right\}\right|\times \delta d.`$ ($`22`$)
Notice that $`[\delta (\mathrm{\Delta }\psi _N)]_J_{\mathrm{}}`$, $`[\delta (\mathrm{\Delta }\psi _N)]_a`$ and $`[\delta (\mathrm{\Delta }\psi _N)]_d`$ do not depend explicitly on the number of orbital revolutions $`N`$. In order to calculate them we need the explicit expressions of the derivatives with respect to $`J_{\mathrm{}}`$

$`{\displaystyle \frac{B}{J_{\mathrm{}}}}`$ $`=`$ $`\left({\displaystyle \frac{𝒢_{\mathrm{}}}{J_{\mathrm{}}}}\right)R^{\mathrm{}}\left[3\mathrm{}d^2a^{\left(\frac{1+2\mathrm{}}{2}\right)}+2\mathrm{}da^{\left(\frac{12\mathrm{}}{2}\right)}3da^{\left(\frac{12\mathrm{}}{2}\right)}\right],`$ ($`23`$)
$`{\displaystyle \frac{D}{J_{\mathrm{}}}}`$ $`=`$ $`\left({\displaystyle \frac{𝒢_{\mathrm{}}}{J_{\mathrm{}}}}\right)R^{\mathrm{}}\left[2a^{\left(\frac{32\mathrm{}}{2}\right)}\right],`$ ($`24`$)
with respect to $`a`$

$`{\displaystyle \frac{A}{a}}`$ $`=`$ $`{\displaystyle \frac{3d}{2\sqrt{a}}},`$ ($`25`$)
$`{\displaystyle \frac{B}{a}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}𝒢_{\mathrm{}}R^{\mathrm{}}\left[3{\displaystyle \frac{(1+2\mathrm{})}{2}}\mathrm{}d^2a^{\left(\frac{3+2\mathrm{}}{2}\right)}+(12\mathrm{})\mathrm{}da^{\left(\frac{1+2\mathrm{}}{2}\right)}3{\displaystyle \frac{(12\mathrm{})}{2}}da^{\left(\frac{1+2\mathrm{}}{2}\right)}\right],`$ ($`26`$)
$`{\displaystyle \frac{C}{a}}`$ $`=`$ $`3\sqrt{a},`$ ($`27`$)
$`{\displaystyle \frac{D}{a}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}𝒢_{\mathrm{}}R^{\mathrm{}}\left[(32\mathrm{})a^{\left(\frac{12\mathrm{}}{2}\right)}\right],`$ ($`28`$)
and those with respect to $`d`$

$`{\displaystyle \frac{A}{d}}`$ $`=`$ $`3\sqrt{a},`$ ($`29`$)
$`{\displaystyle \frac{B}{d}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}𝒢_{\mathrm{}}R^{\mathrm{}}\left[6\mathrm{}da^{\left(\frac{1+2\mathrm{}}{2}\right)}+2\mathrm{}a^{\left(\frac{12\mathrm{}}{2}\right)}3a^{\left(\frac{12\mathrm{}}{2}\right)}\right].`$ ($`30`$)
In order to fix the ideas, let us consider the orbit of the proposed LARES laser–ranged satellite with $`a=12270`$ km, $`i`$=70 deg. Let us assume (Peterson 1997) $`d=5`$ km. With these data we have, for $`\mathrm{}=2`$,

$`[\delta (\mathrm{\Delta }\psi _N)]_{\mathrm{\Delta }T_N}`$ $`=`$ $`7.3\times 10^5\mathrm{s}^1\times {\displaystyle \frac{\delta (\mathrm{\Delta }T_N)}{N}},`$ ($`31`$)
$`[\delta (\mathrm{\Delta }\psi _N)]_{GM}`$ $`=`$ $`7.4\times 10^{14}\mathrm{s}^1\times {\displaystyle \frac{(\mathrm{\Delta }T_N)_{\mathrm{exp}}}{N}},`$ ($`32`$)
$`[\delta (\mathrm{\Delta }\psi _N)]_{J_2}`$ $`=`$ $`8\times 10^{15},`$ ($`33`$)
$`[\delta (\mathrm{\Delta }\psi _N)]_a`$ $`=`$ $`5\times 10^{13}\mathrm{cm}^1\times \delta a,`$ ($`34`$)
$`[\delta (\mathrm{\Delta }\psi _N)]_d`$ $`=`$ $`1\times 10^9\mathrm{cm}^1\times \delta d.`$ ($`35`$)
The accuracy in measuring the difference of the multiples of the orbital periods of the satellites is a crucial factor in obtaining a high precision in $`\mathrm{\Delta }\psi `$. However, it could be possible, in principle, to choose an observational time span covering a very high number of orbital revolutions. It should be noted that $`\delta (\mathrm{\Delta }T_N)`$ accounts for both the measurement errors and the systematical errors induced by various gravitational and non–gravitational aliasing phenomena. The latter ones play a very important role in strongly limiting the possibility of measuring, e.g., the gravitomagnetic clock effect: it amounts to 10-7 s after one orbital revolution for orbits with $`e=i=0`$ while in (Iorio 2001b) it turned out that the systematic errors induced by the present–day level of knowledge of the terrestrial gravitational field are up to 2–3 order of magnitude larger. So, it should not be unrealistic to consider $`\delta (\mathrm{\Delta }T_N)_{\mathrm{sys}}10^410^5`$ s. This would imply $`[\delta (\mathrm{\Delta }\psi _N)]_{\mathrm{\Delta }T_N}\frac{10^910^{10}}{N}`$.

The error due to the difference $`d`$ in the semimajor axes of the satellites, which turns out to be the major limiting factor, cannot be reduced, in principle, by waiting for a sufficiently high number of orbital revolutions because it is independent of $`N`$. A small improvement could be obtained with the use of a larger semimajor axis. A geostationary orbit with $`a=42160`$ km would allow to get $`[\delta (\mathrm{\Delta }\psi _N)]_{\mathrm{\Delta }T_N}=1\times 10^5\mathrm{s}^1\times \frac{\delta (\mathrm{\Delta }T_N)}{N}`$ and $`[\delta (\mathrm{\Delta }\psi _N)]_d=3\times 10^{10}\mathrm{cm}^1\times \delta d.`$ However, in this case, for a fixed time span, we would have at our disposal a smaller $`N`$. In regard to the systematic part of the error $`\delta d`$, it should be noted that there are no secular or long–period perturbations of gravitational origin on the semimajor axis of a satellite. The non–gravitational perturbations could be reduced to a good level by adopting the drag–free technology. In regard to $`\delta d_{\mathrm{exp}}`$, it is important to note that in the present GRACE mission (Davis $`etal`$ 1999), making use of a K/Ka-band intersatellite link that provides dual one–way range measurements, changes in the distance of the two spacecrafts can be established with an accuracy of about 10-2 cm or even better.

It is interesting to note that if we neglect in all calculations $`B`$ and $`D`$ and the related derivatives, i.e., if we neglect the effects of Earth oblateness, it turns out that the numerical results of eqs.($`31`$)-($`35`$) do not change. It is very important because it means that our choice of neglecting the contribution of the non–gravitational perturbations in $`\mathrm{\Delta }n`$ is $`aposteriori`$ correct. Indeed, the perturbing acceleration due to Earth $`J_2`$ on LAGEOS is of the order of 10-1 cm s-2, while the impact of the direct solar radiation pressure, which is the largest non–gravitational perturbation on LAGEOS, amounts to 10-7 cm s-2 (Milani $`etal`$ 1987). Moreover, these conclusions imply that the errors in the inclination $`i`$, which enters the even zonal harmonic perturbations due to the geopotential and the non–gravitational perturbations, can be safely neglected, as done here. The same considerations hold also for the tidal perturbations: suffices it to say that the effects of the 18.6–year and the $`K_1`$ tides, which are the most powerful in perturbing the satellite orbits, on a LAGEOS–type satellite are six orders of magnitude smaller than those due to the static $`J_2`$ even zonal part of geopotential. Also the tiny general relativistic gravitoelectric correction to the orbital period induced by the Schwarzschild part of the metric666For a circular orbit it is given by $`T^{(0)}\times \mathrm{\Theta }_{\mathrm{ge}}=\psi \frac{3\pi \sqrt{GMa}}{c^2}`$ (Mashhoon $`etal`$ 2001)., which depends on $`\sqrt{m_\mathrm{i}/m_\mathrm{g}}`$, can be neglected because for a LAGEOS–type satellite the disturbing acceleration is of the order of 9$`\times 10^8`$ cm s-2 (Milani $`etal`$ 1987).

## 3 The longitude of the ascending node

The longitude of the ascending node $`\mathrm{\Omega }`$ is one of the best accurately measured Keplerian orbital elements of Earth artificial satellites. Then, we wish to examine if it would be possible to use it in order to test the equivalence principle.

Let us recall that there are two kinds of long–period perturbations on the node $`\mathrm{\Omega }`$ of an Earth satellite. First, the static oblateness of Earth induces a secular precession of $`\mathrm{\Omega }`$ through the even zonal harmonics of the geopotential. Second, the time–varying part of Earth gravitational potential induces tidal harmonic perturbations on $`\mathrm{\Omega }`$ (Iorio 2001a). Then, we can pose, by including the solid Earth tidal perturbations

$$\dot{\mathrm{\Omega }}=n\underset{\mathrm{}=2}{}\left(\frac{R}{a}\right)^{\mathrm{}}\left[𝒢_{\mathrm{}}+\underset{m=0}{\overset{\mathrm{}}{}}\left(\frac{H_{\mathrm{}}^m}{R}\right)k_\mathrm{}m^{(0)}A_\mathrm{}m\underset{p=0}{\overset{\mathrm{}}{}}\underset{q=\mathrm{}}{\overset{+\mathrm{}}{}}\frac{dF_{\mathrm{}mp}}{di}\frac{G_{\mathrm{}pq}}{\mathrm{sin}i\sqrt{1e^2}}\mathrm{cos}\gamma _{f\mathrm{}mpq}\right],$$
($`36`$)
where, for $`\mathrm{}=2`$ and $`e=0`$, $`𝒢_2=\frac{3}{2}J_2\mathrm{cos}i`$. By considering a couple of satellites of different compositions freely orbiting along almost identical orbits we could measure the difference of their secular nodal rates $`\mathrm{\Delta }\dot{\mathrm{\Omega }}\dot{\mathrm{\Omega }}^{(2)}\dot{\mathrm{\Omega }}^{(1)}`$. In this case, by posing $`\xi \sqrt{\frac{m_\mathrm{g}}{m_\mathrm{i}}}`$, the violating parameter $`\mathrm{\Delta }\xi =\frac{\eta }{2}`$ can be expressed as

$$\mathrm{\Delta }\xi =\frac{\frac{\mathrm{\Delta }\dot{\mathrm{\Omega }}}{\sqrt{GM}}+𝒜+}{𝒞},$$
($`37`$)
with777Here we neglect the non–gravitational perturbations on the nodes. Contrary to, e.g., the perigees $`\omega `$, the nodes are rather insensitive to such non–geodesic accelerations (Lucchesi 2001; 2002). In regard to the Earth solid tides, they have been neglected because their impact is several orders of magnitude smaller. For example, the amplitude of the nodal rate perturbation induced by the $`K_1`$ tide is five orders of magnitude smaller than that due to the even zonal harmonic coefficient $`J_2`$ of the geopotential for a GPS orbit.

$`𝒜`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}\mathrm{\Delta }𝒢_{\mathrm{}}a^{\left(\frac{3+2\mathrm{}}{2}\right)},`$ ($`38`$)
$``$ $`=`$ $`d{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}𝒢_{\mathrm{}}\left({\displaystyle \frac{3+2\mathrm{}}{2}}\right)a^{\left(\frac{5+2\mathrm{}}{2}\right)},`$ ($`39`$)
$`𝒞`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}𝒢_{\mathrm{}}a^{\left(\frac{3+2\mathrm{}}{2}\right)},`$ ($`40`$)
where $`\mathrm{\Delta }𝒢_{\mathrm{}}=𝒢_{\mathrm{}}^{(1)}𝒢_{\mathrm{}}^{(2)}`$ is the difference in the $`𝒢_{\mathrm{}}`$ functions of the two satellites induced by the inclinations and the eccentricities. The derivatives with respect to $`J_{\mathrm{}}`$, $`a`$ and $`d`$ are

$`{\displaystyle \frac{𝒜}{J_{\mathrm{}}}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}\left({\displaystyle \frac{\mathrm{\Delta }𝒢_{\mathrm{}}}{J_{\mathrm{}}}}\right)a^{\left(\frac{3+2\mathrm{}}{2}\right)},`$ ($`42`$)
$`{\displaystyle \frac{}{J_{\mathrm{}}}}`$ $`=`$ $`d{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}\left({\displaystyle \frac{𝒢_{\mathrm{}}}{J_{\mathrm{}}}}\right)\left({\displaystyle \frac{3+2\mathrm{}}{2}}\right)a^{\left(\frac{5+2\mathrm{}}{2}\right)},`$ ($`43`$)
$`{\displaystyle \frac{𝒞}{J_{\mathrm{}}}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}\left({\displaystyle \frac{𝒢_{\mathrm{}}}{J_{\mathrm{}}}}\right)a^{\left(\frac{3+2\mathrm{}}{2}\right)},`$ ($`44`$)
$`{\displaystyle \frac{𝒜}{a}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}\mathrm{\Delta }𝒢_{\mathrm{}}\left({\displaystyle \frac{3+2\mathrm{}}{2}}\right)a^{\left(\frac{5+2\mathrm{}}{2}\right)},`$ ($`45`$)
$`{\displaystyle \frac{}{a}}`$ $`=`$ $`d{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}𝒢_{\mathrm{}}\left({\displaystyle \frac{3+2\mathrm{}}{2}}\right)\left({\displaystyle \frac{5+2\mathrm{}}{2}}\right)a^{\left(\frac{7+2\mathrm{}}{2}\right)},`$ ($`46`$)
$`{\displaystyle \frac{𝒞}{a}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}𝒢_{\mathrm{}}\left({\displaystyle \frac{3+2\mathrm{}}{2}}\right)a^{\left(\frac{5+2\mathrm{}}{2}\right)},`$ ($`47`$)
$`{\displaystyle \frac{}{d}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{}}R^{\mathrm{}}𝒢_{\mathrm{}}\left({\displaystyle \frac{3+2\mathrm{}}{2}}\right)a^{\left(\frac{5+2\mathrm{}}{2}\right)}.`$ ($`48`$)
For $`\mathrm{}=2`$ and a GPS orbit, by assuming $`a=26578`$ km, $`i=55`$ deg and $`d=5`$ km the errors in $`\mathrm{\Delta }\xi `$ are888In the calculations it turns out that the effect of $`\mathrm{\Delta }𝒢_2`$ in $`𝒜`$ and $`\frac{𝒜}{a}`$, for $`\mathrm{\Delta }i=1`$ deg, can be neglected. On the other hand, with a good quality rocket launcher it is possible to insert two spacecrafts in the same orbital planes up to 10-4 deg, as in the case of the GRACE mission. See on the WEB http://www.csr.utexas.edu/grace/newsletter/2002/august2002.html

$`[\delta (\mathrm{\Delta }\xi )]_{\mathrm{\Delta }\dot{\mathrm{\Omega }}}`$ $`=`$ $`{\displaystyle \frac{\delta (\mathrm{\Delta }\dot{\mathrm{\Omega }})}{\sqrt{GM}(𝒞)}}=2\times 10^8(\mathrm{mas}/\mathrm{yr})^1\times \delta (\mathrm{\Delta }\dot{\mathrm{\Omega }}),`$ ($`49`$)
$`[\delta (\mathrm{\Delta }\xi )]_{GM}`$ $`=`$ $`{\displaystyle \frac{(\mathrm{\Delta }\dot{\mathrm{\Omega }})_{\mathrm{exp}}}{2\sqrt{(GM)^3}(𝒞)}}\times \delta (GM)=3\times 10^{20}(\mathrm{deg}/\mathrm{day})^1\times (\mathrm{\Delta }\dot{\mathrm{\Omega }})_{\mathrm{exp}},`$ ($`50`$)
$`[\delta (\mathrm{\Delta }\xi )]_{J_2}`$ $`=`$ $`\left|\left\{{\displaystyle \frac{\frac{(𝒜+)}{J_2}\times (𝒞)\frac{(𝒞)}{J_2}\times (𝒜+)}{(𝒞)^2}}\right\}\right|\times \delta J_2=2\times 10^9,`$ ($`51`$)
$`[\delta (\mathrm{\Delta }\xi )]_a`$ $`=`$ $`\left|\left\{{\displaystyle \frac{\frac{(𝒜+)}{a}\times (𝒞)\frac{(𝒞)}{a}\times (𝒜+)}{(𝒞)^2}}\right\}\right|\times \delta a=2\times 10^{13}\mathrm{cm}^1\times \delta a,`$ ($`52`$)
$`[\delta (\mathrm{\Delta }\xi )]_d`$ $`=`$ $`\left|\left\{{\displaystyle \frac{\frac{}{d}\times (𝒞+𝒜)}{(𝒞)^2}}\right\}\right|\times \delta d=1\times 10^9\mathrm{cm}^1\times \delta d.`$ ($`53`$)
It can be noticed that the major limiting factor is the term due to the error in the difference of the nodal rates $`\delta (\mathrm{\Delta }\dot{\mathrm{\Omega }})=\delta \dot{\mathrm{\Omega }}^{(1)}+\delta \dot{\mathrm{\Omega }}^{(2)}`$. Indeed, the experimental error in measuring the secular rate of the node is of the order of 1 mas yr-1. The systematic error in $`\dot{\mathrm{\Omega }}`$ due to the uncertainty on $`J_2`$ is, for a GPS satellite, almost 3 mas yr-1. In the case of the error in $`d`$, the same considerations as for the orbital periods hold.

## 4 Conclusions

In this paper we have shown that the comparison of the orbital motions of a pair of artificial satellites of different compositions moving along identical orbits in the gravitational field of Earth in order to test the Weak Equivalence Principle is not competitive with the already performed tests with torsion balances on Earth and the Lunar Laser Ranging technique, and the dedicated space–based missions STEP, GG, $`\mu `$SCOPE and SEE.

We have considered the orbital periods and the secular nodal rates. The analysis of the orbital periods seems to yield more precise measurements. The major limiting factor is represented by the difference in the orbital radiuses induced by the unavoidable orbital injection errors and the related uncertainty. By assuming $`\delta d1`$ cm or less the achievable precision is of the order of $`10^{10}`$$`10^{11}`$.


text:
# The Type Ic Hypernova SN 2003dh/GRB 030329

## 1 INTRODUCTION

Evidence that at least some Gamma Ray Bursts (GRB’s) are connected to Supernovae (SNe) is mounting. After the serendipitous discovery of SN 1998bw in coincidence with GRB980425 (Galama et al. 1998), several other cases of possible SNe in GRB’s have been reported (e.g., Bloom et al. 2002; Garnavich et al. 2003). All of these were however based only on the detection of ‘bumps’ in the GRB afterglows’ light curves, which could be decomposed to look like the light curve of SN 1998bw.

SN 1998bw was no ordinary SN. Its broad spectral features were explained as indicating a very energetic Type Ic SN (arising from the collapse of the bare CO core of a massive star). Because of its high explosion kinetic energy ($`5\times 10^{52}`$ erg in spherical symmetry), and the consequently very broad spectral lines, SN 1998bw was called a ‘hypernova’ (Iwamoto et al. 1998). Other hypernovae have been discovered and analyzed (e.g. SNe 1997ef (Iwamoto et al. 2000; Mazzali, Iwamoto, & Nomoto 2000); SN 2002ap (Mazzali et al. 2002)), but none was associated with a GRB. This may be related to the fact that none of these SNe were either as bright or as powerful as SN 1998bw. So far, the only other positive detection of a SN spectrum in a GRB afterglow is the case of GRB021211/SN 2002lt (della Valle et al. 2003).

Given this intriguing but insufficient evidence, excitement mounted when a very nearby GRB was detected (GRB030329, z=0.1685; Greiner et al. 2003), as a possible SN may be relatively easily observable in the light of the afterglow (AG). Indeed, the detection of broad spectral features characteristic of a supernova was reported by Stanek et al. (2003) and later by Hjorth et al. (2003). The SN (SN 2003dh) appeared to be similar to SN 1998bw. However, spectra at $`1`$ month had changed somewhat, looking more like those of the less energetic hypernova SN 1997ef (Kawabata et al. 2003).

There are significant differences between the SN 2003dh light curves of Matheson et al. (2003) and Hjorth et al. (2003). These may be due to the different observational methods and to major difficulties with subtraction of a) the underlying AG spectrum, which can change with time in an unknown way, and b) the host galaxy background. In particular, the Hjorth et al. (2003) light curve rises much more rapidly, reaches a brighter peak, and then drops much faster. The Matheson et al. (2003) light curve, on the other hand, has a slower rise, and resembles that of SN 1998bw. Unfortunately, neither light curve covers the likely time of peak, $`12`$ – 15 days after the GRB, because this time coincided with full moon. The absolute rest-frame $`V`$ magnitude at peak may have been between $`18.5`$ and $`19.1`$, depending on the dataset used and the estimated extinction. Bloom et al. (2003), however, find that SN 2003dh may have been as bright as $`M_V=19.8\pm 0.4`$ at peak. They used an extinction $`A_V0.2`$ toward SN 1998bw, explaining at least part of the discrepancy. Also, they did not use spectral information in decomposing the light curve.

## 2 Spectral Models

Light curve models alone cannot uniquely constrain the properties of a SN, as models yielding similar light curves may give rise to different synthetic spectra (e.g. Iwamoto et al. 1998). Fitting both light curves and spectra is a much more effective approach.

Unfortunately, in the case of SN 2003dh the exact shape of the light curve is not yet certain. Therefore, we derived fiducial SN spectra by rescaling the early power-law spectrum of the AG to the blue flux of the later spectra, where a SN signature was evident, and subtracting it off. Attributing the entire blue flux to the AG is justified by the fact that because of line-blanketing type Ib/c SNe, like SNe Ia, always show a flux deficiency to the blue of $`3600\mathrm{\AA }`$(e.g. Mazzali et al. 2000). At the third epoch the blue continuum indicative of the AG is very weak, and no subtraction was applied. We fitted three spectra of SN 2003dh using our Monte Carlo code (Mazzali & Lucy 1993; Lucy 1999; Mazzali 2000).

The first spectrum was obtained at the MMT on 2003 April 10 (Matheson et al. 2003). This is $`12`$ days after the GRB, i.e. $`10`$ rest-frame days into the life of the SN, assuming that the SN and the GRB coincided in time. The spectrum is characterized by very broad absorption lines, and it is similar to those of SN 1998bw at comparable epochs. Using the same explosion model as for SN 1998bw (model CO138E50; Table 1), a good match can be obtained for $`\mathrm{log}L=42.83`$ and $`v(\mathrm{ph})=28000`$ km s-1 (Figure 1). The synthetic spectrum has $`V=18.53`$, $`M(\mathrm{Bol})=18.37`$. While the luminosity is lower than in SN 1998bw, the photospheric velocity is comparable (for SN 1998bw we had $`v(\mathrm{ph})=31600`$ km s-1on day 8 and 20700 on day 16). This suggests that SN 2003dh ejected a significant amount of high-velocity material in our direction, which is understandable since we observed the GRB. SN 1998bw may have been even more energetic, but is was probably viewed less on-axis since its associated GRB was much weaker than GRB030329. Note that we cannot confirm the claim that SN 2003dh has higher-velocity ejecta than SN 1998bw. The model used for SN 1997ef (model CO100E18), on the other hand, has much too little mass at high velocities, and it does not yield a broad-lined spectrum similar to the observed one.

The next spectrum we attempted to fit is the MMT 2003 April 24 one (Matheson et al. 2003). The spectrum (rest-frame epoch $`23`$ days) still has broad lines, and it is almost as bright as the 10 April one, so that maximum probably occurred about half way between the two epochs. If model CO138E50 is used, the photosphere falls in the flat inner part of the density distribution (model CO138E50 ‘turns over’ at $`v20000`$ km s-1, which is above the position of the photosphere at this epoch). The rather flat density distribution just above the photosphere leads to very deep but narrow absorption lines, which are not like the observed spectrum. On the other hand, model CO100E18 still has too little mass at these velocities to give any significant spectral features for such large velocities at this relatively advanced epoch (at a similar epoch SN 1997ef had a low $`v(\mathrm{ph})=10000`$ km s-1). The spectrum is therefore in the ‘transition’ phase between SN 1998bw-like and SN 1997ef-like. Clearly, some kind of new model is required. Since it is not clear how such a model should behave, we deal first with the next epoch, which can yield indications about the inner part of the ejecta, and then come back to the April 24 spectrum to verify the results.

Our third spectrum was obtained with Subaru on 2003 May 10 (Kawabata et al. 2003), a rest-frame epoch $`36`$ days, and it resembles that of SN 1997ef at a comparable epoch. Indeed, model CO138E50 yields a narrow-lined spectrum, while a much better synthetic spectrum is obtained using model CO100E18. This model, in fact, turns over at $`v6000`$ km s-1, which is still below the photosphere at this epoch. However, model CO100E18 is too massive at the low velocities near the photosphere, resulting in strong backwarming and a high temperature. We tested different possibilities, and found that a model where the density of CO100E18 is divided by a factor of two throughout (CO100E18/2) gives a good fit using $`\mathrm{log}L=42.43`$ and $`v(\mathrm{ph})=7500`$ km s-1 (Figure 2). The synthetic spectrum has $`V=17.35`$, $`M(\mathrm{Bol})=17.37`$. At a comparable epoch the values for SN 1998bw were $`\mathrm{log}L=42.68`$ and $`v(\mathrm{ph})=7500`$ km s-1, and those for SN 1997ef $`\mathrm{log}L=42.14`$ and $`v(\mathrm{ph})=4900`$ km s-1. Therefore, it appears that while model CO138E50 is adequate for the high velocity part of the ejecta, a less energetic model such as CO100E18/2 is better suited for the low-velocity inner part. This might be expected given the shift in the properties of the spectrum.

The question is how can those two models be merged into one. We have adopted here an empirical approach. Model CO138E50 can be modified below $`v=28000`$ km s-1, since this is below the photosphere of the April 10 spectrum, while model CO100E18/2 can be modified at velocities high enough that the May 10 spectrum is not affected. We tried different boundaries for this, and finally chose $`v=15000`$ km s-1. Therefore, we used model CO100E18/2 for $`v<15000`$ km s-1, model CO138E50 for $`v>25000`$ km s-1, and merged the two models linearly in between. The density structures of the various models are shown in Figure 4. With this new spherically symmetric ‘merged’ model (COMDH) the results for April 10 and May 10 are essentially unchanged. COMDH has a smaller ejected mass than both CO138E50 and CO100E18, $`M_{\mathrm{ej}}`$$`=8M_{}`$, and $`E_{\mathrm{kin}}`$$`=3.8\times 10^{52}`$ erg, similar to CO138E50 as it is dominated by the high-velocity part. Since model COMDH is not computed from the hydrodynamics, but represents only a run of density vs. velocity, multi-dimensional calculations are necessary to see how well COMDH approximates a real multi-dimensional model.

The test for the merged model comes from fitting the April 24 spectrum, when the photosphere falls in the joining region. Using essentially the same parameters as discussed earlier ($`\mathrm{log}L=42.79`$ and $`v(\mathrm{ph})=18000`$ km s-1), we could obtain a reasonably good fit to the observed spectrum (Figure 3). The synthetic spectrum has $`V=18.17`$, $`M(\mathrm{Bol})=18.27`$. Therefore, considering all uncertainties, our ‘merged’ model is probably a fair one-dimensional representation of the density structure of SN 2003dh.

## 3 A Two-Component Light Curve Model

In order to verify that the model we have constructed to reproduce the spectra is indeed valid, it is necessary to compute a light curve with it. Here we can follow the same approach as in the two-component models of Maeda et al. (2003), as this is essentially the same problem. Since we do not know the exact light curve of SN 2003dh, we only tried to reproduce the three points derived from the spectral fitting. These suggest a more rapid rise, a dimmer peak and a faster decline than SN 1998bw, more similar to SN 2002ap.

We computed synthetic bolometric light curves for models CO138E50 and COMDH, using an LTE radiation transfer code and a gray $`\gamma `$-ray transfer code (Iwamoto et al. 2000). The optical opacity was found to be dominated by electron scattering.

The synthetic light curve of CO138E50 with homogeneous mixing of 56Ni out to 40,000 km s-1 is in reasonable agreement with the three points of SN 2003dh inferred from the spectra (Figure 5). The mass of 56Ni ($`0.32M_{}`$) is smaller than in the models for SN 1998bw ($`0.40.7M_{}`$: Nakamura et al. 2001; Iwamoto et al. 1998) and the mixing more extensive. The model light curve of SN 2003dh rises faster and peaks earlier than SN 1998bw.

The synthetic light curve computed with the ‘merged’ model COMDH reproduces the three SN 2003dh points well for a 56Ni mass of 0.35$`M_{}`$. Homogeneous mixing of 56Ni out to 40,000 km s-1 was assumed, as for CO138E50. High velocity 56Ni is a feature that SN 2003dh shares with other hypernovae. It could be attributed to a jet-like asphericity in the explosion (Maeda et al. 2002; Maeda & Nomoto 2003; Woosley & Heger 2003).

The main feature of the ‘merged’ model COMDH is the presence of a high density region at low velocities, as derived from the spectral analysis. The effect of this region should be to slow down the decline of the late-time light curve. In SN 1998bw and other hypernovae this is indeed observed starting at $`50`$ days and becoming more prominent later (Sollerman et al. 2000; Nomoto et al. 2001; Nakamura et al. 2001; Maeda et al. 2003). Figure 5 confirms that the synthetic light curve of model COMDH becomes somewhat brighter than that of CO138E50 at advanced phases. Although the ‘merged’ model is favored because it gives reasonably good fits to both the light curve and the spectra, a definitive extraction of the light curve and further time coverage will be useful to distinguish between the models.

## 4 Discussion

We have shown through spectral models that SN 2003dh was a hypernova. The high ejecta velocities observed in the early spectra require for the outer layers an explosion model similar to that used for SN 1998bw. However, at lower velocities such a model becomes too flat, and it is necessary to replace it with a lower energy one such as CO100, which has a steeper $`\rho (v)`$ at low velocities. From a spherically symmetric calculation, we obtain for SN 2003dh $`E_{\mathrm{kin}}`$$`=3.8\times 10^{52}`$ erg, $`M_{\mathrm{ej}}`$$`=8M_{}`$, and $`M`$(56Ni)$`=0.35M_{}`$.

Woosley & Heger (2003) reached similar results. Their model parameters differ from ours somewhat, but they did not make detailed comparisons with observed data.

Our results suggest that SN 2003dh was intermediate between SNe 1997ef and 1998bw in $`E_{\mathrm{kin}}`$ release and 56Ni production, and that it ejected a comparable, possibly slightly smaller mass. Accounting for the compact remnant, the exploding star may have been a CO core of $`1011M_{}`$, which implies a main-sequence mass of $`3540M_{}`$. SN 2003dh seems to follow the relations between progenitor mass, $`M`$(56Ni) and $`E_{\mathrm{kin}}`$ (Nomoto et al. 2003).

All spherically symmetric models (including CO100 and CO138) are very flat inside, and they all have an inner mass-cut defining a ‘hole’ in the density profile. Evidence for SNe 1997ef, 2002ap and now 2003dh is that such a flat density distribution does not lend itself to successful spectrum synthesis. Actually, synthetic spectra computed for SN 1998bw beyond the epochs published in Iwamoto et al. (1998) are also not perfect, so this may be a common feature. Evidence for a slow-moving, oxygen-dominated inner region in SN 2003dh, to be obtained from the nebular spectra would nicely confirm this picture.

The difficulties one-dimensional explosion models encounter in reproducing the time evolution of hypernovae are almost certainly due to the inadequacy of such models when applied to phenomena that are intrinsically aspherical. Approximate studies suggest that asymmetries in the explosions may affect the light curves (Höflich, Wheeler, & Wang 1999; Woosley & Heger 2003). Therefore, the estimate of the explosion parameters may be subject to errors. In particular, $`E_{\mathrm{kin}}`$ may be overestimated by factors of a few by neglecting the fact that the highest velocity material is only ejected in a narrow cone, as discussed by Maeda et al. (2003), and recently by Woosley & Heger (2003) starting from the collapsar model. Both papers point at a similar scenario. According to Woosley & Heger (2003), if the high velocity ejecta are contained in a cone with opening half-angle $`\theta `$, $`E_{\mathrm{kin}}`$ is reduced by a factor $`(1\mathrm{cos}\theta )`$, and M(56Ni) required to fit the peak luminosity is also reduced by a factor $`(1\mathrm{cos}\theta )/\mathrm{sin}^2\theta `$ (in an extreme case). This is $`0.5`$ if $`\theta `$ is small. The value of $`M`$(56Ni) can be derived from observations in the nebular phase, which is basically independent of geometry. At least for SN 1998bw, these confirmed that $`M`$(56Ni) estimated in spherical symmetry was correct within a factor of 2 or less (Nakamura et al. 2001; Maeda et al. 2002). Early-time spectra are less likely to be affected than the light curve, especially in the case of SN 2003dh, where the high velocity material which dominates most spectral features is probably coming directly towards us. Nebular observations of SN 2003dh are urged.

This work was partly supported by the Japanese Ministry of Education, Science, Culture, Sports, and Technology under Grants-in-Aid 14047206, 14540223, 15204010.


text:
# Chiral perturbation theory for lattice QCD at 𝒪(𝑎²)Talk by O. Bär at Lattice 2003, Tsukuba.

## 1 Introduction

Chiral perturbation theory ($`\chi `$PT) is widely used in analyzing lattice QCD data. Since $`\chi `$PT describes continuum QCD one has to take the continuum limit first before fitting the lattice data to the analytic predictions of $`\chi `$PT. Unfortunately, lattice data for various lattice spacings is often not available in order to perform a reliable continuum extrapolation. Nevertheless, it is common practice to directly fit lattice data obtained for non-vanishing lattice spacing to continuum $`\chi `$PT, thus introducing a systematic error.

This uncertainty can be analytically controlled by including the discretization effects stemming from a non-zero lattice spacing in $`\chi `$PT. For example, the $`𝒪(a^2)`$ taste violations in staggered fermion simulations have been systematically accounted for that way . Fits of MILC simulation data to $`\chi `$PT including $`𝒪(a^2)`$ taste violating terms show a significantly improved $`\chi ^2`$ . In fact, good fits are not possible without these terms.

Similarly, the leading effects linear in $`a`$ for simulations with Wilson fermions have been included in $`\chi `$PT  and the resulting expressions for the pseudo scalar mass and decay constant have already been used in fits to simulation data . However, since the lattice spacings in current unquenched simulations with Wilson fermions are usually not very small, it is important to extend the results in  by including the $`𝒪(a^2)`$ contributions in $`\chi `$PT .

## 2 $`\chi `$PT for Wilson fermions at $`𝓞\mathbf{(}𝒂^\mathrm{𝟐}\mathbf{)}`$

The strategy for including the cut-off artifacts in $`\chi `$PT is essentially a two-step matching to effective theories. We first match the lattice theory to Symanzik’s effective theory , which is a continuum field theory. This effective theory makes the cut-off dependence of the underlying lattice theory explicit in terms of higher dimensional operators in the effective action and the effective operators. In a second step we construct the chiral effective theory for Symanzik’s effective theory, including the extra terms due to the non-zero lattice spacing. The result is a chiral expansion with explicit dependence on the lattice spacing $`a`$.

In order to include $`𝒪(a^2)`$ effects in $`\chi `$PT we need Symanzik’s effective action for lattice QCD with Wilson fermions through $`𝒪(a^2)`$. At this order there are fifteen dimension-6 operators which are compatible with the symmetries of the underlying lattice theory . The way these operators affect the chiral Lagrangian depends on their transformation properties under chiral rotations, and one can basically distinguish three types of operators:

1. the O(4) symmetry breaking term
$`\overline{\psi }\gamma _\mu D_\mu D_\mu D_\mu \psi `$
2. chiral symmetry conserving terms
e.g. $`\overline{\psi }\mathit{}^3\psi `$ and $`(\overline{\psi }\gamma _\mu \psi )(\overline{\psi }\gamma _\mu \psi )`$
3. chiral symmetry breaking terms
e.g. $`(\overline{\psi }\gamma _5\psi )(\overline{\psi }\gamma _5\psi )`$

The first operator, which breaks the full O(4) rotation symmetry, is not forbidden since the underlying lattice theory is invariant only under the hypercubic subgroup of O(4).

The chiral Lagrangian is expanded in powers of momenta $`p^2`$, quark masses $`m`$ and the lattice spacing $`a`$. Generalizing the standard chiral power counting, the leading order (LO) chiral Lagrangian contains the terms of $`𝒪(p^2,m,a)`$, while the terms of $`𝒪(p^4,p^2m,p^2a,m^2,ma,a^2)`$ are of next-to-leading order (NLO). Employing this power counting the first two types of operators given above do not contribute at NLO.

The first operator, for instance, gives rise to terms in the chiral Lagrangian which break the O(4) symmetry and are invariant under the hypercubic subgroup only. These terms necessarily contain four derivatives and, since they are also multiplied by $`a^2`$, are therefore at least of $`𝒪(p^4a^2)`$. These terms are beyond NLO and can be safely neglected.

The second type of operators, the chiral symmetry conserving terms, modify the coefficients in front of already existing terms in the chiral Lagrangian. For example, the LO coefficient $`f^2`$, multiplying the kinetic term $`\text{tr}_\mu \mathrm{\Sigma }\mathrm{\Sigma }_\mu ^{}`$, becomes $`f^2(a^2)=f^2+a^2K+\mathrm{}`$, with some other (unknown) coefficient $`K`$. This leads to the extra term $`a^2K\text{tr}_\mu \mathrm{\Sigma }\mathrm{\Sigma }_\mu ^{}`$ in the chiral Lagrangian, which is of $`𝒪(a^2p^2)`$ and thus beyond NLO. This example can be generalized to all chiral symmetry conserving terms and we conclude that these terms do not affect the NLO chiral Lagrangian.

Only some of the chiral symmetry breaking terms give rise to new operators in the NLO chiral Lagrangian. The way they influence the chiral Lagrangian is obtained by a standard spurion analysis and the final result for the contributions proportional to $`a^2`$ reads

$`[a^2]`$ $`=`$ $`\widehat{a}^2W_6^{^{}}\mathrm{\Sigma }^{}+\mathrm{\Sigma }^2\widehat{a}^2W_7^{}\mathrm{\Sigma }^{}\mathrm{\Sigma }^2`$ (1)
$`\widehat{a}^2W_8^{}\mathrm{\Sigma }^{}\mathrm{\Sigma }^{}+\mathrm{\Sigma }\mathrm{\Sigma }.`$
Here $`\mathrm{\Sigma }`$ is the standard matrix of Nambu-Goldstone fields and the angled brackets denote traces over the flavor indices. The parameter $`\widehat{a}=aW_0`$ absorbs the leading order low-energy constant $`W_0`$ in order to render the new unknown low-energy constants $`W_6^{^{}},W_7^{^{}},W_8^{^{}}`$ dimensionless. Note that only three new operators appear at $`𝒪(a^2)`$ despite the fact that there are many different operators in the underlying Symanzik action at this order.

We are now in the position to compute NLO expressions for physical observables. For example, the $`𝒪(a^2)`$ contribution to the flavor non-singlet pseudo-scalar mass is given by

$$\delta M_{PS}^2[a^2]=\frac{16}{f^2}\widehat{a}^2(N_fW_6^{^{}}+W_8^{^{}}),$$
(2)
where $`N_f`$ denotes the number of flavors. Adding this contribution to the terms of $`𝒪(a)`$ given in  we obtain the complete NLO result through $`𝒪(a^2)`$ for the pseudo-scalar mass.

## 3 Partially quenched QCD

The previous discussion is readily extended to partially quenched lattice QCD with different masses for the sea and valence Wilson fermions. Described by a lattice action with sea, valence and ghost quarks, the Symanzik action through $`𝒪(a^2)`$ for partially quenched lattice QCD is obtained as in the unquenched case, based on locality and the symmetries of the lattice theory. The result for the $`𝒪(a^2)`$ contribution to the partially quenched chiral Lagrangian has the same form as in (1), with the angled brackets now denoting super-traces and the field $`\mathrm{\Sigma }`$ reflecting the larger flavor content of partially quenched $`\chi `$PT.

The discussion of “mixed fermion theories” is slightly more complicated. Mixed theories are a generalization of partially quenched lattice theories. In addition to choosing different quark masses, the lattice Dirac operator is different in the sea and valence sector. The concrete example with the Wilson Dirac operator for the sea quarks and a Dirac operator satisfying the Ginsparg-Wilson relation for the valence quarks has been proposed in .

Such a combination has various advantages. As a consequence of the Ginsparg-Wilson relation the valence sector exhibits an exact chiral symmetry and therefore does not suffer from the shortcomings of Wilson fermions due to their explicit chiral symmetry breaking (additive mass renormalization, cut-off effects linear in $`a`$, etc.). In particular, simulations are possible with valence quark masses much smaller than accessible with Wilson fermions . This should allow numerical simulations deeper in in the chiral regime of partially quenched QCD. Finally, mixed fermion simulations offer a cost-effective compromise towards full unquenched simulations using Ginsparg-Wilson fermions, since mixed simulations only require the computation of correlators in the background of already existing unquenched configurations generated with Wilson fermions.

Deriving the $`𝒪(a^2)`$ corrections in the chiral Lagrangian for the mixed fermion theory with Wilson sea and Ginsparg-Wilson valence quarks goes along the lines outlined above. However, due to the use of Ginsparg-Wilson fermions for the valence quarks there are no $`𝒪(a)`$ terms and fewer $`𝒪(a^2)`$ terms in the valence sector of Symanzik’s effective action for the mixed theory. As before one can show that the O(4) symmetry breaking terms as well as the chiral symmetry conserving operators in Symanzik’s action do not contribute to the NLO chiral Lagrangian. A spurion analysis for the chiral symmetry breaking terms gives the following $`𝒪(a^2)`$ contribution for the chiral Lagrangian:

$`[a^2]`$ $`=`$ $`\widehat{a}^2W_6^{^{}}P_S\mathrm{\Sigma }^{}+\mathrm{\Sigma }P_S^2`$ (3)
$`\widehat{a}^2W_7^{}P_S\mathrm{\Sigma }^{}\mathrm{\Sigma }P_S^2`$
$`\widehat{a}^2W_8^{}P_S\mathrm{\Sigma }^{}P_S\mathrm{\Sigma }^{}+\mathrm{\Sigma }P_S\mathrm{\Sigma }P_S`$
$`\widehat{a}^2W_M\tau _3\mathrm{\Sigma }\tau _3\mathrm{\Sigma }^{}`$
The operators $`P_S`$ and $`P_V`$ denote projection operators to the sea and valence sector, respectively, and $`\tau _3`$ equals $`P_SP_V`$ (as before, $`\widehat{a}=aW_0`$).

The first three terms stem entirely from the Wilson fermion sea sector. They involve sea sector fields only (note the projector $`P_S`$) and they are the same three operators as in (1). The additional operator mixes the sea and valence sector. The reason for the presence of the extra term is the reduced flavor symmetry group of the mixed theory. The use of different Dirac operators for sea and valence quarks does not contain transformations between sea and valence fields and the additional term $`\tau _3\mathrm{\Sigma }\tau _3\mathrm{\Sigma }^{}`$ is not forbidden by the symmetries of the theory.

The presence of this extra term in the chiral Lagrangian does not entail that chiral expressions for observables in the mixed theory depend on more free parameters. For example, taking account of the terms in (3) in the calculation of the flavor non-singlet pseudo-scalar mass one finds that its $`𝒪(a^2)`$ correction actually vanishes. This is not that surprising. Exact chiral symmetry in the massless case for Ginsparg-Wilson fermions implies that the pseudo-scalar meson mass is proportional to the quark mass and vanishes in the chiral limit. Hence any lattice spacing correction to $`M_{PS}^2`$ is suppressed by at least one factor of the valence quark mass $`m_{Val}`$. The largest lattice correction quadratic in the lattice spacing is therefore of $`𝒪(m_{Val}a^2)`$. This example illustrates that certain beneficial properties of Ginsparg-Wilson fermions are preserved even in the presence of a “non-Ginsparg-Wilson” sea sector.