text:
# Is the analysis of flow at the CERN SPS reliable?

## I Introduction

In ultrarelativistic heavy ion collisions, azimuthal distributions of the produced particles with respect to the reaction plane are of particular interest. Their anisotropies, which we refer to as “flow”, result from the interactions between the produced particles. They thus probe the hot stages of the collision and may provide valuable information on the state of the interaction region: thermalized or not, equation of state, etc…  In the recent years, flow has been measured by many experiments at the Brookhaven AGS with Au projectiles  and at the CERN SPS with Pb projectiles . Since the orientation of the reaction plane is not known a priori, flow measurements are usually extracted from two-particle azimuthal correlations. This is based on the idea that azimuthal correlations between two particles are generated by the correlation of the azimuth of each particle with the reaction plane.

The assumption that this is the only source of two particle azimuthal correlations dates back to the early days of the flow . It still underlies the analyses done at ultrarelativistic energies, both at AGS and at SPS (see, however, ). However, there are various sources of direct azimuthal correlations between particles, which do not involve the reaction plane. We have shown in a recent paper  that the quantum correlations due to the HBT effect yield azimuthal correlations. When they are taken into account, the results obtained for the charged pion flow by the NA49 collaboration  change significantly. Here, we investigate other sources of non–flow azimuthal correlations. A first source is the condition that the total transverse momentum of the outgoing particles is zero: this gives a back-to-back correlation between their momenta, which can be evaluated quantitatively. This effect is usually taken into account at lower energies  but has been neglected so far at AGS and SPS. Other correlations are due to resonance decays such as $`\mathrm{\Delta }p\pi `$ and $`\rho \pi \pi `$. Finally, Coulomb and strong interactions between pairs of particles with low relative velocities (which we refer to as “final state interactions”) produce small angle azimuthal correlations.

Most of these correlations are typically of order $`1/N`$, where $`N`$ is the number of particles emitted in the collision. At SPS energies, correlations due to flow are so small that additional correlations become of the same order, and cannot be neglected. The $`1/N`$ dependence also determines the variation of non–flow correlations with the centrality of the collision. Therefore they increase with impact parameter up to very peripheral collisions; on the other hand, they do not vanish for central collisions where the flow is zero by symmetry. They are very important at large impact parameters and should be taken into account, in particular, when studying the centrality dependence of the flow, which has been recently proposed as a sensitive probe of the phase transition to the quark–gluon plasma .

In section II, we briefly recall how the flow can be extracted from two–particle azimuthal correlations. We use the same method and notations as in . The various sources of non–flow correlations are reviewed and their effects are estimated in section III. The experimental data obtained by the NA49 and WA98 collaborations for Pb–Pb collisions at SPS are discussed in section IV. We show in particular that the effect of momentum conservation alone is large enough to reverse the sign of the proton directed flow measured by NA49 , which is used to define the reaction plane. $`\mathrm{\Delta }`$ decays produce correlations which are of the opposite sign, and of the same order of magnitude, although we do not attempt to calculate them accurately. For WA98, final state interactions play an important part, and may explain a large fraction of the observed pion directed flow. Our conclusions are given in section V.

## II Two–particle azimuthal correlations and standard flow analysis

In this section, we define the flow and non–flow contributions to two–particle azimuthal correlations, and we show how azimuthal distributions with respect to the reaction plane can be extracted from these correlations.

We call “flow” the azimuthal correlations between the outgoing particles and the reaction plane. These are best conveniently characterized in terms of the Fourier coefficients $`v_n`$  which we now define. We choose a coordinate system in which the $`x`$ axis is the impact direction, and $`(x,z)`$ the reaction plane. $`\varphi `$ denotes the azimuthal angle with respect to the reaction plane. In this frame, $`v_n`$ can be expressed as a function of the one-particle momentum distribution $`dN_j/d^3𝐩`$ for a particle of type $`j`$ (in this paper, we consider pions and protons):

$$v_n(p_T,y,j)\mathrm{cos}n\varphi =\frac{{\displaystyle _0^{2\pi }}\mathrm{cos}n\varphi {\displaystyle \frac{dN_j}{d^3𝐩}}d\varphi }{{\displaystyle _0^{2\pi }}{\displaystyle \frac{dN_j}{d^3𝐩}}𝑑\varphi }$$
(1)
where the brackets denote an average value over many events with approximately the same impact parameter, and $`p_T`$ and $`y`$ are the transverse momentum and rapidity of the particle. Since the system is symmetric with respect to the reaction plane for spherical nuclei, $`\mathrm{sin}n\varphi `$ vanishes. The goal of the flow analysis is to extract $`v_n`$ from the data. The coefficients $`v_1`$ and $`v_2`$ are usually called directed and elliptic flow, respectively .

Since the reaction plane is not known experimentally, one can only measure relative azimuthal angles between the outgoing particles. In particular, one measures the Fourier coefficients of the relative azimuthal distribution between two species of particles $`j`$ and $`k`$

$$c_n^{\mathrm{measured}}(p_{T1},y_1,j;p_{T2},y_2,k)\mathrm{cos}n(\varphi _1\varphi _2)=\frac{{\displaystyle \mathrm{cos}n(\varphi _1\varphi _2)\frac{dN_{jk}}{d^3𝐩_1d^3𝐩_2}𝑑\varphi _1𝑑\varphi _2}}{{\displaystyle \frac{dN_{jk}}{d^3𝐩_1d^3𝐩_2}𝑑\varphi _1𝑑\varphi _2}}.$$
(2)
The two-particle distribution can generally be expressed as the sum of an uncorrelated distribution and two-particle correlations:

$$\frac{dN_{jk}}{d^3𝐩_\mathrm{𝟏}d^3𝐩_\mathrm{𝟐}}=\frac{dN_j}{d^3𝐩_\mathrm{𝟏}}\frac{dN_k}{d^3𝐩_\mathrm{𝟐}}\left(1+C_{jk}(𝐩_1,𝐩_2)\right)$$
(3)
where $`C_{jk}(𝐩_1,𝐩_2)`$ is the two-particle connected correlation function. Using this decomposition, one can write $`c_n^{\mathrm{measured}}`$ as the sum of two terms:

$$c_n^{\mathrm{measured}}(p_{T1},y_1,j;p_{T2},y_2,k)=c_n^{\mathrm{flow}}(p_{T1},y_1,j;p_{T2},y_2,k)+c_n^{\mathrm{non}\mathrm{flow}}(p_{T1},y_1,j;p_{T2},y_2,k)$$
(4)
where the first term is due to flow:

$$c_n^{\mathrm{flow}}(p_{T1},y_1,j;p_{T2},y_2,k)=v_n(p_{T1},y_1,j)v_n(p_{T2},y_2,k)$$
(5)
while the remaining term comes from direct two-particle correlations:

$$c_n^{\mathrm{non}\mathrm{flow}}(p_{T1},y_1,j;p_{T2},y_2,k)=\frac{{\displaystyle \mathrm{cos}n(\varphi _1\varphi _2)C_{jk}(𝐩_1,𝐩_2)\frac{dN_j}{d^3𝐩_1}\frac{dN_k}{d^3𝐩_2}𝑑\varphi _1𝑑\varphi _2}}{{\displaystyle \frac{dN_{jk}}{d^3𝐩_1d^3𝐩_2}𝑑\varphi _1𝑑\varphi _2}}.$$
(6)
From the flow term (5), one can calculate the Fourier coefficient $`v_n`$, up to a sign (see  for details):

$$v_n(p_T,y,j)=\pm \frac{c_n^{\mathrm{flow}}(p_T,y,j;𝒟)}{\sqrt{c_n^{\mathrm{flow}}(𝒟;𝒟)}}$$
(7)
where $`𝒟`$ denotes the detector used to estimate the reaction plane, $`c_n(p_{T1},y_1,j;𝒟)`$ the average value of $`c_n`$ over $`(p_{T2},y_2,k)`$ in $`𝒟`$, and $`c_n(𝒟;𝒟)`$ the average over both $`(p_{T1},y_1,j)`$ and $`(p_{T2},y_2,k)`$.

Experimental analyses usually assume that the only azimuthal correlation between outgoing particles is due to their correlation with the reaction plane, i.e. they neglect $`c_n^{\mathrm{non}\mathrm{flow}}`$. This means that the published flow data are given by Eq.(7) with $`c_n^{\mathrm{flow}}`$ replaced by the total correlation $`c_n^{\mathrm{measured}}`$; in turn, the total correlation can be reconstructed from the published flow data using Eq.(5).

We shall see in section III that $`c_n^{\mathrm{non}\mathrm{flow}}`$ is typically of order $`1/N`$ where $`N`$ denotes the number of outgoing particles. In a central Pb–Pb collision at SPS, $`N2500`$, so that one expects $`c_n^{\mathrm{non}\mathrm{flow}}4\times 10^4`$. This is certainly a weak effect. However, the total azimuthal correlation is also weak: since the directed and elliptic flow at SPS measured by NA49 are of the order of $`3\%`$ , one obtains from Eq.(5) (since it is assumed that all the measured correlation is from flow) $`c_n^{\mathrm{measured}}9\times 10^4`$, only a factor of 2 larger than non-flow correlations. Therefore one cannot neglect the latter. In order to take them into account, one must estimate the various contributions to $`c_n^{\mathrm{non}\mathrm{flow}}`$, subtract them from $`c_n^{\mathrm{measured}}`$ in (4) to isolate $`c_n^{\mathrm{flow}}`$, and calculate the resulting $`v_n`$ using (7). In , we applied this procedure to the NA49 data on pion flow , taking into account correlations due to the HBT effect. In section IV, we consider in addition the effects of momentum conservation and resonance decays.

## III Sources of non-flow correlations

