text:
# 1 Introduction

## 1 Introduction

The study of the process $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}`$ has recently become possible since LEP now operates at center-of-mass energies above the threshold for on-shell $`\mathrm{Z}`$ boson pair production. In the Standard Model, the process $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}`$ occurs via the NC2 diagrams  shown in Figure 1. The $`\mathrm{Z}`$-pair cross section depends on properties of the $`\mathrm{Z}`$ boson ($`m_\mathrm{Z}`$, $`\mathrm{\Gamma }_\mathrm{Z}`$ and the vector and axial vector coupling of the $`\mathrm{Z}`$ to electrons, $`g_\mathrm{V}^\mathrm{e}`$ and $`g_\mathrm{A}^\mathrm{e}`$) that have been measured with great precision at the $`\mathrm{Z}`$ resonance . The expected $`\mathrm{Z}`$-pair cross section increases from about 0.25 pb at $`\sqrt{s}=183`$ GeV to about 1.0 pb at $`\sqrt{s}=200`$ GeV, but remains more than an order of magnitude smaller than $`\mathrm{W}`$-pair production. In contrast to $`\mathrm{W}`$-pair production, where tree level $`\mathrm{WW}\gamma `$ and $`\mathrm{WW}\mathrm{Z}`$ couplings are important, no $`\mathrm{Z}`$$`\mathrm{Z}`$$`\mathrm{Z}`$ and $`\mathrm{Z}\mathrm{Z}\gamma `$ couplings are expected in the Standard Model. However, physics beyond the Standard Model could lead to effective couplings  which could then be observed as deviations in the measured $`\mathrm{Z}`$-pair cross section from the Standard Model prediction. Such deviations have been proposed in the context of Higgs doublet models  and in low scale gravity theories . In this paper we report on measurements of the NC2 $`\mathrm{Z}`$-pair cross section, including the extrapolation to final states with one or both $`\mathrm{Z}`$ bosons off-shell. These measurements, along with the angular distribution of the observed events, are then used to extract limits on possible $`\mathrm{Z}`$$`\mathrm{Z}`$$`\mathrm{Z}`$ and $`\mathrm{Z}\mathrm{Z}\gamma `$ couplings.

In Section 2 we describe the data sets used and the Monte Carlo simulation of signal and background. In Section 3 we describe the selection of the processes $`\mathrm{Z}\mathrm{Z}\mathrm{}^+\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$, $`\mathrm{Z}\mathrm{Z}\mathrm{}^+\mathrm{}^{}\nu \overline{\nu }`$, $`\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$, $`\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }`$, and $`\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$, where $`\mathrm{}^+\mathrm{}^{}`$ denotes a charged lepton pair of opposite charge and $`\mathrm{q}\overline{\mathrm{q}}`$ any of the five lightest quark-antiquark pairs. We also describe analyses of $`\mathrm{Z}\mathrm{Z}\mathrm{b}\overline{\mathrm{b}}\mathrm{}^+\mathrm{}^{}`$, $`\mathrm{Z}\mathrm{Z}\mathrm{b}\overline{\mathrm{b}}\nu \overline{\nu }`$ and $`\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\mathrm{b}\overline{\mathrm{b}}`$ which use b-tagging methods similar to those used in the OPAL Higgs search . The use of b-tagging improves the separation of the $`\mathrm{Z}`$-pair signal from background and allows us to check the $`\mathrm{b}\overline{\mathrm{b}}`$ content of our $`\mathrm{Z}`$-pair sample for consistency with the Standard Model. The description of the individual selections is followed by a discussion of possible systematic errors (Section 3.6). In Section 4 the selected events are used to measure the $`\mathrm{Z}`$-pair cross section. Then the cross section and angular distribution are compared with the Standard Model predictions and limits on anomalous neutral current triple gauge couplings are derived.

## 2 Data analysis and Monte Carlo

The OPAL detector111OPAL uses a right-handed coordinate system in which the $`z`$ axis is along the electron beam direction and the $`x`$ axis is horizontal. The polar angle, $`\theta `$, is measured with respect to the $`z`$ axis and the azimuthal angle, $`\varphi `$, with respect to the $`x`$ axis., trigger and data acquisition system are described fully elsewhere . Our analyses use approximately $`55\mathrm{pb}^1`$ of data collected at center-of-mass energies between 181–184 GeV and approximately $`178\mathrm{pb}^1`$ collected at center-of-mass energies near 189 GeV. The corresponding luminosity-weighted mean center-of-mass energies are $`182.62\pm 0.05`$ GeV and $`188.63\pm 0.04`$ GeV . The luminosity was measured using small-angle Bhabha scattering events recorded in the silicon-tungsten luminometer  and the theoretical calculation of Reference . The overall error on the luminosity measurement amounts to less than 0.5% and contributes negligibly to our cross-section measurement error.

Selection efficiencies and backgrounds were calculated using Monte Carlo simulations. All events were passed through a simulation  of the OPAL detector and processed as for data. We define the $`\mathrm{Z}`$$`\mathrm{Z}`$ cross section as the contribution to the total four-fermion cross section from the NC2 $`\mathrm{Z}`$-pair diagrams shown in Figure 1. All signal efficiencies given in this paper are with respect to these $`\mathrm{Z}`$-pair processes. Contributions from all other four-fermion final states, including interference with NC2 diagrams, are considered as background. For studies of the signal efficiency we have used grc4f , YFSZZ  and PYTHIA .

Backgrounds are simulated using several different generators. PYTHIA is used to simulate two-fermion final states such as $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}^{}(n\gamma )\mathrm{q}\overline{\mathrm{q}}(n\gamma )`$ and $`\mathrm{e}^+\mathrm{e}^{}\gamma ^{}(n\gamma )\mathrm{q}\overline{\mathrm{q}}(n\gamma )`$, where $`(n\gamma )`$ indicates the generation of one or more initial state photons. HERWIG  and KK2f  are used as checks for these final states. These two-fermion generators include gluon radiation from the quarks which produce $`\mathrm{q}\overline{\mathrm{q}}\mathrm{g}`$, $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ and $`\mathrm{q}\overline{\mathrm{q}}\mathrm{gg}`$ final states. The grc4f generator, with the contribution exclusively due to NC2 diagrams removed, is used to simulate other four-fermion background. KORALW  and EXCALIBUR  are used as checks of the four-fermion background. Multiperipheral (“two-photon”) processes with hadronic final states are simulated by combining events from PYTHIA, for events without electrons scattered into the detector, and HERWIG , for events with electrons scattered into the detector. For the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ final state, TWOGEN  is used to simulate two-photon events with both the electron and positron scattered into the detector. The Vermaseren  generator is used to simulate multiperipheral production of the final states $`\mathrm{e}^+\mathrm{e}^{}\mathrm{}^+\mathrm{}^{}`$.

To avoid background from four-fermion final states mediated by $`\mathrm{Z}\gamma ^{}`$, our selections were optimized to select simulated events with masses, $`m_1`$ and $`m_2`$, that satisfy $`m_1+m_2>170\mathrm{GeV}`$ and $`|m_1m_2|<20\mathrm{GeV}`$. At 189 GeV (183 GeV) more than 90% (80%) of the events produced via the NC2 diagrams are contained in this mass region. Events from the NC2 diagrams dominate in this mass region except for final states containing electron pairs. Backgrounds in these samples from two-photon and electroweak Compton scattering ($`\mathrm{e}\gamma \mathrm{e}\mathrm{Z}`$) processes  are reduced by using electrons detected in the electromagnetic calorimeters with $`|\mathrm{cos}\theta _\mathrm{e}|<0.985`$, where $`\theta _\mathrm{e}`$ is the polar angle of the electron.

## 3 Event selection

In the following subsections we describe event selections which exploit every decay mode of the $`\mathrm{Z}`$ boson. Our selections cover all $`\mathrm{Z}`$$`\mathrm{Z}`$ final states except $`\nu \overline{\nu }\nu \overline{\nu }`$ and $`\tau ^+\tau ^{}\nu \overline{\nu }`$. In hadronic final states, the energy and direction of the jets are determined using reconstructed tracks and calorimeter clusters using the correction for double counting described in Reference . In the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$, $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ and $`\mathrm{q}\overline{\mathrm{q}}\mathrm{b}\overline{\mathrm{b}}`$ analyses four-constraint (4C) and five-constraint (5C) kinematic fits are used. The 4C fit imposes energy and momentum conservation. In the 5C fit the added constraint requires the masses of the two candidate $`\mathrm{Z}`$ bosons to be equal to one another. For final states with one or more $`\mathrm{Z}`$ bosons decaying to tau pairs, the energy and total momentum of the tau leptons are obtained by leaving the reconstructed direction of the four fermions fixed and scaling the energy and momentum of each of the fermions to obtain energy and momentum conservation. The scaled values of the tau momentum and energy are then used in the subsequent steps of the analysis. In the $`\mathrm{q}\overline{\mathrm{q}}\tau ^+\tau ^{}`$ and $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ final states, subsequent kinematic fits are effectively 2C and 3C fits.

### 3.1 Selection of $`𝐙𝐙\mathbf{}\mathbf{}^\mathbf{+}\mathbf{}^{\mathbf{}}\mathbf{}^\mathbf{+}\mathbf{}^{\mathbf{}}`$ events

