text:
# Multiplicity of different hadrons in 𝑒⁺𝑒⁻, pp, and AA collisions

## 1 Introduction

Studying hadron production in $`e^+e^{}`$annihilation and pp scattering has already a long tradition . In recent times particle production has gained renewed interest because it has been argued that in heavy ion collisions the enhanced yield of strange particles (as compared to pp) may serve as "smoking gun” for the conjectured creation of a quark gluon plasma (QGP). Looking carefully at the presently available data of the reaction 158 AGeV Pb \+ Pb one observes indeed this fact  .

The formation of particles with a given mass depends always on the energy which is available in the system which forms them. If this energy is not much higher than the masses of the final state particles one expects that the phase space of the strange particles is much smaller than the one of their lighter non-strange counterparts and therefore they are less frequently produced. This is expected to be the case in two body collisions in an expanding hadron gas with a temperature of T=170 MeV as considered as a good description of the final state of the heavy ion interaction. Contrary, if the energy of the forming system is much larger than the hadron masses, the mass difference should not play an important role and we expect to see (almost) the same yield. This scenario is expected if the hadrons are produced by a decaying QGP.

In this situation it seems to be worthwhile to apply state-of-the -art models in order to check the validity of the above mention conjecture. Here we use the recently developed NEXUS model, whose results have been extensively compared with most of the existing $`e^+e^{}`$, $`pp`$, $`\mu A`$, $`pA`$ and $`AA`$ data , . For our study we compare the results of simulations for 4 different reactions: $`e^+e^{}`$at $`\sqrt{s}`$=91 GeV, $`pp`$ at $`\sqrt{s}`$=17 GeV, central Pb-Pb collisions at $`\sqrt{s}`$=17 AGeV (SPS) and central Au-Au collisions at $`\sqrt{s}`$=200 AGeV (RHIC).

## 2 NEXUS

How can one simulate in a realistic way hadron-hadron or even nucleus-nucleus collisions? A sophisticated approach to such high energy interactions is the so-called Gribov-Regge theory. This is an effective theory, which allows multiple interactions to happen “in parallel”, with phenomenological object called “Pomerons” representing elementary interactions. Using the general rules of field theory, on may express cross sections in terms of a couple of parameters characterizing the Pomeron. Interference terms are crucial, they assure the unitarity of the theory. A big disadvantage of this approach is the fact that cross sections and particle production are not calculated consistently: the fact that energy needs to be shared between many Pomerons in case of multiple scattering is well taken into account when calculating particle production (in particular in Monte Carlo applications), but not for cross section calculations. Another problem is the fact that at high energies one needs a consistent approach to include both soft and hard processes. The latter ones are usually treated in the parton model approach, which only allows to calculate inclusive cross sections.

We recently presented a completely new approach  for hadronic interactions and the initial stage of nuclear collisions, which is able to solve several of the above-mentioned problems. We provide a rigorous treatment of the multiple scattering aspect, such that questions as energy conservation are clearly determined by the rules of field theory, both for cross section and particle production calculations. In both (!) cases, energy is properly shared between the different interactions happening in parallel. This is the most important and new aspect of our approach, which we consider a first necessary step to construct a consistent model for high energy nuclear scattering. The elementary interaction is the sum of the usual soft Pomeron and the so-called semi-hard Pomeron, where the latter one may be obtained from perturbative QCD calculations (parton ladders). To some extend, our approach provides a link between the Gribov-Regge approach and the parton model, we call it “Parton-based Gribov-Regge Theory”. Although the main purpose of neXus is the simulation of nuclear collisions, also electron-positron annihilation and deep inelastic lepton-nucleon scattering have been intensively studied to make sure to have a consistent approach of high energy scattering.

The neXus simulation program is therefore available to simulate lepton-lepton, lepton-hadron, hadron-hadron, hadron-nucleus, and finally nucleus-nucleus collision at high beam energy.

In electron-positron annihilation, starting from the virtual $`Z`$, a time-like cascade of parton develops, which are used to define a kinky string via identifying partons and kinks. Particle production is the calculated using the standard covariant string fragmentation procedures: the string breaks into string segments vis quark-pair production, the segment representing the final hadrons and resonances.

In pp scattering, one has to deal with a complicated multiple scattering process, but at the end the situation is quite similar to electron-positron annihilation in the sense that parton configurations are translated into string language in order to “parameterize” the theoretically unknown hadronization.

Concerning nucleus-nucleus collisions, we have to distinguish primary and secondary interactions. Primary interactions are a straightforward generalization of proton-proton scattering: multiple scattering – parton production – string formation. Whereas in pp only few strings are formed and one may without problem consider independent string fragmentation, this is not realistic for nuclear collisions. This fact is take into account by so-called secondary interactions, going beyond neXus. To get some idea what happens we employ a simple droplet model: based on neXus, we calculate energy densities at a given proper time (say at 1 fm/c). Spatial region with an energy density bigger than some critical value are considered to “quark matter” and we treat the further development macroscopically, assuming Bjorken-like longitudinal expansion and some transverse expansion. Finally, when the energy density drops below the critical density, the droplet decays into hadrons, purely via phase space. This approach proved to provide a very satisfactory description of the CERN heavy ion reactions .

## 3 Results

Before embarking to detailed comparisons and conclusions, we have to make sure that the model describes reasonable well the observables we use here. Here we are interested in the multiplicity in 4$`\pi `$ and therefore it is useful to compare the Nexus results with experiments as well in 4$`\pi `$ .

In fig. 1, we compare the 4$`\pi `$ Nexus results for particle multiplicities (after decay of all not identified resonances) for central events with a compilation of the experimental data. The $`\mathrm{\Omega }`$ yield was obtained by assuming that the ratio of $`\mathrm{\Xi }`$ to $`\mathrm{\Omega }`$ at mid-rapidity obtained by WA97 is the same as the one for $`4\pi `$. On the right hand side, we display mesons and baryons, on the left hand side mesons and anti-baryons. The Nexus results are in good agreement with the extrapolated experimental data over four orders of magnitude. In view of a similar good agreement for $`e^+e^{}`$and pp  it is useful to use Nexus for the question at hand. In fig. 1, we display at the bottom our predictions for central collisions of Au + Au at $`\sqrt{s}`$ = 200 AGeV (here as well non identified resonances have decayed).

The important difference of pp and AA as compared to $`e^+e^{}`$reactions is the distribution of the energy W of the strings, i.e. the energy which is available for the production of particles. For the 4 different reactions it is displayed in fig.2. In the $`e^+e^{}`$reactions at 91 GeV gluon bremsstrahlung rarely produces a quark-anti-quark pair and hence the decaying string has almost always the full energy of 91 GeV. In pp and AA reactions, on the contrary, only a small fraction of the initial energy of the protons is available in the strings. The distribution of the string energies for the reactions at $`\sqrt{s}`$ = 17 GeV has a peak at around 1GeV and falls steeply off towards the maximal value of 17 GeV. The string contains only a small fraction $`W=\sqrt{x_1x_2s}`$ of the available energy $`\sqrt{s}`$, due to the fact that the light cone momentum fractions $`x_1`$ and $`x_2`$ of the partons representing the string ends are effectively given by structure functions which peak at low x. This peak value is small as compared to the threshold for multi strange baryon production (due to strangeness conservation the threshold for $`\mathrm{\Omega }`$ is around 3 GeV).

For RHIC energies the average string energy is already very large as compared to the mass of stable baryons but still smaller than that for $`e^+e^{}`$.

Due to these distributions of string energies, we expect to see differences in the yield of heavy as compared to light hadrons in the different reactions. In fig.3, we display the yields for the different reactions divided by $`m^{(3/2)}(2J+1)`$. In the case of $`AA`$ collisions, only those detected particles are considered which have been created directly from string decay and not via droplets.

In any case, we observe a roughly exponential drop of the yields as a function of the mass: it is always more difficult to produce high mass particles compared to low mass ones. The same exponential form of the mass yield is obtained in this representation in canonical models which describe the mass yield in terms of a temperature and chemical potentials . Looking more closely we observe, however, different slopes of theses exponential functions: $`pp`$ (17 GeV) and $`AA`$ (17 GeV) show a similar behavior and both drop faster than $`e^+e^{}`$(91 GeV) and $`AA`$ (200 GeV). In order to see the effects more clearly, we divide all distributions by the $`e^+e^{}`$results, see fig. 4.

The $`e^+e^{}`$(91 GeV) – by definition – and $`AA`$ (200 GeV) show a horizontal spectrum, whereas $`pp`$ (17 GeV) and $`AA`$ (17 GeV) show both a suppression of heavy particles as compared to $`e^+e^{}`$. At 17 GeV the string masses are small and hence phase space reduces the production of heavy particles, whereas for $`e^+e^{}`$(91 GeV) and $`AA`$ (200 GeV) there is no such limitation.

So far we considered only $`AA`$ collisions with hadrons coming directly from strings. The majority of hadrons which are observed finally is, however, produced from droplets. These droplets decay according to the microcanonical phase space into baryons (the spin 3/2 decuplet and the spin 1/2 octet) and mesons (pseudoscalar and vector octet) respecting all the conserved quantum numbers of the droplet. Their energy is usually large as compared to the hadron masses. In fig. 5

we combine the hadrons produced by droplets and string decays and show the full $`AA`$ results as seen by the detectors, again compared to $`pp`$ and $`e^+e^{}`$, and in all cases divided by the $`e^+e^{}`$results. We observe that all the nuclear results ($`AA`$ (17 GeV) and $`AA`$ (200 GeV)) are almost parallel, whereas the $`pp`$ spectrum drops much faster. If we take $`pp`$ as a reference, we might conclude: we observe indeed an enhancement of heavy (multi-strange) hadrons for nuclear collisions but also for $`e^+e^{}`$, however, in view of the above discussion it is more appropriate to talk about a strangeness suppression due to the limited energy in pp.

This is a very astonishing result in view of the fact that the origin of the hadron production is very different in $`e^+e^{}`$as compared to the hadronic interactions. As said before, the decay of a string in $`e^+e^{}`$is determined by the area law which is close to a sequential decay in a one dimensional longitudinal phase space. The droplet, on the contrary, disintegrates instantaneously and the different decay channels are populated according to the three dimensional phase space. Is it a accident that both mechanism give about the same hadron multiplicity distribution? The answer is clearly yes as can be seen from fig. 6.

The droplet decay pattern depends on the energy density $`ϵ`$ where the decay takes place. Smaller energy density implies larger volumes and therefore decay channels in which many particles are produced are favored due to the additional volume factor for each additional particle. If more particles are produced the heavier ones are more suppressed, as can be seen. Now it happens to happen that at energy densities of $`ϵ_{^{}}1GeV`$, a number which can be derived on theoretical grounds as well as by comparison with the data, $`e^+e^{}`$and hadronic interactions have about the same multiplicity distribution. This has been already observed in thermal models but never understood in physical terms. In the calculation presented here we have used $`ϵ=0.15GeV`$.

Comparing our completely microcanonical results with that of the canonical models , we see that both give about the same mass yield distribution. The strangeness suppression factor which is used in  is a simple parameterization for the fact that strange particles are heavier than their non strange counterparts. It depends on the mass difference and hence is not identical for all hadrons. This is in agreement with data.