Several physical effects yield two–particle correlations, which contribute to the non-flow correlation defined in (6). In this section, we discuss momentum conservation, resonance decays, quantum HBT correlations, and final state interactions, and we estimate their respective contributions to $`c_n^{\mathrm{non}\mathrm{flow}}`$.

### A Transverse momentum conservation

Global momentum conservation gives a back-to-back correlation between the produced particles. Since we are interested in azimuthal correlations, we consider only transverse momentum. Its contribution to $`C_{jk}(𝐩_1,𝐩_2)`$, which we denote by $`C_{jk}^{\mathrm{\Sigma }p_T}(𝐩_1,𝐩_2)`$, is calculated in Appendix A in the case where there is no other correlation between particles. The result is

$$C_{jk}^{\mathrm{\Sigma }p_T}(𝐩_1,𝐩_2)=\frac{2𝐩_{T1}𝐩_{T2}}{p_T^2}$$
(8)
where the sum in the denominator runs over all the particles emitted in the collision. This correlation is clearly of order $`1/N`$. The Fourier coefficient $`c_n^{\mathrm{\Sigma }p_T}`$ is obtained by inserting (8) in Eq.(6). Neglecting azimuthal anisotropies of the one-particle distribution, which are of the order of a few percent , one finds a result which is non-vanishing only for $`n=1`$:

$$c_1^{\mathrm{\Sigma }p_T}(p_{T1},y_1,j;p_{T2},y_2,k)=\frac{p_{T1}p_{T2}}{p_T^2}.$$
(9)
Since $`c_n^{\mathrm{\Sigma }p_T}`$ vanishes for $`n1`$, the effects of momentum conservation must be taken into account only in the measurement of directed flow. Eq.(9) allows to calculate the correlation arising from momentum conservation, as soon as the value of the denominator $`p_T^2`$ is known. This quantity is evaluated in Appendix B in the case of Pb–Pb collisions at SPS. The value (B3) is probably overestimated, which means that we underestimate the correction due to momentum conservation in the numerical calculations of section IV. Eq.(9) shows that the azimuthal correlations between two particles due to momentum conservation are independent of the rapidities $`y_1`$ and $`y_2`$ and of the particle species $`j`$ and $`k`$, but increase linearly with transverse momentum. The effect is therefore stronger for heavier particles, which have larger transverse momenta, as for instance protons. The centrality dependence is also easily determined: since transverse momentum spectra depend weakly on the centrality , $`p_T^2`$ scales like the multiplicity $`N`$, so that $`c_1^{\mathrm{\Sigma }p_T}`$ scales like $`1/N`$. Note, finally, that momentum conservation does not contribute if the window used for the reaction plane determination is symmetric around mid-rapidity .

### B Resonance decays

We now discuss the azimuthal correlation between the decay products of a resonance. It is known that many pions originate from $`\rho \pi \pi `$ decays, and that many nucleons which are excited into $`\mathrm{\Delta }`$ resonances, which then decay into a nucleon and a pion. When pions are used to determine the reaction plane, $`\rho `$ (resp. $`\mathrm{\Delta }`$) decays must be taken into account in the analysis of the pion (resp. proton) flow; when protons are used, $`\mathrm{\Delta }`$ decays must be taken into account in the analysis of the pion flow. Resonances with higher masses are less abundant and will not be considered in this paper; however, they may also give rise to sizeable correlations.

The contribution of a resonance decay, say, $`\mathrm{\Delta }p\pi `$, to $`c_n^{\mathrm{non}\mathrm{flow}}`$ is the product of two factors: the first factor is the average value of $`\mathrm{cos}n(\varphi _1\varphi _2)`$ where $`\varphi _1`$ and $`\varphi _2`$ are the azimuthal angles of a proton and a pion originating from the same $`\mathrm{\Delta }`$; the second factor is the probability that a given proton and a given pion originate from the same $`\mathrm{\Delta }`$. This probability is hard to estimate. There are recent measurements of $`\rho ^0`$  and $`\mathrm{\Delta }^{++}`$  multiplicities, but this is not sufficient: indeed, due to the short resonance lifetimes, the decay products may interact again after they have been produced, in a way which depends on the detailed collision dynamics. Microscopic models of heavy ion collisions (see  for a review) could be used to calculate correlations from resonance decays, but the discussion in this paper will remain at a more qualitative level.

Generally, in a two-particle decay, the azimuthal correlation depends on the relative magnitude of the velocity $`V`$ of the decaying particle and the velocities $`V_1`$ and $`V_2`$ of the products in the rest frame of the decaying particle: in particular, the azimuthal angles of the decay products are equal ($`\varphi _1=\varphi _2`$) in the limit where $`V_1=V_2=0`$, while they are back-to-back ($`\varphi _1=\varphi _2+\pi `$) if $`V=0`$. This is illustrated in Fig.1, which displays the average value of $`\mathrm{cos}(\varphi _1\varphi _2)`$ and $`\mathrm{cos}2(\varphi _1\varphi _2)`$ for the decays $`\mathrm{\Delta }p\pi `$ and $`\rho \pi \pi `$, as functions of the transverse momentum of the decaying resonance.

Assuming that the transverse momentum distribution is exponential in the transverse mass $`m_T=\sqrt{p_T^2+m^2}`$, i.e. $`dN/d^2𝐩_T\mathrm{exp}(m_T/T)`$ (see Eq.(B1)), one can average the correlations displayed in Fig.1 over $`p_T`$. The result is shown in Fig.2, as a function of the inverse slope parameter $`T`$.

Taking for $`\mathrm{\Delta }`$ the same value of $`T`$ as for protons, i.e. $`T300`$ MeV, one sees that the contribution to $`c_1`$ is large and positive, while the contribution to $`c_2`$ is close to zero. In the case of $`\rho `$ decays, assuming that $`T`$ varies linearly with the mass  one may choose $`T270`$ MeV, which results in both $`c_1`$ and $`c_2`$ slightly negative.

In order to obtain the contribution of $`\mathrm{\Delta }`$ decays to $`c_n^{\mathrm{non}\mathrm{flow}}`$, one must multiply the correlation displayed in Fig.(1) by the probability $`p_\mathrm{\Delta }`$ that a given pion and a given proton originate from the same $`\mathrm{\Delta }`$, as explained above. Let $`\alpha `$ denote the fraction of nucleons are excited into $`\mathrm{\Delta }`$ resonances. Then, the probability $`p_\mathrm{\Delta }`$ is $`\alpha /N_\pi `$ where $`N_\pi `$ is the number of pions. From isospin symmetry, one expects this probability to be the same for neutral and charged pions. Since most particles are pions, $`N_\pi N`$. If $`\alpha `$ is of order unity, we obtain finally $`c_n^{\mathrm{\Delta }p\pi }1/N`$. The centrality dependence of this correlation is also in $`1/N`$ if the relative abundance of $`\mathrm{\Delta }`$ is independent of centrality. Similar arguments hold for the correlations due to $`\rho `$ decays.

Note that a way to eliminate correlations from $`\mathrm{\Delta }`$ decays experimentally would be to measure protons and pions in widely separated rapidity windows. Indeed, the maximum rapidity difference between the decay products of a $`\mathrm{\Delta }`$ is $`1.5`$. In the case of the $`\rho \pi \pi `$ decay, the maximum rapidity difference between the outgoing pions is $`3.3`$, so that the corresponding correlation is harder to eliminate.

### C Small angle correlations

Finally, there are correlations between particles with low relative velocities: quantum HBT correlations between identical particles, Coulomb interactions between charged particles, and strong interactions.

Azimuthal correlations between identical particles due to the HBT effect were studied in detail in . The Fourier coefficient $`c_n^{\mathrm{HBT}}(𝐩_\mathrm{𝟏},𝐩_\mathrm{𝟐})`$ is of order unity if the relative momentum $`|𝐩_1𝐩_2|`$ is smaller than $`\mathrm{}/R`$, where $`R`$ is the size of the particle source. It depends weakly on $`n`$. The order of magnitude of $`c_n^{\mathrm{HBT}}`$ integrated over phase space is therefore $`c_n^{\mathrm{HBT}}(\mathrm{}/Rp)^3`$, where $`p`$ denotes a typical momentum. In the case of a high temperature massless pion gas, $`p`$ is of order $`T/c`$, while the particle density $`N/R^3`$ is of order $`(T/\mathrm{}c)^3`$, so that $`c_n^{\mathrm{HBT}}`$ is again of order $`1/N`$ . One would expect the centrality dependence of $`c_n^{\mathrm{HBT}}`$ to follow the same behavior. However, the measured HBT radii vary more slowly than $`N^{1/3}`$ . In collisions with sulphur projectiles, a dependence in $`RN^\alpha `$ was found with $`\alpha 0.2`$ . If the same behavior holds for Pb–Pb collisions, the dependence of azimuthal HBT correlations with centrality will be in $`N^{0.6}`$. Finally, $`c_n^{\mathrm{HBT}}(𝒟)`$ is proportional to the one–particle momentum distribution, hence it is larger at low $`p_T`$.

More generally, all particles with low relative velocities undergo final state interactions, strong and/or electrostatic . One defines the invariant momentum as $`Q_{\mathrm{inv}}=2\mu v_{\mathrm{rel}}`$, where $`\mu =m_1m_2/(m_1+m_2)`$ is the reduced mass and $`v_{\mathrm{rel}}=\sqrt{1m_1^2m_2^2/(P_1P_2)^2}`$ is the relative velocity of the two particles with 4-momenta $`P_1`$ and $`P_2`$. Then the correlations due to strong interactions are sizeable only in a limited $`Q_{\mathrm{inv}}`$ range, independent of the size of the interaction region. For instance, proton–proton correlations, which are dominated by strong interactions, are maximum for $`Q_{\mathrm{inv}}20`$ MeV/c, and negligible for $`Q_{\mathrm{inv}}40`$ MeV/c, for any projectile and target. On the other hand, the range over which Coulomb correlations are significant depends on the system size. For Pb–Pb collisions at SPS, $`p\pi ^+`$ correlations, which are dominated by the Coulomb repulsion, extend up to $`Q_{\mathrm{inv}}50`$ MeV/c . Both Coulomb and strong interactions give rise to small angle azimuthal correlations, similar to those due to the HBT effect. Thus, the resulting contribution to $`c_n^{\mathrm{non}\mathrm{flow}}`$, denoted by $`c_n^{\mathrm{FSI}}`$, is expected to be weakly dependent on the order $`n`$ of the Fourier component, and $`c_n^{\mathrm{FSI}}(p_T,y,𝒟)`$ is roughly proportional to the one–particle distribution.