$`\mathrm{Z}`$-pair events decaying to final states with four charged leptons ($`\mathrm{}^+\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$) produce low multiplicity events with a clear topological signature that is exploited to maximize the selection efficiency. The $`\mathrm{}^+\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$ analysis begins by selecting low multiplicity events (less than 13 tracks or clusters) with visible energy of at least $`0.2\sqrt{s}`$ and at least one track with momentum of 5 GeV or more. Using a cone algorithm, the events are required to have exactly four cones of $`15^{}`$ half angle containing between 1 and 3 tracks. Cones of opposite charge are paired222Two-track cones are assigned the charge of the most energetic track if the momentum of one track exceeds that of the other by a factor of 4. Events with a cone which fails this requirement are rejected. to form $`\mathrm{Z}`$ boson candidates.

Lepton identification is only used to classify events as background or to reduce the number of cone combinations considered by preventing the matching of identified electrons with identified muons. Electrons are identified on the basis of energy deposition in the electromagnetic calorimeter, track curvature and specific ionization in the tracking chambers. Muons are identified using the association between tracks and hits in the hadron calorimeter and muon chambers.

To reduce background from two-photon events with a single scattered electron detected, we eliminate events with forward going electrons (backward going positrons) with the cut $`\mathrm{cos}\theta _\mathrm{e}^{}<0.85`$ ($`\mathrm{cos}\theta _{\mathrm{e}^+}>0.85`$). Here $`\theta _\mathrm{e}^{}`$ ( $`\theta _{\mathrm{e}^+}`$) is the angle of the electron (positron) with respect to the incoming electron beam. Background from partially reconstructed $`\mathrm{q}\overline{\mathrm{q}}(n\gamma )`$ events and two-photon events is reduced by requiring that most of the energy is not concentrated in a single cone, $`E_{\mathrm{vis}}E_{\mathrm{cone}}^{\mathrm{max}}>0.2\sqrt{s}`$. Here $`E_{\mathrm{vis}}`$ is the total visible energy of the event and $`E_{\mathrm{cone}}^{\mathrm{max}}`$ is the energy contained in the most energetic cone.

The invariant masses of the lepton pairs are calculated in three different ways which are motivated by the possibility of having zero, one or two tau pairs in the event. The events are classified according to the number of tau pairs in the event. (i) Events with $`E_{\mathrm{vis}}>0.9\sqrt{s}`$ are treated as $`\mathrm{e}^+\mathrm{e}^{}\mathrm{e}^+\mathrm{e}^{}`$, $`\mathrm{e}^+\mathrm{e}^{}\mu ^+\mu ^{}`$ or $`\mu ^+\mu ^{}\mu ^+\mu ^{}`$ events. We also treat all events with $`|\mathrm{cos}\theta _{\mathrm{miss}}|>0.98`$ ($`\theta _{\mathrm{miss}}`$ is the polar angle associated with the missing momentum in the event) as $`\mathrm{e}^+\mathrm{e}^{}\mathrm{e}^+\mathrm{e}^{}`$, $`\mathrm{e}^+\mathrm{e}^{}\mu ^+\mu ^{}`$ or $`\mu ^+\mu ^{}\mu ^+\mu ^{}`$ events to maintain efficiency for $`\mathrm{Z}`$-pairs with initial state radiation. In these events there are no missing neutrinos and the mass of each pair of lepton cones is evaluated. (ii) Events failing (i) with a cone-pair combination that has energy exceeding $`0.9m_\mathrm{Z}`$ are tried as an $`\mathrm{e}^+\mathrm{e}^{}\tau ^+\tau ^{}`$ or $`\mu ^+\mu ^{}\tau ^+\tau ^{}`$ final state. The mass of the tau-pair system is calculated from the recoil mass of the presumed electron or muon pair. (iii) Any remaining combinations are treated as $`\tau ^+\tau ^{}\tau ^+\tau ^{}`$ final states. The momenta of the tau leptons are determined with the scaling procedure described in the introduction to Section 3 and the invariant masses of the cone pairs are evaluated using the scaled momenta. In any event with more than one valid combination, each combination is tested using the invariant mass cuts listed below.

To reduce the combinatorial background, combinations with pair masses closest to $`m_\mathrm{Z}`$ are selected. In events with one or more combination satisfying $`|m_\mathrm{Z}m_{\mathrm{}\mathrm{}}|>0.1m_\mathrm{Z}`$ the cone-pair combination with the smallest value of $`(m_\mathrm{Z}m_{\mathrm{}\mathrm{}})^2+(m_\mathrm{Z}m_{\mathrm{}^{}\mathrm{}^{}})^2`$ is selected for further analysis. In the other combinations, the combination with the smallest value of $`|m_\mathrm{Z}m_{\mathrm{}\mathrm{}}|`$ or $`|m_\mathrm{Z}m_{\mathrm{}^{}\mathrm{}^{}}|`$ is selected. The final event sample is then chosen with the requirement $`m_{\mathrm{}\mathrm{}}+m_{\mathrm{}^{}\mathrm{}^{}}>160`$ GeV and $`|m_{\mathrm{}\mathrm{}}m_{\mathrm{}^{}\mathrm{}^{}}|<40`$ GeV. The signal detection efficiency, averaged over all $`\mathrm{}^+\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$ final states is given in Table 1 (line a). The efficiency for individual final states range from 30% for $`\tau ^+\tau ^{}\tau ^+\tau ^{}`$ to more than 70% for $`\mu ^+\mu ^{}\mu ^+\mu ^{}`$. The invariant masses of all cone pairs passing one of the selections are shown in Figure 2a. One candidate is found in the 183 GeV data and one candidate is found in the 189 GeV data.

### 3.2 Selection of $`𝐙𝐙\mathbf{}\mathbf{}^\mathbf{+}\mathbf{}^{\mathbf{}}𝝂\overline{𝝂}`$ events

The selection of the $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }`$ and $`\mu ^+\mu ^{}\nu \overline{\nu }`$ final states is based on the OPAL selection of $`\mathrm{W}`$ pairs decaying to leptons . The mass and momentum of the $`\mathrm{Z}`$ boson decaying to $`\nu \overline{\nu }`$ are calculated using the beam energy constraint and the visible decay of the other $`\mathrm{Z}`$ boson to a charged lepton pair. A likelihood selection based on the visible and recoil masses as well as the polar angle of the leptons, is then used to separate signal from background.

The $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }`$ selection starts with OPAL $`\mathrm{W}`$-pair candidates where both charged leptons are classified as electrons. Each event is then divided into two hemispheres using the thrust axis. The highest momentum charged (leading) track is selected from each hemisphere. The sum of the charges of these two tracks is required to be zero. The determination of the visible mass, $`m_{\mathrm{vis}}`$, and the recoil mass, $`m_{\mathrm{recoil}}`$, is based on the energy as measured in the electromagnetic calorimeter and the direction of the leading tracks.

Three variables were chosen for the likelihood selection: $`Q\mathrm{cos}\theta `$, where $`\theta `$ is the angle of the highest momentum charged track and $`Q`$ is its charge, the normalized sum of visible and recoil masses $`(m_{\mathrm{vis}}+m_{\mathrm{recoil}})/\sqrt{s}`$ and the difference of visible and recoil masses, $`m_{\mathrm{vis}}m_{\mathrm{recoil}}`$. The performance of the likelihood is improved with the following preselection: $`25\mathrm{GeV}<\mathrm{m}_{\mathrm{vis}}\mathrm{m}_{\mathrm{recoil}}<15\mathrm{GeV}`$ and $`(\mathrm{m}_{\mathrm{vis}}+\mathrm{m}_{\mathrm{recoil}})/\sqrt{\mathrm{s}}>0.90`$. One event with $`_{\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }}>0.60`$ is selected (see Table 1 (line b) and Figure 3a).

The $`\mu ^+\mu ^{}\nu \overline{\nu }`$ selection starts with the OPAL $`\mathrm{W}`$-pair candidates where both charged leptons are classified as muons. The selection procedure is the same as for the $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }`$ final states except that $`m_{\mathrm{vis}}`$, $`m_{\mathrm{recoil}}`$, and $`E_{\mathrm{vis}}`$ are calculated from the the momentum of the reconstructed tracks of the $`\mathrm{Z}`$ boson decaying to muon pairs. The likelihood preselections $`25\mathrm{GeV}<\mathrm{m}_{\mathrm{vis}}\mathrm{m}_{\mathrm{recoil}}<25\mathrm{GeV}`$ and $`(\mathrm{m}_{\mathrm{vis}}+\mathrm{m}_{\mathrm{recoil}})/\sqrt{\mathrm{s}}>0.90`$ are applied. Two events with $`_{\mu ^+\mu ^{}\nu \overline{\nu }}>0.60`$ are selected (see Table 1 (line c) and Figure 3b).

### 3.3 Selection of $`𝐙𝐙\mathbf{}𝐪\overline{𝐪}\mathbf{}^\mathbf{+}\mathbf{}^{\mathbf{}}`$ events

The lepton pairs in the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ and $`\mathrm{q}\overline{\mathrm{q}}\mu ^+\mu ^{}`$ final states have a distinctive signature making possible selections with high efficiencies and a low background contamination. In the $`\mathrm{q}\overline{\mathrm{q}}\tau ^+\tau ^{}`$ final state, the decay of the tau leptons produces events which are more difficult to identify. The identification of this final state exploits the missing momentum and missing energy carried away by the neutrinos produced in the decay of the tau lepton.

#### 3.3.1 Selection of $`𝐙𝐙\mathbf{}𝐪\overline{𝐪}𝐞^\mathbf{+}𝐞^{\mathbf{}}`$ and $`𝐙𝐙\mathbf{}𝐪\overline{𝐪}𝝁^\mathbf{+}𝝁^{\mathbf{}}`$ events