In view of these results we can only draw the conclusion that at high energies (i.e. if the average string energy is considerably larger than the masses of the hadrons considered, as at RHIC energies) the functional form of the multiplicity as a function of the hadron mass cannot serve as a "smoking gun" for the plasma. A $`e^+e^{}`$reaction, in which certainly no plasma but a string is created, produces the same functional form. Hence the multiplicity distribution cannot tell whether the strings, which are created in the AA reactions as well, have interacted to form a QGP or whether they will disintegrate independently without having it formed. Thus the quest for the QGP is more challenging than thought.

Acknowledgments: We thank Johanna Stachel for fruitful discussions which have initiated this work.


format: text

text:
# Long-Term Variability of the Low-Luminosity Active Galactic Nucleus of M81

## 1 INTRODUCTION

A number of low-luminosity AGNs (LLAGNs), with X-ray luminosities in the range $`10^{4042}`$ erg s-1, have been discovered in nearby normal galaxies with ASCA (e.g. Makishima et al. 1994; Iyomoto et al. 1996; Ishisaki et al. 1996; Iyomoto et al. 1998; Terashima et al. 1998; Ptak et al. 1999). LLAGNs are now considered to be ubiquitous at least in our neighbourhood.

The low luminosities of LLAGNs imply lower black-hole masses and/or lower mass accretion rates as compared to more luminous AGNs. These alternatives may not be distinguished by X-ray spectra, since AGNs exhibit similar X-ray spectra over wide range of luminosities, from LLAGNs through quasars. On the other hand, X-ray variability is expected to provide useful clues to the issue, as a zero-th order approximation, because characteristic time scales of the X-ray emission is considered to correlate positively with the system size and hence the system mass.

During individual observations with ASCA, which typically last for one day, LLAGNs tend to be much less variable than the lower-luminosity end of the Seyfert galaxies (e.g. Ptak et al. 1998). We are therefore tempted to speculate that LLAGNs have longer characteristic time scales of the X-ray variation and hence larger black-hole masses compared to Seyferts. In this case, LLAGNs might be “dead quasars” in our neighbourhood. To examine this possibility, it is important to determine the power spectral density (PSD) of intensity variations of LLAGNs on time scales such as dozens of days or longer. So far, no such attempts have been made, except for a rough estimation by Ishisaki et al. (1996), because only few observations with uneven intervals are usually available for each LLAGN in X-rays.

Using ASCA, we investigate long-term variability of an LLAGN at the nucleus of M81, which is an Sab galaxy at a distance of 3.6 Mpc (Freedman et al. 1994). Its nucleus exhibits weak AGN properties in various wave bands. Optically, it is classified as a LINER (Ho et al. 1997). In radio, the spectrum from the nuclear region is interpreted as a jet emission (Boettcher et al. 1997). The nucleus is also an X-ray source with a 2–10 keV luminosity of $`2\times 10^{40}`$ erg s-1, as confirmed repeatedly by Einstein Observatory (Fabbiano et al. 1988), EXOSAT (Barr et al. 1995), Ginga (Tsuru 1992), ROSAT (Boller et al. 1992), BBXRT (Petre et al. 1993), ASCA (Ishisaki et al. 1996), and BeppoSAX (Pellegrini et al. 2000). These multi-mission results indicate a mild X-ray intensity variation. Ishisaki et al. (1996), using multiple observations with ASCA spanning $``$2.5 years, detected long-term variation as well as a short-term ($`1`$ day) variation on one occasion. Pellegrini et al. (2000) detected a periodic oscillation with a sub-day period. In addition, the spectra of this LLAGN, except that measured with Einstein, are well fitted with a power-law model having a photon index of 1.5–2.0. Therefore, M81 has been considered to possess an LLAGN.

## 2 OBSERVATIONS AND DATA REDUCTION

Owing to the explosion of the supernova SN 1993J (Kohmura et al. 1994) in M81, this particular field has been observed with ASCA 16 times over 5.5 years, as summarized in the observation log of Table 1. It is among the objects most frequently observed with ASCA. The results on SN 1993J have been reported by Kohmura (1994), Kohmura et al. (1994), and Uno (1997). The results on the M81 nucleus of the 1st to 10th observations have been published as Ishisaki et al. (1996).

The 16 observations were carried out with the GIS (Ohashi et al. 1996; Makishima et al. 1996) in the PH mode, while the SIS (Burke et al. 1991) in several different modes as described in Table 1. The M81 nucleus typically gave an intensity of 0.3 c s-1 for the GIS (0.7–10 keV, for each detector) and 0.4 c s-1 for the SIS (0.5–10 keV, for each detector). According to Kohmura et al. (1994), the SIS count rate of SN 1993J was $``$ 25% of that of the nucleus in the 1st observation, and decreased as $`t^{0.7}`$, where $`t`$ is the time after the explosion.

For the first six observations when SN 1993J was bright, we accumulated photons $`1^{}.5`$ around the nucleus to avoid the contamination from SN 1993J. Since SN 1993J exhibit hard spectrum in this epoch (Kohmura et al. 1994), we excluded photons above 4 keV in the following analysis. For the other 10 observations, we accumulated photons $`3^{}`$ around the nucleus, and considered all photons in the analysis. The SN 1993J contribution to the data accumulation region is expected to be $`<5`$%.

For the spectral analysis, we accumulated background spectra using the blank-sky data provided by NASA. The background photons are $`<2`$% of the signal photons. For the timing analysis, we made background light curve using the blank regions of the M81 field.

## 3 SPECTRA AND LIGHT CURVES

Figure 1 shows long-term light curves of the M81 nucleus, where the 2–10 keV flux in each observation is obtained by spectral fitting, in the following manner. The spectra of the M81 nucleus are dominated by a featureless hard continuum without significant absorption or sharp line feature, as Ishisaki et al. (1996) reported, so that we adopted single power-law model. We fixed photon index $`\mathrm{\Gamma }`$ at 1.8 and determined the fluxes in Figure 1, since $`\mathrm{\Gamma }`$ has remained almost constant around 1.8 over 5.5 years. The GIS and SIS give consistent fluxes within systematic errors. The results on the 5th through 10th observations agree with those reported by Ishisaki et al. (1996). Through the 1st to 6th observations conducted during the first 1.5 months, the intensity remained almost constant. On the other hand, over the other 10 observations performed every six months, peak-to-peak intensity variation amounts to a factor of three.

We also analyzed short-term variability of the source intensity as shown in Figure 2, by subdividing each observation into time bins of 1/8-day width. With the 1/8-day bin, typical error due to photon counting statistics becomes $`6\times 10^{13}`$ erg cm-2 s-1, or $``$3% of the average flux.

We summed up the GIS and SIS counts, since the two instruments give consistent results. Considering the integration radius and positional response of individual observations and the spectral information, we converted the count rates into 2–10 keV fluxes. Before the conversion, we subtracted the background and the contamination from SN 1993J, using the average value of each observation. We fitted the light curves assuming a constant count rate, to evaluate the intensity variation through the $`\chi ^2`$ test. The obtained values of $`\chi ^2`$ and degree-of-freedom (d.o.f.) are shown in the left bottom of the figures. Thus, the source exhibits statistically significant short-term variation in the 5th, 6th, 8th, 11th, 12th, 14th and 15th observations, in addition to the 3rd (Serlemitsos et al. 1996) and 9th (Ishisaki et al. 1996) observations, as already reported.

## 4 TIMING ANALYSIS

With the 1/8-day bin width, there are 112 light-curve bins in total, of which the average flux becomes $`2\times 10^{11}`$ erg cm-2 s-1. Root-mean-square (RMS) of the variation becomes $`7\times 10^{12}`$ erg cm-2 s-1, an order of magnitude larger than Poisson error. As already attempted by Ishisaki et al. (1996), we may estimate a gross continuum shape of the M81 PSD using the 112-point light curve, because it includes various sampling intervals: hours, days, months and years.

However, the light curve is subject to an extremely uneven and sparse sampling. To illustrate it, if M81 had been observed continuously during the total observed span (2022 days), we would obtain $`2022\times 8`$=16176 bins with the 1/8-day bin width. In practice, there are only 112 bins. More than 99% of the total span are thus data gaps. Consequently, it is difficult to determine the M81 PSD in general manner, i.e. Fourier transforming the light curve into “periodogram”, and obtaining a PSD as modulus squared of the periodogram.

We have introduced two key methods to estimate the long-term PSD from the sparse light curve (Iyomoto 1999; Iyomoto et al. 1999). The first is “forward method”. Namely, assuming various PSDs, we generate fake light curves and compare them with the actual M81 data to estimate which model PSD can best reproduce the actual M81 data. We utilized Monte-Carlo technique to evaluate the effects of the data sampling gaps. The second is the use of ”structure function”, hereafter SF, which is mathematically equivalent to PSD, but is less affected by data gaps. We convert the light curves into SFs, before comparing the actual data and the simulation, because it is difficult to compare light curves directly.

In §4.1, we calculate the M81 SF. In §4.2, assuming various PSDs, we generate fake SFs via Monte-Carlo technique. In §4.3, we compare the fake and actual SFs to determine break frequency, $`f_\mathrm{b}`$, and slope index, $`\alpha `$, of PSD. In §4.4, we determine the normalization of PSD, $`P_0`$, directly from the RMS of the M81 light curve, and derive confidence regions of the M81 PSD.

### 4.1 Structure Function

The definition of first-order structure function is given by Simonetti et al. (1985) as

$$S(\tau _k)=\frac{1}{N(\tau _k)}\underset{m=0}{\overset{N1n}{}}\left(l_ml_{n+m}\right)^2,$$
(1)
where $`l_n`$ is luminosity at time $`t_n`$. The sum is taken over data pairs satisfying $`t_{n+m}t_m`$=$`\tau _k`$, with $`\tau _k`$ ($`k`$=1,2,…) being the time-lag, and $`N(\tau _k)`$ denotes the number of such pairs. For a stationary random process with equally spaced sampling, SF, auto-correlation-function (ACF), and PSD are equivalent to one another. For example, SF is related to ACF, $`A(\tau _k)`$, as

$$S(\tau _k)=2\sigma ^2[1A(\tau _k)],$$
(2)
where $`\sigma ^2`$ is variance of the process.

Nevertheless, SF is much less sensitive to data gaps compared to ACF and periodogram, and hence to PSD (e.g. Hughes et al. 1992). In astronomy, SF have been used in various research fields including radio (e.g. Hughes et al. 1992), optical (e.g. Paltani et al. 1997, Kawaguchi et al. 1998), and X-ray (e.g. Leighly et al. 1997; Heidt et al. 1998) frequencies.

Figure 3a shows the M81 SFs calculated by eq. (1). In addition to the SF derived from the 1/8-day bin light curve, we superposed small-time-lag portion of SFs derived from 1/4, 1/16, 1/32 and 1/64-day bin light curves, for comparison. The SFs with shorter light-curve bin, especially those with 1/32 and 1/64-day binning, are strongly contaminated by the Poisson noise. On the other hand, it is reasonable to suppose that the Poisson noise becomes negligible for the SFs with light-curve binning longer than 1/8 day, because the 1/4 and 1/8-day bin SFs are similar to one another. It is consistent with the fact that the Poisson noise is $``$1/10 of the RMS amplitude of the M81 variation in the 1/8-day bin light curve. Accordingly, we utilize the 1/8-day binning in this paper.