The centrality dependence differs for Coulomb and strong interactions. For strong interactions, the strength of the correlation at a given $`Q_{\mathrm{inv}}`$ varies like $`1/V`$, where $`V`$ is the volume of the system . This is due to the fact that the range of the interaction is much smaller than the size of the system. Thus the centrality dependence follows a $`1/N`$ behavior. For Coulomb interactions, which are long–ranged, the $`N`$ dependence is weaker. We do not attempt to evaluate this dependence. By analogy with the Coulomb potential, it might vary like $`1/L`$, with $`L`$ the size of the interaction region, i.e. like $`1/N^{1/3}`$.

Both HBT and final state interactions can be eliminated by studying the correlation between particles widely separated in phase space, so that their invariant relative momentum is much larger than, typically, 50 MeV/c.

### D Summary

The properties of the three main types of two–particle correlations studied above are summarized in Table I. Note that while the average over all phase space of the three is of order $`1/N`$, they can be much larger in definite regions of phase space, for example at low relative momentum for HBT correlations. Note also that this catalogue of two-particle correlations may not be exhaustive.

As we have seen, two–particle correlations are generally of order $`1/N`$, much smaller than unity. We therefore assume that the contributions of the various effects can simply be added. If one studies the correlation between identical pions, for instance, there are contributions to $`c_n^{\mathrm{non}\mathrm{flow}}`$ from momentum conservation, denoted by $`c_n^{\mathrm{\Sigma }p_T}`$, from $`\rho \pi \pi `$ decays, denoted by $`c_n^\rho `$, from HBT correlations, denoted by $`c_n^{\mathrm{HBT}}`$, and from final state interactions, denoted by $`c_n^{\mathrm{FSI}}`$. The resulting non-flow correlation is

$$c_n^{\mathrm{non}\mathrm{flow}}=c_n^{\mathrm{\Sigma }p_T}+c_n^\rho +c_n^{\mathrm{HBT}}+c_n^{\mathrm{FSI}}.$$
(10)
|  | Momentum conservation | Resonance decays | Small angle correlations |
| --- | --- | --- | --- |
| Which pairs? | any | $`|y_1y_2|<1.5`$ ($`\mathrm{\Delta }`$) | low relative velocities |
|  |  | $`|y_1y_2|<3.3`$ ($`\rho `$) | identical particles (HBT) |
| $`n`$ dependence | $`n=1`$ only | non-trivial | weak |
| centrality dependence | $`1/N`$ | $`1/N`$ | $`1/N^{0.6}`$ ? (HBT) |
|  |  |  | $`1/N`$ (strong) |
|  |  |  | $`1/N^{1/3}`$ ? (Coulomb) |
| $`p_T`$ range | high $`p_T`$ | non-trivial | low $`p_T`$ |
| under control? | calculated | simulations, microscopic models | measured |

TABLE I: Properties of the various sources of two–particle azimuthal correlations

## IV Application to SPS data

The most detailed flow analyses in Pb–Pb collisions at SPS were performed by the NA49 and WA98 collaborations. NA49 estimates the reaction plane using pions in the forward hemisphere, while WA98 uses protons and fragments close to the target rapidity. In this section, we discuss to what extent non–flow correlations may contribute to the measured azimuthal correlations, which are assumed to be solely due to flow in both experiments. In the case of the NA49 data, we shall explicitly subtract the various non-flow correlations, following the method outlined in section II, and see what the “true” flow might be. Our discussion of WA98 data will remain at a semi-quantitative level due to the complex acceptance of the Plastic Ball detector used for the flow analysis.

### A NA49 Data

The NA49 experiment at CERN measures the directed and elliptic flow of pions and protons in Pb–Pb collisions at 158 GeV per nucleon . Charged pions are used to estimate the reaction plane, in the kinematic window $`4<y<6`$ and $`0<p_T<600`$ MeV/c for directed flow, $`3.5<y<5`$ and $`0<p_T<2000`$ MeV/c for elliptic flow . The azimuthal distribution of identified particles (protons and charged pions) is then measured with respect to this reaction plane. We now show the effect of taking into account the various sources of non–flow correlations discussed in the previous section.

We begin with the directed flow of protons. It is obtained from from Eq.(7), where $`𝒟`$ refers to charged pions in the above mentioned kinematic window. The numerator in this equation corresponds to the correlations between protons and the pions in $`𝒟`$, while the denominator corresponds to the correlations between the pions in $`𝒟`$. The arbitrary sign in Eq.(7) is chosen negative by NA49, in order to obtain a positive $`v_1`$ for protons (see Fig.3), as usually assumed at high energies in the forward rapidity region.

The correction to the denominator in Eq.(7) from non-flow correlations will be discussed below, when we discuss the pion flow. It is a small correction, of a few percent. The numerator of Eq.(7) gets contributions from momentum conservation, $`\mathrm{\Delta }p\pi `$ decays, and final state interactions.

We first discuss momentum conservation. The corresponding azimuthal correlation is given by Eqs.(9) and (B3). Since it is negative, once it is subtracted from the measured correlation (4), this gives a positive correction to the correlations due to flow $`c_1^{\mathrm{flow}}`$. Because of the negative sign in (7), the correction to $`v_1`$ is finally negative. As can be seen in Fig.3, the correction is so large that the sign of the proton directed flow is now negative for almost all $`p_T`$. Now, the sign of the proton directed flow is particularly important since it is used to fix the arbitrary sign in Eq.(7), which is chosen so that $`v_1>0`$ in the forward rapidity region for protons. The surprising result is that the directed flows of pions and protons now have the same sign. This is contrary to the theoretical expectation that nucleons and pions flow in opposite directions, because of pion rescattering , and to other measurements at AGS  and at SPS .

Since the correlation between pions and protons from $`\mathrm{\Delta }p\pi `$ decays is positive (see Fig.2), it may restore the original sign of the proton directed flow $`v_1`$. As explained in the previous section, an accurate computation of this correlation would require a microscopic model. We present in Fig.3 a crude estimate of the effect, based on a simple Monte-Carlo calculation taking into account the acceptance windows of NA49 for protons and pions. We do not take into the possible rescatterings of the decay products. Our assumptions are the following: 50% of the outgoing protons originate from $`\mathrm{\Delta }`$ decays, as suggested by simulations based on the UrQMD model . We take into account the fact that only $`2/3`$ of $`\mathrm{\Delta }`$ decays yield a charged pion. The mass distribution of $`\mathrm{\Delta }`$’s follows a lorentzian, their rapidity distribution is flat between 0 and 5.8 (i.e. between the target and projectile rapidities), and their transverse momentum distribution is exponential in $`m_T`$ (see Eq.(B1)) with an inverse slope parameter $`T=300`$ MeV, as for protons. One sees in Fig.3 that correlations from $`\mathrm{\Delta }`$ decays can indeed restore the positive sign for the proton directed flow, since their contribution is almost as large as that from momentum conservation in absolute value.

In principle, one should also take into account final state interactions. Since the pions used to determine the reaction plane are approximately half $`\pi ^+`$ and half $`\pi ^{}`$, Coulomb effects cancel. Correlations between pions and protons due to strong interactions are small and should not contribute significantly to $`c_n^{\mathrm{measured}}`$.

We now turn to the directed flow of charged pions. Since pions are also used to estimate the reaction plane, pion–pion correlations are involved in both the numerator and the denominator of Eq.(7). Following Eq.(10), we must take into account momentum conservation, $`\rho `$ decays, the HBT effect, and final state interactions. We neglect final state interactions for the same reasons as for protons. According to Fig.2, the correlation from $`\rho `$ decay is small. We also neglect it. Three sets of points are displayed in Fig.4: NA49 data without correction, with correction for HBT correlations (taken from ) and with correction for HBT and momentum conservation. While the correction due to the HBT effect is important at low $`p_T`$, momentum conservation changes the pion flow at high $`p_T`$; the latter effect is less important than for protons, because pions have lower transverse momenta. However, the positive pion flow observed by NA49 at high $`p_T`$ seems to be explained by momentum conservation.

Let us say a few words of elliptic flow. Here, momentum conservation does not contribute, as discussed in the previous section. For protons, $`\mathrm{\Delta }`$ decays probably give a negligible contribution to $`c_2^{\mathrm{measured}}`$, as explained in section III B. For pions, HBT correlations significantly reduce elliptic flow at low $`p_T`$, as shown in . Here, $`\rho `$ decays should be taken into account. As discussed in section III B, they give a negative correlation, i.e. opposite to the correlation due to flow. Therefore, taking $`\rho `$ decays into account would lead to a (slight) increase of the pion $`v_2`$.

Finally, the centrality dependence of the directed and elliptic flow of pions has been studied recently . It is striking to note that both increase continuously as the reaction becomes more peripheral, which is the behavior expected for non–flow correlations. We have shown that at least a fraction of the measured correlation, which is interpreted as flow, is due to non–flow correlations. Even if it is only a small fraction, it will increase with impact parameter and eventually dominate for the most peripheral collisions. It is thus likely that the measurements of flow at large impact parameter are contaminated by non–flow phenomena. However, none of the above mechanisms can explain the observed $`y`$ dependence of directed flow, which increases (in absolute value) up to the projectile rapidity for peripheral collisions.