The selection of $`\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ and $`\mathrm{q}\overline{\mathrm{q}}\mu ^+\mu ^{}`$ final states requires the visible energy of the events to be greater than 90 GeV and at least six reconstructed tracks. Among all tracks with momenta greater than 2 GeV, the highest momentum track is taken as the first lepton candidate and the second-highest momentum track with a charge opposite to the first candidate is taken as the second lepton candidate. Using the Durham  jet algorithm, the event, including the lepton candidates, is forced into four jets and the jet resolution variable that separates the three-jet topology from the four-jet topology, $`y_{34}`$, is required to be greater than $`10^3`$. Excluding the electron or muon candidates and their associated calorimeter clusters, the rest of the event is forced into two jets. The 4C and 5C fits to the two lepton candidates and the two jets are required to converge.333 In the context of this paper, convergence is defined as a fit probability greater than $`10^{10}`$.

In the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ selection no explicit electron identification is used. Electron candidates are selected by requiring the sum of the electromagnetic cluster energies $`E_1+E_2`$ associated to the electrons to be greater than 70 GeV and the momentum of the most energetic electron track to exceed 20 GeV. We also reject the event if the angle between either electron candidate and any other track is less than 5.

In the $`\mathrm{q}\overline{\mathrm{q}}\mu ^+\mu ^{}`$ selection the muons are identified using (i) tracks which match a reconstructed segment in the muon chambers, (ii) tracks which are associated to hits in the hadron calorimeter or muon chambers , or (iii) isolated tracks associated to electromagnetic clusters with reconstructed energy less than 2 GeV. No isolation requirement is imposed on events with both muon tracks passing (i) or (ii). Events with at least one muon identified with criterion (iii) are accepted if both muon candidates in the event have an angle of at least 10 to the nearest track. We require the sum of the momenta of the two leptons to be greater than 70 GeV.

$`\mathrm{Z}`$-pair events are separated from $`\mathrm{Z}\gamma ^{}`$ background by requiring the fitted mass of the 5C fit to be larger than 85 GeV and the invariant masses $`m_{\mathrm{}\mathrm{}}`$ and $`m_{\mathrm{qq}}`$ obtained from the 4C fit to satisfy $`(m_{\mathrm{}\mathrm{}}+m_{\mathrm{qq}})>170`$ GeV and $`|m_{\mathrm{}\mathrm{}}m_{\mathrm{qq}}|<30`$ GeV. Figure 2b (2c) shows the distribution of $`m_{\mathrm{ee}}`$ ($`m_{\mu \mu }`$) and $`m_{\mathrm{qq}}`$ before the cuts on the masses from the 4C and 5C fits.

After all cuts the selection efficiency444 Small amounts of feedthrough from other $`\mathrm{Z}`$$`\mathrm{Z}`$ final states, in this case $`\mathrm{q}\overline{\mathrm{q}}\tau ^+\tau ^{}`$, are counted as signal. for $`\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ signal events is $`(55.2\pm 2.7)`$% at 183 GeV and $`(65.1\pm 2.9)`$% at 189 GeV. The errors on these efficiencies include the systematic errors (see Section 3.6). No candidate is observed at 183 GeV. Six candidate events are found after all cuts in the data taken at 189 GeV. In Table 1 (lines d and e) we give the efficiency, background and observed number of events. The largest source of background after all cuts is from $`\mathrm{Z}\gamma ^{}`$ mediated $`\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ events and from two-photon events.

For the $`\mathrm{q}\overline{\mathrm{q}}\mu ^+\mu ^{}`$ final state the selection efficiency is $`(66.9\pm 2.6)`$% at 183 GeV and $`(72.3\pm 2.8)`$% at 189 GeV. Three events are observed in the 189 GeV data sample and none in the 183 GeV data sample. The errors on these efficiencies include the systematic errors given below. In Table 1 (lines f and g) we give the efficiency, backgrounds and observed number of events. The background after all cuts is expected to come from $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\gamma ^{}\mathrm{q}\overline{\mathrm{q}}\mu ^+\mu ^{}`$ events.

#### 3.3.2 Selection of $`𝐞^\mathbf{+}𝐞^{\mathbf{}}\mathbf{}𝐙𝐙\mathbf{}𝐪\overline{𝐪}𝝉^\mathbf{+}𝝉^{\mathbf{}}`$ events

The $`\mathrm{q}\overline{\mathrm{q}}\tau ^+\tau ^{}`$ final state is selected from a sample of events with track multiplicity greater or equal to six. Events which have been selected as $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ or $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\mu ^+\mu ^{}`$ are excluded from this selection. The tau-lepton candidates are selected using an artificial neural network algorithm which is described in detail in Reference . The tau candidate with the highest neural network output value is taken as the first candidate. The second best candidate is required to have its charge opposite to the first candidate and the highest output value among all remaining candidates. If a second candidate cannot be found the event is rejected.

Motivated by the presence of neutrinos in the final state, the visible energy of the event, $`E_{\mathrm{vis}}`$, is required to exceed 90 GeV and the missing energy $`\sqrt{s}E_{\mathrm{vis}}`$ is required to exceed 15 GeV. In addition, the sum of the momenta of the leading tracks from the tau-lepton decays is required to be less than 70 GeV. Since the direction of the missing momentum in signal events will tend to be along the direction of one of the decaying tau leptons, the angle $`\alpha _{\tau ,\mathrm{miss}}`$ between the missing momentum and a tau-lepton candidate is required to satisfy $`\alpha _{\tau ,\mathrm{miss}}<90^{}`$ for at least one of the two tau candidates.

The two hadronic jets are selected in the same way as in the $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ selection. The initial estimate of the energy and the momenta of the tau candidates is found from the sum of the tracks associated to the tau by the neural network algorithm and all unassociated electromagnetic clusters in a cone with a half angle of 10 around the leading track from the tau decay. A 2C kinematic fit that imposes energy and momentum conservation (see the introduction to Section 3) is required to converge. A 3C kinematic fit, with the additional constraint of the equality of the fermion pair masses is also required to converge.

Using the network output  for each tau lepton, a probability is calculated taking into account the different branching ratios, sensitivities, efficiencies and background levels for 1-prong and 3-prong tau-lepton decays. In the following, we combine the probabilities $`𝒫_1`$ and $`𝒫_2`$ to form a likelihood using

$$=\frac{𝒫_1𝒫_2}{𝒫_1𝒫_2+(1𝒫_1)(1𝒫_2)}.$$
(1)
The likelihood associated with probabilities of the two tau candidates is required to satisfy $`_{\tau \tau }>0.977`$. In addition, the common mass of the 3C fit is required to exceed 85 GeV. Using the 2C fit masses of the tau pair, $`m_{\tau \tau }`$, and the quark pair, $`m_{\mathrm{qq}}`$, as obtained from the kinematic fit, we also require $`m_{\mathrm{qq}}+m_{\tau \tau }>170`$ GeV and $`|m_{\mathrm{qq}}m_{\tau \tau }|<60`$ GeV.

After all cuts the selection efficiency for signal events is found to be $`(22.2\pm 1.7)`$% at 183 GeV and $`(26.9\pm 2.0)`$% at 189 GeV. One candidate event is found in the data at 189 GeV while no candidate is selected at 183 GeV. Figure 2d shows the masses of the candidate events before the invariant mass cuts. In Table 1 (lines h and i) we give the efficiencies, backgrounds and observed number of events.

#### 3.3.3 Selection of $`𝐞^\mathbf{+}𝐞^{\mathbf{}}\mathbf{}𝐙𝐙\mathbf{}𝐛\overline{𝐛}\mathbf{}^\mathbf{+}\mathbf{}^{\mathbf{}}`$

Events with $`\mathrm{b}\overline{\mathrm{b}}`$ final states are selected using the algorithm described in Reference . The b probabilities of the two hadronic jets are combined to form a likelihood, $`_{\mathrm{bb}}`$, according to Equation 1. Because the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{e}^+\mathrm{e}^{}`$ and $`\mathrm{q}\overline{\mathrm{q}}\mu ^+\mu ^{}`$ selections are pure, a relatively loose cut of $`_{\mathrm{bb}}>0.2`$ is used to select the $`\mathrm{b}\overline{\mathrm{b}}\mathrm{e}^+\mathrm{e}^{}`$ and $`\mathrm{b}\overline{\mathrm{b}}\mu ^+\mu ^{}`$ samples. For the selections with electron and muon pairs there are two classes of events since the selected $`\mathrm{b}\overline{\mathrm{b}}\mathrm{}^+\mathrm{}^{}`$ events are a subset of the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$ events. In Table 1 (lines e and g) we give the efficiencies of the b-tagged samples with respect to the expected fraction of $`\mathrm{b}\overline{\mathrm{b}}`$ events. The efficiencies for samples without b-tags are given with respect to the hadronic decays without $`\mathrm{b}\overline{\mathrm{b}}`$ final states.

In the $`\mathrm{b}\overline{\mathrm{b}}\mathrm{e}^+\mathrm{e}^{}`$ selection one candidate is found in the data. The selection efficiency is found to be (47$`\pm `$3)% at 183 GeV and (46$`\pm `$3)% at 189 GeV. 555 The efficiencies given in this section do not include feedthrough from other $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$ final states. The efficiencies given in Table 1 include this feedthrough. In the $`\mathrm{b}\overline{\mathrm{b}}\mu ^+\mu ^{}`$ selection we find no candidate at 183 GeV and one candidate at 189 GeV. The selection efficiency is found to be (50$`\pm `$3)% at 183 GeV and (54$`\pm `$3)% at 189 GeV.

For the $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ selection the $`_{\tau \tau }`$ cut of the $`\mathrm{q}\overline{\mathrm{q}}\tau ^+\tau ^{}`$ selection is loosened and combined with $`_{\mathrm{bb}}`$ as follows. $`_{\tau \tau }`$ and $`_{\mathrm{bb}}`$ are both required to be greater than 0.1. The $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ probability for the event, $`_{\mathrm{bb}\tau \tau }`$, is calculated from Equation 1 with $`_{\tau \tau }`$ and $`_{\mathrm{bb}}`$ as inputs and required to exceed 0.95. After the cut on $`_{\mathrm{bb}\tau \tau }`$, the remaining cuts of the $`\mathrm{q}\overline{\mathrm{q}}\tau ^+\tau ^{}`$ selection are applied, giving a selection efficiency of $`(21\pm 3)`$% at 183 GeV and $`(24\pm 3)`$% at 189 GeV. No candidate event is found in the 183 GeV or 189 GeV data.

