text:
# CERN-TH/2001-288 IFIC/01-59, FTUV/011030 Ginsparg-Wilson fermions: practical aspects and applications

## 1 Introduction

It has been established that lattice Dirac operators $`D`$ can be constructed which are local, do not suffer from the doubling problem and satisfy the Ginsparg–Wilson (GW) relation :

$`\{D,\gamma _5\}=D\gamma _5D\{D^1,\gamma _5\}=\gamma _5.`$ (1)
This relation ensures that Ward identities associated with a standard chiral transformation are satisfied on shell at finite $`a`$ . Furthermore, this relation implies an exact lattice chiral symmetry  and index theorem .

Two apparently disconnected avenues led to operators satisfying the GW relation. The first one was the lattice version of the five dimensional (5D) Dirac fermion Yukawa-coupled to a domain wall (DW), which naturally leads to a chiral fermion in four dimensions, attached to the defect. On a lattice such a construction involves in general two lattice spacings $`a,a_s`$ (in the physical four and the extra dimensions respectively). If the number of sites in the extra dimension, $`N_s`$, is infinite, one still expects some kind of discretized chiral fermion by analogy with the continuum analysis . Indeed it can be shown that in this limit, at finite lattice spacing, the effective action describing the light fermionic degrees of freedom is that of a 4D GW fermion . For the standard DW construction of , the 4D Dirac operator describing the light degrees of freedom at finite $`N_s`$ is

$`aD_{N_s}=1\gamma _5{\displaystyle \frac{(1+\overline{Q})^{N_s}(1\overline{Q})^{N_s}}{(1+\overline{Q})^{N_s}+(1\overline{Q})^{N_s}}},`$
$`\overline{Q}\gamma _5{\displaystyle \frac{a_sA}{2a_sA}}Am_0D_W.`$ (2)
$`m_0`$ is the height of the domain wall and $`D_W`$ is the standard Wilson Dirac operator. Formally

$`\underset{N_s\mathrm{}}{lim}aD_{N_s}1\gamma _5\mathrm{sign}(\overline{\mathrm{Q}})=1\overline{\mathrm{Q}}/\sqrt{\overline{\mathrm{Q}}^2},`$ (3)
which satisfies the GW relation. If also the limit $`a_s0`$ is taken, the overlap operator  is obtained111Note that if $`m_0a1`$, the operator of eq. (1) is not properly normalized. For simplicity, the formulae below assume $`m_0a=1`$ but are easily generalized.:

$`\underset{a_s0}{lim}\underset{N_s\mathrm{}}{lim}aD_{N_s}1\gamma _5\mathrm{sign}(\mathrm{Q}),\mathrm{Q}=\gamma _5\mathrm{A}.`$ (4)
The locality of the operators of eqs. (3) and (4) was shown in .

The second avenue was the search for lattice Dirac operators obtained from a renormalization group (RG) blocking of the continuum one . Starting from a chirally symmetric Dirac operator in the continuum, the effective operator on the lattice is not chirally symmetric in general, but retains a remnant of the original symmetry encoded in the GW relation. A fixed-point (FP) Dirac operator  is invariant under a RG transformation and, thus, not only does it satisfy the GW relation , but it is in principle free of discretization errors. In practice, however, such an operator is very difficult to find and no explicit construction exists.

In this review I concentrate on the uses of GW fermions for lattice QCD. For a review on other topics including chiral gauge theories see .

## 2 Chiral symmetry at finite $`a`$

Given that it is very expensive to simulate GW fermions, it is important to assess the need for an exact chiral symmetry at finite $`a`$ in lattice QCD. There are at least three very good reasons to preserve such a symmetry.

1. The renormalization of fermion operators is enormously simplified by the exact chiral symmetry. For any Dirac operator $`D`$ satisfying the GW relation of eq. (1), we can define the chiral projections

$`\mathrm{\Psi }_{R,L}=\widehat{P}_\pm \mathrm{\Psi },\overline{\mathrm{\Psi }}_{L,R}=\overline{\mathrm{\Psi }}P_\pm `$ (5)
with $`P_\pm (1\pm \gamma _5)/2`$ and $`\widehat{P}_\pm (1\pm \gamma _5(1aD))/2`$, so that the action density splits as

$`\overline{\mathrm{\Psi }}D\mathrm{\Psi }=\overline{\mathrm{\Psi }}_LD\mathrm{\Psi }_L+\overline{\mathrm{\Psi }}_RD\mathrm{\Psi }_R,`$ (6)
implying an exact $`SU(N_f)_L\times SU(N_f)_R`$ flavour symmetry under which the chiral components transform in the standard way:

$`\mathrm{\Psi }_LV_L\mathrm{\Psi }_L,\mathrm{\Psi }_RV_R\mathrm{\Psi }_RV_{L,R}SU(N_f)_{L,R}.`$
Quark masses should be included by adding $`\overline{\mathrm{\Psi }}_Lm\mathrm{\Psi }_R+\overline{\mathrm{\Psi }}_Rm^{}\mathrm{\Psi }_L`$ to the action. The usual spurion analysis implies that operators can mix only if they have the same transformation properties under the exact lattice symmetry

$`\mathrm{\Psi }_LV_L\mathrm{\Psi }_L\mathrm{\Psi }_RV_R\mathrm{\Psi }_RmV_LmV_R^{}.`$ (7)
Several important results follow in a straightforward way from this analysis:

$``$ there is a conserved axial current and no additive mass renormalization;

$``$ two and four-fermion operators relevant to weak matrix elements can be constructed (using these chiral components) in chiral representations of the symmetry group of eq. (7) in such a way that their mixing is the same as in any chirally symmetric continuum regularization.

Recently the perturbative renormalization constants for quark bilinears and four-fermion operators have been computed for both DW and overlap fermions .

2. $`O(a)`$ improvement is almost automatic. A similar spurion analysis as that outlined above shows that there are no invariant operators of dimension $`d=5`$ so the action is $`O(a)`$-improved off shell. Similarly the quark bilinears are easily $`O(a)`$-improved if they are constructed in the same chiral representations as their continuum counterparts .

3. The exact chiral symmetry makes it possible to explore the regime of the light quark masses. Naively one would think that this is anyway very difficult since the pion Compton wavelength would not fit in the lattice volumes accessible to simulations: it would not be possible to ensure the condition $`M_\pi L1`$ resulting in large finite volume effects. Fortunately (quenched) chiral perturbation theory (Q)$`\chi `$PT predicts these effects , including those of topological zero modes, and a finite size scaling analysis in the regime $`M_\pi L1`$, is actually a very useful tool to extract the physical low energy constants. Recently this technique has been applied to quenched QCD using GW fermions  to extract the low energy constant $`\mathrm{\Sigma }`$.

## 3 Implementation of GW fermions

Overlap fermions

The difficulty in the implementation of overlap fermions is in constructing the square root or sign function of eq. (4). Several methods have been proposed to approximate these functions, which are variants of two basic approaches:

$``$ Polynomial approximations (PA), e.g. Chebyshev or Legendre :

$`D_N1\gamma _5QP_N(Q^2)`$

$``$ Rational approximations (RA) :

$`D_N1\gamma _5Q\left(c_0+{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \frac{c_k}{Q^2+q_k}}\right).`$ (8)
In this case, the $`N`$ inversions $`(Q^2+q_k)^1`$ can all be obtained in one go, using a multimass solver conjugate gradient (CG) . The nested CG necessary to invert the resulting operator can be avoided by rewriting the RA as a truncated continued fraction, which in turn can be rewritten as a Gaussian path integral with extra fermionic fields . The problem reduces to inverting a local 5D operator very similar to the DW.

Naively, the cost increase with respect to Wilson fermions is the number, $`n`$, of additional $`Q`$ matrix$`\times `$vector multiplications that are needed to approximate the square root: $`2N`$ for the PA, and twice the number of CG iterations needed to invert the factor $`Q^2+q_1`$, with the minimum pole, for RA. Asymptotically, the convergence of the approximation is exponential in $`n`$ in both cases:

$`D_ND_{ov}<e^{\sqrt{ϵ}n},ϵ=\lambda _{min}(Q^2)/\lambda _{max}(Q^2),`$
but it is controlled by the condition number of the matrix $`Q^2`$. Figure 1 shows the distribution of the four lowest eigenvalues of $`Q^2`$ for three values of $`\beta `$. At large $`\beta `$ the spectrum shows a gap, which shrinks with lowering $`\beta `$ and becomes populated by isolated eigenvalues that can be orders of magnitude smaller. The presence of these eigenvalues results in a very poor convergence of any approximation to the sign function.

DW fermions