### B WA98 Data

We shall now discuss the results obtained by the WA98 experiment at CERN for the directed flow of protons and positive pions . There, the reaction plane is estimated using protons and fragments detected by the Plastic Ball, in the target rapidity region $`0.6<y<0.3`$. One then measures the azimuthal angle of a positive pion or a proton detected in the Plastic Ball with respect to this estimated plane. If the particle is a proton, it is not used in the determination of the reaction plane to avoid autocorrelations. Two points closer to midrapidity are obtained using the tracking system, for which the flow is zero within error bars; we shall not consider these two points here. Contrary to the NA49 analysis where pions were used to estimate the reaction plane, the flow of protons (resp. pions) is now determined from proton–proton (resp. pion–proton) azimuthal correlations. The value of $`v_1`$ is still obtained from Eq.(7), where the arbitary sign is again chosen negative, so that the directed flow of protons be negative, as it should in the backward rapidity region. On the other hand, the directed flow of pions is positive, i.e. pions flow opposite to the protons.

WA98 measures the directed flow $`v_1`$ in the target fragmentation region. Since $`v_1`$ vanishes at midrapidity, it is expected to be larger there than in the more central rapidity region, where NA49 is working. Indeed, the measured values of $`v_1`$ are larger for WA98 than for NA49. For semi-central collisions, $`|v_1|`$ is about $`0.1`$ for both pions and protons. Following Eq.(5), this corresponds to an azimuthal correlation $`|c_1^{\mathrm{measured}}|10^2`$, either between pions and protons or between protons, an order of magnitude larger than for NA49.

The values of $`v_1`$ are measured as a function of centrality. For protons, $`v_1`$ reaches a maximum of $`0.12`$ for an impact parameter $`b8`$ fm, while for pions $`v_1`$ increases with $`b`$ up to very peripheral collisions, where it reaches values as high as $`0.3`$. The latter centrality dependence is reminiscent of the $`1/N`$ behavior of non–flow correlations, which suggests that part of the measured flow is spurious. We thus review the various non–flow correlations listed in section III to evaluate their contributions. We shall see that final state interactions are the only significant ones, and that they may explain the pion flow at large impact parameter.

Momentum conservation is less important than for the NA49 experiment. The reason is that the particles detected in the Plastic Ball have low transverse momenta . For protons, $`p_T`$ is typically 400 MeV, resulting in $`c_1^{\mathrm{\Sigma }p_T}`$ of the order of $`4\times 10^4`$ for the proton–proton correlation, following Eq.(9). This is much smaller than $`c_1^{\mathrm{measured}}10^2`$. For pions, transverse momenta are even smaller. In both cases, the contribution of transverse momentum conservation to $`c_1^{\mathrm{non}\mathrm{flow}}`$ is negligible.

Resonance decays will contribute to the pion flow. The sign of the correlation due to $`\mathrm{\Delta }p\pi `$ depends on the inverse slope parameter $`T`$, as can be seen in Fig.2. This quantity is smaller in the fragmentation region than in the central rapidity region, but probably still larger than 100 MeV. Then, the corresponding pion–proton correlation should be positive. It is thus opposite to the correlation measured by WA98, and cannot explain the observed pion flow.

Finally, let us discuss the small angle correlations due to the HBT effect and final state interactions. We shall first summarize the results obtained by the NA49 experiment for two-particle small angle correlations in the central rapidity region. Then we extrapolate these results to the fragmentation region covered by the WA98 Plastic Ball in order to discuss the correlations measured by WA98.

The NA49 collaboration has measured both proton-proton  and proton-$`\pi ^+`$  small angle correlations for central collisions. Two-proton correlations, resulting from the interplay of the HBT effect and strong and Coulomb interactions, were found to be positive, which means that they are dominated by the attractive strong interaction. They peak at a relative momentum between particles of 20 MeV/c, where they reach $`25\%`$, and are no longer significant as soon as $`|𝐩_1𝐩_2|>40`$ MeV/c. Small angle correlations between positive pions and protons are negative, and typically of the order of $`10\%`$ for invariant momenta smaller than 40 MeV/c, i.e. $`|𝐩_{T\pi }(m_\pi /m_p)𝐩_{Tp}|<40`$ MeV/c and $`y_\pi y_p<40/m_\pi 0.3`$.

Extrapolating these results to the fragmentation region, we can estimate the azimuthal correlation $`c_1^{\mathrm{FSI}}`$ for WA98. In the case of proton–proton correlations, we need to evaluate the fraction of proton pairs having relative momenta less than 40 MeV/c, i.e. with $`|𝐩_{T1}𝐩_{T2}|<\mathrm{\Delta }p=40`$ MeV/c and $`y_1y_2<\mathrm{\Delta }p/m_T0.04`$, with $`m_T1000`$ MeV. Since the acceptance window is typically  $`p_T^{\mathrm{min}}=200<p_T<p_T^{\mathrm{max}}=600`$ MeV/c and $`y^{\mathrm{min}}=0.4<y<y^{\mathrm{max}}=0.1`$, this fraction

$$\frac{\frac{4\pi }{3}\mathrm{\Delta }p^3}{\pi (p_T^{\mathrm{max}2}p_T^{\mathrm{min}2})m_T(y^{\mathrm{max}}y^{\mathrm{min}})}5\times 10^4.$$
(11)
With a correlation strength of 20%, this yields a value of $`c_1^{\mathrm{FSI}}`$ of the order of $`10^4`$, much smaller than the measured correlation $`c_1^{\mathrm{measured}}`$. Thus, the measurement of the proton flow is not influenced by small angle correlations.

In the case of pion–proton correlations, we first note that the velocity windows of pions and protons of the WA98 Plastic Ball overlap , so that small angle correlations between pions and protons should be taken into account. In fact, the rapidity range is almost the same for pions and protons, while the transverse momentum range is roughly scaled by a factor $`m_\pi /m_p`$, so that the acceptance windows roughly coincide in velocity space. Then, for a given pion, the probability that a proton seen in the Plastic Ball has an invariant relative momentum $`|𝐩_{T\pi }(m_\pi /m_p)𝐩_{Tp}|<\mathrm{\Delta }p=40`$ MeV/c is given by (11) with $`\mathrm{\Delta }p`$ replaced by $`(m_p/m_\pi )\mathrm{\Delta }p`$. This gives a factor $`(m_p/m_\pi )^3200`$, compared to Eq.(11), so that the fraction is of the order of $`10^1`$. Using the correlation strength measured by NA49, this gives a correlation $`c_1^{\mathrm{FSI}}(𝒟)`$ of order $`10^2`$, comparable to the measured correlation. On the one hand, this value is overestimated since the acceptance windows for pions and protons do not exactly coincide in velocity space. On the other hand, it is underestimated because the correlations are measured by NA49 for central collisions, and they are even stronger for peripheral collisions. Therefore it is likely that a large fraction of the $`\pi ^+`$ directed flow seen by WA98 is due to the repulsive Coulomb interaction between positive pions and protons.

## V Discussion

We have shown that the assumption that all azimuthal correlations between particles are due to flow is no longer valid at SPS energies, where other sources of correlations become of the same order of magnitude due to the smallness of the flow. All the effects studied above, momentum conservation, resonance decays, HBT correlations, final state interactions, turn out to give important correlations. In particular, we have shown that the effect of momentum conservation is so large that it changes the sign of the proton directed flow measured by NA49: this sign is particularly important since it is used to define the orientation of the reaction plane.

The methods currently used to analyse the flow can be modified in order to take into account these additional correlations and to subtract them from the measured correlations, as explained in section II. Correlations from momentum conservation can be calculated. They are governed by only one unknown quantity, the sum over all particles of squared transverse momenta $`p_T^2`$. Correlations from the decays of short-lived resonances cannot be evaluated so easily. In this paper, we have only given semi-quantitative estimates. More accurate calculations would require a microscopic model of the collision. Finally, correlations which affect particles with low relative velocities, i.e. those due to the HBT effect and to final state (Coulomb, strong) interactions can either be measured independently, or eliminated by measuring correlations between particles in widely separated rapidity and/or transverse velocity windows.

However, our list of correlations may not be exhaustive, and one cannot exclude that other sources exist which are of the same order of magnitude as those studied here. There is no systematic way to separate the flow and the non–flow contribution to the measured azimuthal correlation (4). One could try to use the factorization property of Eq.(5) to isolate the correlation due to flow. However, non-flow correlations may also factorize, as is the case for the correlation due to momentum conservation (9).

A way to circumvent these difficulties would be to use the remarkable property of two–particle correlations, that they vary with the multiplicity $`N`$ like $`1/N`$ (with the possible exception of Coulomb interactions and, to a lesser extent, HBT correlations). This fact has been noted in section III for the various correlations studied, but holds more generally for any extensive system, and would probably still be true for other types of correlations than those listed in section III. This $`1/N`$ behavior determines the centrality dependence of non–flow two–particle correlations. Since the flow is known to vanish for central collisions by symmetry, one could measure (and then subtract) non-flow correlations by measuring them for central collisions and then assuming that they vary like $`1/N`$. Unfortunately, it is impossible to select “true” central collisions (with vanishing impact parameter) experimentally, so that this procedure may not be accurate.

Future heavy ion experiments at RHIC and LHC are likely to be affected in the same way by non–flow correlations. At these energies, directed flow is expected to be weaker than at SPS, maybe undetectable. On the other hand, elliptic flow should be at least of the same magnitude as at SPS. Pions should be used to determine the reaction plane since they are by far the most abundant. As the multiplicity will be higher than at SPS, non-flow correlations will be smaller, since they are proportional to $`1/N`$. However, since the detectors will cover only a limited interval in rapidity, effects of $`\rho `$ decays, HBT correlations and final state interactions may be important and should be taken into account in the measurement of the pion elliptic flow. Hard processes could give an additional contribution to non–flow correlations: there is a strong azimuthal correlation between jet fragments . Clearly, more work is needed, both theoretical and experimental, in order to obtain accurate measurements of flow observables at ultrarelativistic energies. Alternative methods, based on multiparticle correlations, are currently under study .