We also use an alternative jet-based $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ selection  and accept any event which passes either $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ analysis, but events previously selected by another $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$ selection are rejected. The alternative analysis uses a different approach to reconstruct the tau leptons. This event selection consists of a set of preselection cuts and a subsequent multivariate likelihood selection. Events are reconstructed as four jets using the Durham algorithm. Tau-lepton candidates are sought in the four jets using a likelihood technique to separate real tau leptons and fakes in quark jets. The tau and b-tag likelihood values are combined in a $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ likelihood which is maximized to choose the b jets and tau leptons of the event. This $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ likelihood uses tau and b-tag likelihoods and some topological variables as input. Events are accepted if their likelihood exceeds 0.6. In addition, the fitted 2C masses are required to satisfy $`m_{\mathrm{qq}}+m_{\tau \tau }>170`$ GeV and $`|m_{\mathrm{qq}}m_{\tau \tau }|<60`$ GeV.

At $`189\mathrm{GeV}`$ data the alternative selection has an efficiency of $`(30\pm 3)`$% for $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ events. After combining the two selections, the efficiency for $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ events for all cuts except those rejecting events found by the other $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$ selections is 40%. A similar improvement is realized for the $`183\mathrm{GeV}`$ data. The efficiency and backgrounds after rejecting events found by the other selections are given in Table 1 (line j). One exclusive event is selected by the combined $`\mathrm{b}\overline{\mathrm{b}}\tau ^+\tau ^{}`$ analysis at $`189\mathrm{GeV}`$.

### 3.4 Selection of $`𝐙𝐙\mathbf{}𝐪\overline{𝐪}𝝂\overline{𝝂}`$ events

The $`\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }`$ selection is based on the reconstruction of the $`\mathrm{Z}`$ boson decaying to $`\mathrm{q}\overline{\mathrm{q}}`$ which produces somewhat back-to-back jets. The selection uses contained events with a two-jet topology. The beam energy constraint is then used to determine the mass of the $`\mathrm{Z}`$ boson decaying to $`\nu \overline{\nu }`$. The properties of the $`\mathrm{q}\overline{\mathrm{q}}`$ decay and the inferred mass of the $`\nu \overline{\nu }`$ decay are then used in a likelihood analysis to separate signal from background.

Two-jet events are selected by dividing each event into two hemispheres using the plane perpendicular to the thrust axis. The number of charged tracks in each hemisphere is required to be four or more. The polar angles of the energy-momentum vector associated with each hemisphere, $`\theta _{\mathrm{hemi1}}`$ and $`\theta _{\mathrm{hemi2}}`$, are used to calculate the quantity $`\mathrm{cos}\theta _\mathrm{h}=\frac{1}{2}(\mathrm{cos}\theta _{\mathrm{hemi1}}\mathrm{cos}\theta _{\mathrm{hemi2}})`$. Contained events are selected by requiring $`|\mathrm{cos}\theta _\mathrm{h}|<0.80`$. The total energy in the forward detectors and in the forward region of the electromagnetic calorimeter ($`|\mathrm{cos}\theta |>0.95`$) is required to be less than 3 GeV. $`\mathrm{W}`$ decays identified by the OPAL $`\mathrm{W}`$-pair selection are rejected; the likelihood for $`\mathrm{e}^+\mathrm{e}^{}\mathrm{q}\overline{\mathrm{q}}\mathrm{}\nu `$ from Reference , $`_{\mathrm{WW}}`$, is required to be 0.5 or less.

An important background to our selection is $`\mathrm{q}\overline{\mathrm{q}}(n\gamma )`$ events with photons which escape detection. We discriminate against these events by looking for a significant amount of missing transverse momentum, $`p_\mathrm{t}`$. In each event, $`p_\mathrm{t}`$ can be resolved into two components, $`p_{\mathrm{t}i}`$, perpendicular to both the thrust axis and the beam axis and $`p_{\mathrm{t}j}`$, along the thrust axis and perpendicular to the beam axis. Since $`p_{\mathrm{t}i}`$ is based primarily on angular measurements, it is better measured than $`p_{\mathrm{t}j}`$. We approximate $`p_{\mathrm{t}i}`$ as $`p_{\mathrm{t}i}=\frac{1}{2}E_\mathrm{b}\mathrm{sin}\varphi \mathrm{sin}\theta _\mathrm{h}`$. Here $`E_\mathrm{b}=\sqrt{s}/2`$ is the beam energy, $`\varphi `$ is the acoplanarity of the momentum vectors of the two hemispheres and $`\mathrm{sin}\theta _\mathrm{h}=\sqrt{1\mathrm{cos}^2\theta _\mathrm{h}}`$. The resolution on $`p_{\mathrm{t}i}`$, $`\sigma _{p_{\mathrm{t}i}}`$, was parameterized as a function of thrust and $`\mathrm{cos}\theta _\mathrm{h}`$ using data taken at the $`\mathrm{Z}`$ resonance. The variable $`R_{p_{\mathrm{t}i}}=(p_{\mathrm{t}i}p_{\mathrm{t}i}^0)/\sigma _{p_{\mathrm{t}i}}`$ is used as an input to the likelihood. Here $`p_{\mathrm{t}i}^0`$ corresponds to the transverse momentum carried by a photon with half the beam energy which just misses the inner edge of our acceptance ( $`p_{\mathrm{t}i}^0=E_\mathrm{b}\mathrm{sin}(32\text{mrad})/2`$). We also use the variable $`\mathrm{cos}\theta _{\mathrm{miss}}`$, the direction of the missing momentum in the event, to discriminate against the $`\mathrm{q}\overline{\mathrm{q}}(n\gamma )`$ events.

In the final selection of events, we use a likelihood based on the following five variables: (i) the normalized sum of visible and recoil masses $`(m_{\mathrm{vis}}+m_{\mathrm{recoil}})/\sqrt{s}`$, (ii) the difference of visible and recoil masses ($`m_{\mathrm{vis}}m_{\mathrm{recoil}}`$), (iii) $`\mathrm{log}(y_{23})`$, where $`y_{23}`$ is the jet resolution parameter that separates the two-jet topology from the three-jet topology as calculated from the Durham jet algorithm, (iv) $`\mathrm{cos}\theta _{\mathrm{miss}}`$ and (v) $`R_{p_{\mathrm{t}i}}`$. The mass variables are useful for reducing background from $`\mathrm{W}`$-pair production and $`\mathrm{We}\nu `$ final states. The jet resolution parameter is useful in reducing the remaining $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}\nu `$ final states. To improve the performance of the likelihood analysis we use only events with : $`|m_{\mathrm{vis}}m_{\mathrm{recoil}}|<50`$ GeV, $`(m_{\mathrm{vis}}+m_{\mathrm{recoil}})/\sqrt{s}>0.89`$ and $`R_{p_{\mathrm{t}i}}>1.2`$. Events are then selected using $`_{\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }}>0.5`$, where $`_{\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }}`$ is the likelihood for the $`\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }`$ selection. The likelihood distribution of data and Monte Carlo is shown in Figure 3c. For the $`\mathrm{b}\overline{\mathrm{b}}\nu \overline{\nu }`$ selection we require, in addition, the b-tag variable of Reference  to be greater than 0.65.

The efficiencies for the $`\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }`$ selection alone are $`(30.5\pm 2.0)`$% at 183 GeV and $`(33.9\pm 2.3)`$% at 189 GeV. The errors on these efficiencies include the systematic errors discussed below in Section 3.6. The efficiencies after considering the results of the b-tagging, as well as the number of events selected at the two energies are given in Table 1 (lines k and l).

### 3.5 Selection of $`𝐙𝐙\mathbf{}𝐪\overline{𝐪}𝐪\overline{𝐪}`$ events

The fully hadronic channel of the $`\mathrm{Z}`$-pair decay has the largest branching fraction of all channels (about 50%), but suffers from large background from hadronic $`\mathrm{W}`$-pair decays. We apply two different event selections, one of which is based mainly on reconstructed mass information in order to accept all hadronic $`\mathrm{Z}`$-pair decays without flavor requirement, while a second analysis applies a flavor tag in order to select final states involving b quarks, allowing for looser requirements on the reconstructed boson mass.

For both subsamples, hadronic $`\mathrm{Z}^{}/\gamma ^{}\mathrm{q}\overline{\mathrm{q}}`$ events are an important background. We therefore start with a common preselection based on event shape variables which is mainly aimed at reducing this background. We use a likelihood method, described below in Section 3.5.1, in order to choose the most likely jet pairing for $`\mathrm{Z}`$-pair decays.

#### 3.5.1 Preselection and jet pairing

The event selection starts from the inclusive multihadron selection described in Reference . The radiative process $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\gamma \mathrm{q}\overline{\mathrm{q}}\gamma `$ is suppressed by requiring the effective center-of-mass energy after initial state radiation, $`\sqrt{s^{}}`$, to be larger than 150 GeV. $`\sqrt{s^{}}`$ is obtained from a kinematic fit  that allows for one or more radiative photons in the detector or along the beam pipe. The final state particles are then grouped into jets using the Durham algorithm . A four-jet sample is formed by requiring the jet resolution parameter $`y_{34}`$ to be at least $`0.003`$ and each jet to contain at least two charged tracks. In order to suppress $`\mathrm{Z}^{}/\gamma ^{}\mathrm{q}\overline{\mathrm{q}}`$ background, the event shape parameter $`C_{\mathrm{par}}`$ , which is large for spherical events, is required to be greater than 0.25. A 4C kinematic fit using energy and momentum conservation is required to converge. A 5C kinematic fit which forces the two jet pairs to have the same mass is applied in turn to all three possible combinations of the four jets. This fit is required to converge for at least one combination. The efficiencies of these preselection cuts are (86.4 $`\pm `$ 0.5) % and (88.9 $`\pm `$ 0.5) % for signal events at 183 GeV and 189 GeV, respectively.