The cost of DW fermions is less obvious. First, there is a factor of $`N_s`$ more computations involved owing to the extra dimension. The asymptotic convergence to the limit $`N_s\mathrm{}`$ is very similar to the convergence in $`n`$ of the approximations to the overlap: it is exponential in $`N_s`$, with a rate determined by the lowest eigenvalue of the operator $`\overline{Q}^2`$ of eq. (2). The distribution of the low-lying eigenvalues of $`\overline{Q}^2`$ is very similar to that of $`Q^2`$ in Fig. 1 , implying that the convergence in $`N_s`$ also becomes extremely slow at accessible $`\beta `$’s. It is customary to quantify the deviations from the $`N_s=\mathrm{}`$ limit by the anomalous term in the axial Ward identity properly normalized, i.e. the so-called residual mass. The dependence of the residual mass on $`N_s`$ has been extensively studied numerically , and indeed it is found that its exponential decrease with $`N_s`$ is very slow. A recent study showed that the small eigenvalues of $`Q^2`$ are responsible for the residual mass at large $`N_s`$ . An additional cost of DW fermions comes from the fact that the number of CG iterations needed to invert the 5D operator is considerably larger than for the 4D Wilson operator. Although a more detailed comparison is needed, there is at present no evidence that the DW formulation is more efficient than the overlap at some fixed accuracy of the GW relation .

FP Dirac operator

In recent years a big effort has been devoted to find, among the Dirac operators with all the allowed hypercubic couplings, the one that approximates the closest either the FP operator or simply the GW relation. Fig. 2 shows the free spectrum of four such operators presented in refs.  (ordered from left to right and up to down). The spectrum should lie on a circle if the approximation were perfect. These operators are also considerably more costly than Wilson, e.g. a factor $``$ 20 for an accuracy of the GW relation of about $`1\%`$ .

## 4 Improvements

As we have seen, the cost of the implementation of overlap or DW fermions is inversely proportional to the square root of the gap of the operators $`Q^2`$ or $`\overline{Q}^2`$ respectively. Several ideas to improve the convergence in this respect have been proposed and studied:

$``$ A small eigenspace of $`O(10)`$ eigenvalues is treated exactly in the limit $`n,N_s\mathrm{}`$.

In the context of the overlap this means that the exact operator is computed in this eigenspace, while a PA or RA is used only in the space orthogonal to it . The convergence is then controlled by the smallest eigenvalue of the orthogonal subspace, which can be orders of magnitude larger than the smallest eigenvalue of $`Q^2`$. Analogously, in the case of the DW operator the challenge is to construct a local 5D operator such that its effective 4D action is $`\stackrel{~}{D}_{N_s}=D_{N_s}+\mathrm{\Delta }`$, where $`D_{N_s}`$ is the operator in eq. (2), $`\overline{Q}|\mathrm{\Phi }_i=\mu _i|\mathrm{\Phi }_i`$ and

$`a\mathrm{\Delta }\gamma _5{\displaystyle \underset{i=1}{\overset{p}{}}}\left(\mathrm{sign}(\mu _i){\displaystyle \frac{\alpha _i^+\alpha _i^{}}{\alpha _i^++\alpha _i^{}}}\right)|\mathrm{\Phi }_i\mathrm{\Phi }_i|`$
with $`\alpha _i^\pm =1\pm \mu _i`$. The convergence in $`N_s`$ of $`\stackrel{~}{D}_{N_s}`$ is controlled by $`ϵ=|\mu _{p+1}||\mu _1|`$. There is no unique 5D operator satisfying this criterion. At least two explicit methods have been described in refs. . Note that the condition number of the 5D operator might depend on the implementation.

$``$ The gap of $`Q^2`$ and $`\overline{Q}^2`$ is increased by improving the gauge action.

This is a cheap solution in principle. A comparison of different actions in the DW construction has been carried out in . Fig. 3 shows the time dependence of the residual mass for the Wilson action, the one-loop and tadpole-improved Symanzik action and two RG-improved actions: Iwasaki’s and DBW2 . The improvement of the DBW2 action is remarkable. Some understanding of this effect comes from the fact that this action strongly supresses small (and not so small) instantons . The overimprovement of the DBW2 action in particular is very large, resulting in a very slow thermalization of the topological charge . It will be important to clarify if the use of these gauge actions results in a significant improvement also when the method of treating exactly a small eigenspace is implemented.

$``$ Substitute $`D_W`$ by an improved operator in the overlap or DW constructions.

If $`D_W`$ in eq. (2) is substituted by an operator that is closer to satisfy the GW relation, the approximation of the overlap operator or the limit $`N_s\mathrm{}`$ of the DW are much faster, since the size of the gap of the corresponding $`Q^2`$ (or $`\overline{Q}^2`$) is close to $`1`$ . This idea has been implemented in the overlap context with the operators of Fig. 2. The results show that in general the gain in the convergence of the approximations of the sign function does not quite compensate the overhead of using these actions. On the other hand, the overlap operator is more local (at least in those cases where fat links are not involved) and also the scaling properties are improved in the case of the truncated FP actions (see results in ).

Finally, there is room for algorithmic improvement. The authors of  pointed out that, in the context of the RA, the vector updates in the multimass CG solver are costly for a large number of poles. They proposed a new RA which requires less poles than previously used ones and from which they can derive rigorous error bounds. The authors of  proposed a strategy to search for truncated continued fractions that result in better conditioned 5D operators.

## 5 A few physics results

In the past year a number of QCD simulations with GW fermions have been performed. New results on the weak matrix elements using standard DW fermions have been presented. I refer to  for a description of these results. An important step forward this year has been the demonstration that simulations with overlap fermions are also feasible in the quenched approximation. Results on hadron spectroscopy and low energy constants have been presented by several groups using both the standard overlap operator  and also improved ones .

Pseudoscalar and scalar two-point functions have been computed at $`a=0.090.16fm`$ in various physical volumes. Fig. 4 shows the pion mass square versus the bare quark mass in units of $`r_0`$ extracted from the pseudoscalar correlator (full symbols) or the pseudoscalar minus the scalar one (open symbols) 222This combination is expected to be less sensitive to topological zero modes for small quark masses..

Within errors the behaviour is linear and consistent with a vanishing intercept, as expected from the exact chiral symmetry of these actions. Assuming lowest order Q$`\chi `$PT, the low energy constant $`\mathrm{\Sigma }`$ can be extracted from the slope. The needed renormalization factor, $`Z_m^1=Z_S=Z_P`$, has been obtained non-perturbatively . This same renormalization factor can be used to obtain the renormalized $`\mathrm{\Sigma }`$ from a finite-size scaling analysis in a completely different regime . The agreement of these two quantities is very good . Also a value for the strange quark mass and decay constant, $`F_K`$, can be obtained from these measurements. Note that the latter needs no renormalization in GW regularizations. Figs. 5 show the comparison of the overlap results for these quantities with earlier ones obtained with DW fermions  and $`O(a)`$-improved Wilson fermions . These results show no surprises and, although the statistical errors are still large, they seem to indicate rather small $`O(a^2)`$ discretization errors for overlap fermions.

## 6 Conclusions

Having an exact chiral symmetry at finite $`a`$ in lattice QCD permits the exploration of the regime of light quark masses and might turn out to be essential in solving long-standing problems, such as the $`\mathrm{\Delta }I=1/2`$ rule or the prediction of $`ϵ^{}/ϵ`$. Unfortunately, it is also very costly and a simple rule of thumb gives about 2 orders of magnitude overhead with respect to Wilson fermions. Several ideas to improve the situation have been proposed and implemented, but they have not yet resulted in a significant cost improvement in logarithmic terms. In spite of this situation, several quenched simulations have been successfully carried out with GW fermions and show that these actions do not have large discretization errors of $`O(a^2)`$ and can reach a regime of light quark masses which is not accessible to Wilson fermions.


format: text

text:
# A conservative semi-Lagrangian method for oscillation-free computation of advection processes

## 1 Introduction

The semi-Lagrangian method is a category of numerical methods for partial differential equations including an advection (or convection) term, and has been used mainly for atmospheric and geographic problems . In the semi-Lagrangian approach, the advection of a physical quantity is solved as a transportation and interpolation problem on fixed grids. The trajectory of the Lagrangian invariant that will reach a grid point is computed backwardly by the time integration of the advection velocity, and the invariant at the departure point, which is determined by an interpolation within a computational cell, is transported to the grid point. One can use a CFL-free time step in this approach. Since there are some well-established methods for calculating the trajectory such as the Runge-Kutta methods, recent studies on this approach are focused on the construction of interpolation functions.

One of problems of this approach is the difficulty of constructing a completely conservative scheme. There are only a few studies attacking to this problem. In Ref.  Priestley proposed a quasi-conservative semi-Lagrangian method by coupling a flux-corrected transport method and a high-order interpolation. In Ref. , Leonard et al. proposed a conservative, explicit algorithm achieving stable high-CFL computations, by using an integral variable of a Lagrangian invariant as the dependent variable, and showed results for one-dimensional advection at constant velocity by coupling the algorithm with many interpolation schemes. In Ref. , Manson and Wallis modified the QUICKEST scheme  so that the stable high-CFL computation is achieved. Applications of their method to pure advection problems in non-uniform flow fields are shown in Ref. . In Ref.  Lin and Rood discussed the use of readymade Eulerian schemes (the MUSCL  and PPM  schemes) modified to achieve stable computations even with a large time step.