###### Acknowledgements.

We thank the Nuclear Theory group of Brookhaven National Laboratory for its hospitality during the period when this work was done. We thank Terry Awes, Art Poskanzer, Hubertus Schlagheck and Sergei Voloshin for detailed explanations concerning the NA49 and WA98 flow analyses, and Stefan Bass for discussions.

## A Correlations from momentum conservation

We denote by $`𝐩_1\mathrm{}𝐩_N`$ the momenta of the $`N`$ particles emitted in a heavy ion collision and by $`𝐩_{𝐓}^{}{}_{1}{}^{}\mathrm{}𝐩_{𝐓}^{}{}_{N}{}^{}`$ their transverse components. Since we are interested in azimuthal correlations, we consider only transverse momentum conservation: $`𝐩_{𝐓}^{}{}_{1}{}^{}+\mathrm{}+𝐩_{𝐓}^{}{}_{N}{}^{}=\mathrm{𝟎}`$. We assume for simplicity that this is the only correlation between the momenta of the outgoing particles. We want to calculate the distribution of $`𝐩_1,\mathrm{},𝐩_k`$, with $`k<N`$, defined as:

$$f_c(𝐩_1,\mathrm{},𝐩_k)\frac{\left({\displaystyle \underset{i=1}{\overset{k}{}}}f(𝐩_i)\right){\displaystyle \delta ^2(𝐩_{𝐓}^{}{}_{1}{}^{}+\mathrm{}+𝐩_{𝐓}^{}{}_{N}{}^{})\underset{i=k+1}{\overset{N}{}}\left(f(𝐩_i)d^3𝐩_i\right)}}{{\displaystyle \delta ^2(𝐩_{𝐓}^{}{}_{1}{}^{}+\mathrm{}+𝐩_{𝐓}^{}{}_{N}{}^{})\underset{i=1}{\overset{N}{}}\left(f(𝐩_i)d^3𝐩_i\right)}}$$
(A1)
where $`f(𝐩)`$ denotes the single particle normalized transverse momentum distribution. The average transverse momentum is naturally assumed to be zero:

$$𝐩_𝐓𝐩_𝐓f(𝐩)d^3𝐩=\mathrm{𝟎}.$$
(A2)
In order to calculate the distribution $`f_c`$ defined in Eq.(A1), we make use of the following lemma: the sum of $`M`$ uncorrelated momenta $`𝐏_𝐓_{i=1}^M𝐩_{𝐓}^{}{}_{i}{}^{}`$ has a gaussian distribution if $`M`$ is large, according to the central limit theorem:

$`F_M(𝐏_𝐓)`$ $``$ $`{\displaystyle \delta ^2\left(𝐏_𝐓+\underset{i=1}{\overset{M}{}}𝐩_{𝐓}^{}{}_{i}{}^{}\right)\underset{i=1}{\overset{M}{}}\left(f(𝐩_i)d^3𝐩_i\right)}`$ (A3)
$`=`$ $`{\displaystyle \frac{1}{\pi \sigma ^2}}\mathrm{exp}\left({\displaystyle \frac{𝐏_{𝐓}^{}{}_{}{}^{2}}{\sigma ^2}}\right).`$ (A4)
In writing this equation, we have assumed that the transverse momentum distribution is isotropic in the transverse plane, i.e. we have neglected the azimuthal asymmetries (1), which are of the order of a few percent. The width of the gaussian is

$$\sigma ^2=𝐏_{𝐓}^{}{}_{}{}^{2}=Mp_T^2.$$
(A5)
Eqs.(A3) and (A5) allow to write (A1) as

$`f_c(𝐩_1,\mathrm{},𝐩_k)`$ $`=`$ $`\left({\displaystyle \underset{i=1}{\overset{k}{}}}f(𝐩_i)\right){\displaystyle \frac{F_{Nk}\left({\displaystyle \underset{i=1}{\overset{k}{}}}𝐩_{𝐓}^{}{}_{i}{}^{}\right)}{F_N(0)}}`$ (A6)
$`=`$ $`\left({\displaystyle \underset{i=1}{\overset{k}{}}}f(𝐩_i)\right){\displaystyle \frac{N}{Nk}}\mathrm{exp}\left({\displaystyle \frac{({\displaystyle \underset{i=1}{\overset{k}{}}}𝐩_{𝐓}^{}{}_{i}{}^{})^2}{(Nk)𝐩_{𝐓}^{}{}_{}{}^{2}}}\right).`$ (A7)
For $`k=1`$, Eq.(A6) gives the corrected one-particle distribution

$$f_c(𝐩)=f(𝐩)\left(1+\frac{1}{N}\frac{𝐩_{𝐓}^{}{}_{}{}^{2}}{Np_T^2}\right)$$
(A8)
where we have expanded the correction to leading order in $`1/N`$. Similarly, one obtains for k=2

$$f_c(𝐩_1,𝐩_2)=f(𝐩_1)f(𝐩_2)\left(1+\frac{2}{N}\frac{(𝐩_{𝐓}^{}{}_{1}{}^{}+𝐩_{𝐓}^{}{}_{2}{}^{})^2}{Np_T^2}\right).$$
(A9)
The contribution of momentum conservation to the correlation function (3) is thus

$$C^{\mathrm{\Sigma }p_T}(𝐩_1,𝐩_2)=\frac{f_c(𝐩_1,𝐩_2)}{f_c(𝐩_1)f_c(𝐩_2)}1=\frac{2𝐩_{𝐓}^{}{}_{1}{}^{}𝐩_{𝐓}^{}{}_{2}{}^{}}{Np_T^2}$$
(A10)
Similar results were obtained in , Appendix B.

## B Estimate of $`p_T^2`$ at SPS

In a central Pb-Pb collision at SPS, about 680 negatively charged particles are emitted . If we assume that these are mostly $`\pi ^{}`$ and $`K^{}`$ with a ratio $`K^{}/\pi 0.09`$, where $`\pi (\pi ^++\pi ^{})/2`$ , there are some 625 $`\pi ^{}`$ and 55 $`K^{}`$. Taking into account the various kaon-to-$`\pi `$ ratios, as well as the expected isospin symmetry between pion species, this gives a total of 280 kaons and 1875 pions. With these particles come the 416 nucleons (we assume there is no fragment), giving about 2500 particles in all. We shall neglect the contributions to $`p_T^2`$ from other particles as e.g. antiprotons.

In order to calculate the average value of the squared transverse momentum $`p_T^2`$ for a particular particle species, one needs its distribution in transverse momentum and rapidity. The $`p_T`$\- and $`y`$-dependences may factorize, leaving a normalized $`p_T`$ distribution which is parametrized by

$$\frac{dN}{d^2𝐩_𝐓}=\frac{e^{m/T}}{2\pi T(m+T)}\mathrm{exp}\left(\frac{m_T}{T}\right).$$
(B1)
where $`m_T=\sqrt{m^2+p_T^2}`$ is the transverse mass, while the inverse slope parameter $`T`$ depends on the particle species and the energy of the collision. This parametrization can be used for the pions, nucleons and kaons at SPS, with $`T_\pi =180`$ MeV, $`T_N=300`$ MeV and $`T_K=215`$ MeV respectively . These numerical values fit transverse momentum spectra in the central rapidity region. Away from mid-rapidity, transverse momenta are somewhat smaller, so that we overestimate $`\mathrm{\Sigma }p_T^2`$. The spectrum (B1) yields an average $`p_T^2`$ given by

$$p_T^2=2mT\frac{1+3\frac{T}{m}+3\frac{T^2}{m^2}}{1+\frac{T}{m}}.$$
(B2)
Substituting in Eq.(B2) the values of the inverse slope parameter $`T`$ given above, and multiplying the resulting average squared transverse momentum by the multiplicity of the particle species, one obtains for a central collision at SPS.

$$\underset{j}{}p_T^2930\text{GeV}^2.$$
(B3)
As mentioned above, this is somewhat overestimated since we have neglected the decrease of inverse slopes away from mid-rapidity.

For a non-central collision, we assume that the sum of $`p_T^2`$ scales with the multiplicity, i.e. we neglect the centrality dependence of inverse slopes and particle ratios. Therefore, in the case of semi-central collisions used by NA49 for flow studies, with only 40 to 55% of the maximum multiplicity , the value (B3) must be divided by 2.


format: text

text:
# Mathematical Software: Past, Present, and FutureDedicated to John R. Rice on occasion of his 65th birthday. Contribution of the National Institute of Standards and Technology, not subject to copyright. Mention of commercial products in this paper does not imply recommendation or endorsement by NIST. Conversely, omission of a product’s name does not imply unsuitability for use. Author’s address: Mathematical and Computational Sciences Division, Information Technology Laboratory, NIST, Stop 8910, 100 Bureau Drive, Gaithersburg, MD 20899-8910, USA

## 1 Introduction

### 1.1 The study of mathematical software

The field of mathematical software is concerned with the science and engineering of solving mathematical problems with computers. The primary focus is the development of general-purpose software tools applicable to problems in a variety of disciplines. There are a large number of facets to this work, including the following.

* the development and analysis of algorithms for standard mathematical problems which occur in a wide variety of applications
* the practical implementation of mathematical algorithms on computing devices, including study of interactions with particular hardware and software systems
* the environment for the construction of mathematical software, such as computer arithmetic systems, languages, and related software development tools
* software design for mathematical computation systems, including user interfaces
* testing and evaluation of mathematical software, including methodologies, tools, testbeds, and studies of particular systems
* issues related to the dissemination and maintenance of mathematical software