In order to determine which pair of jets comes from each $`\mathrm{Z}`$, we calculate a likelihood function using the mass obtained from the 5C fit, the corresponding fit probability, and the difference between the two di-jet masses obtained from the 4C fit. In the YFSZZ simulation of $`\mathrm{Z}`$-pair decays the correct jet pairing has the highest likelihood output in (86.8 $`\pm `$ 0.5) % of the events. This fraction rises to (93.8 $`\pm `$ 0.5) % for the events after the final selection.

#### 3.5.2 Likelihood for the inclusive $`𝐙𝐙\mathbf{}𝐪\overline{𝐪}𝐪\overline{𝐪}`$ event selection

We use a likelihood selection with eight input variables for the selection of $`\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ events. The first variable is the jet pairing likelihood described above. Excluding the jet pairing with the largest difference between the two di-jet masses as obtained from the 4C fit, we identify among the remaining two possible pairings the one for which the 5C-fit mass is closer to the $`\mathrm{W}`$ mass. We use the difference between this 5C-fit mass and the $`\mathrm{W}`$ mass in order to discriminate against hadronic $`\mathrm{W}`$-pair events. Two variables that are sensitive to unobserved particles along the beam direction are the fitted center-of-mass energy and the sum of the cosines of the polar angles of the four jets. In order to discriminate against $`\mathrm{Z}^{}/\gamma ^{}`$ events, we use the difference between the largest and smallest jet energies after the 4C fit, and the angular variable $`j_{\mathrm{ang}}=E_4(1\mathrm{cos}\theta _{12}\mathrm{cos}\theta _{13}\mathrm{cos}\theta _{23})/\sqrt{s}`$, where $`E_4`$ is the smallest of the four jet energies, and the $`\theta _{ij}`$ are the opening angles between jets $`i`$ and $`j`$, with the jets ordered by energy. Finally we calculate from the momenta configuration of the four jets the effective matrix element for the QCD processes $`\mathrm{Z}^{}/\gamma ^{}\mathrm{q}\overline{\mathrm{q}}\mathrm{gg}`$ and $`\mathrm{Z}^{}/\gamma ^{}\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ as defined in Reference , and the matrix element for the process $`\mathrm{WW}\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ from Reference .

The distribution of the likelihood function calculated from these eight variables is shown in Figure 3d. In order to maximize the significance of the measured cross section, assuming a 10% relative systematic error on the background, we place a cut on the likelihood at 0.65. This cut leads to an efficiency of $`(33.0\pm 1.8)`$% (189 GeV) relative to all fully hadronic NC2 $`\mathrm{Z}`$-pair final states. In the data, 52 events are selected. We expect a total of 57.0 events from Standard Model processes, of which 17.6 originate from the $`\mathrm{Z}`$-pair signal, while 27.0 events are expected from hadronic $`\mathrm{W}`$-pair decays and 12.4 events from hadronic two-fermion processes. At 183 GeV the selection efficiency for fully hadronic $`\mathrm{Z}`$-pair decays is $`(25.3\pm 1.3)`$%. In the data 8 events are selected. The Standard Model expectation is 1.7 signal events and 7.2 background events.

#### 3.5.3 Likelihood for $`𝐙𝐙\mathbf{}𝐪\overline{𝐪}𝐛\overline{𝐛}`$ event selection

Jets originating from b-quarks are selected using the same b-tagging algorithm used in the $`\mathrm{b}\overline{\mathrm{b}}\mathrm{}^+\mathrm{}^{}`$ and $`\mathrm{b}\overline{\mathrm{b}}\nu \overline{\nu }`$ selections. We evaluate the probability for each of the four jets to originate from a primary b quark, and use the two highest probabilities as input variables for a likelihood to select $`\mathrm{Z}\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}\mathrm{b}\overline{\mathrm{b}}`$ events. In addition, we use the parameters $`y_{34}`$, $`C_{\mathrm{par}}`$, the difference between the largest and smallest jet energies and the output of the jet pairing likelihood. We also use the fit probabilities of a 5C kinematic fit which constrains one boson mass to the $`\mathrm{Z}`$ mass, and the probability of a 6C fit which forces both masses to be equal to the $`\mathrm{W}`$ mass.

Figure 3e shows the distribution of the likelihood function calculated from these eight variables for the preselected events. The signal likelihood is required to be larger than 0.80 for both 183 GeV and 189 GeV data. After the likelihood selection, we perform an additional cut on the mass obtained from the 5C-mass fit in the most likely jet pairing, which is required to be larger than 86 GeV. This choice of the cuts on the likelihood and on the 5C-fit mass was made by maximizing the expected statistical significance of signal over background. The final efficiency is $`(35.1\pm 2.4)`$% for the 183 GeV data and $`(38.6\pm 2.7)`$% for the 189 GeV data. The observed number of events, expected signal and background are given in Table 1 (line o).

#### 3.5.4 Combination

To account for overlap between the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ and $`\mathrm{q}\overline{\mathrm{q}}\mathrm{b}\overline{\mathrm{b}}`$ selections we divide the data into three logical classes, exclusive $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ events (Table 1 (line m)), exclusive $`\mathrm{q}\overline{\mathrm{q}}\mathrm{b}\overline{\mathrm{b}}`$ (line n) events and the overlapping region (Table 1 (line o)). In Table 1 (line m) we give the efficiency relative to all fully hadronic final states without b-quarks. The other efficiencies are given relative to fully hadronic final states with b-quarks.

### 3.6 Selection systematic errors

Systematic errors have only a modest effect on our final result because of the large statistical error associated with the small $`\mathrm{Z}`$-pair cross section.

Detector effects can best be studied by comparing calibration data taken at the $`\mathrm{Z}`$ resonance with a simulation of the same process. These comparisons are important for final states such as $`\mathrm{}^+\mathrm{}^{}\nu \overline{\nu }`$, $`\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }`$, and $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$, where the tight cuts are needed to separate signal and background. In these cases, we add additional smearing to the energy and momentum of the simulated events to match data and simulation. We then apply the same smearing to the signal and background Monte Carlos and then correct our efficiencies and background accordingly. The full difference is used as the systematic error in these cases. At $`\sqrt{s}=189`$ GeV, these differences give relative systematic errors on the efficiency of 2.5% for the $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }`$ final state, 5.2% for the $`\mu ^+\mu ^{}\nu \overline{\nu }`$ final state and 3.8% for $`\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }`$ final state. In the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ inclusive analysis, where uncertainties on the reconstructed angles are also important, similar studies lead to a relative detector systematic error of 6%.

Detector systematic errors for the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$ and $`\mathrm{}^+\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$ selections without $`\tau `$-pairs in the final state are small because of the good separation of signal and background. In the $`\mathrm{}^+\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$ final state the largest effect (3%) is from modeling of the multiplicity requirement which is important for final states containing $`\tau `$-pairs. In the $`\mathrm{q}\overline{\mathrm{q}}\tau ^+\tau ^{}`$ selection, the systematic uncertainties in the efficiencies were determined by overlaying hadronic and tau decays taken from $`\mathrm{Z}`$ resonance data giving a contribution to the systematic error of 6.0%. When the final states are combined to determine the cross section and limits on the anomalous triple gauge couplings, we assume a common relative systematic error of 3% on the efficiencies.

Another important detector effect comes from the simulation of the variables used by the OPAL b-tag which is discussed in Reference . We allow for a common 5% error in the efficiency of the b-tag, consistent with our studies on $`\mathrm{Z}`$ resonance data and Monte Carlo.

These detector effects were propagated through to our background errors.

In each channel the signal and background Monte Carlo generators have been compared against alternative generators. In almost all cases the observed differences are consistent within the finite Monte Carlo statistics and the systematic error has been assigned accordingly. One notable exception is in the $`\mathrm{q}\overline{\mathrm{q}}\mathrm{b}\overline{\mathrm{b}}`$ background where differences in the PYTHIA, grc4f and EXCALIBUR simulation of the $`\mathrm{W}`$-pair background are as large as 20%. In this case the full difference has been assigned as the systematic error.

In the $`\mathrm{}^+\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$ channels, which has a large background from two-photon events, we have compared the number of selected events at an early stage of the analysis with the Monte Carlo prediction and based our background systematic error on the level of agreement. This results in 20% systematic error.

## 4 Results

We combine the information from all of the analyses reported above using a maximum likelihood fit to determine the production cross section for $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}`$. The information which was used in the fit, as well as the Standard Model prediction for $`\mathrm{Z}`$-pair production, is summarized in Table 1. For each channel the table gives the number of events observed, $`n_{\mathrm{obs}}`$, the Standard Model prediction for all events, $`n_{\mathrm{SM}}`$, the expected signal, $`n_{\mathrm{ZZ}}`$, the expected background, $`n_{\mathrm{back}}`$, the efficiency $`ϵ_{\mathrm{chan}}`$, and the integrated luminosity, $`L_{\mathrm{int}}`$. $`B_{\mathrm{ZZ}}`$ is the branching ratio of $`\mathrm{Z}`$-pairs to the given final state, calculated from $`\mathrm{Z}`$ resonance data . In the table we give the overlap between the b-tag and non-b-tag analyses. Possible overlap between $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ and $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$ has been studied, and found to be an order of magnitude smaller than the overlap of $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ and $`\mathrm{q}\overline{\mathrm{q}}\mathrm{b}\overline{\mathrm{b}}`$ and has therefore been ignored.