Recently, Tanaka et al. proposed a way to solve this difficulty. In Ref. , they modified the CIP method, which is a non-conservative semi-Lagrangian method using the Hermite cubic interpolation function , into a conservative one by employing the cell-integrated quantity of a Lagrangian invariant as an additional dependent variable. A 4th- or 2nd-order polynomial is used as an interpolation function, and is constructed using a physical quantity, its first spatial derivative (only for the 4th-order one), and its cell-integrated quantity. The conservation of total mass is achieved by updating the cell-integral quantity in a conservation sense. The CIP method and its conservative variant using the 2nd-order polynomial (the CIP-CSL2 method) are briefly reviewed in the next section.

The aim of this paper is to propose a variant of the CIP-CSL2 method by employing the hybrid-cubic-rational (HCR) interpolation function introduced recently by the author . The hybrid function is constructed with a combination of the cubic polynomial and a rational function, and is proposed to obtain oscillation-free results having higher resolution than those given by the rational method  being a variant of the CIP. In the results using the CIP-CSL2 method one can sometimes observe numerical oscillation caused by the use of the classical polynomial . As is well known, the numerical oscillation becomes a serious problem for a certain examples such as the propagation of shock waves and the multiphase flow problems including large density jump. In Ref.  the use of the rational function employed in the rational method is suggested to overcome the numerical oscillation. The author expects that combining the hybrid function with the concept of the CIP-CSL can derive an oscillation-free, conservative semi-Lagrangian method that has higher resolution than that by the rational method. In Sec. 3 we show that the hybrid method can easily be modified into a conservative one, and in Sec. 4 we demonstrate the accuracy and stability of the present method. In this paper, we discuss only the case of $`\mathrm{CFL}\mathrm{number}1`$; however, we believe that the present method is applicable to a high-CFL condition by incorporating the concept discussed in Ref. .

## 2 CIP method and its conservative variants

The CIP method  is a semi-Lagrangian solver for an advection equation,

$$\frac{f}{t}+u\frac{f}{x}=0,$$
(1)
where $`f`$ is a dependent variable and $`u`$ is the advection velocity. The fact that the solution of this equation can be represented as

$$f(x,t+\mathrm{\Delta }t)=f(X(x,t),t)$$
(2)
allows us to solve Eq. (1) as an interpolation problem. Where $`X`$ is the trajectory of fluid particle that locates at $`x`$ at the time $`t+\mathrm{\Delta }t`$,

$$X(x,t)=x+_t^{t\mathrm{\Delta }t}u(X(x,\tau ),\tau )𝑑\tau .$$
(3)
In the CIP method, the dependent variable $`f`$ is interpolated by the Hermite cubic function,

$$C(X)=f_i+d_i(Xx_i)+C2_i(Xx_i)^2+C3_i(Xx_i)^3,$$

constructed using both $`f_i`$ and its derivative $`d_i`$ ( $`=f/x|_i`$) defined at each grid points $`x_i`$, where $`i=1,2,\mathrm{},N`$, $`N`$ is the number of grids, and $`C2_i`$ and $`C3_i`$ are determined from the condition of continuity . Equation (1) and its spatial derivative, i.e.,

$$\frac{f(x,t+\mathrm{\Delta }t)}{x}=\frac{X(x,t)}{x}\frac{f(X(x,t),t)}{X},$$
(4)
are used to update $`f_i`$ and $`d_i`$, respectively. In the case where we assume the first-order accuracy in time for example, Eqs. (2) and (4), respectively, are reduced to

$$f(x,t+\mathrm{\Delta }t)=f(xu(x,t)\mathrm{\Delta }t,t),$$
(5)
$$\frac{f(x,t+\mathrm{\Delta }t)}{x}=\left[1\frac{u(x,t)}{x}\mathrm{\Delta }t\right]\frac{f(xu(x,t)\mathrm{\Delta }t,t)}{X},$$
(6)
where $`\mathrm{\Delta }t`$ is the time interval during a computational step. In the case of $`u_i<0`$, those equations are represented using the Hermite cubic function as

$$f_i^{n+1}=f_i^n+d_i^n\xi +C2_i^n\xi ^2+C3_i^n\xi ^3,$$
(7)
$$d_i^{n+1}=d_i^n+2C2_i^n\xi +3C3_i^n\xi ^2,$$
(8)
where

$$C2_i^n=(d_{i+1}^n+2d_i^n3S_{i+1/2}^n)/h,$$
(9)
$$C3_i^n=(d_{i+1}^n+d_i^n2S_{i+1/2}^n)/h^2,$$
(10)
$$S_{i+1/2}^n=(f_{i+1}^nf_i^n)/h,$$

$$\xi =u_i\mathrm{\Delta }t,$$

$`h`$ is the grid width, and the superscript $`n`$ denotes the time $`t=n\mathrm{\Delta }t`$.

Basically, the CIP-CSL2 method is a solver for the equation of continuity,

$$\frac{f}{t}+\frac{uf}{x}=0.$$
(11)
In this method, a quadratic function of

$$Q_i(X)=f_i+2q1_i(Xx_i)+3q2_i(Xx_i)^2$$
(12)
is used to interpolate the physical quantity $`f`$, where the coefficients $`q1_i`$ and $`q2_i`$ are determined, for example in the case of $`u_i<0`$, from the constrains of

$$Q_i(x_{i+1})=f_{i+1}\mathrm{and}_{x_i}^{x_{i+1}}Q_i(x)𝑑x/h=\rho _{i+1/2},$$

as

$$q1_i=(f_{i+1}+2f_i3\rho _{i+1/2})/h,$$
(13)
$$q2_i=(f_{i+1}+f_i2\rho _{i+1/2})/h^2,$$
(14)
where $`\rho `$ is the cell-integrated average of $`f`$ within a computational cell, defined at the cell centers, and is used as a dependent variable in this method (We should note that the definition of $`\rho `$ is slightly different from that in Ref. ). $`\rho `$ is updated so that the total mass is conserved. The mass flux flowing out from the cell \[$`x_i`$, $`x_{i+1}`$\] through the point $`x_i`$ during $`\mathrm{\Delta }t`$ can be calculated by

$$\mathrm{\Delta }\rho _ih=_{x_i}^{x_iu_i\mathrm{\Delta }t}Q_i(x)𝑑x=f_i\xi +q1_i\xi ^2+q2_i\xi ^3,$$

where we assume again the first-order accuracy in time. Using this quantity, $`\rho `$ is updated as follows:

$$\rho _{i+1/2}^{n+1}=\rho _{i+1/2}^n+\mathrm{\Delta }\rho _{i+1}^n\mathrm{\Delta }\rho _i^n.$$
(15)
The use of this formula allows us the complete conservation of the total mass, i.e., the total sum of $`\rho `$.

The quantity $`f`$ is updated by a time splitting technique . We solve Eq. (11) by splitting into two phases as

$$\frac{f}{t}+u\frac{f}{x}=0,$$

$$\frac{f}{t}=\frac{u}{x}f.$$

The former one shows the advection of $`f`$ and, therefore, is solved as an interpolation problem by using Eq. (12) as

$$f_i^{n+1}=f_i^n+2q1_i^n\xi +3q2_i^n\xi ^2.$$
(16)
The latter one represents the change of $`f`$ due to the compression or expansion and is solved in general by a finite difference technique .

As was pointed out already , the coefficients of the first- and second-order terms in Eq. (16) correspond to those in Eq. (8) by replacing $`f`$ and $`\rho `$ in Eq. (16) with $`d`$ and $`S`$, respectively. Furthermore, Eq. (15) also can be rewritten into a form like the conventional CIP. By introducing a variable $`D`$ defined as

$$D_i=_{x_0}^{x_i}f(x)𝑑x\left(=\underset{j=0}{\overset{i1}{}}_{x_j}^{x_{j+1}}Q_j(x)𝑑x\right),$$

one obtains

$$D_i^{n+1}=D_i^n+\mathrm{\Delta }\rho _i^nh=D_i^n+f_i^n\xi +q1_i^n\xi ^2+q2_i^n\xi ^3,$$
(17)
and

$$\rho _{i+1/2}^n=(D_{i+1}^nD_i^n)/h,$$

where we assume $`\mathrm{\Delta }\rho _0^n=0`$, namely no transfer of mass occur at the boundary $`x_0`$. Those equations reproduce Eq. (15), and replacements of $`D`$ and $`f`$ in Eq. (17) with $`f`$ and $`d`$, respectively, yield Eq. (7) (Note that the former replacement reduces $`\rho `$ to $`S`$). This result suggests us that the variants of the CIP method such as the rational method  and the HCR method  can be modified to a conservative one only by those replacements.

## 3 A conservative method by the hybrid-cubic-rational interpolation