In 1977 John Rice aptly characterized the need for specialized study in this area with the following observation .

> Many sophisticated scientists produce naive software just as many sophisticated computer programmers produce naive science.

Tremendous progress has been made in the mathematical software field in the past 25 years. Yet, there continues to be a wide range of quality in existing software, in both the research and commercial domains. Good mathematical software results from the application of certain principles, methodologies, and practices derived from both applied mathematics and computer science. The study of these principles and practices is central to the field of mathematical software. To this end, typical software engineering practices, while beneficial to the production of mathematical software systems, are not sufficient. Mathematical software operates in the milieu of scientific computing, which has a number of characteristics that distinguish it from other areas. Among these are the following.

* Floating-point arithmetic. Most scientific computations are performed with floating-point arithmetic. Consequently, rounding errors occur in most arithmetic operations. Mathematical algorithms, therefore, must not only be correct in a strict mathematical sense, but they must control the accumulation of round-off, avoid catastrophic loss of significance from the subtraction of like quantities, and avoid unnecessary overflows and underflows. Such problems are sometimes unavoidable; software systems must be designed so that they do not fail when these anomalies do occur.
* Approximations. Floating-point arithmetic certainly implies approximation at a very fine level. However, more substantial approximations occur in mathematical computation. Infinite series are truncated, difficult-to-compute functions are approximated by polynomials, derivatives are approximated by differences, integrals are approximated by finite sums, curved domains are approximated by polygonal ones. The combined effect of such approximations on the final result can be quite difficult to assess. Analysis must be used to show that the correct solution is obtained as the approximations are made more precise (i.e., that the algorithm is convergent). However, good software must do more. It must provide mechanisms for a user to assess the quality of the result, and to alert the user when the result is suspect. Well engineered software can use such metrics to automatically control the level of approximations, optimally adapting the algorithm to the situation at hand.
* Infinite processes. Many mathematical computations consist of applying some infinite process that obtains the desired result only in the limit. Such processes must be truncated for practical use. Considerable research efforts have been involved in finding iterations that converge quickly. Deciding when to stop is always a difficult problem of practical concern to software developers. Good software must employ techniques that detect divergence or too slow convergence and take appropriate action.

Coupled with these fundamental mathematical challenges are practical concerns about portability. How can developers produce software with reliable, reproducible behavior when it must run in very different environments, with different types of processor architectures, arithmetic systems, memory hierarchies, operating systems, and language processors? Such questions are critical in the study of mathematical software.

### 1.2 The contributions of John R. Rice

At this conference, we are celebrating John Rice’s long and influential research career. John has made fundamental contributions to the areas of approximation theory, numerical analysis, mathematical software, and computer science. In the area of mathematical software, his technical contributions have had three overriding themes.

1. Architecture of scientific software systems. John has participated in the design and development of a variety of widely distributed mathematical software systems . In the course of this work he pioneered a number of design concepts which have influenced many systems. Among these are polyalgorithms , meta-algorithms , and software parts .
2. Raising the level of abstraction. Software users are more efficient when they can express their computational needs in the language of their technical field. For applied mathematics, abstractions are based upon concepts of the calculus, not simple arithmetic operations encapsulated in programming languages like Fortran. John Rice’s work in high-level components and languages , intelligent interfaces , and problem-solving environments  have served to push abstractions to higher and higher levels.
3. Understanding software via experimentation. Understanding the behavior of software is necessary in order to make practical decisions regarding algorithm selection . John has often stressed the importance of the use of experimentation in such evaluations. The many small engineering decisions made in the course of translating an abstract algorithm into a working computer program can have an enormous impact on its performance characteristics. John has devoted much time to developing testing and evaluation methodology , and applying it to particular situations . Indeed, one of the principle applications of the ELLPACK system  and its successors has been to the performance evaluation of software for partial differential equations.

John’s contributions to the field of mathematical software have been voluminous and far-reaching. (In this paper I have only cited a few examples of his many writings on this subject.) In the remainder of this paper I will enumerate some of the major events in math software, pointing out some of John’s key contributions along the way. I will then describe several current issues facing the field and make several hazy predictions of the future.

## 2 Mathematical Software Past

### 2.1 Beginnings

The earliest applications of electronic computers were in science and engineering, for which mathematical computation played a central role. Programming was a very difficult chore, done without modern aids like high-level languages, compilers and debuggers. The first publication of a piece of mathematical software in a research journal probably occurred in 1949, when Mathematical Tables and Other Aids to Computation printed a UNIVAC code for the solution of Laplace’s equation written in machine language . Such codes were very difficult to produce, and the need for reuse of software was recognized very early on. In 1951, Wilkes, Wheeler and Gill presented one of the earliest program libraries, which was developed for the EDSAC111The EDSAC (Electronic Delay Storage Automatic Computer) was built in the late 1940s at the Mathematical Laboratory of the University of Cambridge. It was operational from 1949 until 1958. .

By the 1960s, the introduction of high-level programming languages, e.g. Algol and Fortran, had greatly eased the task of producing reusable mathematical software. The use of such languages was not without controversy, of course. Compiled code was not quite as efficient as hand-tuned assembly code, but most people were willing to accept this in light of the great savings in programmer time. Also, the subprogram structure provided by these languages provided a simple framework for the construction and maintenance of libraries of utilities.

In 1960, the Association for Computing Machinery (ACM) began a new editorial department in the Communications of the ACM (CACM) devoted to the publishing of algorithms. Edited by J. H. Wegstein of the National Bureau of Standards (NBS), this section printed the code of contributed Algol procedures (most such codes were quite short). Also, remarks on and certifications of previously published codes were solicited. The first such contribution was a code for numerical quadrature submitted by R. J. Herbold of NBS . Each algorithm was given a number, and the set of algorithms later became known as the Collected Algorithms of the ACM (CALGO).

Computer manufacturers also began to develop libraries for their users. The most prominent of these was probably the IBM Scientific Software Package (SSP). A number of laboratories, such as Bell Labs, Boeing, Harwell, and Monsanto, began the development of math software libraries for their internal use. Several organizations, such as SHARE, the IBM user’s group, began to collect such utilities for redistribution.

Of course, subroutine libraries were not the only focus of researchers in this new field. Some were imagining ways in which these new powerful computers could be used to transform the way in which applied mathematics was practiced. Many of these ideas were discussed at the Symposium on Interactive Systems for Experimental Applied Mathematics held in in Washington, D.C. in August 1967 . The vision there is remarkably clear; many of the participants reported on developments in technologies which would only finally begin to be realized in the 1980s and 1990s. At Purdue, for example, John Rice, Saul Rosen and colleagues designed NAPSS (Numerical Analysis Problem-Solving System), an interactive mathematical problem-solving system which would accept input akin to normal mathematical notation , and would employ a variety of heuristics to automate numerical analysis. Unfortunately, the resources necessary for such an ambitious system exceeded even the supercomputers of the day (like the CDC 6400), and a fully functional system was never realized.

### 2.2 A community emerges

Perhaps the first event that provided a real sense of community for researchers interested in the production and dissemination of reusable mathematical software was the Mathematical Software Symposium held at Purdue University in April 1970. John Rice organized the symposium222The organizing committee included Robert Ashenhurst, Charles Lawson, M. Stuart Lynn, and Joseph Traub, which was sponsored by ACM and the Office of Naval Research, and the proceedings were published as a book by Academic Press . Included in the proceedings were 23 papers, four descriptions of selected mathematical software, and more than 40 pages of introductory material prepared by Rice.

One of the recommendations from the Symposium was for the establishment of a journal that would publish papers related to mathematical software. John Rice vigorously pursued this possibility333A committee that included Wayne Cowell, Lloyd Fosdick, Tom Hull, M. Stuart Lynn, and Joseph Traub worked with him.. After considerable negotiations with ACM and the Society for Industrial and Applied Mathematics (SIAM), ACM agreed to publish the new journal. Papers from an NSF-sponsored conference were used to provide articles to seed the journal. Mathematical Software II was held at Purdue in May 1974. There were 225 attendees, with 82 papers presented. The best of those papers make up the majority of the first volume of the ACM Transactions on Mathematical Software (TOMS) which published its first issue in 1975 with John Rice as Editor-in-Chief. John continued in that position until 1993.

TOMS was chartered not only to publish traditional research papers, but also algorithms (with included code which would be refereed), certifications, translations, and remarks on previously published algorithms. The Algorithms section of CACM was moved to TOMS, and hence TOMS algorithms were numbered beginning at 493. One of the important features of the new journal was the establishment of a reliable Algorithms Distribution Service for CALGO. The distribution, on magnetic tape, was performed on a cost-recovery basis by IMSL, Inc. Obtaining software in machine-readable form was much more useful than reading code on paper. This also allowed TOMS to adopt the policy of not printing the code of algorithms in the pages of its journal, thus saving much in production costs.

A third conference organized by John Rice, Mathematical Software III, was held at University of Wisconsin in 1977 . By the end of the decade, mathematical software had emerged as a viable research area with an active community to support it. After publishing 25 volumes, TOMS remains a vibrant outlet for the work of this community .

### 2.3 Software emerges

Another important activity in the 1970s were the numerous efforts to develop carefully constructed, systematized collections of mathematical software. One of the first of these was the NATS project, the National Activity to Test Software, which was conceived in 1970. A joint venture of Argonne National Laboratory, Stanford University, and the University of Texas at Austin, NATS was designed to study problems in producing, certifying, distributing, and maintaining quality numerical software. A key part of this effort was the production of two Fortran software packages, EISPACK  for eigenvalue problems, and FUNPACK  for special functions.