Thanks to the various sampling intervals involved in our data, we can calculate the SF over a wide frequency range. SFs with time-lags shorter than one day come from within individual observations. SFs with time-lags between two days and 50 days are owing to various combinations among the 1st to the 6th observations. SFs with time-lags longer than six months are provided by combinations of observations which are separated by multiples of six month.

For time-lags shorter than one day, the M81 SF is smooth. Relatively large values of $`N(\tau _k)`$ may account for the smoothness. On the contrary, in the larger time-lag region, the calculated SF exhibits a large scatter, due to small $`N(\tau _k)`$. Therefore, hereafter, we average the SF into appropriate bins in time-lag, as shown in Figure 3b. This is justifiable because we are interested in the PSD continuum, rather than its fine structure.

The binned M81 SF seems to keep increasing towards the longer end of the time coverage. In other words, it is likely that the characteristic time scale of the M81 nucleus is longer than the total observed span of 2022 days. However, to make this inferences more secure, we need a more quantitative evaluation. Even the binned SF is subject to strong sampling bias, which depends on the sampling window. Consequently, it is difficult to directly convert the M81 SF into a PSD. Similarly, it is difficult to directly compare SFs suffering from different sampling windows. Instead, we here take a forward-method approach to find out a set of PSDs that are consistent with the observed SF. That is, we start from an assumed PSD, produce a number of Monte-Carlo light curves based on the assumed PSD and employing the actual sampling window, and convert each fake light curve into a fake SF. Then, we can study how the ensemble of SFs derived from a single PSD behaves statistically under the particular sampling window, and whether the actual M81 SF can be considered as a member of this ensemble or not. This is carried out in the following subsections.

### 4.2 Monte-Carlo Simulation

In order to qualitatively study the characteristic time scale of the variation and to estimate its lower limit, we hypothesized that the M81 PSD is expressed by a convex broken power-law as