In Ref. , the author proposed a numerical solver for the advection equation by employing both the cubic  and rational  functions. In the method, the following combination of those functions is used as an interpolation function:

$$F(X)=\alpha R(X)+(1\alpha )C(X),$$
(18)
where $`R`$ and $`C`$ show the rational and the cubic functions, respectively, and $`\alpha `$ denotes the weighting parameter whose range is $`\alpha [0,1]`$. Equation (18) is reduced to the rational function for $`\alpha =1`$ and to the cubic one for $`\alpha =0`$. The value of $`\alpha `$ is determined theoretically so as to be the minimum of the values to which the convexity-preserving condition  is satisfied. Where the convexity-preserving condition is expressed as that $`F_{xx}(X)>0`$ for a concave data of $`d_i<S_{i+1/2}<d_{i+1}`$ or $`F_{xx}(X)<0`$ for a convex data of $`d_i>S_{i+1/2}>d_{i+1}`$ is satisfied for the region of $`X[x_i,x_{i+1}]`$. The resulting formulae are

$$F(k)=f_i+d_ihk+(G1_i+G2_i)k^2,$$
(19)
$$\frac{F(k)}{x}=d_i+[G1_i\frac{Q_i+D_i}{D_i}+2G2_i+(1\alpha )(Q_iD_i)]\frac{k}{h},$$
(20)
where

$$G1_i=\alpha P_i^2/D_i,G2_i=(1\alpha )(2P_iD_i),$$

$$D_i=Q_i+(P_iQ_i)k,$$

$$\alpha =\frac{M_i(M_i2)}{M_i(M_i2)+1},$$

$$M_i=\mathrm{max}[2,\mathrm{max}(\frac{Q_i}{P_i},\frac{P_i}{Q_i})]$$

with

$$P_i=(S_{i+1/2}d_i)h,Q_i=(d_{i+1}S_{i+1/2})h,$$

and

$$k=\xi /h=u_i\mathrm{\Delta }t/h.$$

For the case of $`u_i>0`$, we need replacements of $`i+1i1`$ and $`hh`$. This interpolation function provides oscillation-free results, unlike the cubic one, and higher resolution of solution than that by the conventional rational function . Same as the conventional CIP, this method is constructed using a physical quantity $`f_i`$ and its first spatial derivative $`d_i`$; thus this method may be modified into a conservative method only by the replacements as discussed in the last section. Though the explicit definition and use of the cell-integrated average $`\rho `$ may be more convenient in an actual application, we now use $`D`$, instead of $`\rho `$, and adapt the replacements of