EISPACK, which first appeared in 1972, was based upon algorithms published in the 1960s in Numerische Mathematik and later collected by Wilkinson and Reinsch in the Handbook for Automatic Computation . Although the core of EISPACK was largely a Fortran translation of these existing Algol codes, the project was enormously influential. It set a new standard for quality transportable mathematical software, rigorously tested in a wide variety of computing environments. Its success inspired the development of many systematized collections, or “PACKs”, in other areas: LINPACK for linear systems , FISHPAK for separable elliptic problems , DeBoor’s B-spline interpolation package , MINPACK for nonlinear systems , DEPAC for ordinary differential equations , Fullerton’s function library FNLIB , Swarztrauber’s FFTPACK for fast Fourier transforms , and QUADPACK for numerical quadrature .

Of all the early “PACKs”, LINPACK undoubtably saw the most widespread use. One of the keys to LINPACK’s success was the decision to base its coding on the newly proposed Basic Linear Algebra Subprograms (BLAS) . The BLAS performed elementary vector operations, such as norms, dot products, scaling, and vector sums. The innermost loops in LINPACK’s column-oriented algorithms occurred inside the BLAS. This allowed optimization of the whole package by simply optimizing the BLAS. This approach proved quite successful, and many machine-specific versions of the BLAS were developed and supported by computer manufacturers.

The 1970s also saw great advances in software for ordinary differential equations (ODEs). Gear’s code DIFSUB  provided a well designed framework for automatic integration of both stiff and non-stiff problems using linear multistep methods. Shampine and Gordon’s ODE  did the same for Runge-Kutta methods. Many subsequent packages were built using these basic designs. Other influential packages included COLSYS for two-point boundary-value problems  and DASSL444DASSL won the 1991 Wilkinson Prize for Numerical Software for differential algebraic systems .

A number of high quality multi-purpose libraries also had their start in the 1970s. In 1970, six British computing centers began an effort to develop a library for their ICL 1906A/S computers. The next year Mark 1 of the Nottingham Algorithm’s Group (NAG) library was released. Implementations for other systems followed, and in 1976 a not-for-profit company, Numerical Algorithms Group, Ltd., was formed to continue development and distribution. The NAG effort continues today . The first commercial math library effort was also begun in 1970 with the incorporation in Houston, Texas of IMSL, Inc. by Charles W. Johnson and Edward Battiste. By the time the company showed its first profit in 1976, there were 430 library subscribers; the IMSL library remains a viable commercial product . Bell Laboratories also developed a library, PORT, whose single-source approach to portability influenced many subsequent efforts .

The development of ELLPACK, a system for elliptic boundary-value problems, also began in the mid 1970s. This effort, which was led by John Rice, was a cooperative project of Purdue University, the University of Texas at Austin, Yale University, and others. In ELLPACK, the solution process was partitioned into a number of distinct phases (domain processing, discretization, indexing, linear system solution, and output), and the interfaces between these phases were carefully defined. This design allowed the development of a large library of components which could be easily composed to build algorithms for solving particular problems. ELLPACK also proved to be an excellent testbed for the evaluation of software for elliptic problems. To ease the use of the system, John Rice designed a very-high-level language to describe the problem to be solved, and to select the components to be used to solve it. The system first became fully operational in 1978 . Many of the basic concepts in ELLPACK’s design, such as high-level user interfaces and plug-and-play software parts technology, are in common use in modern problem-solving environments.

The development of mathematical software in the 1970s and early 1980s is described in detail in the book Sources and Development of Mathematical Software edited by Wayne Cowell .

### 2.4 Increased access

By the beginning of the 1980s a substantial collection of mathematical software, mostly in the form of Fortran subprograms, was available for use. The user base for this software had grown substantially, and with it came a new problem: how to locate that needed software component. The National Bureau of Standards (now NIST) developed an extensive catalog of such software. Their Guide to Available Mathematical Software (GAMS), based upon a detailed tree-structured problem-oriented classification system , allowed readers to see which components of which libraries and packages, both public domain and commercial, were available to solve each problem. The catalog remains available today as an online resource .

Another barrier to the widespread use of software developed by the research community was simply the process of obtaining the code. One had to locate the author, request a magnetic tape, and attempt to decipher its format. Useful software was often lost to the community when an author changed institutions and there was no longer support for distributing it. In 1985, Jack Dongarra, then at Argonne National Labs, and Eric Grosse at Bell Labs, started a software repository they called netlib , which pioneered the use of computer networks in software distribution. Software could be obtained automatically by return email after sending requests to an address whose email was processed by a Unix daemon. The ready availability of such software changed the way in which many researchers worked. Many more made routine use of high quality software, and many others were freed of the necessity of maintaining their own private repositories. Now Web-accessible and supported by the University of Tennessee at Knoxville and Bell Labs, with mirrors worldwide, netlib remains the premier repository of software developed by the mathematical software community .

The 1980s also saw the first commercial success for general-purpose interactive systems for mathematics. A system for matrix computations developed as a teaching aid during the period 1977-84 by Cleve Moler at the University of New Mexico, was commercialized as a tool for control system engineers. Today MATLAB is a very popular system for scientific computing . The overall structure of modern interactive mathematics systems was greatly influenced by the system Mathematica developed by Stephen Wolfram in 1988 . Mathematica was the first commercial system to integrate symbolic, numerical, and graphical capabilities into a single package. The growing availability of personal computers and workstations was an important factor in the success of these systems. With these tools, the use of mathematical software was beginning to expand to those with little experience in numerical methods or even programming.

### 2.5 New architectures

The 1980s also brought vector and parallel computers into widespread use, and with them additional challenges to the design of mathematical software. Vector processor vendors developed specialized math libraries tuned for their systems, mainly containing software for linear systems and FFTs. These solutions emerged because the performance of linear algebra software such as LINPACK was disappointing on vector register architectures like the Cray and Convex. The main reason for this was the fact that moving data from memory to vector registers was very costly, and that LINPACK’s column-oriented algorithms, based on the BLAS, necessitated more data movement than was really necessary.

In 1984 John Rice hosted a workshop at Purdue (“ParVec Workshop Number 4”) in which a variety of schemes for developing portable high-performance software for vector parallel systems were proposed . Jack Dongarra and Sven Hammarling proposed the development of new classes of BLAS: Level 2 BLAS for matrix-vector operations, and Level 3 BLAS for matrix-matrix operations. Encapsulating $`O(n^2)`$ and $`O(n^3)`$ operations, respectively, as fundamental operations would provide much more opportunity to optimize core operations on different processors. These new BLAS  would provide the basis for a major new linear algebra package released in 1992. LAPACK , which included the functionality of both EISPACK and LINPACK, used block-oriented algorithms in which the fundamental operations were now matrix operations encapsulated in the Level 3 BLAS. These have proven to be highly efficient on modern vector processors and symmetric multiprocessors. Every major computer manufacturer now supports tuned BLAS for their systems and incorporates LAPACK in their math library. Community efforts are currently underway to extend the BLAS in new directions, such as sparse matrix operations

In the late 1980s multiprocessor systems of widely differing design were becoming routinely available, and with them a host of new programming models, supported by specialized message-passing primitives. Developing portable software for the class of distributed memory (multiple instruction multiple data, or MIMD) systems became a new challenge. The PVM system developed in 1991 provided a useful abstraction for parallel programming and was very widely adopted . Its implementation on many parallel machines demonstrated the usefulness and feasibility of a common message-passing infrastructure. This led to a grass root message-passing standardization effort. The resulting Message Passing Interface (MPI) transformed the landscape for distributed parallel computing .

One of the first portable math software libraries for distributed architectures was ScaLAPACK, a distributed memory counterpart of LAPACK linear system solvers . This package became the core of several multi-purpose distributed memory math software libraries which first appeared in the 1990s. Among these are the NAG Parallel Library , IBM’s PESSL , and the European PINEAPL effort .

The increasing complexity of scientific software systems being developed in the 1990s led to an interest in new software architectures. Object-oriented approaches to the development of mathematical software began to be seriously considered. The notorious inefficiencies of pure object-oriented design, and the lack of language standardization made such pursuits difficult. Nevertheless, approaches that allowed many of the advantages of object-oriented design without sacrificing efficiency were developed. LAPACK++, a subset of the linear systems solvers in LAPACK written in C++, was one of the first such successful packages . Today object-oriented approaches are routinely used in scientific computing.

### 2.6 Expanding vision

By the 1990s, rapidly increasing computer power was leading to new visions for the future of mathematical software systems. During that period, for example, John Rice and colleagues led in the establishment of a new community of researchers interested in exploiting the promise of expert systems for numerical computing. In a series of conferences held at Purdue , the use of AI approaches for such tasks as algorithm selection, automatic programming, and process management were explored.

By this time, computation had become an essential ingredient in the practice of science and engineering. Interest in computational science as a new field of study was beginning, and interdisciplinary programs for training its practitioners were being established in many universities. John Rice and others began to develop a new vision for mathematical software systems to support computational science research. These systems, called problem-solving environments , would provide natural graphical user interfaces in which scientists describe their problems using the vocabulary of their native discipline. They would provide access to rich libraries of problem-solving components enabling Web-based parallel and distributed computation. Users would be able to interact with ongoing computations, to easily visualize results, to manage a large database of experimental results, and to ask advice of expert advisory systems. Many small-scale special-purpose systems now under development and use can be classified as problem-solving environments, and research groups throughout the world are working on infrastructure necessary for the routine construction and use of such systems. Work at Purdue on parallel ELLPACK , WebELLPACK , and PYTHIA  are serving to address issues in PSE design. Examples of current work in network-based scientific computing are Netsolve , the NEOS optimization server , and the computational grid .

