text:
# First Observation of the Exclusive Decays Λ_𝑐⁺→Λ𝜋⁺𝜋⁺𝜋⁻𝜋⁰ and Λ_𝑐⁺→Λ𝜔𝜋⁺

CLEO Collaboration

(October 18, 2002)

## Abstract

Using data collected by the CLEO III detector at CESR, we report on the first observation of the decays $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$. The branching fractions are measured relative to the $`pK^{}\pi ^+`$ mode; we find $`(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)/(\mathrm{\Lambda }_c^+pK^{}\pi ^+)=0.36\pm `$0.09(stat.)$`\pm `$0.09(syst.) and $`(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+)/(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$ = 0.24$`\pm `$0.06(stat.)$`\pm `$0.06(syst.). We also observe a small excess of events in $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$, $`\eta \pi ^+\pi ^{}\pi ^0`$ from which we find $`(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+)/(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$ = 0.41$`\pm `$0.17$`\pm `$0.10 with a corresponding limit of $`<`$0.65% at 90% confidence level. This rate is consistent with previous measurements of $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$, $`\eta \gamma \gamma `$. We find that the combined branching fractions for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ comprise (87$`\pm `$13)% of the total $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ rate. We set an upper limit on $`(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)/(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$ for non-$`\omega `$ and non-$`\eta `$ modes of 0.13 at 90% confidence level.

preprint:   CLNS 02/1800  CLEO 02-13

Over the last decade large data samples have been used to measure charm meson and charm baryon branching fractions to exclusive final states. Exclusive final states for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }(n)\pi `$ decays, with n=1,2,3 have been measured lampi\_modes , but modes with four (or more) pions have yet to be observed. The large rates for $`D^+\overline{K^0}\pi ^+\pi ^+\pi ^{}\pi ^0`$ ($`5.4_{1.4}^{+3.0})\%`$PDG2000  and $`D^0K^{}\pi ^+\pi ^+\pi ^{}\pi ^0`$ (($`4.0\pm 0.4)\%`$PDG2000  suggest non-negligible rates for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$. In this paper, we report the first measurement of the branching fractions for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$, $`\omega \pi ^+\pi ^{}\pi ^0`$. There is also some indication that a portion of the $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ occurs through $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$, $`\eta \pi ^+\pi ^{}\pi ^0`$. We find that the resonant $`\omega `$ and $`\eta `$ sub-modes nearly account for most (if not all) of the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ rate. Throughout this paper, charged conjugate states are implied unless otherwise specified.

The data for this analysis were collected using the CLEO III detector located at the Cornell Electron Storage Ring (CESR). The data sample consists of 2.7 fb-1 on or near the $`\mathrm{{\rm Y}}`$(4S), 1.4 fb-1 on or near the $`\mathrm{{\rm Y}}`$(3S), and 0.02 fb-1 on or near the $`\mathrm{{\rm Y}}`$(1S). The CLEO III detector is an upgrade of the CLEO II/II.V detector CLEOdetector  which includes a four-layer double-sided silicon vertex detector, a tracking drift chamber dr3  consisting of 16 inner axial layers and 31 outer small-angle stereo layers, and a ring-imaging Cerenkov detector (RICH) rich . The RICH is located outside the tracking chambers and uses LiF radiators to produce Cherenkov photons from incident charged tracks. The Cerenkov cone expands in a 20-cm thick $`\mathrm{N}_2`$-filled expansion volume and are then photo-converted in a gas mixture of Methane and TEA and detected by a multi-wire chamber with cathode pad readout. The electromagnetic calorimeter consists of 7800 CsI crystals which are used to measure the position and energy of photons and to aid in lepton identification. These detector components are immersed in a 1.5 Tesla solenoidal magnetic field that is uniform to within 0.2% over the tracking volume. Outside the solenoid there is a muon system covering 85% of the solid angle.

Signal candidates are formed using charged tracks and photons. We only use charged tracks which have a high quality track fit and have $`|\mathrm{cos}\theta |<0.925`$, where $`\theta `$ is the polar angle of the track with respect to the $`e^{}`$ direction. Charged particles (except those used in forming $`\mathrm{\Lambda }`$ candidates) are required to have a distance of closest approach (to the interaction point) in the $`r\varphi `$ plane of 5 mm and 5 cm along the $`\widehat{z}`$ axis (the direction of the $`e^{}`$ beam, apart from the small crossing angle of the $`e^+e^{}`$ beams). Photons used in forming $`\pi ^0`$ candidates must be isolated from charged tracks and have a shower shape that is consistent with the shape of an electromagnetic shower. The isolation criterion requires that no charged track points to any of the crystals which are used to make up the electromagnetic shower.

Particle identification is performed using information on the specific ionization ($`dE/dx`$) measured in the drift chamber and particle likelihoods obtained from the RICH. For $`dE/dx`$, we quantify the consistency of the measured $`dE/dx`$ with $`p`$, $`K`$ and $`\pi `$ hypotheses by defining a significance, $`S_i=((dE/dx)_i^{obs}dE/dx_i)/\sigma _i(dE/dx)`$, where $`(dE/dx)_i^{obs}`$ and $`dE/dx_i`$ are the measured ionization (corrected for path length) and expected (momentum-dependent) mean value of $`dE/dx`$ for particle type $`i`$ ($`i=p,K,\pi `$) respectively, and $`\sigma _i(dE/dx)`$ is the resolution on the measured $`dE/dx`$. For the RICH, we use the measured track momentum and a given particle hypothesis to search for Cerenkov photons whose measured positions are consistent (within 3 standard deviations) with expected positions. For each particle hypothesis, $`i`$, we compute a likelihood of the form:

$$L_i=\underset{j=1}{\overset{N_\gamma ^{obs}}{}}[G(\theta _j^{obs}|\theta _i^{exp})+B].$$

Here, $`G`$ is a Gaussian-like function which expresses the probability of observing the $`j^{th}`$ photon at an angle $`\theta _j^{obs}`$ with respect to the expected Cerenkov angle $`\theta _i^{exp}`$ for particle type $`i`$. The second term, $`B`$, is the background probability distribution which is a flat function in $`\theta _i^{exp}\theta _j^{obs}`$. Using the likelihoods, we compute differences in negative-log likelihoods between two particle-type hypotheses: proton-kaon, proton-pion, and kaon-pion (i.e., $`(K\pi )=2\mathrm{log}(L_K/L_\pi )`$). Defined in this manner, true kaons will peak at negative values of $``$ and true pions will peak at positive values of $``$.

Proton candidates are identified by requiring either $`|S_p|<3`$, $`|S_\pi |>3`$ and $`|S_K|>3`$, or at least three reconstructed Cerenkov photons with $`(p\pi )<5`$ and $`(pK)<0`$. Kaon candidates are only required to have either $`|S_K|<3`$ and $`|S_\pi |>3`$, or at least three reconstructed Cerenkov photons with $`(K\pi )<0`$. Because pions are the majority particle in hadronic events, only loose particle identification requirements are imposed. Pion candidates are required to have $`|S_\pi |<4`$ or at least three reconstructed Cerenkov photons with $`(K\pi )>5`$ and $`(p\pi )>5`$. By construction, the list of pion, kaon and proton candidates are not exclusive of each other, but in forming a given $`\mathrm{\Lambda }_c^+`$ candidate, all daughter tracks must be distinct from one another.

Candidate $`\mathrm{\Lambda }`$ baryons are constructed by pairing oppositely charged proton and pion candidates and performing a constrained vertex fit. To suppress combinatorial background we require $`\mathrm{\Lambda }`$ candidates to have a decay length of at least 1 cm. The invariant mass distribution of $`\mathrm{\Lambda }`$ candidates is shown in Fig. 1(a). The distribution is fit using two Gaussians (with the same mean) to describe the signal and a first order polynomial to describe the background. The relative fraction of events in the Gaussian components is consistent with simulation. The signal (sideband) region includes candidates which have a mass in the range from 1111-1121 MeV/$`c^2`$ (1100-1105 or 1130-1135 MeV/$`c^2`$). The $`\mathrm{\Lambda }`$ sidebands are used to check for a possible selection bias. Candidate $`\pi ^0`$ mesons are formed by combining pairs of photons and requiring that the pull is less than 2. The pull is defined as the difference in the $`\gamma \gamma `$ invariant mass from the known $`\pi ^0`$ mass, divided by the uncertainty in the $`\gamma \gamma `$ mass, where the uncertainty on the invariant mass is computed on a combination-by-combination basis. The pull distribution for the candidate $`\pi ^0`$’s in the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ sample are shown in Fig. 1(b). Superimposed is the result of a fit to the sum of a Gaussian and a linear background. The fit has a mean of zero and a width of 1.1 units which is consistent with simulation. To reduce combinatorial background, we require that the energy of each photon is larger than 50 MeV and that the normalized energy asymmetry of the two photons $`|E_1E_2|/(E_1+E_2)<0.8`$.

Candidates for the decay $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ are formed by combining $`\mathrm{\Lambda }`$ candidates with one neutral and three charged pions. The combinatoric background is suppressed by requiring $`\mathrm{\Lambda }_c^+`$ candidates to have a total momentum in excess of 3.5 GeV/$`c`$. Kinematically, this removes the contribution of $`\mathrm{\Lambda }_c^+`$ baryons from $`B`$ meson decay, leaving only the contribution from charm fragmentation. To improve the signal-to-noise ratio, we require $`\mathrm{\Lambda }`$, $`\pi ^0`$ and $`\pi ^\pm `$ candidates to have momenta in excess of 500 MeV/c, 250 MeV/c, and 100 MeV/c respectively.

Candidates for the $`\mathrm{\Lambda }_c^+pK^{}\pi ^+`$ decay (reference mode) are formed by combining proton candidates with charged kaon and charged pion candidates. As with the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ candidates, we require that the $`\mathrm{\Lambda }_c^+pK^{}\pi ^+`$ candidates have total momentum larger than 3.5 GeV/$`c`$. The charge of the kaon is required to be opposite to that of the proton. As the decay length of the $`\mathrm{\Lambda }_c^+`$ is too small to be observed in the CLEO detector, we require that all three tracks are consistent with coming from the $`e^+e^{}`$ interaction point. We also require that the proton have a total momentum larger than 500 MeV/$`c`$, while the kaon and pion are required to have momenta larger than 300 MeV/$`c`$.

Figure 2(a) shows the invariant mass distribution for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ candidates for data after all selection criteria. A clear excess of events is observed at the $`\mathrm{\Lambda }_c^+`$ mass. The distribution is fit to the sum of a Gaussian whose width is fixed to the expected width from simulation (8.2 MeV/$`c^2`$) and a linear background from which we extract a signal yield of 49.5$`\pm `$11.9 events at a mass of 2287$`\pm `$1.3 (stat.) MeV/$`c^2`$. Also shown in the figure (shaded histogram) are combinations where we pair the $`\pi ^+\pi ^+\pi ^{}\pi ^0`$ with a $`\overline{\mathrm{\Lambda }}`$. No $`\mathrm{\Lambda }_c^+`$ signal should be present in this mode, and the data show no significant excess in this mass range. The data have also been separated into $`\mathrm{\Lambda }_c^+`$ and $`\overline{\mathrm{\Lambda }_c^+}`$ in Fig. 2(b) and (c), respectively. They are fit using the same functional form as for the entire signal. The extracted yields are 21.3$`\pm `$8.1 $`\mathrm{\Lambda }_c^+`$ and 28.3$`\pm `$8.5 $`\overline{\mathrm{\Lambda }_c^+}`$ events. We therefore confirm that $`\mathrm{\Lambda }_c^+`$ and $`\overline{\mathrm{\Lambda }_c^+}`$ appear in approximately equal numbers.

We search for substructure in the $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ final state by forming $`\pi ^+\pi ^{}\pi ^0`$ combinations (2 per event) among the 4 pions using $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ combinations which have an invariant mass between 2270 and 2300 MeV/$`c^2`$. The $`\pi ^+\pi ^{}\pi ^0`$ invariant mass distributions are shown in Fig. 3 for the $`\mathrm{\Lambda }_c^+`$ signal region (2270-2300 MeV/$`c^2`$) and the $`\mathrm{\Lambda }_c^+`$ sideband regions (2245-2260 and 2320-2335 MeV/$`c^2`$). A strong peak is observed near $`M_\omega =782`$ MeV/$`c^2`$ and a smaller excess at $`M_\eta =547`$ MeV/$`c^2`$. The distribution is fit to the sum of two Gaussians, one for the $`\omega `$ and one for the $`\eta `$, and a second-order polynomial for the background. The signal shapes used in the fit are determined from simulation. The fitted masses are 552$`\pm `$2 MeV/$`c^2`$ and 778$`\pm `$3 MeV/$`c^2`$, consistent with the known masses of the $`\eta `$ and $`\omega `$ mesons. The probability that the enhancement at the $`\eta `$ mass is caused by a background fluctuation is about 4%. The $`\omega `$ signal is approximately 6 standard deviations above the expected fitted background. We also observe three narrow peaks between 825 and 900 MeV/$`c^2`$. Each of these peaks correspond to $``$1.5-2 standard deviations above the background prediction, depending on the assumed Gaussian width. Since there are no known particles at these masses, we assume they are statistical fluctuations.

The signal yields in the $`\omega `$ ($`\eta `$) modes are determined by requiring the $`\pi ^+\pi ^{}\pi ^0`$ invariant mass (either combination) to fall within $`\pm `$30 MeV/$`c^2`$ ($`\pm `$14 MeV/$`c^2`$) of their known masses and fitting for a signal in the $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ invariant mass spectrum. To ensure non-overlapping sub-samples, we exclude $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ candidates from the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ sample. The signal regions for the $`\omega `$ and $`\eta `$ correspond to about $`\pm `$3.5 times the Gaussian width obtained from simulations of these decay modes. The $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ invariant mass spectra for the $`\omega `$ and $`\eta `$ modes are shown in Fig. 4(a) and (b), respectively. The shaded regions show the corresponding distributions when we select combinations which fall into the $`\omega `$ sidebands (690-720 MeV/$`c^2`$ or 843-873 MeV/$`c^2`$) and $`\eta `$ sidebands (500-513 MeV/$`c^2`$ or 580-593 MeV/$`c^2`$) respectively. The distributions are fit to the sum of a Gaussian and a linear background. The widths of the Gaussians are constrained to 8.4 MeV/$`c^2`$ for the $`\omega `$ and 6.6 MeV/$`c^2`$ for the $`\eta `$. These widths are determined from Monte Carlo studies of these processes, which includes the Breit-Wigner width of the $`\omega `$ and detector resolution effects. The fitted yields are 31.9$`\pm `$7.1 $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ candidates and 10.8$`\pm `$4.1 $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ candidates. Due to the limited statistics in the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ mode, we also quote an upper limit of 16.8 events at the 90% confidence level. Thus these two modes account for 87% of the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ candidates.

The contribution from the $`\omega `$ and $`\eta `$ are combined and shown in Fig. 5(a). The spectrum is fit to the sum of a Gaussian and a linear background shape. The width is constrained to the event-weighted average of the widths obtained from the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ simulations (7.9 MeV/$`c^2`$). The number of fitted candidates is 42.5$`\pm `$8.2, which agrees with the sum from the individual fits. Events which do not have a $`\pi ^+\pi ^{}\pi ^0`$ with invariant mass in the $`\eta `$ or $`\omega `$ region are shown in Fig. 5(b). We fit the spectrum assuming a Gaussian signal shape and a first order polynomial for the background shape. The mean of the Gaussian is fixed to the world-average value (2285 MeV/$`c^2`$) and the width is constrained to the value found from a simulation of non-resonant $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ (9.8 MeV/$`c^2`$). The fitted yield of $`5.0_{7.7}^{+8.4}`$ events is consistent with the expected difference of 7.0 events as well as zero. This corresponds to an upper limit of 17.7 events at the 90% confidence level.

The fitted $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ signals are corrected for possible contributions from non-resonant production by subtracting off the expected non-resonant contribution, as determined from simulation. The number of non-resonant events is taken to be 7.0 events as discussed above. This number must be corrected for the contribution from $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$, because some of the signal lies outside the $`\omega `$ mass window. The $`\eta `$ width is negligible compared to the mass resolution, and therefore no such correction is necessary for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$. From simulation, we find that (3$`\pm `$2)% of the $`\omega `$ signal events fall outside the $`\omega `$ mass window. We therefore estimate that 6.0 of the 49.5 events are from sources other than $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$. We assume that these 6 events are from non-resonant sources, and using simulation we find that 16% (4%) are expected to have a $`\pi ^+\pi ^{}\pi ^0`$ combination in the $`\omega `$ ($`\eta `$) mass window. The corrected yields are therefore 31.0$`\pm `$7.1 and 10.6$`\pm `$4.1 for the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ modes. The 3% loss of $`\omega `$ events due to the mass window definition are accounted for in the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ efficiency. We take the systematic uncertainty on the non-resonant subtraction to be 100% of the correction.

We have investigated the possibility that the $`\omega \pi `$ system may come from higher mass resonance decays, such as $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\rho ^+,\rho ^+\omega \pi ^+`$, which may contribute to this final state due to its large width. We compared the $`\omega \pi ^+`$ invariant mass distribution in data to a non-resonant $`\omega \pi ^+`$ simulation and the expectation from $`\rho ^+\omega \pi ^+`$. No peaks are observed in data, and a Kolmogorov-Smirnov test of the shapes give a consistency of 80% with respect to a pure non-resonant $`\mathrm{\Lambda }\omega \pi ^+`$ final state and a 10% consistency with pure $`\mathrm{\Lambda }\rho ^+`$. Thus, the data are more compatible with a pure non-resonant $`\omega \pi ^+`$ contribution, although a $`\rho ^+`$ contribution cannot be ruled out.

The relative branching fraction $`(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)/(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$ is:

$$\frac{(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)}{(\mathrm{\Lambda }_c^+pK^{}\pi ^+)}=\frac{N(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)}{N(\mathrm{\Lambda }_c^+pK^{}\pi ^+)}\frac{ϵ(\mathrm{\Lambda }_c^+pK^{}\pi ^+)}{ϵ(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)},$$

where $`N(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)`$ and $`N(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$ are the number of reconstructed $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ and $`\mathrm{\Lambda }_c^+pK^{}\pi ^+`$ signal events and $`ϵ(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)`$ and $`ϵ(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$ are the relative efficiencies for reconstructing these modes respectively.

The efficiencies for the $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ and the $`pK\pi `$ final states are obtained using a Monte Carlo simulation of these decays (see Table 1). Events are simulated using the Lund/Jetset model lund  for the production mechanism ($`e^+e^{}q\overline{q}`$) and are allowed to decay according to the known lifetime and decay branching fractions PDG2000 . The response of the detector to the particles in the generated events is simulated using a GEANT-based simulation geant  of the detector. The fully simulated events are then passed through the same analysis chain as data. For the $`pK\pi `$ final state, our simulation includes the resonant sub-modes $`pK^{}`$ and $`\mathrm{\Delta }^{++}K`$ as well as a non-resonant $`pK\pi `$ contribution PDG2000 . The efficiency for the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ final state is taken to be the branching fraction weighted efficiency for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$, $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$, and non-resonant $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$. Table 1 summarizes the inputs for the computation of the branching fractions.

The branching fractions for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$, $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ relative to $`\mathrm{\Lambda }_c^+pK^{}\pi ^+`$ are measured to be:

$$\frac{(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0)}{(\mathrm{\Lambda }_c^+pK^{}\pi ^+)}=0.36\pm 0.09\pm 0.09$$

$$\frac{(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+,\omega \pi ^+\pi ^{}\pi ^0)}{(\mathrm{\Lambda }_c^+pK^{}\pi ^+)}=0.21\pm 0.05\pm 0.05$$

$$\frac{(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+,\eta \pi ^+\pi ^{}\pi ^0)}{(\mathrm{\Lambda }_c^+pK^{}\pi ^+)}=0.093\pm 0.038\pm 0.022(<0.15\mathrm{at}90\%\mathrm{CL}),$$

where the first uncertainty is statistical and the second is systematic. The systematic uncertainties are discussed below.

These measurements are the first observations of $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ modes, which are (or almost) entirely through $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ and $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$. This is the first observation of the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ decay mode. The $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ decay has been previously observed in $`\eta \gamma \gamma `$ decays lametapi , but this is the first indication in the $`\eta \pi ^+\pi ^{}\pi ^0`$ mode. Our measurement compares well with those results which yielded $`(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+)/(\mathrm{\Lambda }_c^+pK^{}\pi ^+)=0.35\pm 0.05`$(stat.)$`\pm `$0.06(syst.).

As a systematic check of our procedure, we have also measured the relative branching fraction for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+`$ and find $`(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+)/(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$=0.17$`\pm `$0.02(stat.). This is consistent with the current world average of 0.18$`\pm `$0.03 PDG2000 .

Significant systematic uncertainties arise from sources which affect the $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ mode but not the $`pK^{}\pi ^+`$ mode (or vice versa). Systematic uncertainties include estimation of reconstruction efficiencies of tracks, $`\mathrm{\Lambda }`$ and $`\pi ^0`$ reconstruction, particle identification selection, background estimation and signal estimation.

Studies of data and simulation show that the charged particle tracking simulation reproduces data within 2%. A 2% systematic uncertainty per track introduces a 6% relative systematic uncertainty. The systematic uncertainty in the $`\pi ^0`$ reconstruction efficiency is estimated by comparing the differential yield for $`\mathrm{{\rm Y}}(4S)\pi ^0X`$ in CLEO III data to CLEO II data. The difference is found to be 8% for our selection criteria. When combined with the 5% systematic uncertainty on the CLEO II $`\pi ^0`$ efficiency, we assign a 10% uncertainty to the $`\pi ^0`$ reconstruction efficiency. The systematic uncertainty in the $`\mathrm{\Lambda }`$ reconstruction efficiency is estimated by comparing our measured value for $`(\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+)/(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$ \[0.17$`\pm `$0.02(stat.)\] to the world average value (0.18$`\pm `$0.03). While these results are consistent with one another, we take the uncertainty to be the the quadrature sum of the difference in the two results (6%) and the relative uncertainty in the world average value (17%). The uncertainty in the $`\mathrm{\Lambda }`$ reconstruction efficiency is therefore estimated to be 18%. The systematic uncertainty in the particle identification was investigated using a sample of protons from $`\mathrm{\Lambda }p\pi ^{}`$ and kaons from $`D^+D^0(K^{}\pi ^+)\pi ^+`$. The studies indicate that the systematic uncertainty is less than 3% for protons and less than 1% for kaons. The uncertainty in the background is estimated by using alternate background shapes derived from (a) the $`\overline{\mathrm{\Lambda }}\pi ^+\pi ^+\pi ^{}\pi ^0`$ invariant mass distribution in data and (b) the $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ invariant mass distribution using the $`\mathrm{\Lambda }`$ sidebands. We also varied the number of bins used in the fit to the background; these sources lead to an uncertainty of 6% on the signal yield.

Additional uncertainty in the signal estimation may arise from an imperfect understanding of the signal resolution and efficiency. The resolution uncertainty is estimated by varying the widths by one standard deviation ($`\pm `$ 1 MeV/$`c^2`$) and refitting the data. The uncertainty in the signal yield is found to be 5%. For the efficiency uncertainty, we consider one source which only affects the full $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ branching fraction measurement, and one which affects the resonant $`\omega `$ and $`\eta `$ contributions. For the full $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ sample, we use a branching fraction weighted efficiency. To evaluate the systematic uncertainty, we shift the number of $`\omega `$ events by $`\pm `$7.6 events and the corresponding number of $`\eta `$ events by $``$7.6 events and compute the corresponding changes in the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ branching fraction. We find a shift of $`\pm `$3% in the $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ branching fraction. For the resonant modes, we also include a 3% systematic which corresponds to a 100% relative uncertainty on the subtraction for possible non-resonant contributions which have a $`\pi ^+\pi ^{}\pi ^0`$ invariant mass in either the $`\omega `$ or $`\eta `$ mass region. A compilation of the sources of systematic uncertainty is presented in Table 2. The total systematic uncertainty is estimated to be 13%.

These measurements indicate that ($`87\pm 13`$)% of the four-pion state is feed-down through resonant $`\omega `$ and $`\eta `$ production. The branching fractions (after dividing through by the $`(\omega \pi ^+\pi ^{}\pi ^0)`$ and $`(\eta \pi ^+\pi ^{}\pi ^0)`$) are shown in Table 3. The common systematic uncertainty in the branching fraction for $`\mathrm{\Lambda }_c^+pK^{}\pi ^+`$ has not been included in the last column in Table 3. The $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\eta \pi ^+`$ mode has been measured previously in the $`\eta \gamma \gamma `$ mode lametapi , and the result presented here is consistent with those findings. We find that the branching fraction for $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\omega \pi ^+`$ is (0.24$`\pm `$0.06$`\pm `$0.06) of $`(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$. For the $`\eta `$ mode, we find a rate relative to $`(\mathrm{\Lambda }_c^+pK^{}\pi ^+)`$ of (0.41$`\pm `$0.17$`\pm `$0.10) with an upper limit of 0.65 at 90% CL. The non-$`\omega ,\eta `$ signal constitutes $`5.0_{7.7}^{+8.4}`$ events which corresponds to an upper limit of 17.7 events at 90% confidence level. This corresponds to an upper limit on the branching fraction relative to $`(\mathrm{\Lambda }_c^+pK^{}\pi ^+`$) of 0.13 at 90% confidence level.

The fact that the $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ rates are competitive with $`\mathrm{\Lambda }\pi ^+`$ ((0.9$`\pm `$0.3)%), $`\mathrm{\Lambda }\pi ^+\pi ^0`$ ((3.6$`\pm `$1.3)%) and $`\mathrm{\Lambda }\pi ^+\pi ^{}\pi ^+`$ ((3.3$`\pm `$1.0)%) is not surprising, as charm mesons exhibit similar branching fractions for $`DKn(\pi )`$, where $`n=14`$. In fact, it is interesting to note that the exclusive rate for $`D^0K^{}\omega \pi `$ ((2.7$`\pm `$0.5)%) constitutes about 68% of the rate for $`D^0K^{}\pi ^+\pi ^+\pi ^{}\pi ^0`$ ((4.0$`\pm `$0.4)%). This is quite comparable to what we find here for the ratio of the exclusive $`\mathrm{\Lambda }\omega \pi ^+`$, $`\omega \pi ^+\pi ^{}\pi ^0`$ rate to the total $`\mathrm{\Lambda }\pi ^+\pi ^+\pi ^{}\pi ^0`$ rate (59%).

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. M. Selen thanks the PFF program of the NSF and the Research Corporation, and A.H. Mahmood thanks the Texas Advanced Research Program. This work was supported by the National Science Foundation, and the U.S. Department of Energy.


format: text

text:
# Emergence of a complex and stable network in a model ecosystem with extinction and mutation

## 1 Introduction

The relationship between the complexity and stability of an ecosystem has been one of the most fascinating topics in theoretical biology for decades (Pimm 1991). In the 1950s and 1960s the proposition that highly complex communities are more stable than simple ones was widely supported (MacArther, 1955; Elton, 1958). However, this view was challenged by theorists in the 1970s, who discussed the stability of a community of species interacting randomly (Gardner & Ashby, 1970; May, 1972, 1974; Roberts, 1974; Gilpin & Case, 1976; Taylor, 1988a). The most significant result was given by May (1972), who considered a large-dimensional ecological equation with an $`n`$-dimensional random interaction matrix $`\{a_{ij}\}`$ whose diagonal elements are $`1`$ and whose off-diagonal elements are assigned as Gaussian random numbers with mean $`0`$ and variance $`\alpha ^2`$ (with probability $`C`$) or $`0`$ (with probability $`1C`$). He found a kind of phase transition from stability to instability in which the equilibrium solution corresponding to the coexistence of all species becomes unstable if $`\alpha >(nC)^{1/2}`$. He concluded that an ecological system cannot be stable if it is complex. After May’s work and the clear conclusions he reached, a number of mathematical biologists attempted to explain the discrepancy between the observed complexity of ecosystems in nature and the results of these mathematical studies. Their works can be classified into two groups with regarded to the methodology of assembling a complex and stable ecosystem.

The first pioneering idea was formulated by Tregonning & Roberts (1979). They prepared a $`50\times 50`$ random interspecies interaction matrix and checked the feasibility of the equilibrium solution of a population equation. Here, an equilibrium solution is termed feasible if all species have positive populations. In the case that the solution under consideration was not feasible, they removed that species whose equilibrium population was most negative. Repeatedly applying this procedure, they acquired a stable equilibrium of positive population with considerably higher diversity than that predicted by May’s theorem \[see also Roberts & Tregonning (1980)\]. Let us call their approach the eliminating approach. Although this approach can produce a complex and stable ecosystem, the resulting diversity is in fact much less than the initial diversity ($`=50`$). Moreover, the biological significance of their modeling was not very clear because negative population is never achieved in nature and because the origin of the random interactions was not discussed. \[See also Tokita & Yasutomi (1999) for a paleontological explanation of the formation of such random interspecies interactions.\] Furthermore, according to the analysis of extinction dynamics (Tokita & Yasutomi, 1999), which is a type of population dynamics similar to the elimination process, the resulting diversity is independent of the initial diversity. This result suggests that when the interactions are assigned randomly, a system with arbitrarily large diversity cannot emerge through such an elimination process. However it is still worthy of attention that Tregonning & Roberts (1979) reported that the most often-observed interspecies interaction of the survived subsystem was a prey-predator relationship, because it was recently demonstrated that coexistence of a number of species can be achieved by an ecological model with such prey-predator relationships even if the interactions are random (Chawanya & Tokita, 2002).

Other authors have used a second method, which is known as community assembly (Pimm 1991). In this method, a model community is assembled through a sequence of invasions of species. This type of invasion is believed to reflect the natural tendency of species to arrive sequentially rather than simultaneously. The method further can be classified into two subgroups. The first trailblazing one is represented by that of Post & Pimm (1983). They prepare a pool of species ex ante and then add them sequentially. The properties of the newly introduced species are then determined through a sampling process that is completely independent of the properties of the species existing in the community at the time of introduction. A considerable number of studies have followed this line (Robinson & Valentine, 1979; Yodzis, 1988; Taylor, 1988b; Drake, 1990; Case, 1990; Law & Morton, 1996; Happel & Stadler, 1998). These models are realistic in the case that we consider systems in which new species come from outside of the existing community. We thus refer to this as the invasion approach.

The second subgroup concerns a longer-term evolutionary process (Colwell & Winkler, 1984; Ginzburg et al., 1988; Happel & Stadler, 1998). In such models, a new species is introduced as a mutant of an existing species. More precisely, parameters characterizing a new species are generated by making appropriate changes to the parameters of one of the existing species, rather than by sampling the parameter values independently from a separate predetermined ensemble. We refer to this method as the mutation approach.

Although these works of community assembly succeed in producing stable communities whose numbers of species are larger than that predicted by May’s theorem, the coexistence of hundreds of species has not been demonstrated. On the other hand, it is worthy of mention that in a recent mutation approach (Drossel et al., 2001), considerably higher diversity has been realized than that by other mutation approaches. They focused on influence of prey-predator interactions to the evolution of food web structure. Consequently, mutualisms and other interspecies relationships were not considered, and “connectance” (the percentage of non-zero values of interactions) was low. In contrast, in the present study, we consider not only the prey-predator type but also high connectance of interactions of mutualisms, competitions and parasitisms.

In the present model, species are defined by interactions matrices, following the preceding studies mentioned above. Such definition is, however, still a controversial paradigm because interaction coefficients are not real biological traits, like body size, but rather summaries of these traits filtered through some interaction rules. We neverthless adopt the interaction matrices because one of the motivations here is to clarify the relationship between complexity and stability of large ecosystems, i.e. the paradox of ecology, raised by May (1972). For that purpose, analyses of multidimensional trait space would be essential, which is, in general, not a simple task for the real biological traits. Here we only cite some approaches on community assembly that does deal in real traits (Brown & Vincent, 1992; Sasaki 1997; Geritz, et al., 1998).

In this paper, we observe the behavior of a system based on the replicator equation in which a species is eliminated when its population becomes too small, and mutants or invaders are added periodically. We consider the cases of three different rules governing the generation of new species. We call these invasive, global and local. The invasive rule can be compared with the model of the invasion approach, and the global rule with the global mutation approach. The local rule is introduced here for the first time.

We will show that the local rule exclusively allows for the assembly of a very large and complex ecosystem which contains hundreds of species strongly interacting with each other. In addition, we observe that the diversity and mutualism increase together. The key point in our model is the introduction of a mutant species whose relationships with the dominant species (with large populations) are virtually the same as those of their parent species, but who can, nevertheless, interact differently to minor species (with small populations) as well as to those species which will come in the future. The fitness of such a mutant is almost identical to that of its parent when it is introduced, and hence it often behaves as a temporally neutral mutant.

We will show that the minor mutants play a key role in the emergence of a complex ecosystem. This idea can be discussed in connection with neutral molecular evolution (Kimura, 1983), which stresses the importance of the accumulation of nondirectional mutations. Our theory is closely related to the ecological neutral evolution theory for the abundance and diversity of species in tropical rain forests and coral reefs proposed by Hubbell (1979, 1997, 2001). The size distribution of the relative abundance of species obtained in our model is in good agreement with a well-known distribution obtained from detailed observations of real ecosystems.

We will also discuss the resistance to invasion of assembled ecosystems. We find that an ecosystem assembled under the local rule has perfect resistance to invasion if we keep it free from invaders which are assembled according to the invasive rules for a sufficiently long time. Contrastingly, an ecosystem assembled under the global rule cannot resist invaders.

This paper is organized as follows. In the next section we observe the dynamics of the number of surviving species in our model. In the third section, we discuss the resistance to invasion of the assembled ecosystems. The final section is devoted to discussion of the relationship between the generalized Lotka-Volterra equations and the replicator equations.

## 2 The replicator equation with extinction and mutation

### 2.1 Model

Here, we investigate a model based on the system of ordinary differential equations

$$\frac{\mathrm{d}x_i(t)}{\mathrm{d}t}=x_i(t)\left(\underset{j=1}{\overset{N_I}{}}a_{ij}x_j(t)\underset{j,k=1}{\overset{N_I}{}}a_{jk}x_j(t)x_k(t)\right),$$
(1)
called the  replicator equations (RE)(Hofbauer & Sigmund, 1998), on a $`N_I`$-dimensional simplex:

$$\underset{i=1}{\overset{N_I}{}}x_i(t)=1(0x_i(t)1).$$
(2)
The variable $`x_i`$ denotes the population (or relative abundance) of the species $`i`$, and $`N_I`$ denotes the initial number of species, that is, the initial value of the diversity. The $`(i,j)`$-th element of the matrix $`A=(a_{ij})`$ determines the effect of species $`j`$ on the growth rate of species $`i`$.

It is known that the $`N_I`$-dimensional RE is topologically equivalent to an $`N_I1`$ dimensional generalized Lotka-Volterra equation (GLVE) (Hofbauer & Sigmund, 1998). We use the RE here, simply because they are more tractable for a numerical simulation of dynamics of a large ecosystem with high diversity and complex interactions. A comparison between the GLVE and RE in the context of the dynamics of diversity is given in the last section.

It should be noted that the RE may possess heteroclinic orbit even in low dimensions ($`N_I4`$) (May & Leonard 1975; Chawanya 1995 & 1997). When a heteroclinic orbit approaches a saddle, where some species are ‘extinct’, their population take extremely small values but they never exactly reach zero, because the orbit is bound in the interior of the simplex (2). In the vicinity of such a saddle, these population have such small values that they cause underflow in naive numerical calculations. However, if we continue a calculation in such a region using the more sophisticated calculation method of Chawanya (1995, 1997), we often observe some of these species start to revive and cause the orbit to eventually leave for another saddle, in particular for a system with high diversity and complex interactions. This transition among saddles continues cyclically or chaotically. The exponential approach of a population to zero and its subsequent revival to order $`O(1)`$ plays an important role in heteroclinic orbits. However, in the real world, such small population cannot be realized.

Considering this problem, we introduce a parameter $`\delta (1)`$, the extinction threshold, into the dynamics described by (1) and (2) : At each discrete time step, the population $`x_k`$ is set to zero if this quantity becomes less than $`\delta `$. The population of the surviving species {$`x_i`$}$`(ik)`$ are then renormalized to satisfy $`_{ik}x_i=1`$. The diversity $`N(t)`$ decreases through the above described process. We refer to this as the extinction dynamics (ED) model (Tokita & Yasutomi, 1999). With $`\delta `$, a stochastic effect on extinctions is introduced into the fully deterministic RE, and hence, $`\delta `$ plays a role analogous to that of random genetic drift. The introduction of $`\delta `$ constitutes a finite size effect on the total population into RE, because $`\delta `$ represents a minimum unit of reproduction for each species, and its reciprocal $`1/\delta `$ represents to the permissible population size. That is, the system explicitly possesses an energy constraint and is free from the problem of population explosion often observed in the simulation of a large-scale GLVE pointed out by Taylor (1988b).

In contrast to the typical ED simulation with a large $`N_I`$, in the present simulation, we start from only one species (i.e. $`N_I=1`$), following Ginzburg et al.(1988). The initial intraspecies interactions is set to $`v(<0)`$, which is one of the parameters in the present model. A simulation proceeds by first computing the ED for a fixed number ($`T`$) of time steps. This is followed by the introduction of a new species. This procedure is then repeated many times. We set $`T=900`$ in our simulations. Hereafter, we refer to $`T`$ time steps as one period, and we use $`\tau `$ to count the number of periods. At the end of each period, we select one parent species, $`j`$, among the existing species \[$`1jN(\tau )`$\] with a sampling probability proportional to the population $`x_j`$. Then we produce a new species $`k[=N(\tau )+1)]`$ from the parent species $`j`$, copying $`\{a_{ji}\}`$ and $`\{a_{ij}\}`$ to $`\{a_{ki}\}`$ and $`\{a_{ik}\}`$ for all $`i`$, with some variations applied in the manner described below. We then set the new species’ population as $`x_k=\zeta `$ ($`\zeta >\delta `$) and normalize all $`x_i`$ in order to maintain the condition $`x_i=1`$. Thus completing the introduction of the new species, we then return to the ED, and the entire process is repeated.

As we discussed above, we consider three different rules for the variations of the interactions $`\{a_{ki}\}`$ and $`\{a_{ik}\}`$ for the new species: invasive, global and local. In the case of the invasive rule, we add a new species whose interaction coefficients are assigned randomly and are independent of all existing species (including its parent). In the case of the global rule, we produce a mutant whose interaction coefficients are all only slightly different from those of its parent species. In the case of the local rule, we introduce a mutant that is identical to its parent species, apart from its interaction with one particular surviving species. While the invasive (Robinson & Valentine, 1979; Yodzis, 1988; Happel & Stadler, 1998) and global (Ginzburg et al., 1988; Happel & Stadler, 1998; Drossel et al., 2001) rules have been employed in previous studies, the local rule is examined for the first time in the present study. In terms of theoretical population genetics, the invasive rule corresponds to the “house-of-cards” approximation (Turelli, 1983). The details of these rules are as follows.

* Invasive: We assign the values of $`\{a_{ki}\}`$ and $`\{a_{ik}\}`$ for all $`i(k)`$ as Gaussian random numbers (mean 0 and variance 1). This implies that the mutant species is independent of its parent species.
* global: Quantities $`\{\upsilon _{ki}\}`$ and $`\{\upsilon _{ik}\}`$ representing Gaussian random noise with mean $`0`$ and variance $`1`$ are independently added to $`\{a_{ki}\}`$ and $`\{a_{ik}\}`$ for all $`i`$ according to

$$(\begin{array}{c}a_{ki}\gamma a_{ki}+(1\gamma )\upsilon _{ki}\\ a_{ik}\gamma a_{ik}+(1\gamma )\upsilon _{ik}\end{array}$$
(3)
for all $`i`$, where $`\gamma `$ ($`0<\gamma <1`$) denotes the strength of the noise. In our simulations, this value was set to $`0.9`$.
* local: We select a species $`l`$ whose population is non-zero and change only the relationship between the mutant species and species $`l(k)`$. Thus, the effects of a mutation are concentrated on only two elements, $`a_{kl}`$ and $`a_{lk}`$. We assigned the values of these two elements using Gaussian random numbers with mean $`0`$ and variance $`1`$.

In all cases, the intraspecies interaction $`a_{ii}`$ takes the value $`v`$ for each species $`i`$. Finally, we introduce a threshold $`\xi `$ ($`\xi >\zeta >\delta `$), and we consider a species to be unobservable if its population $`x_i`$ is smaller than $`\xi `$. Since $`\xi `$ is larger than $`\zeta `$, a mutant species is not observable until its population increases and becomes greater than $`\xi `$. The function $`S(\tau )`$ denotes the number of the observable species existing at the end of the $`\tau `$th period. We regard $`S(\tau )`$ as the diversity in this model. It should be noted that, in general, there appear much more unobservable species than observable species although the population of the latter is much larger than that of the former.

### 2.2 Results

Figure 1 (a) displays the development of $`\overline{S}(\tau )`$ for each of the three rules, where $`\overline{S}(\tau )`$ is defined as $`S(\tau )`$ averaged over independent simulations with different series of random numbers. We set the initial conditions as $`x_1=1`$ \[that is, $`N(0)=S(0)=1`$\] and $`a_{11}=v`$. This figure shows that diversity grows most under local rule. The diversity $`\overline{S}(\tau )`$ under the local rule increases to approximately $`50`$ at $`\tau =10,000`$, and as can be inferred from this plot, it continues to increase beyond this value as more mutants are added. We averaged the diversity over $`50`$ simulations and found that this averaged value obeys $`\overline{S}(\tau )\tau ^{0.38}`$ \[Figure 1(b)\], which should be compared with the square root law $`S(\tau )\sqrt{\tau }`$ in the coinfection model (May & Nowak, 1995).

Figure 2 (a) displays the development of the average fitness under the three rules. As in the case of the diversity, it is seen here that only the local rule allows for the average fitness to continue to increase. Since the average fitness $`\overline{f}=_{j,k=1}^{N(\tau )}a_{jk}x_j(\tau )x_k(\tau )`$ is a kind of average of the interaction coefficients, we can regard this quantity as an index of mutualism. The results here therefore indicate that the diversity and the level of the mutualism increase together. It is interesting that the average fitness does not display power-law behavior, but rather increases much more slowly: $`\overline{f}c\mathrm{log}(\tau )(c0.25)`$ (Figure 2(b)). This is analogous to the slow increase of the diversity observed in a GLVE model of super-infection in a host-parasite system (May & Nowak, 1994).

Figure 3(a) plots the development of $`S(\tau )`$, where the largest number of species observed for this case is $`196`$ at $`\tau =144,700`$. In the plots of Figures 3 (a) and (b), only one sample for the development of each $`S(\tau )`$ and $`\overline{f}`$ is used. By contrast, the plots in Figure 1 represent averages over $`50`$ samples. From Figure 3(a), we wee that the diversity does not increase monotonically, and repeated avalanches of extinction followed by rapid recoveries are observed. Figure 3(b) gives an enlarged view of the time evolution plotted in Figure 3(a) from $`\tau =7,500`$ to $`7,900`$ in order to clarify the relationship between $`S(\tau )`$ and $`\overline{f}`$. Here we observe a clear negative correlation; that is, the average fitness drastically increases during an avalanche of extinction. This behavior is typically observed in extinction dynamics (Tokita & Yasutomi, 1999). However, the average fitness does not increase during the time that the diversity increases. This implies that extinction plays a more important role than speciation in the increase of the average fitness. This is because large extinction events swiftly exterminate less mutualistic species and produce a more mutualistic network of species. The refined network then prepares for a rapid recovery of the diversity in the aftermath of the extinction.

In the case that $`S(\tau )`$ is defined with just one sample, its fluctuation suggests simultaneous occurence of extinction and speciation of multiple species. On the other hand, $`\overline{f}(\tau )`$ develops relatively smoothly, because the definition of $`\overline{f}`$ is a kind of average of the interactions $`\{a_{ij}\}`$ over $`\{x_i\}`$, and $`x_i`$ at the time that the $`i`$th species appears or extincts is very small ($`\delta `$ for extinction and $`\zeta `$ for speciation).

Figure 4 plots the dependence of $`\overline{S}`$ on the ratio of the value of the diagonal elements, $`v`$, and the variance of the Gaussian random numbers ($`=1`$) that are used to assign the interaction coefficients of the mutant in the case of the local rule. We see that $`\overline{S}`$ depends only weakly on $`v`$. We should note that May (1972) showed that the internal equilibrium point becomes unstable when this value becomes less than $`1`$ if the system is assembled randomly at one time, rather than through the periodic introduction of mutants, as done in our study. The strength of the interactions between species affects the stability of the ecosystem if it is assembled at one time, but it does not do so when the ecosystem gradually develops for a long time on the evolutionary time scale.

Why does the diversity $`S(\tau )`$ increase so drastically only in the case of the local rule? To answer this question, let us consider Figure 6, which gives an example of our numerical integration. This figure shows that a few species enter the existing ecosystem more rapidly than exponentially and virtually simultaneously. This kind of rapid population growth results from the fact that the relationship between these mutant species is more mutualistic than those between the dominant species. When the relationship between mutants $`p`$ and $`q`$ is more mutualistic than those between the dominant species, and this relationship therefore makes these mutants “fitter” than the average species, their population increase. The growth of the population of species $`p`$ thus increases the fitness of species $`q`$, thereby enhancing its growth rate, and hence increasing the fitness of $`p`$. These unobservable species constitute a new mutualistic unobservable subsystem and eventually merge into the observable existing ecosystem. The implication of this behavior is that the formation of a mutualistic network is more effective as a invasion strategy than competition or even exploitation under the pressure of natural selection, because a mutualistic invasion has the abovementioned positive feedback effect, while the exploitation involves a negative feedback (as the more one species preys on another, the less abundant this prey becomes).

In order to form such a subsystem, mutants must wait for the arrival of partners. In order to avoid extinction while waiting, such a mutant must have virtually the same fitness as its parent species. More precisely, they must have the same interaction coefficients with respect to the dominant species as their parents. If this is not the case, they will independently invade the existing ecosystem, replacing their parent species, or they will simply be repulsed. Of course, the invasive and global rules rarely produce such “sleeping” mutants, because in general in these cases all interactions of a given mutant will be different from those of their parents. Figure 6 gives an example of our numerical integration for the case of the global rule. Here a new species invades independently, replacing an existing species, and we can observe no faster-than-exponential population growth. Only the local rule can produce such mutants, and this is the reason why only it can yield a remarkable increase in the diversity.

Since a subsystem of the type described above has a higher degree of mutualism than the observable system, we expect that the ecosystem will become more mutualistic as the diversity increases. Since the average fitness $`\overline{f}=_{j,k=1}^{N(t)}a_{jk}x_j(t)x_k(t)`$ is a kind of average of the interaction coefficients, we can regard this average as an index of mutualism. As seen in Figure 1, it continues to increase only under the local rule. Figure 8 exhibits the distribution of $`\{a_{ij}\}`$ for an ecosystem that has emerged under the local rule. It should be noted that these values are ‘naturally selected’ from Gaussian random numbers with mean $`0`$ and variance $`1`$. As seen in the figure, the average of the distribution shifts to the positive area. This suggests that the relationships between species tend to become more and more mutualistic through the struggle for survival.

Let us consider in detail the process of speciation by introducing a $`2S(\tau )`$-dimensional trait space, where each species $`i`$ is characterized by the point $`(a_{i1},a_{i2},\mathrm{},a_{i,S(\tau )},a_{1i},a_{2i},\mathrm{},a_{S(\tau ),i})`$. In order to visualize the distribution of the species in this high-dimensional space, we project this point onto a point in a two-dimensional space:

$$\stackrel{}{b}^{(i)}=(b_1^{(i)},b_2^{(i)})\frac{1}{S(\tau )}\left(\sqrt{\underset{k}{}a_{ki}^2}_,\sqrt{\underset{j}{}a_{ij}^2}\right)_.$$
(4)
If the ecological status of two species $`k`$ and $`l`$ are similar, $`\stackrel{}{b}^{(k)}`$ and $`\stackrel{}{b}^{(l)}`$ are close to each other.

Figure 8 displays the sets $`\{\stackrel{}{b}^{(i)}\}`$ for each species with the interaction matrices resulting from invasive, global and local rules. For the sake of contrast, we also display the points $`\stackrel{}{r}^{(i)}`$, defined in the same way as the $`\stackrel{}{b}^{(i)}`$, for a $`50\times 50`$ matrix whose off-diagonal elements are given as Gaussian random numbers with mean $`0`$ and variance $`1`$. We can see that the sets $`\{\stackrel{}{b}^{(i)}\}`$ for these three interaction matrices differ greatly from $`\{\stackrel{}{r}^{(i)}\}`$. This suggests that it is quite inappropriate to use a random matrix to approximate an interaction matrix of a real ecosystem. Among these three sets $`\{\stackrel{}{b}^{(i)}\}`$, that obtained from the local rule in particular is shifted in relation to $`\{\stackrel{}{r}^{(i)}\}`$ more mutualistic region in the trait space. This is suggested also by Figure 8. Moreover, several groups of species are observed in this case, indicating that the speciation processes under the local rule do not result in a simple random drift or a random diffusion in the trait space; that is, the groups of species that appear in this case are fixed by the pressure of natural selection.

Figure 9 displays the relationship between species in the case of the local rule. It is observed that this ecosystem starts from a single species, and the diversity increases through speciation. Figure 9(a) displays the history of the speciation. The black lines indicate a genealogy of species existing at $`\tau =50,000`$, and the gray lines at $`\tau =100,000`$ from their ancestors at $`\tau =0`$ and $`\tau =50,000`$, respectively. Lines representing the extinct species before $`\tau =50,000`$ and $`\tau =100,000`$ are not plotted, although there are in fact more branches for the extinct species. At $`\tau =50,000`$ there exist 58 species, and they are separated into two large groups which branches at a very early stage.

Figure 9(b) displays the points $`\stackrel{}{b}^{(i)}`$ for each species. The set $`\{\stackrel{}{b}^{(i)}\}`$ at $`\tau =50,000`$ is plotted with black dots and those at $`\tau =100,000`$ as circles. We observe two large groups in this figure: Groups A and B contain 39 and 19 species at $`\tau =50,000`$, respectively. Those species on the same branch of the genealogy are plotted near each other in the 2-dimensional space. We also observe several subgroups in each group: Group A contains four subgroups and group B contains two. These groups also reflect the structure of genealogy, as those species in the same subgroup appear on the same sub-branch of the genealogy. If we think of each dot and circle as a ‘species’, then each subgroup can be thought of as a ‘genus’ and an entire group is a ‘family’. The two families A and B still exist at $`\tau =100,000`$, shifting toward the upper right. All species at $`\tau =50,000`$ are extinct, and those descendants of only three species survive at $`\tau =100,000`$. Each arrow denotes a radiation from one of these three to a group or a subgroup at $`\tau =100,000`$. Descendants of one species of family A radiated into $`32`$ species. This family separated into two genera, and even these genera formed subsubgroups. There are $`13`$ species of descendants that emerged from two species of family B: one of these genera contains $`4`$ species and the other 9 species. The latter genus is separated into $`4`$ sub-genera.

Finally, we note the interesting correspondence of a result obtained in the present study to a well-known ecological law of population abundance: As seen from Figure 10, the dominance-diversity distribution is $`S`$-shaped, and its middle part displays a power law behavior. Such behavior is observed in many species-rich communities (Hubbell, 1997; Magurran, 1988). This suggests that the present model shares some characteristics with populations in nature. Let us here note that our data were not obtained through any averaging, but from only one simulation. This is significant because in order to obtain such a distribution using one sample, it is required that the sample produce a sufficiently large number of species. This has been realized for the first time in the present study.

In the present study, we have mentioned “invasion type” mutants with respect to species-specific interaction by the local rule. However we can also think of ”mutant type” mutants with respect to species-specific interaction, which have the parameter deviated from the parent’s value by a random quantity. We have executed similar simulations by such a mutative local rule and found a similar diversification like the local rule. Under the mutative local rule, the diversity $`S(\tau )`$ increases little slower but is fitted by a power function like under the local rule.

From the above discussion, we conclude that an ecological system which is composed of many hierarchically structured species strongly interacting with each other can emerge through an evolutionary speciation process by the local rule.

## 3 Resistance to Invasion

### 3.1 Method

As mentioned in the first section, Elton (1958) suggested that complex communities are more stable than simple ones. In this context, as pointed out by Case (1990, 1991), Elton’s definition of the ‘stability’ is not the asymptotic stability of the assembled matrix but the resistance to invasion of communities. We thus can understand his hypothesis as an assertion that species-rich ecosystems are more resistant to invasion by exotics than are species-poor ones.

Case (1990, 1991) constructed models to check Elton’s hypothesis from this point of view. He randomly assembled ecosystems and selected stable and feasible ones. Then he tested their resistance to invasion by adding new species generated by a rule very similar to our invasive rule. He concluded that a more diverse and more strongly connected ecosystem has greater resistance to invasion.

### 3.2 Results

We carried out a similar test on ecosystems formed under the global and local rules and compared the resistance to invasion of these ecosystems. There is an important difference between these two rules and the invasive rule. Under the global and local rules, new species are mutants of the existing species. However, under the invasive rule these are independent of the existing species, invading from without. In our tests, we added mutants generated by the local and global rules for only the first $`R`$ periods, and then we let invaders generated by the invasive rule attack the ecosystem for 200 periods. $`R`$ denotes a kind of maturity index of the ecosystem.

The results of our experiment are depicted in Figure 11, which shows the dependence of the resulting diversity $`S(\tau )`$ on $`R`$. As seen there, when $`R`$ is small, the diversity is drastically reduced under either rule. However, when $`R`$ is larger than $`200`$, ecosystems developed under the local rule maintain their diversity, resisting invasion. Contrastingly, those developed under the global rule exhibit a reduction in their diversity for any value of $`R`$, as a result of destructive invasions. This result suggests that an ecosystem developed under the local rule has strong resistance to invasion if it is sufficiently mature.

The difference between the local and global rules results from the difference in their levels of mutualisms. As we have seen in the previous section, an ecosystem evolving under the local rule is more mutualistic than one evolving under the global rule. We should note that the average fitness $`\overline{f}`$ here again turns out to play an important role as an index indicating the resistance to invasion . If an ecosystem has an average fitness $`\overline{f}`$, the fitness $`f_I`$ of an invader must be larger than $`\overline{f}`$ in order for the invader to succeed. Since $`f_I`$ can be approximated by a Gaussian random number with mean $`0`$ and variance $`1`$, for an ecosystem with sufficiently large diversity, the probability for a successful invasion is given by the error function,

$$\mathrm{Prob}.(f_I>\overline{f})=\mathrm{Erfc}(\overline{f})=\frac{1}{\sqrt{2\pi }}_{\overline{f}}^{\mathrm{}}\mathrm{exp}\left(\frac{t^2}{2}\right)dt,$$
(5)
which decreases exponentially as $`\overline{f}`$ increases. From Figure 1(b), it is possible to estimate the probability (5) by using the value of $`\overline{f}`$ for $`R200`$ at which the difference between the local and global becomes larger in Figure 11. The average fitness for the local rule, $`\overline{f}(\tau =200)1`$, gives $`\mathrm{Erfc}(1)0.16`$, while that for the global rule, $`\overline{f}(\tau =200)0.25`$, gives $`\mathrm{Erfc}(0.25)0.4`$. Figure 1(b) indicates that $`\overline{f}`$ is a monotonically increasing function of $`R`$ for the local rule, and hence we expect that a random invasion becomes almost impossible against a fully-matured ecosystem (i.e., one of large $`R`$) under the local rule. That is, only temporally neutral mutants can influence or create extinction events in an existing ecosystem evolved under the local rule.

## 4 Discussion

### Generalized Lotka-Volterra equations vs. replicator equations

The generalized Lotka-Volterra equations (GLVE) are more commonly used than the replicator equations (RE) in the field of mathematical biology. However, it is known that an $`N`$-dimensional RE is mathematically identical to the ($`N1`$)-dimensional GLVE (Hofbauer & Sigmund, 1998)

$$\frac{\mathrm{d}y_i}{\mathrm{d}t}=y_i(r_i+\underset{j=1}{\overset{N1}{}}b_{ij}y_j)(i=1,2,\mathrm{},N1),$$
(6)
where the $`i`$th species’ population $`y_i`$ is defined in terms of the population in the RE (1), $`x_i`$ by

$$y_i=\frac{x_i}{x_N}.$$
(7)
The value of $`y_i`$ can take any non-negative real number. The interactions $`\{b_{ij}\}`$ and the intrinsic growth (or decay) rate $`r_i`$ here are defined in terms of interactions in the RE by

$`b_{ij}`$ $`=`$ $`a_{ij}a_{Nj}^{}{}_{,}{}^{}`$ (8)
$`r_i`$ $`=`$ $`a_{iN}a_{NN}=a_{iN}+v_.`$ (9)
Any orbit of the GLVE can therefore be mapped onto an equivalent orbit of the RE, and any behavior exhibited by these two systems can be studied using either, at least in principle.

Here we explain why we nevertheless study the RE rather than the GLVE. As described in the first section, Taylor (1988b) pointed out that most of the dynamics resulted in explosions of populations as proportion of species with $`r_i>0`$ increased or proportion of positive $`b_{ij}`$ increased. Moreover, the populations of species with low fitness often take extremely small values that cause underflow in naive numerical computational scheme. Consequently, numerical analysis of the GLVE often encounters both divergence and underflow, in particular for a system with high diversity and complex interactions. On the other hand, in the RE, population explosions are avoided by definition, because any orbit of the population $`\{x_i\}`$ is bound in the simplex $`_{i=1}^Nx_i=1`$. Moreover, in the present model, the problem of underflow is also avoided because the heteroclinic orbits, the cause of underflow, are excluded by the introduction of the extinction threshold (Tokita & Yasutomi 1999). These are the reasons we use the RE rather than the GLVE.

Let us note that the meaning of the interaction coefficients $`a_{ij}`$ in the RE is different from that of the $`b_{ij}`$ in the GLVE: even if the $`a_{ij}`$ and $`a_{ji}`$ take only positive values (mutualistic relationships), the $`b_{ij}`$ and $`b_{ji}`$ can take negative values (competitive or exploitative relationships) depending on the values of $`a_{ij}`$ and $`a_{Nj}`$ in Eq. (8). As the growth rate of the species $`i`$ is defined by the fitness $`_ja_{ij}x_j`$ deducted by its average $`_{j,k}a_{jk}x_jx_k`$ in the definition of the replicator equations (1), the ecological significance of the developed ecosystem in the present study should be discussed by the matrix $`B=(b_{ij})`$ in the GLVE than $`A=(a_{ij})`$ in the RE. We therefore study the nature of matrices $`B`$ transformed from $`A`$ by Eq. (8). Here we should note that there are $`N`$ ways of transformation from $`A`$ to $`B`$ because of the arbitrary option of the species $`N`$ in Eq. (7). To characterize a $`L\times L`$ matrix $`M=\{m_{ij}\}`$, here we define “antisymmetricity” as

$$D_{AS}(M)\frac{(m_{ij}m_{ji})^2_M}{4m_{ij}^2_M}_.$$
(10)
The bracket $`(\mathrm{})_M\frac{1}{L(L1)}_{i=1}^L_{ji}^L(\mathrm{})`$ denotes the average over all off-diagonal elements of $`M`$. For the RE with interaction $`A`$ at $`\tau `$, $`L`$ becomes $`S(\tau )`$, and $`L=S(\tau )1`$ for $`B`$ of the corresponding GLVE. In general, the antisymmetricity $`D_{AS}(M)`$ becomes $`0`$, $`1/2`$ and $`1`$ when $`M`$ is symmetric ($`m_{ij}=m_{ji}`$), randomly asymmetric ($`m_{ij}m_{ji}`$) and antisymmetric ($`m_{ij}=m_{ji}`$), respectively. Figure 12 plots the distribution of $`D_{AS}(B)`$ at $`\tau =4,000,170,000\text{ and }700,000`$, where $`S(4,000)=45,S(170,000)=125\text{ and }S(700,000)=216`$. The value of the antisymmetry $`D_{AS}(A)`$ of the RE at $`\tau =4,000,170,000`$ is also indicated. From the figure, the interaction matrix $`A`$ of the RE turns out to approach the symmetric point ($`D_{AS}=0`$; mutualism). On the other hand, the matrix $`B`$ tends to approach the intermediate region ($`D_{AS}0.3`$) between symmetric and asymmetric. This implies that the system evolved under the local rule has complex ecological interspecies interactions including mutualism, competition, predation and parasitism.

From the figure 8 and 12, we conclude that the frequency of the interspecies interactions of the emerged system can be approximated by a symmetric Gauss distribution with mean $`m(>0)`$ and variance $`1`$

$`P(a_{ij})`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\mathrm{exp}({\displaystyle \frac{(a_{ij}m)^2}{2}}),(\text{for}ij\text{and}a_{ij}=a_{ji})`$ (11)
$`a_{ii}`$ $`=`$ $`v(<0)`$ (12)
in the infinite limit of the number of species ($`N\mathrm{}`$), where $`m>0`$ is the level of mutualism. This type of the replicator equation, in general, has a number of saturated fixed points, which is proved by the symmetry of the interaction (Hofbauer & Sigmund, 1998).

Finally, let us see how {$`r_i`$} and {$`b_{ij}`$} in the corresponding GLVE mutate in the local rule. First, values $`r_k`$, {$`b_{ik}`$} and {$`b_{ki}`$} for all $`i`$ of a new species $`k`$ are copied from $`r_j`$, {$`b_{ij}`$} and {$`b_{ji}`$} of a parent species $`j`$, respectively. Second, by the definition of the local rule and the transformation (8) and (9), $`r_k`$ or $`b_{lk}`$ (and $`b_{kl}`$) for arbitrary $`l`$ is replaced by a random quantity. Note that $`r_k`$ is replaced by a random number with mean $`v`$ and variance $`1`$ while $`b_{lk}`$ ($`b_{kl}`$) with mean $`0`$ and variance $`1`$. From the transformation (8) and the distribution (11) and (12) of $`a_{ij}`$, the distribution of $`b_{ij}`$ becomes a Gauss distribution with mean $`0`$ and, therefore, is not mutualistic. Minor mutualist species in the RE, therefore, are not necessarily mutualists in the GLVE, and do not necessarily lead population explosion. On the other hand, $`r_i`$ in Eq. (9) is biased by $`v(>0)`$ and the mean value of $`r_i`$ increases over time because mean value of $`a_{iN}`$ increases as indicated by Fig. 8. Accordingly, in terms of the GLVE, proportion of producers increases as time goes on in the local rule. As pointed out by Taylor (1988b), this may trigger population explosion. Another chance for explosion is implied in the transformation (7) when the $`N`$-th “base” species goes extinct ($`x_N(t)0`$), which breaks down the simulation of the corresponding GLVE. We stress here that even if the extinction of one “base” species causes the explosion in the GLVE, the evolution in the RE proceeds successfully and there are still other equivalent GLVE systems by transformations by other $`N1`$ “base” species as $`y_i=x_i/x_M`$ ($`M=1,2,\mathrm{},N1`$). The RE, therefore, traces a possible path of evolution avoiding the breakdown of the simulation which may occur in the GLVE. This methodological advantage of the RE is related to the circumstance that the dimension (the degree of freedom) of the RE is one higher than the GLVE and the total population is conserved (Eq. (2)). This suggests that such a conservative quantity, e.g. resource limitations, would be essential in the modeling of multispecies evolution using the GLVE.

The authors would like to thank Yoh Iwasa, Takashi Ikegami and Tsuyoshi Chawanya for their helpful comments on the manuscript and for their encouragement. This work was partially supported by the Japan Society for the Promotion of Science, a Grant-in-Aid from the Ministry of Education, Science, Sports and Culture of Japan (No. 12740235, 13831007 and 14740232), and the Suntory Foundation. Some of the numerical calculations were carried out on a Fujitsu VPP500/40 at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo and on the Machikaneyama PC Cluster System in the Condensed Matter Theory Group, Graduate School of Science, Osaka University (URL: http://wwwacty.phys.sci.osaka-u.ac.jp/~mhill).