$$\{\begin{array}{cccccccccccccccccccc}fD,& & & & & & & & & & & & & & & & & & & \\ df,& & & & & & & & & & & & & & & & & & & \end{array}$$
(21)
for demonstration. Namely, the subroutine “$`\mathrm{HCR}(f,d)`$” of the hybrid method whose parameters are $`f`$ and $`d`$ is used as a function of $`D`$ and $`f`$, i.e., $`\mathrm{HCR}(D,f)`$. The initial condition is set as follows. At first, $`f`$ is initialized using an arbitrary function $`G(x)`$, which expresses the initial profile of $`f`$, as

$$f_i^0=G(x_i).$$

Next, the initial value of $`D`$ is calculated using a recurrence procedure of

$$D_0^0=0,$$

$$D_i^0=D_{i1}^0+h(f_i^0+f_{i1}^0)/2,\mathrm{for}i=1,2,\mathrm{},N.$$
(22)
Here the last term of Eq. (22) means the approximated integration of $`f`$ over a cell between $`x_{i1}`$ and $`x_i`$.

In the next section we show some numerical results using those procedures.

## 4 Numerical experiments

### Example 1:

The first example is the linear propagation of square waves. In this example we use $`u=1`$, $`h=1`$ and

$$G(x_i)=\{\begin{array}{cccccccccccccccccccc}1,& \mathrm{for}\mathrm{\hspace{0.33em}\hspace{0.33em}13}i21,& & & & & & & & & & & & & & & & & & \\ 1,& \mathrm{for}\mathrm{\hspace{0.33em}\hspace{0.33em}40}i48,& & & & & & & & & & & & & & & & & & \\ 0,& \mathrm{elsewhere}.& & & & & & & & & & & & & & & & & & \end{array}$$

Thus, the width of the square waves is $`9h`$. Figures 1 and 2 show the results using $`\mathrm{CFL}=0.2`$ at $`n=200`$ and $`n=2000`$, respectively. Here we show results using four types of method, the HCR, CIP (the case of $`\alpha =0`$), rational ($`\alpha =1`$), and modified rational methods. In the modified rational method , a switching technique, which is represented as

$$\gamma =\{\begin{array}{cccccccccccccccccccc}1,& d_id_{i+1}<0,& & & & & & & & & & & & & & & & & & \\ 0,& \mathrm{otherwise},& & & & & & & & & & & & & & & & & & \end{array}$$

is added to the conventional rational method to select the interpolation function. The rational function is adopted in the case of $`\gamma =1`$ and the cubic one in the case of $`\gamma =0`$. It was pointed out that this technique sometimes breaks the preservation of the convexity of solution while this improves the dissipation property of the rational method . In the present example, this technique is adapted by using $`f`$, not by $`d`$, because of the replacements. In the figures, we see clearly that the HCR method provides most accurate results, which are less diffusive than those by the conventional rational method and less oscillatory than those by the CIP and modified rational methods. In Table 1, we show the calculated values at $`n=200`$ around the point $`i=12`$ at which the left discontinuity of the negative pulse locates. The calculated data by the CIP and modified rational methods are positive near the left side of the discontinuity and become minimums at $`i=14`$, which are smaller than $`1`$. On the contrary, the data by the HCR and conventional rational methods decrease monotonously as $`i`$ increases, and are not less than $`1`$. Those results mean that only the HCR and conventional rational methods can keep the convexity of solution.

For comparison, we shall show results without the replacement of Eq. (21). Figure 3 shows the results at $`n=2000`$ given by using the same initial condition as the previous one. In this case, we set the initial value of the first derivative as $`d_i^0=0`$ at all points. As was pointed out in Ref. , the CIP-CSL2 and CIP methods (the CIP with and without the replacements, respectively) give quite similar results each other. However, another methods provide obviously different results by whether the replacements are employed or not. It is interesting to point out that the result by the HCR method with the replacements is less diffusive than that by without the replacements, while we cannot explain this result for the moment.

### Example 2:

Next, we solve the problem used in Ref.  to compare the performance of the conventional rational method with those of the MUSCL  and PPM  schemes (See Fig. 1 in the paper). In this example, we use $`h=1`$ and $`\mathrm{CFL}=0.2`$. The initial value of $`f`$ is set to

$$G(x_i=i)=\{\begin{array}{cccccccccccccccccccc}(i20)/11,& \mathrm{for}\mathrm{\hspace{0.33em}\hspace{0.33em}20}i<31,& & & & & & & & & & & & & & & & & & \\ 1(i31)/20,& \mathrm{for}\mathrm{\hspace{0.33em}\hspace{0.33em}31}i<41,& & & & & & & & & & & & & & & & & & \\ 1/2,& \mathrm{for}\mathrm{\hspace{0.33em}\hspace{0.33em}41}i<60,& & & & & & & & & & & & & & & & & & \\ 1,& \mathrm{for}\mathrm{\hspace{0.33em}\hspace{0.33em}60}i<80,& & & & & & & & & & & & & & & & & & \\ 0,& \mathrm{elsewhere},& & & & & & & & & & & & & & & & & & \end{array}$$

and the velocity, assumed as constant in time and space, is set to $`u=1`$. Figures 4 and 5 show the results at $`n=440`$, given by the four types of methods, with and without the replacements, respectively. In Ref.  it was concluded that the conventional rational method gives better representation of the sharp corner at the top of the triangular wave than those by the MUSCL and PPM methods. From the present result, however, we know that the HCR method and its conservative version achieve still better representation of that than the conventional rational method. The calculated maximum values of $`f`$ at the corner are 0.935 (conservative HCR), 0.937 (HCR), 0.916 (conservative rational), and 0.923 (rational). Furthermore, we know that the conservative HCR method gives almost equivalent (or slightly better) resolution of the square wave to that by the PPM, which gives the best one among those by the schemes discussed in Ref. .

### Example 3:

Now we show results to a nonlinear problem in which the inviscid Burgers’ equation,

$$\frac{u}{t}+u\frac{u}{x}=0,$$
(23)
is selected as the governing equation . This equation can be solved directly by the conservative methods discussed in this paper. For the calculation of $`u`$, we can use the procedures for $`f`$ because Eq. (23) is an advection equation (here $`f`$ means that after the replacements). For the calculation of $`D`$ ($`=u𝑑x`$), we need to rewrite Eq. (23) into a conservation form as

$$\frac{u}{t}+\frac{(u^2/2)}{x}=0.$$
(24)
From this formula we found that the tranportation velocity $`u`$ used in the procedures for $`D`$ should be replaced with $`u/2`$ .

The initial condition and the grid width are set to

$$G(x)=0.5+0.4\mathrm{cos}(2\pi x/100)$$

and $`h=1`$. Figure 6 shows results at $`t=100`$ with $`\mathrm{\Delta }t=0.1`$. In this figure, the solid line shows the result by the first-order upwind method using the conservation formula (24) with $`h=1/5`$ and $`\mathrm{\Delta }t=0.1/5`$. Same as the result by the CIP-CSL2, that by the conservative HCR method represents the correct position of the shock wave. Figure 6(c) shows the result by the HCR method without the replacements. This result is obviously incorrect. Those results prove that the HCR method has gained complete conservation by the replacements. Figure 7 shows the maximum and minimum values of $`u`$ as a function of time. In the result by the CIP-CSL2 we can observe unphysical strong oscillation due to the numerical dispersion of a method using a classical polynomial interpolation; however, such the oscillation cannot be seen in the result by the present method. This result allows us to consider that the present method is applicable to a shock problem even without any artificial viscosity, and is oscillation free even for the nonlinear problem.

### Example 4:

Now we shall try an experimental study which uses further replacements. We already know that the CIP method provides an almost equivalent result whether the replacements (21) are adapted or not. If the variables $`D`$ and $`f`$ are replaced with those integrated once more, what result is obtained? Figure 8 shows results to Example 2 obtained by replacing $`D_i^0`$ and $`f_i^0`$ with

$$E_i^0_{x_0}^{x_i}D^0𝑑x\mathrm{and}D_i^0$$

and by updating those variables with the HCR and conventional methods. The parameters are same as those used in Example 2. $`f_i`$ shown in Fig. 8 is given by $`f_i=(E_{i+1}2E_i+E_{i1})/h^2`$. Despite the adaptation of two times of the replacements, the CIP method provides an almost equivalent result to those given without or with one-time replacements. However, one can observe some differences in the results using another methods from the previous ones. The numerical oscillation around the discontinuities of the square wave appears even in the results obtained by using the HCR and the conventional rational methods. In the result by the modified rational method, the oscillation can be observed not only at the top of the square wave but also at the bottom of it.

## 5 An additional comment

It is known that when no artificial viscosity is adapted, a numerical method for an advection equation based on the non-conservative formula cannot give correct propagation speed of shock in the region where the direction of the advection velocity changes spatially. In Ref.  the use of high-order polynomial that covers two cells are discussed for solving this problem. Now we present an alternative method to overcome this difficulty.

Here we use again the Burgers’ equation (23) as the governing equation. Introducing $`Uu+C`$, where $`C`$ is a constant, reduces Eq. (23) to

$$\frac{U}{t}+U\frac{U}{x}C\frac{U}{x}=0.$$

This is also the Burgers’ equation, but has a linear advection term (the last term of left hand side). We can solve this equation by splitting into two stages. First, we solve $`U/t+UU/x=0`$ and update $`U`$ during the time interval $`\mathrm{\Delta }t`$. Next, we shift the position of $`U`$ by $`C\mathrm{\Delta }t`$. If $`C`$ is determined so that $`U_i>0`$ or $`U_i<0`$ for all $`i`$, we may be able to calculate the position of shock correctly. In the case where we set $`C`$ so that the condition of $`m\left|C\mathrm{\Delta }t\right|=h`$ (where $`m`$ is a positive integer) is satisfied, the shift step can be solved only by one substitution like $`U_{i+1}=U_i`$ per $`m`$ computational steps. This procedure is convenient since we need no additional interpolation function or difference equation; however, further discussions should be required for the practical application.

## 6 Conclusion

In this paper we have derived a conservative variant of the hybrid-cubic-rational method for an advection equation. It was proved by some numerical experiments that the hybrid method is oscillation free even when it is modified to a conservative one. Interestingly, the conservative hybrid method gives less-diffusive representation of the discontinuities of the square wave than that by the non-conservative one. We should try to make this result clear theoretically.

We expect that combining the present method with the hybrid interpolation-extrapolation method, which realizes the discontinuous representation of the density interface , yields a powerful tool for a semi-Lagrangian computation of multiphase flows like those treated in Refs. .


text:
# The Deconfinement Transition in 𝑆𝑈(4) Lattice Gauge Theory hep-lat/0110054, TIFR/TH/01-37

## 1 INTRODUCTION

The nature of the phase transition(s) to Quark-Gluon Plasma, which we hope to see in the experiments of RHIC and later of LHC, and the physics driving it, have always been of great interest. While the real world has presumably two very light ($`u`$,$`d`$) flavours of quarks and one somewhat heavier ($`s`$) flavour, both analytical and numerical methods in lattice QCD begin from the limiting cases of either massless or infinitely massive quarks. One talks of the chiral symmetry restoring phase transition and the deconfinement phase transition in these two cases respectively and has suitable order parameters to investigate them. For quarks with $`N`$ colours and $`N_f`$ massless flavours, these transitions are related to spontaneous breaking of a global $`Z(N)`$ and $`SU(N_f)\times SU(N_f)`$ chiral symmetry. Since these symmetries are broken explicitly to various extents in the real world, which of them is more relevant is a priori not clear. The low masses of the light flavours suggest chiral symmetry to be the dominant one. However, it is seen in numerical simulations that the energy density shows a large change at the chiral transition and even the order parameter for the deconfinement transition,

$$L(\stackrel{}{x})=\frac{1}{N}\mathrm{tr}\underset{t=1}{\overset{N_t}{}}U_4(\stackrel{}{x},t),$$

also rises to nonzero values there. These apparently mysterious observations can be explained using large $`N`$ arguments, if the deconfinement transition for $`N4`$ is of second order. $`SU(4)`$ is clearly the simplest case to check this out.

Numerical simulations of $`SU(4)`$ theory at finite temperatures have been done in the past and recently as well. All of them used the Wilson action, or the more general mixed action:

$$S=\underset{P}{}\left[\beta \left(1\frac{\mathrm{Re}\mathrm{tr}U_P}{N}\right)+\beta _A\left(1\frac{\mathrm{tr}_\mathrm{A}U_P}{N}\right)\right],$$
(1)
A well known problem in the simulations with these actions, especially for large $`N`$, is the presence of a bulk transition which is a lattice artifact. The phase diagram of the mixed action,

in its coupling plane is as shown in Fig. 1. The solid lines in it show first order bulk transitions lines. The dotted line after the end point D is drawn to suggest the impact D may have on the Wilson axis ($`\beta _A`$ = 0). For $`N4`$, D is expected to be where E is shown, causing a first order bulk transition for the usual Wilson action. In order to avoid it, simulations at negative $`\beta _A`$  and for larger  $`N_t`$ = 6 were made. They lead to a first order deconfinement phase transition for $`SU(4)`$.

From our extensive studies of the deconfinement phase transition for the action above but for the $`SU(2)`$ theory, we know that bulk transitions affect the order and location of the deconfinement transition in subtle and inexplicable ways, even leading to apparent qualitative violations of universality. Universality was restored  in that case only after eliminating the bulk transitions associated with the $`Z(2)`$ vortices and $`Z(2)`$ monopoles by adding large chemical potentials for them. It seems natural to expect that the bulk transitions for $`N>2`$ can also be cleaned off by suppressing the corresponding $`Z(N)`$ objects. We pursue this idea here for $`SU(4)`$ to investigate its deconfinement phase transition.

## 2 SIMULATIONS AND RESULTS

Generalizing the idea of positive plaquette models in the literature, we use the action

$$S=\beta \underset{P}{}\left(1\frac{\mathrm{tr}U_P}{N}\right)\theta (\frac{\pi }{N}|\alpha |),$$
(2)
where $`\pi <\alpha \pi `$ is the phase of tr $`U_P`$. By adding the adjoint term of eq. (2) to the action (3), one sees that the phase diagram of the resultant mixed should not have the bulk lines AC or BC and hence the endpoint D or E.

We have simulated the above action on $`N_s^3\times N_t`$ lattices for $`N_s`$ =8, 10, 12, 15, 16 and $`N_t`$ = 4, 5, 6, 8 using a 15-hit Metropolis et al. algorithm. The calculations were done on a cluster of pentiums. Typically short runs to look for points of rapid variations in $`|L|`$ were followed by long runs (a few million sweeps) to determine the susceptibility $`\chi _{|L|}`$ using the histogramming technique. Usual finite size scaling techniques were used to determine the order of the transition and its exponents.

In simulations on $`N_s^3\times 4`$ lattices, $`N_s=`$ 8, 10, 12, one sees hot and cold starts converge quickly at couplings a little away from the transition point on its both sides but a clear co-existence of states is visible for all lattices at the transition point. The tunneling frequency goes down with spatial volume. The histograms of $`|L|`$ show peaks which become narrower with increasing volume and the gap between them remains unchanged. These classic signs of a first order phase transition are confirmed by a quantitative analysis of the linear growth of $`\chi _{|L|}`$ with volume, as seen in Fig. 2. The horizontal lines in each case are predictions obtained by scaling the $`N_s=8`$ results linearly with volume.

For larger $`N_t`$, we used many longer runs in the region of strong variation of $`|L|`$ to obtain the susceptibility directly and used the histogramming only for the finer determination of the critical coupling. Our results for $`|L|`$ as a function of $`\beta `$ clearly show the expected shift for a deconfinement phase transition for $`N_t`$ = 5, 6 and 8. This is evident in the $`\beta _c`$ determinations from the $`\chi _{|L|}^{\mathrm{max}}`$, as seen in Table 1 for $`N_t`$ = 6 for two different spatial volumes. Again, using the peak height for the smaller spatial volume, the $`\chi ^{\mathrm{max}}`$ on the bigger lattice can be predicted, assuming a linear growth with volume. The prediction listed in the table can be seen to be in very good agreement with the direct Monte Carlo determination. Along with the shifts in $`\beta `$, it confirms that the same physical deconfinement phase transition is being simulated on these lattices thus approaching the continuum limit of $`a0`$ in a progressive manner by keeping the transition temperature $`T_c`$ constant in physical units.

## 3 LATENT HEAT

While the results above for the deconfinement order parameter $`|L|`$ and the corresponding susceptibility, $`\chi _{|L|}`$, are indicative of a first order deconfinement phase transition, one needs to make sure that they indeed are not due to a coincident first order bulk transition. Apart from the characteristic (logarithmic) shift of the transition point with $`N_t`$, seen above, the latent heat of a first order deconfinement phase transition should also remain constant as $`N_t\mathrm{}`$. Requiring the pressure to be continuous at the deconfinement phase transition, the latent heat can be obtained from two different observables $`\mathrm{\Delta }_1\mathrm{\Delta }(ϵ3p)/T_c^4`$, and $`\mathrm{\Delta }_2\mathrm{\Delta }(ϵ+p)/T_c^4`$, where

$`\mathrm{\Delta }_1=48N_t^4a{\displaystyle \frac{g^2}{a}}\mathrm{\Delta }P,`$
$`\mathrm{\Delta }_2=32N_t^4{\displaystyle \frac{C(g^2)}{g^2}}(\mathrm{\Delta }P_t\mathrm{\Delta }P_s),`$ (3)
and $`C(g^2)=(10.2366g^2+O(g^4))`$ for $`SU(4)`$. $`\mathrm{\Delta }`$ denotes discontinuities across the transition in the respective variables. In order to obtain the $`\mathrm{\Delta }P`$, $`\mathrm{\Delta }P_s`$ and $`\mathrm{\Delta }P_t`$, the minimum of the histogram $`N(|L|)`$ was used to separate the two phases in each case. From eq. (3), it is clear that the plaquette discontinuity $`\mathrm{\Delta }PN_t^4`$ in order to obtain the same latent heat in physical units, as $`N_t\mathrm{}`$. Indeed, its decrease with $`N_t`$ was seen to be consistent with expectations for $`N_t5`$. Furthermore, both estimates must agree in this limit, as the neglected cut-off corrections become then insignificant. Table 2 suggests this to be the case, leading to an estimate which is $`\frac{2}{3}`$ of the ideal gas energy density at $`T_c`$ and agrees with earlier results.

## 4 SUMMARY

The linear growth of $`\chi _{|L|}^{\mathrm{max}}`$ with volume for $`N_t`$ = 4, suggests a first order deconfinement phase transition for $`SU(4)`$. Various indicators, such as, histograms and evolutions, are in accord with this. Increasing $`N_t`$ to 5, 6 and 8 shows the expected shift of the deconfinement transition, with a growth in $`\chi _{|L|}^{\mathrm{max}}`$ that is consistent with being linear in volume. The plaquette discontinuity $`\mathrm{\Delta }P`$ decreases as the fourth power of $`N_t`$, indicating both a lack of a bulk transition and a first order deconfinement phase transition. The large estimated latent heat, being about 2/3 of the corresponding ideal gas energy density, suggests the deconfinement transition to grow stronger in nature as the number of colours in increased.


text:
# On Āryabhaṭa’s Planetary Constants

## 1 Introduction

An old problem in the history of Indian science is whether ideas at the basis of Āryabhaṭa’s astronomy were borrowed from outside or were part of India’s own tradition. This problem was first raised in the context of the now discredited thesis that sound observational astronomy did not exist in India prior to India’s encounter with the West. Thus in a recent paper,1 Abhyankar argues that “Āryabhaṭa’s values of bhagaṇas were probably derived from the Babylonian planetary data.” But Abhyankar makes contradictory assertions in the paper, suggesting at one place that Āryabhaṭa had his own observations and at another place that he copied numbers without understanding, making a huge mistake in the process.

In support of his theory, Abhyankar claims that Āryabhaṭa used the Babylonian value of 44528 synodic months in 3600 years as his starting point. But this value is already a part of the Śatapatha altar astronomy reconciling lunar and solar years in a 95-year yuga. In this ritual, an altar is built to an area that is taken to represent the nakṣatra or the lunar year in tithis and the next design is the same shape but to a larger area (solar year in tithis), but since this second design is too large, the altar construction continues in a sequence of 95 years. It appears that satisfactory reconciliation by adding intercalary months to the lunar year of 360 tithis amounted to subtracting a certain number of tithis from the 372 tithis of the solar year, whose most likely value was 89 tithis in 95 years.2

The areas of the altars increase from $`7\frac{1}{2}`$ to $`101\frac{1}{2}`$ in the 95 long sequence in increments of one. The average size of the altar is therefore $`54\frac{1}{2}`$, implying that the average difference between the lunar and the solar year is taken to be one unit with $`54\frac{1}{2}`$ which is about $`6.60`$ tithis for the lunar year of 360 tithis. This is approximately correct.

Considering a correction of 89 tithis in 95 years, the corrected length of the year is $`37289/95=371.06316`$ tithis. Since each lunation occurs in 30 tithis, the number of lunations in 3600 years is 44527.579. In a Mahāyuga, this amounts to 53,433,095. In fact, the number chosen by Āryabhaṭa (row 1 in Table 1) is closer to this number rather than the Babylonian number of 53,433,600.

Table 1 presents the Babylonian numbers given by Abhyankar together with the Āryabhaṭa constants related to the synodic lunar months and the revolutions of the lunar node, the lunar apogee, and that of the planets. The so-called Babylonian numbers are not actually from any Babylonian text but were computed by Abhyankar using the rule of three on various Babylonian constants.

Table 1: Revolutions in one Mahāyuga

| Type | Babylonian | Ārybhaṭa |
| --- | --- | --- |
| Synodic lunar months | 53,433,600 | 52,433,336 |
| Lunar node | -232,616 | -232,352 |
| Lunar apogee | 486,216 | 488,219 |
| Mercury | 17,937,000 | 17,937,020 |
| Venus | 7,022,344 | 7,022,388 |
| Mars | 2,296,900 | 2,296,824 |
| Jupiter | 364,216 | 364,224 |
| Saturn | 146,716 | 146,564 |

We see that no numbers match. How does one then make the case that Āryabhaṭa obtained his numbers from a Babylonian text? Abhyankar says that these numbers are different because of his (Āryabhaṭa’s) own observations “which are more accurate.” But if Āryabhaṭa had his own observations, why did he have to “copy” Babylonian constants, and end up not using them, anyway?

Certain numbers have great discrepancy, such as those of the lunar apogee, which Abhyankar suggests was due to a “wrong reading of 6 by 8” implying–in opposition to his earlier view in the same paper that Āryabhaṭa also had his own observations– that Āryabhaṭa did not possess his own data and that he simply copied numbers from some manual brought from Babylon!

The Āryabhaṭa numbers are also more accurate that Western numbers as in the work of Ptolemy.3 Given all this, there is no credible case to accept the theory of borrowing of these numbers from Babylon.

Abhyankar further suggests that Āryabhaṭa may have borrowed from Babylon the two central features of his system: (i) the concept of the Mahāyuga, and (ii) mean superconjunction of all planets at some remote epoch in time. In fact, Abhyankar repeats here an old theory of Pingree4 and van der Waerden5 about a transmission from Babylon of these two central ideas. In this paper, we show that these ideas were already present in the pre-Siddhāntic astronomy and, therefore, a contrived connection with Babylonian tables is unnecessary.

## 2 The Indic tradition of yugas and superconjunctions

In the altar ritual of the Brāhmaṇas,6 equivalences by number connected the altar area to the length of the year. The 5-year yuga is described in the Vedāṅga Jyotiṣa, where only the motions of the sun and the moon are considered. The Śatapatha Brāhmaṇa describes the 95-year cycle to harmonize the solar and the lunar years. The Śatapatha Brāhmaṇa also describes an asymmetric circuit for the sun7, which the Greeks speak about only around 400 BC.

Specifically, we find mention of the nominal year of 372 tithis, the nakṣatra year of 324 tithis, and a solar year of 371 tithis. The fact that a further correction was required in 95 years indicates that these figures were in themselves considered to be approximate.

In the altar ritual, the primal person is made to an area of $`7\frac{1}{2}`$ puruṣas, when a puruṣa is also equated with 360 years leading to another cycle of 2700 years. This is the Saptarṣi cycle which was taken to start and end with a superconjunction.

The Śatapatha Brāhmaṇa 10.4.2.23-24 describes that the Ṛgveda has 432,000 syllables, the Yajurveda has 288,000 and the Sāmaveda has 144,000 syllables. This indicates that larger yugas in proportion of 3:2:1 were known at the time of the conceptualization of the Saṃhitās.

Since the nominal size of the Ṛgveda was considered to be 432,000 syllables (SB 10.4.2.23) we are led to the theory of a much larger yuga of that extent in years since the Ṛgveda represented the universe symbolically.

Elsewhere, I show8 how the Vedāṅga Jyotiṣa serves as a coordinate system for the sun and the moon in terms of the 27 nakṣatras. Such a coordinate system implies a calculation where whole cycles are subtracted from large numbers. Such modular arithmetic appears to lie at the basis of the idea of a superconjunction. Traditionally, the Vedāṅga Jyotiṣa has been dated to around 1350 BC, but a new paper by Narahari Achar9 argues for a much earlier date of 1800 BC.

Van der Waerden10 has argued that a primitive epicycle theory was known to the Greeks by the time of Plato. He argued such a theory might have been known in the wider Indo-European world by early first millennium BC. With new ideas about the pre-history of the Indo-European world emerging, it is possible to push this to an earlier millennium. An old theory may be the source which led to the development of very different epicycle models in Greece and India.

The existence of an independent tradition of observation of planets and a theory thereof as suggested by our analysis of the Śatapatha Brāhmaṇa helps explain the puzzle why the classical Indian astronomy of the Siddhānta period uses many constants that are different from those of the Greeks.

## 3 More on the Great Year

Since the yuga in the Vedic and the Brāhmaṇa periods is so clearly obtained from an attempt to harmonize the solar and the lunar years, it appears that the consideration of the periods of the planets was the basis of the creation of an even longer yuga.

There is no reason to assume that the periods of the five planets were unknown during the Brāhmaṇa age. I have argued that the astronomical numbers in the organization of the Ṛgveda indicate with high probability the knowledge of these periods in the Ṛgvedic era itself.11

Given these periods, and the various yugas related to the reconciliation of the lunar and the solar years, we can see how the least common multiple of these periods will define a still larger yuga.

The Mahābhārata and the Purāṇas speak of the kalpa, the day of Brahmā, which is 4,320 million years long. The night is of equal length, and 360 such days and nights constitute a “year” of Brahmā, and his life is 100 such years long. The largest cycle is 311,040,000 million years long at the end of which the world is absorbed within Brahman, until another cycle of creation. A return to the initial conditions (implying a superconjunction) is inherent in such a conception. Since the Indians and the Persians were in continuing cultural contact, it is certain that this old tradition became a part of the heritage of the Persians. This explains how we come across the idea of the World-Year of 360,000 years in the work of Abū Ma’shar, who also mentioned a planetary conjunction in February 3102 BC.

The theory of the transmission of the Great Year of 432,000 years, devised by Berossos, a priest in a Babylonian temple, to India in about 300 BC, was advanced by Pingree.12 But we see this number being used in relation to the Great Year in the Śatapatha Brāhmaṇa itself, a long time before Berossos.13

The idea of superconjunction seems to be at the basis of the cyclic calendar systems in India. The Śatapatha Brāhmaṇa speaks of a marriage between the Seven Sages, the stars of the Ursa Major, and the Kṛttikās; this is elaborated in the Purāṇas where it is stated that the ṛṣis remain for a hundred years in each nakṣatra. In other words, during the earliest times in India there existed a centennial calendar with a cycle of 2,700 years. Called the Saptarṣi calendar, it is still in use in several parts of India. Its current beginning is taken to be 3076 BE.

The usage of this calendar more than 2000 years ago is confirmed by the notices of the Greek historians Pliny and Arrian who suggest that, during the Mauryan times, the Indian calendar began in 6676 BC. It seems quite certain that this was the Saptarṣi calendar with a beginning which starts 3600 years earlier than the current Saptarṣi calendar.

The existence of a real cyclic calendar shows that the idea of superconjunction was a part of the Indic tradition much before the time of Berossos. This idea was used elsewhere as well but, given the paucity of sources, it is not possible to trace a definite place of origin for it.

## 4 Conclusions

More than thirty years ago, Roger Billard showed14 the falsity of the 19th century notion that India did not have observational astronomy. His analysis of the Siddhāntic and the practical karaṇa texts demonstrated that these texts provide a set of elements from which the planetary positions for future times can be computed. The first step in these computations is the determination of the mean longitudes which are assumed to be linear functions of time. Three more functions, the vernal equinox, the lunar node and the lunar apogee are also defined.

Billard investigated these linear functions for the five planets, two for the sun (including the vernal equinox) and three for the moon. He checked these calculations against the values derived from modern theory and he found that the texts provide very accurate values for the epochs when they were written. Since the Siddhānta and the karaṇa models are not accurate, beyond these epochs deviations build up. In other words, Billard refuted the theory that there was no tradition of observational astronomy in India. But Billard’s book is not easily available in India, which is why the earlier theory has continued to do rounds in Indian literature.

Āryabhaṭa’s constants are more accurate than the one’s available in the West at that time. He took old Indic notions of the Great Yuga and of cyclic time (implying superconjunction) and created a very original and novel siddhānta. He presented the rotation information with respect to the sun which means that his system was heliocentric to a certain extent.15 Furthermore, he considered the earth to be rotating on its own axis. Since we don’t see such an advanced system amongst the Babylonians prior to the time of Āryabhaṭa, it is not reasonable to look outside of the Indic tradition or Āryabhaṭa himself for the data on which these ideas were based.

## Notes and References

1. K.D. Abhyankar, “Babylonian source of Āryabhaṭa’s planetary constants.” Indian Journal of History of Science 35: 185-188, 2000.
2. S. Kak, The Astronomical Code of the Ṛgveda. New Delhi: Munshiram Manoharlal, 2000, page 87.
3. K.S. Shukla and K.V. Sarma, Āryabhaṭīya of Āryabhaṭa. New Delhi: Indian National Science Academy, 1976.
4. D. Pingree, “Astronomy and astrology in India and Iran.” Isis 54: 229-246, 1963.
5. B.L. van der Waerden, “The great year in Greek, Persian and Hindu astronomy.” Archive for History of Exact Sciences 18: 359-384, 1978.
6. S. Kak, The Astronomical Code of the Ṛgveda. New Delhi: Munshiram Manoharlal, 2000;
S. Kak, “Birth and early development of Indian astronomy.” In Astronomy Across Cultures: The History of Non-Western Astronomy, H. Selin (ed.), pp. 303-340. Dordrecht: Kluwer Academic Publishers, 2000.
7. S. Kak, “The sun’s orbit in the Brāhmaṇas.” Indian Journal of History of Science 33: 175-191, 1998.
8. S. Kak, “The astronomy of the age of geometric altars.” Quarterly Journal of the Royal Astronomical Society 36: 385-396, 1995.
9. B.N. Narahari Achar, “A case for revising the date of Vedāṅga Jyotiṣa.” Indian Journal of History of Science 35: 173-183, 2000.
10. B.L. van der Waerden, “The earliest form of the epicycle theory.” Journal for the History of Astronomy 5: 175-185, 1974.
11. S. Kak, The Astronomical Code of the Ṛgveda. New Delhi: Munshiram Manoharlal, 2000.
12. See note 4, above.
13. B.N. Narahari Achar, “On the astronomical basis of the data of Śatapatha Brāhmaṇa: a re-examination of Dikshit’s theory.” Indian Journal of History of Science 35: 1-19, 2000.
14. Roger Billard, L’astronomie Indienne. Paris: Publications de l’ecole francaise d’extreme-orient, 1971;
B.L. van der Waerden, “Two treatises on Indian astronomy.” Journal for History of Astronomy 11: 50-58, 1980.
15. H. Thurston, Early Astronomy. New York: Springer-Verlag, 1994, page 188.


text:
# Stellar Populations in Star Clusters: The Rôle Played by Stochastic Effects

## 1. Introduction

When photometric properties of star clusters are compared to the predictions of population synthesis models, it is common to find that the scatter among the observed colors and the difference between these colors and the model predictions, are larger than allowed by the observational errors. This is particularly true for intermediate age clusters in optical-IR colors, like $`VK`$, but to a lesser extent the scatter is also large in $`BV`$. Fig 1 illustrates this point.

In a recent paper, Maraston et al. (2001) show that the population synthesis models of Maraston (1998) reproduce quite well the $`VK`$ colors of the cluster population observed in NGC 7252 at an age between 400 and 600 Myr. In these models, built on the fuel consumption theorem (Renzini & Buzzoni 1986), the contribution of AGB stars to the bolometric light in the synthetic population has been increased over that expected from stellar evolution theory, in such a way as to reproduce the empirical values determined by Frogel et al. (1990) for a sample of LMC clusters as a function of SWB class. Even though this is a valid procedure, in this paper I show that the optical and optical-IR colors of the intermediate age clusters of the LMC and NGC 7252 are consistent with the colors expected from models (Bruzual & Charlot 2001; Bruzual 2000, 2001) based on current stellar evolution theory, if stochastic fluctuations in the number of stars in sparsely populated evolutionary stages are properly included into the models. In what follows I use the Montecarlo technique pioneered by Barbaro & Bertelli (1977), Chiosi et al. (1988), Girardi et al. (1995), Santos & Frogel (1997), and Cerviño, Luridiana & Castander (2000), to study the role that stochastic effects on the initial mass function (IMF) have on star cluster colors. I follow the Santos & Frogel (1997) approach to study the LMC and NGC 7252 cluster samples.

## 2. Method

The derivation below follows closely Santos & Frogel (1997). The IMF

$$\mathrm{\Phi }(m)=dN/dm=Cm^{(1+x)},$$
(1)
normalized as usual,

$$C=\frac{x}{m_l^xm_u^x},$$
(2)
obeys, $`\mathrm{\Phi }(m)0`$, and

$$_{m_l}^{m_u}\mathrm{\Phi }(m^{})𝑑m^{}=1.$$
(3)
$`\mathrm{\Phi }(m)`$ can be interpreted as a probability distribution function which gives the probability that a random mass $`m^{}`$ is in the range between $`m`$ and $`m+dm`$. $`\mathrm{\Phi }(m)`$ can be transformed into another probability distribution function $`g(N)`$, such that the probability of occurrence of the random variable $`N^{}`$ within $`dN`$ and the probability of occurrence of the random variable $`m^{}`$ within $`dm`$ is the same,

$$|\mathrm{\Phi }(m)dm|=|g(N)dN|.$$
(4)
where $`N`$ is a single-valued function of $`m`$. From (1),

$$N(m)=_{m_l}^m\mathrm{\Phi }(m^{})𝑑m^{},$$
(5)
is a cumulative distribution function which gives the probability that the mass $`m^{}`$ is less or equal to $`m`$, and using (4), it follows that

$$g(N)=1,0N1.$$
(6)
$`g(N)`$ is thus a uniform distribution for which any value is equally likely in the interval $`0N1`$. If we sample $`N`$ using a random number generator (Press et al. 1992), from (5) we can obtain $`m`$ as a function of $`N`$,

$$m=[(1N)m_l^x+Nm_u^x]^{\frac{1}{x}}.$$
(7)
For each star of mass $`m`$ generated in this way, we obtain its observational properties from the log $`T_{eff}`$ and log $`L`$ corresponding to this mass in the isochrone at the age of interest. We repeat the procedure until the cluster mass, i.e. the sum of $`m`$ for all the stars generated, including dead stars, reaches the desired value. Adding the contribution of each star to the flux in different photometric bands, we obtain the cluster magnitude and colors in these bands. Throughout this paper I use $`x=1.35`$ (Salpeter 1955), and for the lower and upper mass limits of star formation, I assume $`m_l=0.09`$ and $`m_u=125`$ M, respectively.

## 3. Results

Figs 2 to 6 show the results of the Montecarlo simulations. For reasons of space the details of each figure given in the caption are not repeated in the text. From Figs 2 and 3 it is apparent that the fluctuations in the colors become larger as the cluster mass decreases. The number of stars in the $`1\times 10^3`$ M cluster may be unrealistically small, and some evolutionary stages are not sampled. There are not enough high mass MS stars in the $`1\times 10^3`$ M and $`1\times 10^4`$ M clusters to make them as blue in $`UB`$ at early ages as the model computed with the analytical IMF (equivalent to infinite number of stars). In the $`1\times 10^5`$ M case the upper MS is well sampled and both models are equally blue in $`UB`$. Fig 3 shows that the lower mass cluster models at early ages are bluer in $`VK`$ than the analytic model because the lower MS is over populated with respect to the short lived phases with redden this color. For the $`1\times 10^5`$ M cluster the agreement is good at early ages. The expected fluctuations in $`VK`$ for a $`1\times 10^4`$ M cluster amount to almost 2 magnitudes, in close agreement with the range of colors observed at a given age, and is considerably broader than for the $`1\times 10^4`$ M case. Fig 4 shows an enlargement of the diagrams for the $`1\times 10^4`$ M star cluster. At the intermediate ages the models are redder in $`VK`$ than the analytic IMF model because of a larger number of AGB stars, which appear naturally as a consequence of stochastic fluctuations in the IMF. The regions in these planes occupied by the observations and the cluster simulations overlap quite well. This comparison seems to favor a not very large mass for these clusters, if stochastic fluctuations are responsible of the color variations. This conclusion is supported by Fig 5 in which I plot $`K`$ vs $`VK`$ at the distance of the LMC. The region covered by the $`1\times 10^4`$ M clusters matches the region covered by the observations very well. This is not the case for the $`1\times 10^5`$ M clusters, which show a much lower dispersion in color for the same $`K`$ magnitude than the observations and the lower mass simulations.

## 4. Conclusions

Even though the Montecarlo technique used in this paper has been used successfully before to explain cluster colors, I have revisited it in this paper to show that there is no need to introduce ad-hoc assumptions into population synthesis models (Maraston 1998, Maraston et a. 2001), which represent a departure from our current understanding of stellar evolution theory, in order to explain the observed range of values of cluster colors and magnitudes. In this paper I have briefly summarized the results of simple simulations which show clearly how the range of colors observed in intermediate age star clusters can be understood on the basis of current stellar evolution theory, if we take into account properly the expected variation in the number of stars occupying sparsely populated evolutionary stages due to stochastic fluctuations in the IMF.

Our preliminary conclusion that cluster masses in the range around $`1\times 10^4`$ M are preferred, as well as the dependence on this mass of the fluctuations in magnitude and colors will be explored in detail by Bruzual & Charlot (2001).

The interested reader may request results for other quantities equally sensitive to stochastic fluctuations, like the number of ionizing photons from a young stellar population, not shown here for lack of space (see also Cerviño, Luridiana & Castander 2000).

## 5. References

## References

Barbaro, G. & Bertelli, G. 1977, A&A, 54,243

Bruzual A., G. 2000, in Proceedings of the XI Canary Islands Winter School of Astrophysics on “Galaxies at High Redshift”, eds. I. Pérez-Fournon and F. Sánchez, Cambridge Contemporary Astrophysics, in press, (BC2000, astro-ph/0011094)

Bruzual A., G. 2001 in Proceedings of the Euroconference “The Evolution of Galaxies. I - Observational Clues”,eds. G. Stasińska and J. Vílchez, Dordrecht: Kluwer Academic Publishers, in press

Bruzual A., G. & Charlot, S. 2001, in preparation

Cerviño, M., Luridiana, V. & Castander, F.J. 2000, A&A, 360, L5

Chiosi, C., Bertelli, G., & Bressan, A. 1988, A&A, 196, 84

Fagotto, F., Bressan, A., Bertelli, G., & Chiosi, C. 1994a, A&AS, 100, 647

Fagotto, F., Bressan, A., Bertelli, G., & Chiosi, C. 1994b, A&AS, 104, 365

Fagotto, F., Bressan, A., Bertelli, G., & Chiosi, C. 1994c, A&AS, 105, 29

Frogel, J.A., Mould, J., & Blanco, V.M. 1990, ApJ, 352, 96

Girardi, L., Chiosi, C., Bertelli, G., & Bressan, A. 1995, A&A, 298, 87

Lejeune, T., Cuisinier, F., & Buser, R. 1998, A&AS, 130, 65

Maraston, C. 1998, MNRAS, 300, 872

Maraston, C., Kissler-Patig, M., Brodie, J.P., Barmby, P., & Huchra, J.P. 2001, A&A, 370, 176

Persson, S.E., Aaronson, M., Cohen, J.G., Frogel, J.A. & Mathews, K. 1983, ApJ, 266, 105

Press, W.H., Teukolsky, S.A., Vetterling, W.T., & Flannery, B.P. 1992, Numerical Recipes in FORTRAN, 2nd ed., Cambridge: Cambridge Univ. Press

Renzini, A., & Buzzoni, A. 1986, in “Spectral Evolution of Galaxies”, eds. C. Chiosi and A. Renzini, Dordrecht: Reidel, p. 195

Salpeter, E.E. 1955, ApJ, 121, 161

Santos, Jr, J.F.C. & Frogel, J.A. 1997, ApJ, 479, 764

van den Bergh, S. 1981, A&AS, 46, 79