The vision of scientific computation in heterogeneous distributed environments places stringent requirements on the portability and interoperation of scientific software that are extremely difficult to achieve . Such needs have sparked interest in the use of common virtual environments such as Java555Java is a trademark of Sun Microsystems. for computational science and engineering. The Java language and its environment (the Java Virtual Machine), which has become available on nearly every computing platform, provides a fixed floating-point model, threads, remote execution, standard GUIs, and other facilities within a simple object-oriented programming language. While these are the main facilities necessary for the construction of problem-solving environments, there has been some reluctance to adopt Java within the scientific community due to concerns about efficiency and the lack of several programming conveniences important to scientists and engineers . Community efforts such as the Java Grande Forum are seeking to improve this situation .

Virtual environments do not necessarily solve the problem of performance portability, since virtual machine instructions must be mapped on to local computer hardware for execution. Modern computing hardware is extremely complex, characterized by multiple processing units, vector pipes, register farms, several levels of cache (with increasing access times), local memory, remote memory, and disk storage. Getting the highest performance possible requires that the programmer take into account all the special properties of the system in use. This leads to extremely complex software even for the simplest of tasks. Matrix multiplication can turn into a 10,000-line polyalgorithm. Recent approaches have provided new hope for overcoming this software development nightmare. Clint Whaley and Jack Dongarra have recently developed a system, called ATLAS, for Automatically Tuned Linear Algebra Subprograms . ATLAS generates highly efficient BLAS for a given architecture using an experimental approach. By running many hundreds of tests, ATLAS determines the most efficient way to implement a given operation. The result is consistently on par with, and often exceeding, code which takes expert programmers weeks to develop. Matteo Frigo and Stephen Johnson have taken a similar approach in the computation of fast Fourier transform . For FFTs, hardware also interacts with the prime factorization of the sequence length $`n`$ to add further complication. FFTW, the Fastest Fourier Transform in the West, uses heuristics and experimentation to develop a just-in-time strategy for fast computation for a given $`n`$ on a given processor.666FFTW won the 1999 Wilkinson Prize for numerical software.

## 3  Mathematical Software Present

In this section I point out a variety of meta-issues that face mathematical software researchers today.

### 3.1 Mass-market software

Until recently mathematical software was produced mostly by experts in numerical analysis as a byproduct of their research in algorithms. Users of this software also were fairly sophisticated, with some experience in numerical algorithm development themselves. They had an appreciation of the limitations of numerical algorithms, and the necessity of careful verification of results, even when using software developed by experts.

Today’s community of mathematical software developers and users is much larger, and much more diverse. The great demand for mathematical computations has made mass-marketed mathematical software profitable. Commercially supported mathematical and statistical software is now widely available, with high-level interfaces that allow use by non-programmers. Such users often do not have the background necessary to recognize the difference between a difficult problem and a routine one. The mathematical landscape is still littered with pitfalls, and these users may be too trusting of the results produced by the scientific software systems that they use. Programmers who add mathematical and statistical capabilities to commercial software systems are no longer experts in numerical analysis. They may be content to code up a formula from a book without giving thought to its numerical properties. The problem may be even more severe in systems that are not overtly mathematical in nature. Mathematical computations are increasingly being done in embedded devices, coded by programmers whose mathematical sophistication may be suspect.

As a result, in spite of tremendous progress in numerical methods and software, many users of modern mathematical software are at risk. There is now a desperate need for numerical analysts to develop and apply methodologies for the validation of mathematical and statistical software. Techniques, tools, reference data, and reference software are needed to support critical evaluations of mathematical software by developers and users . Unfortunately, there is little interest and support within the research community for such activities.

Those software developers who seek advice regarding numerical software production are likely to look to popular sources like Numerical Recipes . Books like this provide a reasonably good introduction to numerical methods, and the programs they include provide good examples of the basic techniques. Programs like these are often incorporated wholesale into applications, in spite of the fact that they are typically less efficient, robust, and reliable than state-of-the-art mathematical software. The mathematical software community needs more popularizers who can bring the message of good numerical software design to those in other fields.

### 3.2 Tower of Babel

For many years there was one language for scientific computing: Fortran. This greatly simplified the development and reuse of mathematical software components. Today we are faced with a plethora of programming languages in use for scientific computing. Though officially obsolete, Fortran 77 is still the language of many. Good compilers are now available for Fortran 90, and many have been extended to support Fortran 95, the current Fortran standard, although their adoption by programmers has been slow in coming. The C language has proven much more popular, for which excellent compilers are now available. Most GAMS users who cannot find the software they seek are looking for C procedures. C++, the object-oriented extension to C, is the choice for a growing number of new mathematical modeling projects. Unfortunately, C++ has not, until recently, had an agreed-to standard, and, as a result, developing portable software has been difficult. Java, the popular network-aware object-oriented programming language developed by Sun, is being seriously considered by many, although its performance and language features leave much to be desired. The fact that Java is now being widely taught in universities insures its future. Finally, many software components are being developed in the very-high-level languages used in specialized systems; MATLAB is the primary example.

We are clearly facing a transition in computer languages for science and engineering computation. Numerical analysts no longer have much influence on the choice of language of those doing numerical computing. Language choices are more often made based on other considerations, such as the need for convenient and portable graphical user interfaces, visualization tools, and other critical system services. While such services are largely unavailable to Fortran programmers wishing to develop portable systems, they are conveniently at hand in C, C++, and Java. The increased portability afforded by the widespread availability of Java Virtual Machines on Windows, Unix, and Apple platforms, has made Java a very attractive option. While mixed-language programming is possible, and does provide the ability to reuse legacy Fortran software, this option is not popular among users. It adds complexity to the software project, while making the code more difficult to transport.

Unfortunately, there is almost no support for the migration of the existing base of Fortran mathematical software components to other languages. As a result, this well-engineered software is being increasingly bypassed in favor of inferior home-grown solutions.

### 3.3 The risks of self-publishing

The rise of the Internet has greatly eased the exchange of information among researchers. It is simple and convenient for research groups to develop a Web page to distribute software and documentation to potential users. While this has led to increased access to research software, it places the long-term maintenance of the output of the research community in jeopardy. Project Web pages on departmental servers are not permanent fixtures. Nevertheless, many researchers are using such facilities in place of submitting their software to more permanent archives such as netlib or the Collected Algorithms of the ACM. There is a danger that much of the currently available expertise embedded in such software will be lost to future researchers.

## 4 Mathematical Software Future

In this section, I offer a few predictions regarding future mathematical software research and usage.

Prediction: Within five years ACM will cease print publication of TOMS.

Subscriptions to TOMS have been dropping at the rate of about 5 percent per year for some time. Other ACM journals, and indeed most other mathematics and computer science journals, are experiencing the same phenomenon. One of the reasons for this is the proliferation of specialized journals, principally developed by commercial publishers.

For some time, ACM has been considering mechanisms for maintaining their publication program as a viable service to the community. The solution to this problem is found in the ACM Digital Library (ACMDL) which premiered in 1998 . The ACMDL provides its subscribers with online access to all ACM journal articles and conference proceedings published since 1985 at a subscription fee which is less than the cost of three printed journals. Currently this accounts for more than 350,000 pages of text. Acceptance of the ACMDL by members and subscribers has been overwhelming, providing ACM with the additonal revenue to begin the work of extending the ACMDL holdings to include all material published by ACM since its inception in 1947. At the same time, the success of the ACMDL has contributed to a further 25 percent drop in print subscriptions in 1999. If present trends continue, printed versions of ACM journals will be no longer be sustainable in five years time. Instead, they will be superceded by their electronic counterparts.

The ACMDL will serve to blur the distinctions between individual ACM journals. The concept of a journal will be replaced by that of an input stream to the ACMDL controlled by a certification authority, i.e., a board of editors supported by volunteer referees. In such an environment, it will be much easier (and much less financially risky) for ACM to initiate new refereed input streams, and to phase out those which have become less productive. Rather than subscribe to individual journals, ACMDL subscribers will have access to an individually tailored notification service which will alert them to the availability of new articles in their areas of interest.

Publications in the ACMDL will not be restricted to articles with a severe page limit. Extended appendices will be easily accommodated, as will other artifacts such as software, audio, video, etc.

Prediction: Users will no longer install mathematical software on their workstations.

The need for instantaneous distribution and use of mathematical software components in heterogenous network environments will put new pressures for software portability. A key element of the solution will be standardized virtual environments in which software can execute. Java is an example of such an environment. Its widespread availability also provides a new model for software distribution. Complex conglomerations of source code will no longer need to be explicitly downloaded and installed on the local systems in advance of their use. Instead pre-compiled bytecodes for the virtual machine will be able to be downloaded from sites of developers or vendors on demand. This also provides a solution to the problem of distributing patches and updates to software. Rather than purchasing an entire library, software users will have the option of subscribing to a service, paying only for the portions of the library that they actually use.

Another new model for software reuse in a network environment is based on a remote execution paradigm. In this case, problem-solving services are made available to users over the network. When a problem need be solved, a message containing a high-level specification of the problem is sent to the service provider, who provides both the software and the execution cycles needed to solve it. This model is probably more appropriate for access to large scale systems, like finite-element modeling packages.

Prediction: The percentage of people directly using math software libraries will decrease.

The wide availability of problem-solving environments (PSEs) for various domains will bring computational capabilities to an even wider audience than today. These users will make use of the services of the PSE, blissfully unaware of the complex system, involving software libraries, expert systems, and remote execution, which are being marshalled on their behalf.

However, if this vision is to be realized, a new class of software designers must be trained. They must be well-versed in numerical analysis, mathematical algorithms, modern software design, and high-performance computing and communications. Additional research in mathematical software must be performed to provide new methods for improving the robustness and adaptibility of mathematical software systems, and to address new problem areas. And finally, new methods for assessing the correctness and reliability of complex mathematical software systems must be devised and deployed.

Mathematical software is still a vital and vibrant research area that will increase in importance in the coming decades. We are grateful to John Rice for his vision and leadership in getting us here.