The cross section at each energy is determined with a maximum likelihood fit using a Poisson probability density convolved with Gaussians to describe the uncertainties on efficiencies and backgrounds. The expected number of events in each channel, $`\mu _\mathrm{e}`$, as function of the $`\mathrm{Z}`$-pair cross section, $`\sigma _{\mathrm{ZZ}}`$, is given by

$$\mu _\mathrm{e}=\sigma _{\mathrm{ZZ}}L_{\mathrm{int}}ϵ_{\mathrm{chan}}B_{\mathrm{ZZ}}+n_{\mathrm{back}}.$$
(2)
The efficiencies, $`ϵ_{\mathrm{chan}}`$, include the effects of off-shell $`\mathrm{Z}`$ bosons that are produced outside of our kinematic acceptance. Our main result, the NC2 $`\mathrm{Z}`$-pair cross sections obtained from the fits, is

$$\begin{array}{cccc}\sigma _{\mathrm{ZZ}}(183\mathrm{GeV})\hfill & =& 0.12_{0.18}^{+0.20}(\text{stat.})_{0.02}^{+0.03}(\text{syst.})\hfill & \mathrm{pb}\hfill \\ & & & \\ \sigma _{\mathrm{ZZ}}(189\mathrm{GeV})\hfill & =& 0.80_{0.13}^{+0.14}(\text{stat.})_{0.05}^{+0.06}(\text{syst.})\hfill & \mathrm{pb}.\hfill \end{array}$$

The 183 GeV result corresponds to a 95% C.L. upper limit (normalized to the region $`\sigma _{\mathrm{ZZ}}>0`$) of $`0.55\mathrm{pb}`$. The comparison of these measurements with the YFSZZ prediction, using the coupling of $`\mathrm{Z}`$ bosons to electrons measured at the $`\mathrm{Z}`$ resonance, is shown in Figure 4. The results are consistent with the YFSZZ prediction and with the measurements presented in Reference .

Our results assume that, apart from the backgrounds discussed above, only $`\mathrm{Z}`$$`\mathrm{Z}`$ production contributes inside our kinematic region. Possible effects of Higgs boson production have been ignored.

In order to check the results for consistency with the expected fraction of $`\mathrm{b}\overline{\mathrm{b}}`$ final states, we perform a second fit where the branching ratio of the $`\mathrm{Z}`$ boson to $`\mathrm{b}\overline{\mathrm{b}}`$ is a free parameter. The relative branching ratios of the $`\mathrm{Z}`$ to other fermion pairs are fixed to their measured values. The resulting branching ratio and cross section for the 189 GeV data are $`\text{Br(}\mathrm{Z}\mathrm{b}\overline{\mathrm{b}})=0.21_{0.06}^{+0.07}(\mathrm{stat}.)\pm 0.01(\mathrm{syst}.)`$ and $`\sigma _{\mathrm{ZZ}}=0.75_{0.14}^{+0.15}(\text{stat.})_{0.05}^{+0.06}(\text{syst.})\text{pb}`$. The measured branching ratio can be compared with the world average as measured at the $`\mathrm{Z}`$ resonance of $`\text{Br(}\mathrm{Z}\mathrm{b}\overline{\mathrm{b}})=0.1516\pm 0.0009`$. The cross section measured without constraining the branching ratio is also consistent with our main result and the YFSZZ prediction. The smaller error on the main result cross section illustrates the advantage of classifying the hadronic systems as $`\mathrm{b}\overline{\mathrm{b}}`$ or non-$`\mathrm{b}\overline{\mathrm{b}}`$. The 183 GeV data sample is too small to extract a meaningful value of the $`\mathrm{Z}\mathrm{b}\overline{\mathrm{b}}`$ branching ratio.

Limits on anomalous triple gauge couplings were set using the total cross section and the $`|\mathrm{cos}\theta _\mathrm{Z}|`$ distribution of our data. Here $`\theta _\mathrm{Z}`$ is the polar angle of the $`\mathrm{Z}`$ bosons produced. In this study we varied the real and imaginary parts of the $`\mathrm{Z}`$$`\mathrm{Z}`$$`\mathrm{Z}`$ and $`\mathrm{Z}\mathrm{Z}\gamma `$ anomalous couplings parameterized by the form factors $`f_4^{\mathrm{ZZZ}}`$, $`f_5^{\mathrm{ZZZ}}`$, $`f_4^{\mathrm{ZZ}\gamma }`$ and $`f_5^{\mathrm{ZZ}\gamma }`$ as defined in Reference  and implemented in the YFSZZ Monte Carlo. The real and imaginary parts of each coupling were varied separately with all others fixed to zero.

For this study we consider the effect of the anomalous couplings on the total cross section at $`\sqrt{s}=183`$ GeV and on the cross section in four bins of $`|\mathrm{cos}\theta _\mathrm{Z}|`$ at $`\sqrt{s}=189`$ GeV. The selection efficiencies for all final states are parameterized as function of anomalous couplings. At $`\sqrt{s}=189`$ GeV the parameterization is done separately for each bin in $`|\mathrm{cos}\theta _\mathrm{Z}|`$. For values of the anomalous couplings larger than unity, much of the production of the final state fermions occurs at $`|\mathrm{cos}\theta |1`$ where the efficiency for most channels is reduced by a factor of $`0.5`$. An uncertainty of 10%, dominated by Monte Carlo statistical errors, is assigned to the correction we apply to these efficiencies.

The 95% C.L. limits on the anomalous couplings obtained from the maximum likelihood fit are given in Table 2. With the exception of the $`Re\{f_5^{\mathrm{ZZZ}}\}`$ coupling, the limits are insensitive to the sign and complex phase of the couplings.

## 5 Conclusion

The production cross section of $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}`$ has been measured using the final states $`\mathrm{}^+\mathrm{}^{}\mathrm{}^+\mathrm{}^{}`$, $`\mathrm{}^+\mathrm{}^{}\nu \overline{\nu }`$, $`\mathrm{q}\overline{\mathrm{q}}\mathrm{}^+\mathrm{}^{}`$, $`\mathrm{q}\overline{\mathrm{q}}\nu \overline{\nu }`$, and $`\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$. The number of observed events, the background expectation from Monte Carlo and the calculated efficiencies have been combined to measure the production cross section of the process $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}`$. Our measured cross sections include the effects of background and efficiency uncertainties.

We have determined the cross section for $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Z}\mathrm{Z}`$ separately at average center-of-mass energies of $`182.62\pm 0.05`$ GeV and $`188.63\pm 0.04`$ GeV. The NC2 $`\mathrm{Z}`$-pair cross sections were determined to be

$$\begin{array}{cccc}\sigma _{\mathrm{ZZ}}(183\mathrm{GeV})\hfill & =& 0.12_{0.18}^{+0.20}(\text{stat.})_{0.02}^{+0.03}(\text{syst.})\hfill & \mathrm{pb}\hfill \\ & & & \\ \sigma _{\mathrm{ZZ}}(189\mathrm{GeV})\hfill & =& 0.80_{0.13}^{+0.14}(\text{stat.})_{0.05}^{+0.06}(\text{syst.})\hfill & \mathrm{pb}.\hfill \end{array}$$

At the lower center-of-mass energy, the 95% C.L. upper limit on the cross section is 0.55 pb. The measurements at both energies are consistent with the Standard Model expectations. No evidence is found for anomalous neutral current triple gauge couplings. The 95% confidence level limits are listed in Table 2.

## Acknowledgements

We particularly wish to thank the SL Division for the efficient operation of the LEP accelerator at all energies and for their continuing close cooperation with our experimental group. We thank our colleagues from CEA, DAPNIA/SPP, CE-Saclay for their efforts over the years on the time-of-flight and trigger systems which we continue to use. In addition to the support staff at our own institutions we are pleased to acknowledge the
Department of Energy, USA,
National Science Foundation, USA,
Particle Physics and Astronomy Research Council, UK,
Natural Sciences and Engineering Research Council, Canada,
Israel Science Foundation, administered by the Israel Academy of Science and Humanities,
Minerva Gesellschaft,
Benoziyo Center for High Energy Physics,
Japanese Ministry of Education, Science and Culture (the Monbusho) and a grant under the Monbusho International Science Research Program,
Japanese Society for the Promotion of Science (JSPS),
German Israeli Bi-national Science Foundation (GIF),
Bundesministerium für Bildung, Wissenschaft, Forschung und Technologie, Germany,
National Research Council of Canada,
Research Corporation, USA,
Hungarian Foundation for Scientific Research, OTKA T-029328, T023793 and OTKA F-023259.


format: text

text:
# Reversing type II migration: resonance trapping of a lighter giant protoplanet

## 1 Introduction

In the past few years a number of extrasolar giant planets have been discovered around nearby solar–type stars. These objects masses range from $`0.17`$ $`M_J`$ to $`11`$ $`M_J`$ (where $`M_J`$ is Jupiter’s mass) and their orbital semi-major axis range from $`0.038`$ AU to $`3.3`$ AU (Marcy, Cochran & Mayor, 1999). Although many uncertainties remain about planet formation, it is now commonly accepted that planets have formed in and from protoplanetary discs. Necessarily, there must be some time interval over which a giant planet and the surrounding gaseous disc material coexist. The planet and the disc exchange angular momentum through tidal interactions which generally make the planet lose angular momentum This mechanism is called migration. It can roughly be divided in two regimes:

* If the planet mass is small enough, the disc response is linear. The migration rate, in that regime, is proportional to the planet and disc masses, is independent of the viscosity and weakly dependent of the disc surface density and temperature profiles. This is the so-called type I migration (Ward, 1997).
* When the protoplanet mass is above a certain threshold, the torques acting locally on the surrounding disc material open a gap (Papaloizou & Lin, 1984), whose width and depth are controlled by the balance between the tidal torques, which tend to open the gap, and the viscous torques which tend to close it. The disc response is significantly non-linear, and most of the protoplanet Lindblad resonances fall in the gap and therefore cannot contribute to the planet-disc angular momentum exchange. The migration rate slows down dramatically compared to type I migration. Furthermore, the tidal truncation splits the disc into two parts and the planet is locked to the disc viscous evolution (Nelson et al. 2000). This is the type II migration, which describes the orbital evolution of giant protoplanets.

In this letter we consider the coupled evolution of a system of giant protoplanets consisting of two non-accreting cores with masses $`1`$ $`M_J`$ and $`0.29`$ $`M_J`$, which we are going to call from now on respectively “Jupiter” and “Saturn”. Attempts have already been made to describe the behaviour of a system of planets embedded in a disc. Melita & Woolfson (1996) and Haghighipour (1999) considered an embedded Jupiter and Saturn system orbiting a solar mass star, and showed how resonance trapping would affect their evolution. However the dissipative force in these works was due to the dynamical friction with a uniform density interplanetary medium, hence type II migration effects were not taken into account. Resonance trapping of planetesimals by a fixed orbit Jupiter sized protoplanet has also been investigated by Beaugé et al. (1994), and shown to be able to build up a single planetary core with orbital characteristics close to Saturn’s ones. Kley (2000) studied the orbital evolution of two maximally accreting giant cores embedded in a minimal mass protosolar disc, and showed how the migration of the inner core could be halted by the presence of the outer one, and how the eccentricity of the inner core is pumped up by the outer one.

## 2 Results

### 2.1 Numerical codes description

In order to investigate the long-term behaviour of the embedded Jupiter and Saturn system, we have used two independent hydrocodes, which have been described elsewhere in full detail (Nelson et al. 2000). These two codes are fixed Eulerian grid based codes, one of them is NIRVANA (Ziegler & Yorke, 1997) and the other one has been written by one of us (FM). Both have been endowed with the fast advection FARGO algorithm (Masset, 2000), and can run either with this algorithm or with a standard advection algorithm. They gave very similar results. They consist of a pure N-body kernel based on either a fourth (NIRVANA) or fifth order adaptative timestep Runge-Kutta solver (sufficient for the short time-scales involved in this dissipative problem) embedded in a hydrocode which provides a tidal interaction with a 2D non self-gravitating gaseous disc. The simulations are performed in the non-inertial non-rotating frame centered on the primary. The grid outer boundary does not allow inflow nor outflow and is chosen sufficiently far from the planets in order for the spiral density waves that they launch to be damped before they reach it, while the grid inner boundary only allows outflow (inwards), so that the disc material can be accreted on to the primary. Failing to do so may lead to overestimate the inner disc density and artificially favours an outwards migration. In the following our length unit is $`5.2`$ AU, the mass unit is one solar mass, and the time unit is the initial orbital period of Jupiter (the actual period may vary as Jupiter migrates). The disc aspect ratio $`H/R`$ is uniform and constant. In the run presented here the grid resolution adopted is $`N_r=122`$ and $`N_\theta =300`$ with a geometric spacing of the interzone radii — such that all the zones are “as square as possible”, i.e. $`N_r\mathrm{log}(1+2\pi /N_\theta )=\mathrm{log}(R_{max}/R_{min})`$ —. The grid outer boundary is at $`R_{max}=5`$ and its inner boundary is at $`R_{min}=0.4`$. The geometric spacing is the most natural one since the disc thickness scales as $`r`$. On the other hand, a constant spacing leads to an oversampling of the outer disc and an undersampling of the inner one, and therefore is likely to favour an inwards migration.

### 2.2 Initial setup

The cores we consider are embedded in a gaseous minimal mass protosolar nebula around a unit mass central object, and we assume they start their evolution with semi-major axis $`a_j=1`$ for Jupiter and $`a_s=2`$ for Saturn. The disc surface density is uniform and corresponds to two Jupiter masses inside Jupiter’s orbit. The effective viscosity $`\nu `$, the nature of which remains unclear and is usually thought to be due to turbulence generated by MHD instabilities (Balbus & Hawley 1991), is assumed to be uniform through the disc and corresponds to a value of $`\alpha 610^3`$ in the vicinity of Jupiter’s orbit. The disc aspect ratio is $`H/r=0.04`$.

The mass of Jupiter is sufficient to open a deep gap and hence it settles in a type II migration (Nelson et al. 2000), whereas Saturn is unable to fully empty its coorbital region because: (i) its mass is smaller; (ii) The planet is in a regime known as the inertial limit (Ward & Hourigan, 1989) where the inwards migration speed is so high that it makes the planet pass through what would be the gap inner edge before it had time to actually open it.

Therefore Saturn does not clear a deep gap initially, and its migration rate is typical of type I migration, since all its Lindblad resonances can still contribute to the angular momentum exchange with the disc.

### 2.3 Run results

We present in fig. 1 the central star–planet distance curves as a function of time. We see how initially Jupiter migrates as if it was the only planet in the disc (see test run). In the meantime, Saturn starts a much faster migration (the obvious initial acceleration of its migration will be discussed elsewhere), and reaches the 1:2 resonance with Jupiter at time $`t110`$. The eccentricities at that time are small (see fig. 2), and in particular Saturn’s eccentricity is much smaller than the eccentricity threshold below which the capture into resonance is certain if the “adiabatic” condition on the migration rate is satisfied (Malhotra 1993): $`|\underset{s}{\overset{.}{a}}|/(a_s\mathrm{\Omega }_s)0.5j(j+1)\mu _Je_s`$ for the $`j`$:$`j+1`$ resonance, where $`\mu _J`$ is the mass ratio of Jupiter to the central object, and where $`e_s`$ is Saturn’s eccentricity. This condition is not satisfied when Saturn reaches the 1:2 resonance, and it passes through.

The planets then obtain higher eccentricities, and Saturn’s migration rate is reduced. Saturn’s eccentricity increases again rapidly as it passes through the 3:5 resonance with Jupiter at $`t220`$. Eventually the adiabatic condition on the migration rate is satisfied for the 2:3 resonance and Saturn’s eccentricity is still below the corresponding critical threshold, so it gets trapped into the 2:3 resonance with Jupiter (both $`e`$ and $`e^{}`$ resonances, since the two critical angles $`\varphi =3\lambda _s2\lambda _j\stackrel{~}{\omega }_s`$ and $`\varphi ^{}=3\lambda _s2\lambda _j\stackrel{~}{\omega }_j`$ librate, where $`\lambda `$ is the mean longitude and $`\stackrel{~}{\omega }`$ the longitude of perihelion). At that time both planets steadily migrate outwards.

### 2.4 Interpretation

We define the system of interest as the system composed of the two planets. This resonance locked system interacts with the inner disc through torques proportional to $`M_J^2`$, at Jupiter’s inner Lindblad resonances (ILR), whereas it interacts with the outer disc through torques proportional to $`M_S^2`$ at Saturn’s outer Lindblad resonances (OLR), as indicated on fig. 3.

It can be seen that Saturn’s ILR fall in Jupiter’s gap and Jupiter’s OLR fall in Saturn’s gap so their effect is weakened compared to the situation where Jupiter is alone. As $`M_J^2/M_S^210`$, the torque imbalance does not favour an inwards migration as strongly as in a one planet case, and may even lead to a positive differential Lindblad torque on the two planet system. Actually one can estimate what the maximum mass ratio of the outer planet to the inner one should be to get a migration reversal, if one neglects the Inner Lindblad torque on the outer planet and the Outer Lindblad torque on the inner planet. The Inner Lindblad torque on the inner planet reads as:

$$T_{ILR}=C_{ILR}\mu _J^2\mathrm{\Sigma }_0a_J^2(a_J\mathrm{\Omega }_J)^2h^3$$
(1)
where $`C_{ILR}`$ is a dimensionless coefficient which is a sizable fraction of unity (Ward, 1997), and where $`h^{}`$ is the disc aspect ratio. There is a similar formula for the Outer Lindblad torque on the outer planet (obtained by substituting the $`ILR`$ and $`J`$ indices in Eq. (1) respectively with $`OLR`$ and $`S`$). The resulting torque imbalance will be positive if: $`T_{ILR}>T_{OLR}`$, which reads here as:

$$\frac{\mu _S}{\mu _J}<\left(\frac{C_{ILR}}{C_{OLR}}\right)^{1/2}\left(\frac{2}{3}\right)^{1/3}$$
(2)
If we assume that $`C_{ILR}=C_{OLR}`$ then we get: $`\mu _S/\mu _J<0.87`$, whereas if we make the conservative assumption that $`C_{ILR}=\frac{1}{2}C_{OLR}`$, we have: $`\mu _S/\mu _J<0.62`$. This threshold is much bigger than the actual ratio, therefore if the common gap is deep enough to shut off Jupiter’s OLR torques (and Saturn’s ILR torques) then the net Lindblad torque on the two planet system is positive. As the two planet system proceeds outwards in the disc, it does not act on the gas as a snow-plough, but rather it allows the material from the outer disc to travel across the common gap and eventually feed the inner disc. We can find the gap “permeability” condition by requiring that the rate of angular momentum change of the ring of material lying immediately outside Saturn’s gap that is required to expand accordingly to Saturn’s orbit (snow-plough effect) is greater than the torque available from Saturn (at most the sum of its outer Lindblad torques, in which case we need to assume that the waves excited at their OLR are damped locally). In our case, this turns out not to be the case and most of the outer disc material flows through the common gap to the inner disc. We find that in all our runs it is possible to check that the rate of mass flow through the common gap (see fig. 4) can be expressed as :