$$P(f)=\{\begin{array}{cc}P_0\hfill & (f<f_\mathrm{b})\hfill \\ P_0(f/f_\mathrm{b})^\alpha \hfill & (f>f_\mathrm{b})\hfill \end{array},$$
(3)
as a start point of the forward method. Here $`f`$ denotes the frequency, $`P_0`$ is a constant, $`f_\mathrm{b}`$ is a characteristic frequency called “break frequency”, and $`\alpha `$ is a slope index. It is an empirical PSD of AGNs: white at low frequency and red at high frequency (e.g. Pounds & McHardy 1988; Edelson & Nandra 1998; Hayashida et al. 1998). In this paper, we use this formula as a template to describe our result on M81, and try to estimate the values of $`f_\mathrm{b}`$, $`\alpha `$ and $`P_0`$. We examine in § 4.3 whether this hypothesis is appropriate or not. We regard $`f_\mathrm{b}`$ and $`\alpha `$ as free parameters of the forward method. The value of $`P_0`$ is adjusted in § 4.4. The trial values of 1/$`f_\mathrm{b}`$ are $`10^{1.0}`$, $`10^{1.2}`$, $`10^{1.4}`$, …, $`10^{4.8}`$ days, and those of $`\alpha `$ are $`1.10`$, $`1.15`$, $`1.20`$, …, $`2.10`$. We consider Fourier periods from 1/$`f_1=65536`$ days to 1/$`f_2=1/4`$ day, so that there are 218(=262144) bins in the PSD. We adopted the 1/$`f_1`$ that is $`32`$ times longer than the total observed span of the actual data, considering that 1/$`f_1`$ of the actual PSD is infinite, and that the Fourier components longer than the total observed span also contribute to the observed light curve.

Assuming a particular set of $`f_\mathrm{b}`$ and $`\alpha `$, we generate fake SFs in the following steps.
(i) Using Monte-Carlo simulation, we generate a set, namely 218(=262144), of random numbers, uniformly distributed between 0 and 2$`\pi `$, and use them for the phase of each Fourier component of the assumed PSD.
(ii) By a Fourier transform of this Monte-Carlo generated periodogram, we obtain a light curve of 219(=524288) bins, with 1/8-day bin width and 65536-days total span.
(iii) From the 524288-point light curve, we select 112 data points by applying the M81 sampling window.
(iv) We normalize the 112-point light curve to have an RMS equal to that of the M81 light curve.
(v) We calculate SF from the 112-point light curve.
If we use another set of random numbers and go through the above steps, we can generate a different light curve from the same PSD. In this way, we generate 1000 fake SFs for each set of $`f_\mathrm{b}`$ and $`\alpha `$.

Figure 3c provides examples of fake SFs that were produced from a single PSD having $`1/f_\mathrm{b}=10^{4.0}`$ day and $`\alpha =1.4`$; the SFs with $`1/f_\mathrm{b}=10^{4.0}`$ day represents those having $`1/f_\mathrm{b}`$ enough longer than the total observed span. On large-time-lag range, the fake SF scatters as largely as the M81 SF due to the sparse sampling. Even on the time-lag range smaller than one day, there are significant scatters between individual fake SFs, arising from the sparse sampling, while each SF is smooth in such time-lag range. These are explained by mutual dependences among the SFs at different time-lags. Namely, the smoothness in the small time-lag range is not only because $`N(\tau _k)`$ is large, but because there are mutual dependences among the SFs at different time-lags.

Using ensemble properties of the 1000 fake SFs, we can reduce the effect of the sparse sampling. The ensemble average becomes smooth, as examplified in Figure 3d. Moreover, using standard deviation of these fake SFs, we can estimate how much the SF scatter under the sparse sampling. In Figure 3d, we show the standard deviation as error bars. Thus, we compare each ensemble-averaged fake SF with the actual SF considering the standard deviation, to examine whether the particular PSD assumed is appropriate or not.

It is worth pointing out that the 1/$`f_1`$ of the assumed PSD really have an effect on the fake SFs. Figure 3e shows the ensemble averages of fake SFs with 1/$`f_1`$ of 211, 212, …, 218 days. The fake SFs having 1/$`f_1>2^{16}`$ days well reproduce “drops” around one day in the actual M81 SF. We confirmed that similar dependence on 1/$`f_1`$ is seen in the fake SFs with shorter characteristic time scale, 1/$`f_\mathrm{b}`$, such as 10 days, 100 days and 1000 days. These effect of 1/$`f_1`$ depends on the sampling window. Therefore, it is important to select 1/$`f_1`$ that is enough longer than the total observed span of the actual light curve. In this paper, we adopted 1/$`f_1`$ of 216 day.

### 4.3 Comparison between Fake and Actual Data

In order to quantitatively compare the actual and fake SFs, we utilized the standard $`\chi ^2`$ evaluation technique. We defined $`\chi ^2`$ as

$$\chi ^2=\underset{\tau _k}{}\left\{\frac{S_{\tau _k}^{\mathrm{M81}}\overline{S_{\tau _k}^{\mathrm{fake}}}}{\sigma _{\tau _k}^{\mathrm{fake}}}\right\}^2,$$
(4)
where $`S_{\tau _k}^{\mathrm{M81}}`$ is the M81 SF at a time-lag $`\tau _k`$, while $`\overline{S_{\tau _k}^{\mathrm{fake}}}`$ and $`\sigma _{\tau _k}^{\mathrm{fake}}`$ are the ensemble average and standard deviation among the 1000 fake SFs. In this paper, we use the term $`\sigma _{\tau _k}^{\mathrm{fake}}`$ to refer to the data scattering due to the sampling window, and not to the Poisson statistical nature of the data. We regard this $`\chi ^2`$ as a measure of relative goodness of the assumed PSD. We then calculated $`\chi ^2`$ against various $`f_\mathrm{b}`$ and $`\alpha `$ and obtained the minimum $`\chi ^2`$, namely $`\chi _{\mathrm{min}}^2`$ of 25, at $`1/f_\mathrm{b}=10^{4.6}`$ days and $`\alpha =1.40`$. We can regard these particular set of parameters as the best estimate for $`f_\mathrm{b}`$ and $`\alpha `$ of the M81 PSD, as long as the PSD of eq. (3) is assumed. We should notice that, however, we give weight to the lower-limit to 1/$`f_\mathrm{b}`$ rather than the best-fit value, because the obtained best-fit value is much longer than the total observed span, and in such 1/$`f_\mathrm{b}`$ range, the $`\chi _{\mathrm{min}}^2`$ values at different 1/$`f_\mathrm{b}`$ are almost similar to one another, as seen in Figure 5a.

The next step is to estimate 68.3%, 90% and 99% confidence regions. However, we cannot apply general relation between confidence levels and the $`\mathrm{\Delta }\chi ^2`$/d.o.f. values, such as $`\mathrm{\Delta }\chi ^2`$/d.o.f. = 2.71 for the single parameter 90% confidence limits. It is because the d.o.f. is ambiguous in the present case, due to the mutual dependences among the SFs at different time-lags. This effect is particularly prominent in the small-time-lag range, as mentioned in $`\mathrm{\S }`$ 4.2.

Instead, we estimate the $`\mathrm{\Delta }\chi ^2`$ values that correspond to 68.3%, 90% and 99% confidence boundaries, using again the Monte-Carlo simulation. Employing the M81 best parameters obtained above, we generated another set of 5000 fake light curves and hence fake SFs; we would like to emphasize that the $`\mathrm{\Delta }\chi ^2`$ distribution, which is given below, is almost independent of the assumed PSD parameters. Of course, the random numbers of these fake light curves are different from those used in $`\mathrm{\S }`$ 4.2. For each of the newly produced fake SF, we determined the optimum $`f_\mathrm{b}`$ and $`\alpha `$, in exactly the same way as performed for the M81 data. We also calculated $`\chi ^2`$ associated with each simulation at the assumed PSD parameters themselves: $`1/f_\mathrm{b}=10^{4.6}`$ days and $`\alpha =1.40`$. We denote it $`\chi _{\mathrm{true}}^2`$. Note that the optimum parameters do not always coincide with the assumed parameters, and that $`\chi _{\mathrm{min}}^2`$, calculated for the optimum parameters, is equal to or smaller than $`\chi _{\mathrm{true}}^2`$. For each fake SF, we then calculated

$$\mathrm{\Delta }\chi ^2\chi _{\mathrm{true}}^2\chi _{\mathrm{min}}^2$$
(5)
which can be regarded as a measure of statistical fluctuation associated with eq. (4). Figure 4a shows the $`\mathrm{\Delta }\chi ^2`$ distribution of the 5000 simulations. As can be seen there, 4500 out of 5000 simulations have $`\mathrm{\Delta }\chi ^223`$. Hence, we regard $`\mathrm{\Delta }\chi ^2`$ = 23 as definition of the 90% confidence boundary. Similarly, $`\mathrm{\Delta }\chi ^2`$ = 3 and $`\mathrm{\Delta }\chi ^2`$ = 55 provides the 68.3% and 99% confidence boundary, respectively.

Now that we have determined the confidence limits associated with the $`\chi ^2`$ of eq. (4), we return to the M81 SF and the original set of 1000 fake SFs, and examine how much errors are allowed. The two-parameter confidence ranges determined in this way is shown in Figure 5a on the plane of 1/$`f_\mathrm{b}`$ and $`\alpha `$. The 90% confidence regions of the characteristic time scale and the PSD index are determined to be 1/$`f_\mathrm{b}800`$ days and $`\alpha =1.4\pm 0.2`$, respectively. Our detailed analysis has revealed longer 1/$`f_\mathrm{b}`$ compared to the rough estimation of Ishisaki et al. (1996): they considered RMS of a short-term light curve (the 9th observation) and a long-term light curve (the 1st through 10th observations) and assumed $`\alpha =1.40`$ to obtain $`f_\mathrm{b}`$=(1–3)$`\times 10^6`$ Hz (1/$`f_\mathrm{b}`$ = 4–12 days).

Now we can also examine whether the hypothesized shape of PSD is appropriate or not, using the $`\chi _{\mathrm{true}}^2`$ distribution of the 5000 simulations. Generally, a reasonable fit gives a $`\chi ^2`$ value similar to d.o.f. However, as mentioned in § 4.2, d.o.f. is not well determined in the present case. Instead, we compared the $`\chi _{\mathrm{min}}^2`$ value of M81 data with the $`\chi _{\mathrm{true}}^2`$ distribution of fake data, as a d.o.f. independent estimation of the absolute goodness. The $`\chi _{\mathrm{true}}^2`$ distribution of the fake light curves is shown in Figure 4b; again, the distribution is nearly independent of the assumed PSD parameters of the 5000 fake SFs. Because the $`\chi _{\mathrm{min}}^2`$ value of the actual M81 data, shown as an arrow, lies within 90% regions of the distribution, we conclude that the assumed two-component PSD is appropriate within the limitation of the current observation. Similarly, the single power-law PSD also turned out to be appropriate.

We admit that, however, it dose not mean that we can exclude other shapes of PSD which we have not yet examined. Especially, the low frequency portion of the M81 SF has such large scatters that other types of PSD might explain the observed data equally well, or even better. For example, it may have a negative or positive index in low frequencies, while we assumed zero index (white noise). A more complex PSD, for example three-component PSD seen in low-mass X-ray binaries (e.g. Wijnands et al. 1999) and Galactic black holes (e.b. Miyamoto et al. 1994), which consists of a red-noise part at low frequency, a white-noise part at intermediate frequency, and another red-noise part at high frequency, might explain the M81 SF as well.

### 4.4 Power Spectral Density

The remaining task is to determine $`P_0`$, the PSD normalization. By definition, the variation power, $`_{f_1}^{f_2}P(f)𝑑f`$, is equivalent to RMS of the light curve with bin width of 1/$`f_2`$ over the total observed span of 1/$`f_1`$. Consequently, we can relate $`P_0`$ to the RMS of the light curve, when we specify $`f_\mathrm{b}`$ and $`\alpha `$.

In § 4.2, step (iv) determines $`P_0`$ that corresponds to each 112-point light curve. If there were no data gap, all light curves that are generated from the same $`f_\mathrm{b}`$ and $`\alpha `$ would have the same $`P_0`$. In practice, however, $`P_0`$ depends on the random numbers used in the Monte-Carlo simulation. We accordingly made a histogram of 1000 values of $`P_0`$ assuming $`1/f_\mathrm{b}=10^{4.6}`$ days and $`\alpha =1.40`$, and determined 90% confidence region of $`P_0`$. In reference to the best-estimate $`P_0^{\mathrm{best}}`$, the 90% confidence region in $`P_0`$ becomes (0.8–1.2) $`P_0^{\mathrm{best}}`$. Taking into account the uncertainty associated with $`P_0`$, as well as $`f_\mathrm{b}`$ and $`\alpha `$ (Figure 5a), we have calculated a range of allowed PSDs with 68.3%, 90% and 99% confidences. In Figure 5b, we plot the PSD after dividing it by the square of its average flux; the plot is called normalized PSD (NPSD) and convenient to compare random variation of different objects (Miyamoto et al. 1994).

In connection with the normalization of the fake light curves, we would like to stress one point. In Figure 3, the difference in the assumed $`f_\mathrm{b}`$ is best reflected in the absolute level of SF over the small-time-lag region; shorter values of 1/$`f_\mathrm{b}`$ yield larger values of SF. However this does not mean that the information on longer time scales is unimportant. In fact, the clear separation among SFs with different PSDs occurs because we have normalized the RMS of each fake light curve. For the normalization to work properly, in turn, we need to know the overall data variance. When 1/$`f_\mathrm{b}`$ is longer than a few times the total span, the averaged SFs with different values of 1/$`f_\mathrm{b}`$ are no longer discriminated from one another, as shown in Figure 5a.

## 5 DISCUSSION

### 5.1 Comparison with Other Objects

Let us compare the estimated PSD of the M81 nucleus with those of other objects. We utilize the plot invented by Hayashida et al. (1998): namely, a (frequency) $`\times `$ (NPSD) plot, or $`fP(f)`$ plot. As Hayashida et al. (1998) pointed out, the $`fP(f)`$ plot displays relative contribution of the power at each frequency decade, like the $`\nu f_\nu `$ plot of energy spectrum. As illustrated in Figure 6, the power concentrates around $`f_\mathrm{b}`$, so that $`f_\mathrm{b}`$ is an important parameter characterized each system.

So far, Ptak et al. (1998) and Awaki et al. (1999) compared the high-frequency information of LLAGNs with those of the other objects and suggested that LLAGNs may have long time scale and/or large black-hole mass. However, there was an implicit or explicit assumption in their analysis, that the individual AGNs have the same PSD shape and same total area of NPSD on the $`fP(f)`$ plot. There was no concrete evidence of the assumption, and it should be examined by measuring the low frequency portion of the PSD, as we have performed. Moreover, the slope index of PSD differs object to object ($`1<\alpha <2`$). Consequently, their result strongly depend on the frequency where the NPSD amplitude was measured. Therefore, it is important to consider the PSD over wide frequency band.

In Figure 6, we show the M81 NPSD, derived from Figure 5b, together with those of Seyfert galaxies and a Galactic black-hole candidate taken from Hayashida et al. (1998) and Pounds et al. (1988). The figure reveals two facts. First, the power peak of M81 locates at a very low frequency for its low luminosity. Second, the total area of the M81 NPSD is similar to those of the other objects. Consequently, the observed lack of significant short-term variability in M81 is not because the total variation power is small, but because the power concentrates on low frequencies. This suggests that the other LLAGNs with small NPSD amplitude at high frequency also have large power on low frequencies.

### 5.2 Black-Hole Mass

Finally, we estimate the black-hole mass, $`M_{\mathrm{BH}}`$, of the M81 nucleus. Conventionally, characteristic variation time scales of AGNs have been translated into $`M_{\mathrm{BH}}`$ assuming a linear proportionality between the time scale and mass, and often utilizing Cygnus X-1 as a standard, although this estimation is only appropriate if the physics of accretion are similar between Cygnus X-1 and M81. This allows us to estimate $`M_{\mathrm{BH}}`$ in the M81 nucleus to be $`>4.0\times 10^7M_{}`$.

So far, $`M_{\mathrm{BH}}`$ of the M81 nucleus has been estimated in various ways. From measurements of the velocity of the broad-line region or narrow-line region, several authors estimated the mass within these regions of the M81 nucleus, although their results have large scatter: $`1\times 10^6M_{}`$$`1.8\times 10^7M_{}`$ (Peimbert & Torres-Peimbert 1981), $`(38)\times 10^5M_{}`$ (Filippenko & Sargent 1988), and $`7\times 10^5M_{}`$$`3\times 10^6M_{}`$ (Ho et al. 1996). From stellar velocity dispersion measurements, Keel (1989) obtained an upper limit of $`2\times 10^7M_{}`$, and Bower et al. (1996) obtained a lower limit of $`3\times 10^6M_{}`$. Our result falls on the relatively higher side of these previous estimates, though still consistent considering large errors associated with various estimates. Our estimate reinforces the view that LLAGNs have low X-ray luminosities because they have very low accretion rates, rather than because they have less massive black holes than the more luminous AGNs. This further encourages us to interpret LLAGNs as a relic of the past luminous AGNs.

## Acknowledgments

Thanks are due to Professor Hajime Inoue for permission to use the ASCA data of 15th and 16th observations and Professor Kazuhisa Mitsuda for helpful discussions on several points in the paper. This research was supported in part by a grant from Japan Society for the Promotion of Science.

References

Awaki H., Sakano M., Terashima Y., Hayashida K., 1999, AN, 320, 305

Barr P., Giommi P., Wamsteker W., Gilmozzi R., Mushotzky R., 1985, BAAS, 17, 608

Boettcher M., Reuter H.-P., Lesch H., 1997 A&A 326 33

Boller Th., Meurs E.J.A., Brinkmann W., Fink H., Zimmermann U., Adorf H.-M., 1992, A& A, 261, 57

Bower G.A., Wilson A.S., Heckman T.M., Richstone D.O., 1996, in The Physics of Liners in view of recent observations. ASP Conference Series Vol.103, p.163 edited by M. Eracleous

Burke B.E., Mountain R.W., Harrioson D.C., Bautz M.W., Doty J.P., Ricker G.R., Daniels P.J. 1991, IEEE Trans., ED-38 p1069

Edelson R., Nandra J., 1999, ApJ, 514, 682

Fabbiano G., 1988, ApJ, 325, 544

Filippenko A.V., Sargent W.L.W., 1988, ApJ, 324, 134

Freedman W.L., Hughes S.M., Madore B.F., Mould J.R., Lee M.G., Stetson P., Kennicutt R.C., Turner A. et al., 1994, ApJ, 427, 628

Hayashida K., Miyamoto S., Kitamoto S., Negoro H., Inoue H., 1998, ApJ, 500, 642

Heidt J., Wagner S.J., 1998, A&A, 329, 853

Ho L.C., Filippenko A.V., Sargent W.L.W., 1996, ApJ, 462, 183

Ho L.C., Filippenko A.V., Sargent W.L.W., 1997, ApJ, 487, 568

Hughes P.A., Aller H.D., Aller M.F., 1992, ApJ, 396, 469

Ishisaki Y., Makishima K., Inoue H., Iyomoto N., Kohmura Y., Mitsuda K., Mushotzky R.F., Petre R. et al., 1996, PASJ, 48, 237

Iyomoto N., 1999 PhD Thesis, University of Tokyo

Iyomoto N., Makishima K., 1999, AN, 320, 300

Iyomoto N., Makishima K., Fukazawa Y., Tashiro M., Ishisaki Y., Nakai N., Taniguchi Y., 1996, PASJ, 48, 231

Iyomoto N., Makishima K., Matsushita K., Fukazawa Y., Tashiro M., Ohashi T., 1998, ApJ, 503, 168

Kawaguchi T., Mineshige S., Umemura M., Turner E.L., 1998, ApJ, 504

Keel W., 1989, AJ, 98, 195

Kohmura Y., 1994, PhD Thesis, University of Tokyo

Kohmura Y., Inoue H., Aoki T., Ishida M., Itoh M., Kotani T., Tanaka Y., Ishisaki Y. et al., 1994, PASJ, 46, 157

Leighly K.M., O’brien P.T., 1997, ApJL, 481, L15

Makishima K., Fujimoto R., Ishisaki Y., Kii T., Loewenstein M., Mushotzky R., Serlemitsos P., Sonobe T. et al., 1994, PASJ, 46, L77

Makishima K. Tashiro M., Ebisawa K., Ezawa H., Fukazawa Y., Gunji S., Hirayama M., Idesawa E. et al. 1996 PASJ, 48, 171

Miyamoto S., Kitamoto S., Iga S., Hayashida K., Terada K., 1994, ApJ, 435, 398

Ohashi T., Ebisawa K., Fukazawa Y., Hiyoshi K., Horii M., Ikebe Y., Ikeda H., Inoue H. et al., 1996, PASJ, 48, 157

Paltani S., Courvoisier T.J.-L., Blecha A., Bratschi P., 1997, A&A, 327, 539

Peimbert M., Torres-Peimbert S., 1981, ApJ, 245, 845

Pellegrini S., Cappi M., Bassani L., Malaguti G., Palumbo G.G.C., Persic M., 2000, A&A, 353, 447

Petre R., Mushotzky R.F., Serlemitsos P.J., Jahoda K., Marshall F.E., 1993, ApJ, 418, 644

Pounds K.A., McHardy I.M., 1988 in PHYSICS OF NUETRON STARS AND BLACK HOLES p.285 edited by Tanaka Y. Universal Academy Press

Ptak A., Serlemitsos P., Yaqoob T., Mushotzky R., 1999, ApJS, 120, 179

Ptak A., Yaqoob T., Mushotzky R.F., Serlemitsos P., Greffiths R., 1998, ApJL, 501, L37

Serlemitsos P., Ptak A., Yaqoob T., 1996 in The Physics of Liners in view of recent observations. ASP Conference Series Vol.103, p.70 edited by M. Eracleous

Simonetti J.H., Cordes J.M., Heeschen D.S., 1985, ApJ, 296, 46

Terashima Y., Ptak A., Fujimoto R., Itoh M., Kunieda H., Makishima K., Serlemitsos P.J., 1998, ApJ, 496, 210

Tsuru T., 1992, PhD Thesis, University of Tokyo

Uno S., Mitsuda K., Inoue H., Takahashi T., Makino F., Makishima K., Ishisaki Y., Kohmura Y. et al., 1998, IAU Symopsium No. 188 P.245

Wijnands R., van der Klis M., 1999, ApJ, 522, 965


text:
# Partial synchronicity and the (max,+) semiring

##

Abstract.  In this paper we illustrate how non-stochastic (max,+) techniques can be used to describe partial synchronization in a Discrete Event Dynamical System. Our work uses results from the spectral theory of dioids and analyses (max,+) equations describing various synchronization rules in a simple network. The network in question is a transport network consisting of two routes joined at a single point, and our Discrete Events are the departure times of transport units along these routes. We calculate the maximum frequency of circulation of these units as a function of the synchronization parameter. These functions allow us further to determine the waiting times on various routes, and here we find critical parameters (dependent on the fixed travel times on each route) which dictate the overall behavoiur. We give explicit equations for these parameters and state the rules which enable optimal performance in the network (corresponding to minimum waiting time).

## 1 Introduction

Considerable advances have been made in recent years in the algebraic description of Discrete Event Dynamical Systems (DEDS). Chief among these is the approach using the so-called (max,+) semiring , with corresponding development of algorithms for spectral calculations  and concrete applications .

In this paper we illustrate how non-stochastic (max,+) techniques can be used to describe partial synchronization (in a manner to be defined below) in a DEDS. Our work uses results from the spectral theory of dioids and analyses (max,+) equations describing various synchronization rules in a simple network.

The paper is structured as follows. In Section 2 we review necessary mathematical properties of the (max,+) semiring. Section 3 introduces the network on which our DEDS is based, and in Section 4 we describe the different (max,+) equations we will solve. The resulting solutions lead, in Section 5, to critical phenomena in the network, and we calculate which among the (max,+) equations give optimal results. A numerical example illustrates our results in Section 6 we conclude in Section 7.

## 2 (max,+) semiring

The (max,+) semiring (or dioid) $`_{\mathrm{max}}`$ is the set $`\{\mathrm{}\}`$ with the two operations $`,`$ defined by $`ab=\mathrm{max}(a,b)`$ and $`ab=a+b`$. Note that both operations possess identity elements, $`\mathrm{}`$ for $``$ and 0 for $``$, but while $``$ is invertible, $``$ is idempotent: $`aa=a`$. In the remainder of this paper we will often omit the explicit multiplication symbols $`\times ,`$, as is the convention. Our convention will be that if an expression contains any explicit symbol from $`_{\mathrm{max}}`$ we will assume all such hidden symbols are also from $`_{\mathrm{max}}`$: $`a+bc`$ means $`a+(b\times c)`$, $`abc`$ means $`a(bc)`$, and $`(ab)(c+d)`$ means $`(ab)(c+d)`$.

These operations can be extended to matrices $`A,B_{\mathrm{max}}^{n\times n}`$ by defining $`(AB)_{ij}=A_{ij}B_{ij}`$ and $`(AB)_{ij}=_k(A_{ik}B_{kj}).`$ As this may seem unusual to the novice reader, a small example is in order.

$`\left(\begin{array}{cc}3& 1\\ 0& 5\end{array}\right)\left(\begin{array}{c}4\\ 0\end{array}\right)`$ (5)
$`=\left(\begin{array}{c}(34)(10)\\ (04)(50)\end{array}\right)=\left(\begin{array}{c}7\\ 5\end{array}\right).`$ (10)
We now look at solutions to the (max,+) spectral problem $`AX=\lambda X`$.
Theorem 1. Let $`A_{\mathrm{max}}^{n\times n}`$. There exists a maximal eigenvalue

$`\lambda ={\displaystyle \underset{k=1}{\overset{n}{}}}{\displaystyle \underset{i_1,i_2,\mathrm{},i_k}{}}(A_{i_1i_2}A_{i_2i_3}\mathrm{}A_{i_ki_1})/k.`$ (11)
If we extend our notation to $`a^n\stackrel{\mathrm{def}}{=}aa\mathrm{}a=n\times a`$ and $`\mathrm{Tr}(A)=_{k=1}^nA_{kk}`$ equation (11) reads

$`\lambda ={\displaystyle \underset{k=1}{\overset{n}{}}}(\mathrm{Tr}(A^k))^{1/k}.`$ (12)
In the standard way, we associate with $`A`$ a weighted digraph (precedence graph) $`𝒢(A)`$ with nodes $`𝒩=\{1,2,\mathrm{},n\}`$ and edges $`=\{(i,j)|A_{ji}\mathrm{}\}`$. The weight of edge $`(i,j)`$ is simply $`A_{ji}`$. Denoting by $`ij`$ the existence of a path from $`i`$ to $`j`$ we define $`A`$ to be irreducible (and $`𝒢(A)`$ to be strongly connected) iff $`ij1i,jn`$.
Theorem 2. If $`A_{\mathrm{max}}^{n\times n}`$ is irreducible, then $`\lambda `$ as given by (11) is unique.
Note that from a graph-theoretic point of view, $`\lambda `$ is merely the maximum cycle mean of $`𝒢(A)`$.

Our DEDS will be described by equations in $`_{\mathrm{max}}`$: Specifically, denoting the time of the $`k`$th occurence of event $`i`$ by $`x_i(k)`$, we are interested in equations of the form

$`x_i(k)=f(x_j(kl))`$ (13)
where we have a finite family of $`1i,jn`$ events and $`1lk1`$. For linear equations in $`_{\mathrm{max}}`$, (13) becomes

$`x_i(k)={\displaystyle \underset{j,l}{}}(f_{ij}^lx_j(kl))`$ (14)
and we will show in Section 4 how this can be written as a matrix equation with unitary retard

$`x_i(k)={\displaystyle \underset{j}{}}(A_{ij}x_j(k1)).`$ (15)
We are interested in the asymptotic properties of this system as $`k\mathrm{}`$ so we study $`lim_k\mathrm{}A^k`$.
Theorem 3. If $`A_{\mathrm{max}}^{n\times n}`$ is irreducible, there exists integers $`M,c(A)`$ such that

$`A^{kc(A)}=(\lambda )^{c(A)}A^kkM.`$ (16)
We call $`c(A)`$ the cyclicity of $`A`$ (see  for details of how to calculate this integer).

## 3 Model

We consider an elementary model of two tourbuses with partial synchronization. Each tourbus goes on a fixed circular route $`R_1`$ and $`R_2`$, stopping every so often to drop off and pick up passengers. The tour buses meet at one location, a downtown station S, where passengers can pass from one bus to the other one. The physical network consists of two loops meeting at S, as shown in Figure 1.

We assume fixed travel times $`T_1`$ and $`T_2`$ for $`R_1`$ and $`R_2`$ respectively, and without loss of generality let $`T_1T_2`$. The simplest model is to describe this system using the (max,+) equations $`X(k+1)=AX(k)`$, with

$`A=\left(\begin{array}{cc}a_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right),X(k)=\left(\begin{array}{c}x_1(k)\\ x_2(k)\end{array}\right)`$ (21)
and $`x_i(k)`$ is the departure time of bus $`i`$ from station $`S`$ along route $`R_i`$. We consider the following extreme cases.
CASE M1 (no synchronization):

We have

$`A=\left(\begin{array}{cc}T_1& \mathrm{}\\ \mathrm{}& T_2\end{array}\right)`$ (24)
so that the graph $`𝒢(A)`$ is disconnected into two parts each having its own cycle mean $`\lambda _1=T_1`$ and $`\lambda _2=T_2`$ as shown in Figure 2 (the nodes here correspond to departure times  of buses).

CASE M2 (complete synchronization):

We have

$`A=\left(\begin{array}{cc}T_1& T_2\\ T_1& T_2\end{array}\right)`$ (27)
so that the graph $`𝒢(A)`$ is strongly connected. There are 3 circuits in this graph with cycle means $`T_1`$, $`T_2`$, and $`\frac{T_1+T_2}{2}`$, so the maximum cycle mean gives an eigenvalue of $`\lambda =T_2`$ (see Figure 3).

## 4 Partial Synchronicity

For the case of complete synchronization in section 3, it is clear that the large value of $`\lambda =T_2`$ slows down the overall system. We further suppose in this section that we cannot alter the travel time $`T_2`$, but that we have a large amount of additional stock (buses) which we can add to the network.

We add to route $`R_2`$ a large number $`m`$ of buses such that the interval between departure times on this route becomes a very small number $`\delta `$ (which we can think of as $`T_2/m`$ if required). Denote by $`y_1(k),y_2(k),\mathrm{},y_m(k)`$ the departure times of the $`m`$ new buses, which operate according to the rules

$`y_i(k)=`$ $`\delta y_{i1}(k1)`$
$`y_0(k)=`$ $`\delta y_m(k1).`$ (28)
For this route we do not produce a timetable, as we assume to a first approximation there is virtually no waiting for a bus.

For this system, it would seem there should be no interdependence between $`x_1(k)`$ and $`x_2(k)`$. However, what of a passenger who wants to make a non-stop round trip $`R_1R_2R_1`$ – would some synchronization rule benefit such a passenger? To investigate further we set $`T_2=nT_1+r`$ and look at the $`_{\mathrm{max}}`$ equation $`x_1(k)=f(x_1(kl))`$ for different values of $`l`$ and linear $`f`$.
CASE P1:

$`x_1(k)=T_1x_1(k1)T_1T_2x_1(k(l+1))`$ (29)
Using $`n`$ artificial states (i.e. states that do not in fact correspond to the departure of any physical bus) $`z_1,z_2,\mathrm{},z_l`$, we can rewrite (29) as

$`x_1(k)`$ $`=`$ $`T_1x_1(k1)rT_1z_l(k1)`$
$`z_i(k)`$ $`=`$ $`T_1^{n/l}z_{i1}(k1),\text{}i=1,\mathrm{},l`$ (30)
$`z_0(k)`$ $`\stackrel{\mathrm{def}}{=}`$ $`x_1(k).`$
At the expense of increasing the dimension, our system of equations now has unitary retard as in (15) with column vector

$`X(k)=\left(\begin{array}{c}x_1(k)\\ z_l(k)\\ z_{l1}(k)\\ \mathrm{}\\ z_1(k)\end{array}\right)`$ (36)
and matrix

$`A=\left(\begin{array}{cccccc}T_1& T_1+r& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& nT_1/l& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& nT_1/l& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ nT_1/l& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right).`$ (42)
The full precedence graph $`𝒢(A)`$ for the system $`x_1,z_1,\mathrm{},z_l,x_2,y_1,\mathrm{},y_m`$ is shown in Figure 4.

The disconnected parts of the graph $`𝒢(A)`$ have maximum cycle means $`\lambda _2=\delta `$ and

$`\lambda _1=\{\begin{array}{cc}T_1& \text{if }l>n\\ \frac{(n+1)T_1+r}{l+1}& \text{otherwise}\end{array}`$ (45)
since $`T_1r=T_1`$ (by definition of $`r`$) and, for $`l>n`$ (equivalently $`l=l(n+1)`$),

$`{\displaystyle \frac{(n+1)T_1+r}{l+1}}<{\displaystyle \frac{(n+1)T_1+T_1}{l+1}}={\displaystyle \frac{(n+2)T_1}{l+1}}T_1.`$ (46)
Note that $`l>n`$ effectively corresponds to no synchronization, while for the choice $`ln`$, $`(nl+1)`$ out of every $`(n+1)`$ buses are synchronized. The minimum possible $`\lambda _1`$ is $`\lambda _1=T_1`$ given by $`l>n`$.
CASE P2:

$`x_1(k)=T_1x_1(k1)rT_1^{l+1}x_1(k(l+1))`$ (47)
As in the previous case P1, we can rewrite (47) as

$`x_1(k)`$ $`=`$ $`T_1x_1(k1)rT_1z_l(k1)`$
$`z_i(k)`$ $`=`$ $`T_1z_{i1}(k1),\text{}i=1,\mathrm{},l`$ (48)
$`z_0(k)`$ $`\stackrel{\mathrm{def}}{=}`$ $`x_1(k).`$
and the system is described by $`X(k+1)=AX(k)`$ with $`X(k)`$ as in (36) and

$`A=\left(\begin{array}{cccccc}T_1& T_1+r& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& T_1& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& T_1& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ T_1& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right).`$ (54)
From the precedence graph shown in Figure 5 we have $`\lambda _2=\delta `$ and

$`\lambda _1={\displaystyle \frac{(l+1)T_1+r}{l+1}}`$ (55)
since obviously

$`{\displaystyle \frac{(l+1)T_1+r}{l+1}}T_1`$ (56)
for all values of $`l`$.

This system corresponds to one out of every $`(l+1)`$ buses being synchronized. The limiting case is

$`\underset{l\mathrm{}}{lim}{\displaystyle \frac{(l+1)T_1+r}{l+1}}=T_1`$ (57)
so the minimum value of $`\lambda _1`$ occurs at $`l=\mathrm{}`$.

## 5 Optimal Solutions

In both cases P1 and P2 of Section 4 the minimum value of $`\lambda _1`$ is $`T_1`$: In P1 it is given by any choice $`l>n`$ while in P2 it is given by the extreme value $`l=\mathrm{}`$. These values both correspond to no synchronization, and this is the result we would expect, that the system runs faster if the buses do not wait for each other.

If our goal is to choose the synchronization rule which maximizes the speed of the system, the problem is solved. In this Section we look at the separate problem of minimizing the waiting time of passengers that make a circular journey over all or part of the network (by circular we mean a journey that begins and ends at the same point). We itemize such journeys $`J_i`$ as follows:

$`J_1`$ $`:`$ $`R_1R_1`$
$`J_2`$ $`:`$ $`R_1R_2R_1`$
$`J_3`$ $`:`$ $`R_2R_1R_2`$
$`J_4`$ $`:`$ $`R_2R_2`$
For cases P1, P2 we look at regions $`ln`$ and $`l>n`$ to give synchronization rules $`S_i`$ as follows:

$`S_1`$ $`:`$ $`\text{Equation (}\text{29}\text{)},ln`$
$`S_2`$ $`:`$ $`\text{Equation (}\text{29}\text{)},l>n`$
$`S_3`$ $`:`$ $`\text{Equation (}\text{47}\text{)},ln`$
$`S_4`$ $`:`$ $`\text{Equation (}\text{47}\text{)},l>n`$
For $`i,j\{1,2,3,4\}`$ we pair each journey $`J_i`$ with rules $`S_j`$ to give 16 different models. Let $`w_{ij}`$ be the total waiting time on journey $`J_i`$ with rule $`S_j`$. These waiting times can be broken into two components, the waiting time before boarding the bus $`w_{ij}^b`$ and the waiting time in transit $`w_{ij}^t`$, which can only occur for journeys $`J_2`$ and $`J_3`$. From our calculation of $`\lambda _1`$ and $`\lambda _2`$ in Section 4 we determine $`w_{ij}^b,w_{ij}^t`$ and $`w_{ij}=w_{ij}^b+w_{ij}^t`$ to be the values given by equations (69)–(79) in the Appendix.

We now make some observations on the solutions given by equation (79). Our goal is to minimize the waiting times, or more specifically the average waiting time per passenger as a function of the synchronization parameter $`l`$: If $`n_i`$ people take journey $`J_i`$ and $`N=_in_i`$ then we want to find the minimum of

$`\overline{W_i}\stackrel{\mathrm{def}}{=}{\displaystyle \frac{_jn_jw_{ji}}{N}}`$ (58)
Note that

* $`w_{i1}`$ is minimized by choosing the largest possible value of $`l`$, i.e. $`l=n`$.
* $`w_{i2}=lim_l\mathrm{}w_{i4}`$
* similarly $`w_{i3}`$ is minimized by choosing $`l=n`$.

Finally examine $`w_{i4}`$. $`w_{14}`$ and $`w_{34}`$ are both clearly minimized by $`l=\mathrm{}`$. To minimize $`w_{24}`$ we view it as a continuous function of $`l`$ and solve

$`{\displaystyle \frac{dw_{24}}{dl}}={\displaystyle \frac{d}{dl}}{\displaystyle \frac{(3l2n+1)T_1+(2n2l+3)r}{2(l+1)}}=0.`$ (59)
$`l`$ drops out of equation (59), meaning the existence of an extremum does not depend on $`l`$ but rather on the solution of $`(2n+2)T_1=(2n+5)r`$. Thus we have a critical point

$`r_c=\left({\displaystyle \frac{2n+2}{2n+5}}\right)T_1.`$ (60)
The slope of the function $`w_{24}(l)`$ can now be written in terms of the critical remainder $`r_c`$ as

$`{\displaystyle \frac{dw_{24}}{dl}}={\displaystyle \frac{(2n+5)(r_cr)}{(l+1)^2}}.`$ (61)
Hence $`w_{24}(l)`$ is

* strictly increasing for $`r<r_c`$: minimum is obtained by choosing lowest value of $`l`$, i.e. $`l=n+1`$
* strictly decreasing for $`r>r_c`$: minimum is obtained by choosing highest value of $`l`$, i.e. $`l=\mathrm{}`$
* constant for $`r=r_c`$: choose any value of $`l`$

In equation (84) we rewrite $`W`$ with the optimal values of $`l`$ chosen as indicated above: for $`w_{i4}`$ we insert in the matrix $`W^m(W^{\mathrm{min}})`$ the values at $`l=n+1`$, since the values at $`l=\mathrm{}`$ correspond to column 2, $`w_{i2}`$.

Now that we have chosen $`l`$ we should choose the synchronization rule $`S_j`$ that minimizes $`\overline{W_i}`$. In rows 1, 3 and 4 of $`W^m`$ it is clear that the minimum values are in column 2 corresponding to $`S_2`$. Let us look at row 2, $`w_{2j}^m`$. In Figure 6 we plot $`w_{21}^m,w_{22}^m`$ and $`w_{24}^m`$ as functions of $`r`$.

$`r_s=\left({\displaystyle \frac{2n+2}{2n+3}}\right)T_1`$ (62)
along with equation (60) defines three different regimes as follows:

* $`r>r_s:w_{22}^m<w_{21}^m<w_{24}^m`$
* $`r_c<r<r_s:w_{21}^m<w_{22}^m<w_{24}^m`$
* $`r<r_c:w_{21}^m<w_{24}^m<w_{22}^m`$

Further we note the importance of the critical value $`r_s`$ is that for $`r>r_s`$ rule $`S_2`$ gives shortest waiting times for all journeys $`J_i`$, and in this regime $`\overline{W_2}`$ is minimal irrespective of the values $`n_i`$ (see equation (58)). For $`r<r_s`$ the values of $`n_i`$ dictate which $`\overline{W_i}`$ is minimal: In particular, if

$`n_2>{\displaystyle \frac{r}{(2n+3)(r_sr)}}(3n_1+n_3)`$ (63)
then $`\overline{W_1}`$ is minimal. We now have a rule that governs optimum performance in our network: If $`r>r_s`$ implement rule $`S_2`$, otherwise monitor passenger numbers $`n_i`$ and if at any stage equation (63) holds, implement rule $`S_1`$.

Observe further from equation (62) that $`4T_1/5r_s<T_1`$, since $`1n<\mathrm{}`$. This means that even for a random choice of $`T_1`$ and $`T_2`$, the chances that $`r`$ is less than $`r_s`$ are greater than 80%:

$`P(r<r_s|\text{ arbitrary }T_1,T_2)0.8`$ (64)
## 6 Example

Let $`T_1=3`$ units of time and $`n=2`$ so that $`r_s=18/7`$ as given by equation (62).

To examine firstly the regime $`r>r_s`$, let $`r=8/3`$, given by $`T_2=26/3`$ units of time. Under rule $`S_2`$ we calculate $`w_{22}^m=33/18`$ with actual departure times $`x_1(k)=0,3,6,9,12,15,\mathrm{}`$ Under rule $`S_3`$ we have $`w_{23}^m=35/18`$ with departure times $`x_1(k)=0,3,6,35/3,44/3,53/3,\mathrm{}`$ Thus rule $`S_2`$ gives optimal performance on $`J_2`$.

Secondly, let $`r=7/3`$ given by $`T_2=25/3`$, so that we are in the regime $`r<r_s`$. Rule $`S_2`$ gives $`w_{22}^m=39/18`$ with departure times $`x_1(k)=`$ 0, 3, 6, 9, 12, 15,$`\mathrm{}`$ Under rule $`S_3`$ we have $`w_{23}^m=34/18`$ with departure times $`x_1(k)=`$ 0, 3, 6, 34/3, 43/3, 52/3,$`\mathrm{}`$ Here optimal performance on $`J_2`$ is given by rule $`S_3`$.

Note in both regimes of this example the departure times $`x_1(k)`$ under $`S_2`$ are identical, the reason for the differences in $`w_{22}^m`$ lying in the transfer waiting time.

## 7 Conclusion

In this paper we have solved the (max,+) equations for a specific model with two routes joining at one point. We have shown how partial synchronicity can be described by these equations by introducing a synchronization parameter $`l`$. Using results from (max,+) spectral theory we have calculated the maximum frequency (minimum cycle mean) of buses on such routes. These numbers allow us to calculate waiting times on various routes: We observe new critical behaviour which depends on the relation between $`T_1`$ and $`T_2`$, the (fixed) travel times on each route. Finally we state the rules for optimal performance in this network.

Further work is envisaged in looking at more general (max,+) equations, in particular a system with $`m_1`$ buses on $`R_1`$ and $`m_2`$ on $`R_2`$. Another direction in which this work could possibly be extended is to have further routes $`R_3,R_4,\mathrm{}`$

## Appendix

$`W^b=\left(\begin{array}{cccc}\frac{(n+1)T_1+r}{2(l+1)}& \frac{T_1}{2}& \frac{(l+1)T_1+r}{2(l+1)}& \frac{(l+1)T_1+r}{2(l+1)}\\ \frac{(n+1)T_1+r}{2(l+1)}& \frac{T_1}{2}& \frac{(l+1)T_1+r}{2(l+1)}& \frac{(l+1)T_1+r}{2(l+1)}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)`$ (69)
$`W^t=\left(\begin{array}{cccc}\frac{(nl)T_1+r}{l+1}& 0& \frac{r}{l+1}& \frac{r}{l+1}\\ 0& T_1r& \frac{(nl)r}{(l+1)}& \frac{T_1+(T_1r)[l(n+1)]}{(l+1)}\\ \frac{(n+1)T_1+r}{2(l+1)}& \frac{T_1}{2}& \frac{(l+1)T_1+r}{2(l+1)}& \frac{(l+1)T_1+r}{2(l+1)}\\ 0& 0& 0& 0\end{array}\right)`$ (74)
$`W=W^b+W^t=\left(\begin{array}{cccc}\frac{(3n2l+1)T_1+3r}{2(l+1)}& \frac{T_1}{2}& \frac{(l+1)T_1+3r}{2(l+1)}& \frac{(l+1)T_1+3r}{2(l+1)}\\ \frac{(n+1)T_1+r}{2(l+1)}& \frac{3T_1}{2}r& \frac{(l+1)T_1+(2n2l+1)r}{2(l+1)}& \frac{(3l2n+1)T_1+(2n2l+3)r}{2(l+1)}\\ \frac{(n+1)T_1+r}{2(l+1)}& \frac{T_1}{2}& \frac{(l+1)T_1+r}{2(l+1)}& \frac{(l+1)T_1+r}{2(l+1)}\\ 0& 0& 0& 0\end{array}\right)`$ (79)
$`W^m(W^{\mathrm{min}})=\left(\begin{array}{cccc}\frac{T_1}{2}+\frac{3r}{2(n+1)}& \frac{T_1}{2}& \frac{T_1}{2}+\frac{3r}{2(n+1)}& \frac{T_1}{2}+\frac{r}{2(n+2)}\\ \frac{T_1}{2}+\frac{r}{2(n+1)}& \frac{3T_1}{2}r& \frac{T_1}{2}+\frac{r}{2(n+1)}& \frac{T_1}{2}+\frac{2T_1+r}{2(n+2)}\\ \frac{T_1}{2}+\frac{r}{2(n+1)}& \frac{T_1}{2}& \frac{T_1}{2}+\frac{r}{2(n+1)}& \frac{T_1}{2}+\frac{r}{2(n+2)}\\ 0& 0& 0& 0\end{array}\right)`$ (84)
Michael Mc Gettrick
Department of Information Technology
National University of Ireland
Galway
Ireland
e-mail: michael.mcgettrick@nuigalway.ie


text:
# 1 Dynamical model for the 𝑝𝑑→³He(𝜋𝜋)⁰ reaction in terms of sequential 𝑝𝑝→𝑑𝜋⁺ and 𝜋⁺𝑛→(𝜋𝜋)⁰𝑝 processes. There is an analogous contribution from intermediate neutral pions.

TSL/ISV-2000-0233

Production of isoscalar pion pairs in the
$`𝒑𝒅\mathbf{}^\mathrm{𝟑}`$He$`𝝅𝝅`$ reaction near threshold

Göran Fäldt$`^{a,}`$111Electronic address: faldt@tsl.uu.se, Anders Gårdestig$`^{a,}`$222Electronic address: grdstg@uclapp.physics.ucla.edu$`^,`$333Current address: Department of Physics & Astronomy, UCLA, Box 951547,
xxxLos Angeles, CA 90095-1547, USA, Colin Wilkin$`^{b,}`$444Electronic address: cw@hep.ucl.ac.uk

a Division of Nuclear Physics, Box 535, 751 21 Uppsala, Sweden

b Physics and Astronomy Department, UCL, London, WC1E 6BT, UK

## Abstract

The production near threshold of isoscalar pion pairs in the $`pd^3`$He$`(\pi \pi )^0`$ reaction is estimated in a two-step model which successfully describes the production of $`\eta `$, $`\omega `$ and $`\eta ^{}`$ mesons. A virtual pion beam, generated through an $`NNd\pi `$ reaction on one of the nucleons in the deuteron, produces a second pion via a $`\pi N\pi \pi N`$ reaction on the other nucleon. Using the same scale factor as for heavy meson production, the model reproduces the total $`\pi ^0\pi ^0`$ production rate determined at an excess energy of 37 MeV. There are some indications in the data for a suppression of events with low $`\pi \pi `$ masses, as in the $`\pi ^{}p\pi ^0\pi ^0n`$ reaction, and this is confirmed within the model. The model suggests that a significant fraction of the charged pion production in the $`pd^3`$He$`\pi ^+\pi ^{}`$ reaction at $`Q=70`$ MeV might be associated with isoscalar pion pairs, though this does not explain the strong dependence observed on the $`\pi ^+\pi ^{}`$ relative momentum angle.

Keywords: double pion production

PACS: 25.10.+s, 25.40.Ve, 25.60.Dz Corresponding author:
Colin Wilkin,
Physics & Astronomy Dept., UCL, Gower St., London WC1E 6BT.

The study of neutral two-pion production through the $`pd^3`$He$`X^0`$ reaction has a long history. At excess energies $`Q`$ (the c.m. kinetic energy in the final state) around 200-300 MeV, sharp structure is seen at missing masses of about 310 MeV/c2 . The absence of any significant strength in the $`pd^3`$H$`X^+`$ channel at low $`m_X`$ means that the effect is associated with isospin-zero pion-pion pairs which, because of the available energy, must be dominantly in $`s`$-waves. Although a quantitative explanation of the ABC enhancement has not yet been provided for this reaction, a similar effect in $`npdX^0`$ has been shown to originate from the excitation of two $`\mathrm{\Delta }`$-isobars . The prominent ABC peaks in the $`dd^4`$He$`X^0`$ case have also been shown to be due to double pion $`p`$-wave production .

The experimental picture changes dramatically at lower energies. For values of $`Q`$ around 70-90 MeV, the angular distributions observed by the MOMO group  for the exclusive $`pd^3`$He$`\pi ^+\pi ^{}`$ reaction suggest strongly that the $`\pi ^+\pi ^{}`$ spectrum is mainly $`p`$-wave in nature, and hence has isospin-one. The measurement of the $`\pi ^0\pi ^0/\pi ^+\pi ^{}`$ charge ratio at CELSIUS at an excess energy (with respect to the $`\pi ^0\pi ^0`$ threshold) of $`Q=37`$ MeV  shows that there is significant $`I=1`$ production even at this much smaller $`Q`$. The isoscalar $`\pi ^0\pi ^0`$ spectrum, determined either by direct measurement  or through the subtraction of exclusive 3He$`\pi ^+\pi ^{}`$ data from an inclusive $`pd^3`$He$`X^0`$ measurement , shows that there is no $`s`$-wave ABC enhancement at low $`Q`$. On the contrary, there are rather indications that the $`s`$-wave cross section is actually suppressed at low $`\pi ^0\pi ^0`$ masses as compared to phase space . It is the aim of the present paper to demonstrate that near-threshold isoscalar two-pion production in the $`pd^3`$He$`\pi \pi `$ can be described in terms of sequential single-pion production.

The large momentum transfers required to produce heavy mesons, such as the $`\eta `$ or $`\omega `$, through the $`pd^3`$He$`X^0`$ reaction mean that two-step processes which minimise the momentum mismatch in the nuclear wave functions can provide the dominant driving force. In one such model, a pion produced on one of the nucleons in the target deuteron is converted into the observed heavy meson through an interaction on the second of the target nucleons . Apart from an ad hoc overall normalisation factor $`N2.4`$, which may reflect the retention of only bound intermediate deuteron states in the calculation, this approach describes well the threshold amplitudes for producing $`\eta `$, $`\omega `$ and $`\eta ^{}`$, though the experimental $`\varphi `$ yield is a little too high . This success may be attributed to the fact that the intermediate pion in the diagram is close to its mass shell. We wish to apply the same model to isoscalar $`\pi \pi `$ production by introducing rather a final $`\pi N(\pi \pi )^0N`$ process, as in Fig. 1. For definiteness, we consider $`\pi ^0\pi ^0`$ production; estimates of charged pion production in the $`I=0`$ channel will then follow from isospin invariance, after correcting for the pion mass difference.

The principal difference with the earlier work  is that the low mass $`\pi ^0\pi ^0`$ system is in a $`0^+`$ state and so there is a parity change at the $`\pi ^+n\pi ^0\pi ^0p`$ vertex. Parameterising this amplitude in terms of two-component Pauli spinors $`u_{p(n)}`$ as

$$M(\pi ^+n\pi ^0\pi ^0p)=a(m_{\pi \pi },W_{\pi N})u_p^{}𝝈𝒑_\pi u_n,$$
(1)
the corresponding differential cross section for $`s`$-wave production is

$$\text{d}\sigma (\pi ^+n\pi ^0\pi ^0p)=\text{d}\sigma (\pi ^{}p\pi ^0\pi ^0n)=\frac{1}{64\pi ^3}\frac{pp^{}}{W_{\pi N}^2}|a(m_{\pi \pi },W_{\pi N})|^2k_\pi ^{}\text{d}m_{\pi \pi }.$$
(2)
Here $`p`$ and $`p^{}`$ are the incident and final nucleon momenta in the overall c.m. system where the total energy is $`W_{\pi N}`$. In the $`\pi \pi `$ rest frame, $`k_\pi ^{}`$ is the relative momentum, which is related to the $`\pi \pi `$ invariant mass through $`m_{\pi \pi }=2\sqrt{k_\pi ^{\mathrm{\hspace{0.17em}2}}+m_\pi ^2}`$. We shall neglect the angular dependence of the amplitudes in the present work.

Because of the nature of the two-step process in Fig. 1, only small Fermi momenta are required. Working to first order in these momenta, as in , we find that the amplitudes are proportional to the complex form factors

$$S_{\alpha \beta }(𝐖,𝐕)=(2\pi )^3_0^{\mathrm{}}\text{d}te^{it\mathrm{\Delta }E_0}\psi _\alpha ^{}(t𝐖)\phi _\beta (t𝐕).$$
(3)
These involve integrals over configuration–space deuteron $`(\phi _\beta )`$ and 3He $`(\psi _\alpha )`$ wave functions, where $`\alpha ,\beta =(0,2)`$ represent nuclear $`S`$– and $`D`$–state components. The energy mismatch $`\mathrm{\Delta }E_0`$ between the intermediate and external energies for zero Fermi momenta is generally small for near-threshold heavy meson production in this model.

The relativistic relative velocity vectors $`𝐕`$ and $`𝐖`$,

$`𝐕`$ $`=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{1}{E_\pi (\frac{2}{3}𝐩_{\pi \pi }\frac{1}{2}𝐩_d)}}𝐩_{\pi \pi }{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{E_\pi (\frac{2}{3}𝐩_{\pi \pi }\frac{1}{2}𝐩_d)}}+{\displaystyle \frac{1}{E_n(\frac{1}{2}𝐩_d)}}\right]𝐩_d,`$ (4)
$`𝐖`$ $`=`$ $`{\displaystyle \frac{2}{3}}\left[{\displaystyle \frac{1}{E_\pi (\frac{2}{3}𝐩_{\pi \pi }\frac{1}{2}𝐩_d)}}+{\displaystyle \frac{1}{E_d(\frac{2}{3}𝐩_{\pi \pi })}}\right]𝐩_{\pi \pi }+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{E_\pi (\frac{2}{3}𝐩_{\pi \pi }\frac{1}{2}𝐩_d)}}𝐩_d,`$
where $`𝐩_{\pi \pi }`$ is the total $`\pi \pi `$ momentum vector and $`𝐩_d`$ that of the initial deuteron in the overall c.m. frame. The component of $`𝐕`$ along $`𝐩_d`$ must be subjected to a Lorentz contraction . The relativistic energies $`E_i`$ are evaluated at the values of the momenta indicated.

For zero Fermi momenta, the $`ppd\pi ^+`$ amplitudes should be evaluated in the forward direction for threshold heavy meson production. The forward direction assumption is also very good even away from threshold provided that $`p_dp_{\pi \pi }`$, as it is in cases under investigation. There are two $`ppd\pi ^+`$ amplitudes in the forward direction but, at the energies required here, the helicity-zero completely dominates over the helicity-one . Keeping then only the dominant amplitude $`A`$, the c.m. differential cross section is

$$\frac{\text{d}\sigma }{\text{d}\mathrm{\Omega }}(ppd\pi ^+)=\frac{1}{128\pi ^2}\frac{p_\pi }{p_pW_{pp}^2}|A|^2.$$
(5)
The evaluation of the unpolarised differential cross section in this model is similar to that for the production of single heavy mesons  and leads to

$`\text{d}\sigma (pd^3\text{He}\pi ^0\pi ^0)={\displaystyle \frac{p_{\pi \pi }}{pW_{pd}^2m_p^2E_\pi (\frac{2}{3}𝐩_{\pi \pi }\frac{1}{2}𝐩_d)^2}}{\displaystyle \frac{9}{2^{21}\pi ^{10}}}N\left\{|S_a|^2+|S_b|^2\right\}`$
$`\times |A|^2|a(m_{\pi \pi },W_{\pi N})|^2k_\pi ^{}p_\pi ^2\text{d}m_{\pi \pi }\text{d}\mathrm{\Omega }_{\text{He}},`$ (6)
where an isospin factor of $`\frac{9}{4}`$ has been included to account for the $`\pi ^0`$-exchange term in Fig. 1. The form factor combinations required are

$$S_a=S_{00}S_{20}\sqrt{2},S_b=S_{02}S_{22}\sqrt{2}.$$
(7)
In order to describe $`(\eta ,\omega ,\eta ^{},\varphi )`$ production, it was found necessary to multiply the analogous prediction by a normalisation factor $`N=2.4`$ .

Single pion production in pion-nucleon collisions has been measured in many charge states near threshold and parameterisations given for the total cross sections as functions of the beam energy . The charge dependence indicates that the cross section is dominated by $`I=0`$ pion pairs for $`Q<100`$ MeV. This is consistent with the smallness of the anisotropy in the angular distribution of the $`\pi ^+\pi ^{}`$ relative momentum for $`\pi ^{}p\pi ^+\pi ^{}n`$, which arises from $`s`$-$`p`$, and hence $`I=0/I=1`$ interference . Data on the $`\pi ^{}p\pi ^0\pi ^0n`$ reaction in the $`Q50100`$ MeV region show clear evidence for the suppression of events with low $`m_{\pi \pi }`$ . This is also seen for $`\pi ^{}p\pi ^+\pi ^{}n`$ but not $`\pi ^+p\pi ^+\pi ^+n`$, where only $`I=2`$ pion pairs are produced . In the dynamical model of the Valencia group , this shift towards higher $`m_{\pi \pi }`$ is due to an accidental cancellation between two contributions, one of which involves the double pion $`p`$-wave decay of the Roper resonance $`N^{}(1440)N\pi ^0\pi ^0`$.

Although the production of isospin-two $`\pi \pi `$ pairs is small , it has to be subtracted from the $`\pi ^0\pi ^0`$ data to get the $`I=0`$ rate required in Eq. (6) This subtraction is model-dependent and, for this purpose, we have used the predictions of the Valencia model , which describes reasonably well the shape of the experimental data . A global fit to their $`I=0`$ predictions, renormalised slightly to agree with the overall $`\pi ^0\pi ^0`$ amplitude analysis of Lowe and Burkhardt , gives

$$\frac{1}{64\pi ^3}|a(m_{\pi \pi },Q)|^2=(1.0920.0211Q+0.00015Q^2)$$
(8)
$$+(4.18+0.0075Q0.00098Q^2)x+(47.650.935Q+0.00743Q^2)x^2\mu \text{b/MeV}^2,$$

which is valid up to $`Q^{}=Q(1+m_\pi /m_{\text{He}})^2/(1+m_\pi /m_p)^2100`$ MeV. Given that $`x=m_{\pi \pi }/m_\pi 2`$, this illustrates the suppression of the matrix element at low $`m_{\pi \pi }`$, a feature which becomes even more pronounced at higher $`Q`$.

Our predictions for the $`pd^3`$He$`(\pi \pi )^0`$ total cross section divided by $`Q^2`$, obtained using the same value of $`N=2.4`$ which gave good agreement for heavy meson production, are to be found in Fig. 2. The steady increase with $`Q`$ is mainly a reflection of the energy dependence of the $`ppd\pi ^+`$ and $`\pi ^{}p\pi \pi n`$ amplitudes; the average form factor changes comparatively little. The solid curve passes close to the CELSISUS $`\pi ^0\pi ^0`$ point , but the broken one is significantly too low, indicating the presence of some $`I=1`$ $`\pi ^+\pi ^{}`$ production. The IUCF point was obtained at $`Q=0.67`$ MeV , and hence can be assumed to be purely $`s`$-wave, though it is heavily influenced by Coulomb effects. The comparison of our total cross section predictions with the MOMO $`\pi ^+\pi ^{}`$ point at $`Q=70`$ MeV  would suggest that it is mainly $`I=0`$ pairs which are being produced, though this is at variance with the strong dependence observed on the angle of the $`\pi ^+\pi ^{}`$ relative momentum.

The $`I=0`$ $`m_{\pi \pi }`$ distributions expected at the CELSIUS energy are illustrated in Fig. 3 and these demonstrate the shift to higher masses as compared to phase space, which is apparent in the $`\pi N\pi \pi N`$ input. These experimental data  have insufficient statistics to draw definitive conclusions on the shape of the spectrum. It should be noted that the CELSIUS integrated cross section points shown in Fig. 2 are mainly determined by the higher statistics of their inclusive measurement, which was carried out simultaneously . On the other hand, a low $`m_{\pi \pi }`$ suppression in the 70 MeV MOMO data is clear in their high statistics exclusive $`\pi ^+\pi ^{}`$ production results shown in Fig. 4.

The MOMO data  show a strong dependence upon the angle $`\theta _{\pi \pi p}`$ between the relative $`\pi \pi `$ momentum and that of the beam direction. Taken together with the suppression of events at low $`m_{\pi \pi }`$, this suggests the production of $`I=1`$, $`\mathrm{}=1`$ $`\pi \pi `$ pairs with spin projection $`m=\pm 1`$ along the beam direction. Such an interpretation is backed by the group’s preliminary data on the inclusive $`pd^3`$He$`X^0`$ reaction , which indicate a $`\pi ^0\pi ^0`$ production rate less than half of that predicted in Fig. 2.

The production of $`I=0`$, $`\mathrm{}=0`$ $`\pi \pi `$ pairs would give no dependence upon $`\theta _{\pi \pi p}`$, though an interference with an $`I=0`$, $`\mathrm{}=2`$ contribution could lead to such a variation. However, there is no sign of any effect of this kind in $`\pi ^{}p\pi ^+\pi ^{}n`$ . Since $`\pi ^+\pi ^{}`$ $`p`$-waves are so small in $`\pi ^{}p\pi ^+\pi ^{}n`$ near threshold , any simple extension of our model to include $`p`$-wave production cannot lead to $`\pi \pi `$ $`p`$-wave dominance.

For two-pion production near threshold, the intermediate pion in Fig. 1 gets closer to its mass shell when the 3He emerges along the direction of the initial proton beam and this increases the magnitude of the average form factor. The two-step model therefore predicts that, for low $`m_{\pi \pi }`$, the dipion should be produced preferentially in the backward hemisphere. This effect will, of course, disappear at high masses because the situation then approaches one of near-threshold kinematics.

We have shown that the gross features found in the production of isoscalar pion pairs in the $`pd^3`$He$`(\pi \pi )^0`$ reaction near threshold can be understood in terms of the creation of an intermediate virtual pion beam, which in turn produces a second meson. The only way that such a model could generate an ABC peak in the $`Q250`$ MeV region is if this were already present in the $`\pi ^{}p\pi ^0\pi ^0n`$ input. The parameterisation of the results of the Valencia model in Eq. (8) corresponds to a parabola in $`m_{\pi \pi }`$, whose minimum moves to higher values as $`Q`$ increases. This is due to the enhanced importance of the Roper contribution and may leave space at low masses for an ABC effect at higher $`Q`$. The question may soon be resolved, because data on $`\pi ^{}p\pi ^0\pi ^0n`$ at $`p_\pi =750`$ MeV/c are currently being analysed by the Crystal Ball collaboration .

Further theoretical work is needed to include a more detailed description of the $`\pi N\pi \pi N`$ input, though any improvement will, inevitably, be rather model-dependent.

It is a pleasure to thank R. Jahn for extensive discussions regarding the data of Ref.  and M. Andersson, P.E. Tégner, and K. Wilhelmsen Rolander for those pertaining to Ref. . We are greatly indebted to M. Vincente Vacas for giving us numerical values of the calculations reported in  and providing the $`I=0`$ projections. This work has been supported by the Royal Swedish Academy of Sciences within the framework of the European Science Exchange Programme. One of the authors (CW) is grateful to the The Svedberg Laboratory for its generous hospitality.