$$\stackrel{.}{M}3\pi \nu \mathrm{\Sigma }_0+2\pi r_s\underset{s}{\overset{.}{r}}\mathrm{\Sigma }_0$$
(3)
with a reasonable precision ($`1020`$ %). Furthermore we have performed many “restart runs” which consist in restarting a run once Jupiter and Saturn are locked into resonance, and then by varying one parameter at one time, e.g. the viscosity or the aspect ratio (which changes the Lindblad torques and therefore the migration rate). The mass flux through the gap rapidly switches (in a few tens of orbits) to a new value after the restart, so that Eq. (3) remains fulfilled.

From the considerations above we can conclude that the presence of Saturn unlocks Jupiter from the disc evolution : the two planet system evolution (outwards) and the disc viscous evolution (inwards) are basically decoupled. This decoupling and the corresponding mass flow through the common gap has two consequences :

* A refilling of the inner disc, which is too depleted for the torques at Jupiter’s ILR to have any sizable effect in the one planet case (the inner disc is accreted on the primary on its short viscous time-scale and maintaining its surface density at not too low a value implies a permanent flow of material from the outer disc to the inner).
* The angular momentum lost by the material which flows from the outer disc to the inner one is gained by the planets. The exchange of angular momentum between a planet and a gas fluid element occurs during a “close encounter” between these two, the one planet version of which corresponds to the angular momentum exchange at each end of a horseshoe orbit of the fluid element. The resulting torque is the so-called coorbital corotation torque (Goldreich & Tremaine 1979, Ward 1991 and 1992). To the best of our knowledge, an analytical evaluation of the corotation torque in the case of a non-vanishing net mass flow through the orbit (either due to a viscous accretion on to the primary or to a radial migration or both) has not been performed yet. Obviously even the one planet case deserves a large amount of work on this specific topic, therefore the estimate of the corotation torque in this two planet problem is far beyond the scope of this paper. We will just comment that the corotation torque in our case might not be negligible compared to the differential Lindblad torque at some stage.

## 3 Discussions

We have performed a series of restart runs (see section 3) in order to check for a variety of behaviours.

### 3.1 Differential Lindblad torque sign

The one sided Lindblad torque has been shown to be proportional to $`h^3`$ (Ward 1997). We have performed two restart runs ($`h^{}=0.040.03`$ and $`h^{}=0.040.05`$) in order to check that the migration rate variation is consistent with this dependence. This is indeed the case. We note in passing that the migration rate varies as $`h^3`$, and not as $`h^2`$ as it would be the case in a one planet problem, since the Outer/Inner Lindblad torque asymmetry does not vanish as the disk thickness tends to zero (the OLRs would pile-up at Saturn’s orbit, whereas the ILRs would pile-up at Jupiter’s orbit). These results confirm that the behaviour we observe occurs mainly due to the differential Lindblad torque and shows as well that this latter quantity is positive, as expected from Eq. (2).

### 3.2 $`\alpha `$-viscosity vs. uniform viscosity

So far we have only considered a uniform viscosity. Switching to a uniform-$`\alpha `$ viscosity of the form $`\nu =\alpha c_sH`$ makes $`\nu `$ scale here as $`r^{1/2}`$, so the viscosity at the outer edge of the common gap is higher, whereas it is smaller in the inner disc. This has the following effect, which plays in favour of enhancing the migration reversal mechanism: the viscous time-scale of the inner disc is higher and therefore its surface density increases accordingly, since the material brought through the gap piles-up in the inner disc for a longer time before being accreted on the primary. This has been checked with a restart run.

### 3.3 Accretion on to the planets

The cores considered above do not accrete gas from the disc. One can wonder what would be the effects of accretion. We have performed a number of restart runs in order to investigate the effect of accretion on the mechanism presented here. We have only considered accretion on to Jupiter, as it is likely that the accretion rate on Saturn can be regarded as being negligible (i.e. its mass doubling time is much longer than the timescale of the outwards migration, see e.g. Pollack et al. 1996). The prescription we used to model accretion on to Jupiter consists in removing a proportion of the material which lies in the inner Roche lobe (i.e. a sphere with a radius of half the Hill radius). The amount which is removed in one timestep is calculated from the half emptying time of the inner Roche lobe $`\tau _{1/2}`$. We have performed four different restart runs, corresponding to the following values of $`\tau _{1/2}`$: $`\tau _{1/2}=T_0`$ (maximally accreting core, see Kley 1999), $`\tau _{1/2}=3T_0`$, $`\tau _{1/2}=10T_0`$ and $`\tau _{1/2}=30T_0`$, where $`T_0=2\pi /\mathrm{\Omega }_J`$ is Jupiter’s orbital time. In each of these cases, turning on accretion had no impact on the system migration rate, at least in the early stages: in the first case, the mass doubling time for Jupiter is relatively short, and when Jupiter’s mass is significantly larger than its initial mass some additional effects, which will be presented in much greater detail elsewhere, affect the migration rate which then differ from the non-accreting case.

### 3.4 Smoothing

The smoothing parameter of the potential can have a dramatic impact on Saturn’s initial migration rate. This rate is controlled by a subtle balance between outer disc and inner disc torques. In the case of Saturn, all the Lindblad resonances play a role, since there is no gap. Many prescriptions for the smoothing are unable to give trustworthy results for the balance between the outer and inner torques since, depending on the prescription, these two quantities are affected in a different way. On the other hand Jupiter’s migration rate is much more robust, since the presence of the gap prevents high-$`m`$ Lindblad resonances playing a role in the migration, which is therefore controlled only by remote, low $`m`$ resonances and thus almost insensitive to the smoothing parameter. For this reason we have adopted an approach which involves choosing a smoothing prescription which endows Saturn with a migration velocity of the order of magnitude of the linear analytical predictions (type I migration), which is needed to give correct results for the capture into resonance. Once Saturn is trapped into resonance with Jupiter, it is dynamically slaved by the latter and the system evolution is only very weakly affected by the exact value of the outer disc torque exerted on Saturn. We have found that using either of the two prescriptions below satisfactorily preserves the analytical torque imbalance on Saturn and therefore gives it a type I migration rate:

* The potential of a planet acting on the disc is smoothed over the length $`\epsilon =0.4R_H`$ where $`R_H`$ is the Hill radius of the planet under consideration, whereas the potential of the disc acting back on the planet is smoothed over $`\epsilon ^{}=\sqrt{H^2+d^2}`$ where $`H`$ and $`d`$ are respectively the local disc thickness and zone diagonal. Since $`\epsilon ^{}\epsilon `$ the action-reaction law is not fulfilled and the numerical biases which arise favour an inwards migration, as can be easily checked.
* The potential of a planet acting on the disc and the potential of the disc acting on the same planet are smoothed over $`\epsilon =0.4R_H`$. This prescription does fulfill the action-reaction law.

In both these two cases, as in any other which gives Saturn a type I migration rate, including runs performed with a uniform radial spacing, the migration gets reversed. The run presented here corresponds to the first prescription.

### 3.5 Impact of mass ratio and Long-term behaviour

One can wonder about the size of the interval of “Saturn”’s mass which causes the migration to be significantly slowed down or reversed. If “Saturn” is not massive enough it will not significantly affect Jupiter’s evolution (the common “gap” will be too full on Saturn’s side, and therefore Jupiter’s OLR torques will not be shut off), whereas if it is too massive, the torque imbalance will be negative again. Work is in progress to accurately determine which range of parameters leads to a migration reversal. It should be noted that the results presented here depend on the artificial initial conditions. We have performed other runs in which Saturn is initially very close either to the $`1:2`$ or $`3:5`$ resonance, and it turns out that neither of these resonances is able to struggle against the strong Lindblad torques on Saturn: no resonance angle can be found which provides a resonant torque on Saturn which counteracts the tide. Therefore a trapping into the $`2:3`$ resonance is the most likely outcome when the system is still embedded in a massive disc, whatever the initial conditions: catching-up of “Saturn” or in-situ assembling from smaller, type I migrating bodies.

The long-term behaviour of the system is twofold:

* The system is locked into resonance as long as :

+ The two-planet system can adjust its resonance angle in order to prevent the planets being “pushed” towards each other by the Lindblad torques exerted by the disk on each of them. In all our runs we have never observed this behaviour. Now, given the small eccentricities involved here, and given the fact that the adiabatic criterion threshold increases as $`j(j+1)`$, the most probable outcome is that Saturn would then be captured in the next order resonance, that is to say $`3:4`$, and all the physics exposed in this paper would still be valid (presence of a common gap, sharing of the coorbital material by the two planets, mass-weighted torque imbalance, etc.)
    + The planets are not pulled apart by any other torques. Now we have mentioned the possibly important role of the coorbital corotation torque in this problem, which may be sufficient to move the planets apart at some stage, in which case we may ultimately get a low eccentricity double giant planet system when the disc disappears. This will be presented in greater detail elsewhere.
* If the planets happen to be locked into resonance at the time that the gas effects become negligible, then the system is likely to be unstable (we mentioned already that at least two angles librate simultaneously, which strongly suggests a possible chaotic behaviour; see also Kley 2000), and the most likely outcome is that one planet will be ejected whereas the other planet will end up on an eccentric orbit. This could account for the observed eccentricities of the extrasolar planets which are not orbiting close to their host star, i.e. which have not migrated all the way to the star.

## 4 Acknowledgements

We wish to thank J.C.B. Papaloizou, R.P. Nelson, C. Terquem, J.D. Larwood, A.A. Christou and an anonymous referee for useful comments and criticism. This work was partially supported (for F.M.) by the research network “Accretion onto black holes, compact stars and protostars” funded by the European Commission under contract number ERBFMRX-CT98-0195, and additionally supported (for M.S.) by funding from a PPARC research studentship. Computational resources of the Grand HPC consortium were available and are gratefully acknowledged. We thank Udo Ziegler for making a FORTRAN version of his code NIRVANA publicly available.