text:
# 1 Introduction

## 1 Introduction

Electron-positron annihilation into hadrons at high energies provides a precise means to test Quantum Chromodynamics (QCD), the theory of strong interactions. Due to the purely leptonic initial state, there are many experimental quantities for which the long-distance (non-perturbative) effects are modest, especially when compared to similar quantities in hadron-hadron collisions or deep inelastic scattering. These quantities, for instance the total cross section and jet-related correlations, can be calculated in perturbative QCD as a function of a single parameter, the strong coupling strength, $`\alpha _s`$. Therefore, many QCD tests based on measurements of $`\alpha _s`$ have been carried out at LEP, the Large Electron-Positron collider at CERN (see, for example, ).

A key ingredient of QCD is the underlying gauge group. The simultaneous measurement of the strong coupling and the eigenvalues of the quadratic Casimir operator of the underlying gauge theory, the $`C_F`$ and $`C_A`$ colour factors, provides a more general and comprehensive test of QCD than measurements of $`\alpha _s`$ alone. The possible existence of a light gluino111The exclusion of gluinos with mass between 3–5 GeV has been debated in the literature recently . influences both $`\alpha _s`$ and the measured value of the colour factors (or, assuming SU(3) dynamics, the number of light fermionic degrees of freedom $`N_f`$). A simultaneous fit of these parameters to data provides a means to investigate whether the data favour or exclude the additional degrees of freedom.

Previous analyses of LEP data have been performed to probe the underlying gauge structure of the theory . The leading order perturbative prediction was fitted to distributions of four-jet angular correlations (normalised to the total number of the selected four-jet events), leaving the ratios of the colour factors as free parameters. Due to the normalization, the leading order prediction of these distributions is independent of the strong coupling. However, the independence of the leading order prediction from the strong coupling does not imply that the radiative corrections are negligible. In fact, in Ref.  the value of $`T_R/C_F`$ (see the Appendix) was estimated to increase about 25% if next-to-leading order (NLO) theoretical predictions are used instead of the leading order (LO) ones. It is therefore desirable to explicitly verify the effect of the next-to-leading order corrections on these normalised angular distributions.

Beyond leading order in the theoretical predictions, the strong coupling and the colour factors are interdependent. On the one hand, even the normalised next-to-leading order predictions for the angular distributions depend upon $`\alpha _s`$ , while on the other hand, $`\alpha _s`$ itself depends on the colour factors through its running. In this sense, beyond leading order, the only consistent way to determine the three free parameters is a simultaneous measurement of the strong coupling and the colour factors.

Predictions beyond leading order are made possible by theoretical developments achieved in the last few years. For multi-jet rates as well as numerous event shape distributions with perturbative expansions starting at O($`\alpha _s`$), matched next-to-leading order and next-to-leading logarithmic approximations provide very precise descriptions of the data over a wide range of the available kinematic region . Also, recently, the next-to-leading order predictions for the distributions of four-jet angular correlations have been calculated . To make use of these developments, we perform a simultaneous fit of the strong coupling and ratios of colour factors using next-to-leading order predictions of four-jet rates, differential two-jet rates and four-jet angular correlations. The data were collected by the OPAL Collaboration at LEP.

The outline of the paper is as follows: In Section 2, we present the observables used in the analysis and describe their best available perturbative predictions. The analysis procedure is explained in detail and the fit results are presented in Section 3. Section 4 contains the discussion of the systematic checks which were performed and the corresponding systematic errors. We collect our results and compare them to the results of previous studies in Section 5. Section 6 contains our conclusions.

## 2 Observables

To perform a simultaneous measurement of the strong coupling and the colour factors we use multi-jet related variables to gain sensitivity to $`\alpha _s`$ and four-jet angular correlations to gain sensitivity to the colour factors.

There are many ways of defining jets. For the present analysis we use the Durham scheme . Starting by defining each particle to be an individual jet, a resolution variable $`y_{ij}`$ is calculated for each pair of jets $`i`$ and $`j`$:

$$y_{ij}=\frac{2\mathrm{min}(E_i^2,E_j^2)(1\mathrm{cos}\theta _{ij})}{E_{\mathrm{vis}}^2},$$
(1)
where $`E_i`$ and $`E_j`$ are the energies of jets $`i`$ and $`j`$, $`\theta _{ij}`$ is the angle between them and $`E_{\mathrm{vis}}`$ is the sum of the energies of the visible particles in the event or of the partons in a theoretical calculation. If the smallest value of $`y_{ij}`$ is smaller then a predefined value $`y_{\mathrm{cut}}`$, the pair is replaced by a pseudo-jet with four momentum $`p_{ij}^\mu =p_i^\mu +p_j^\mu `$ and the clustering starts again with momenta $`p_i^\mu `$ and $`p_j^\mu `$ dropped and $`p_{ij}^\mu `$ added to the final state. Clustering ends when the smallest value of $`y_{ij}`$ is larger than $`y_{\mathrm{cut}}`$.

In our analysis we use the differential two-jet rate,

$$D_2(y_{23})\frac{1}{\sigma _{\mathrm{tot}}}\frac{\mathrm{d}\sigma }{\mathrm{d}y_{23}},$$
(2)
where $`y_{23}`$ is the $`y_{\mathrm{cut}}`$ value for which the two- and three-jet configurations are separated in a given event, and the four-jet rate,

$$R_4(y_{\mathrm{cut}})\frac{\sigma _{4\mathrm{jet}}(y_{\mathrm{cut}})}{\sigma _{\mathrm{tot}}},$$
(3)
to gain sensitivity to $`\alpha _s`$. For theory, $`\sigma _{\mathrm{tot}}`$ is the total hadronic cross section at O($`\alpha _s`$) accuracy. For experiment, $`\sigma _{\mathrm{tot}}`$ is the total visible hadronic cross section. The $`D_2`$ distribution has been extensively used at LEP to determine $`\alpha _s`$ . Use of $`R_4`$ allows us to obtain sensitivity to $`\alpha _s`$ using a four-jet based variable. For these observables both next-to-leading order (NLO)  and next-to-leading logarithmic (NLL)  perturbative predictions are known. The former provide a fixed order approximation valid for large values of the resolution variable, while the latter provide an all order summed but approximate cross section valid for small values of the resolution variable. The explicit NLO and NLL formulae and their matched expressions used in our fits are given in the Appendix. Using those formulae, we obtain various theoretical predictions for the four-jet rates ($`R_4`$) and the differential two-jet rates ($`D_2`$) which are shown in Fig. 1. In these plots, the theoretical predictions with the exception of the fitted curves were obtained using the world average of the strong coupling $`\alpha _s(M_Z)=0.118`$  and standard QCD values for the colour factors, $`C_F=4/3`$ and $`C_A=3`$. The theoretical predictions in Fig. 1 are shown in comparison to OPAL data corrected to the parton level. The manner in which the data are corrected to the parton level is explained in Section 3.3. It is seen that the theoretical predictions provide a good description of the data. For the theoretical plots in Fig. 1, we used the $`R`$-matched expression for $`R_4`$ and $`\mathrm{ln}R`$-matched expression for $`D_2`$ (see the Appendix).

The other class of observables we employ in our measurements is the four-jet angular correlations. The angular variables are

* the Bengtsson-Zerwas angle : $`\chi _{BZ}=\mathrm{}[(\stackrel{}{p}_1\times \stackrel{}{p}_2),(\stackrel{}{p}_3\times \stackrel{}{p}_4)]`$;
* the modified Nachtmann-Reiter angle : $`\mathrm{\Theta }_{NR}=\mathrm{}[(\stackrel{}{p}_1\stackrel{}{p}_2),(\stackrel{}{p}_3\stackrel{}{p}_4)]`$;
* the Körner-Schierholtz-Willrodt angle :
$`\mathrm{\Phi }_{KSW}=\frac{1}{2}\left(\mathrm{}[(\stackrel{}{p}_1\times \stackrel{}{p}_4),(\stackrel{}{p}_2\times \stackrel{}{p}_3)]+\mathrm{}[(\stackrel{}{p}_1\times \stackrel{}{p}_3),(\stackrel{}{p}_2\times \stackrel{}{p}_4)]\right)`$;
* the angle between the two lowest energy jets : $`\alpha _{3,4}=\mathrm{}[\stackrel{}{p}_3,\stackrel{}{p}_4]`$,

where the $`\stackrel{}{p}_i`$ denote the three-momenta of the energy-ordered jets $`(E_1>E_2>E_3>E_4)`$. We defined the jets using the Durham clustering and selected four-jet events at $`y_{\mathrm{cut}}=0.008`$. In order to have a high statistics four-jet sample, $`y_{\mathrm{cut}}`$ should be chosen in the range from $`10^2`$ to $`10^3`$. However, we cannot choose $`y_{\mathrm{cut}}`$ to be much smaller than 0.01 for technical reasons (see Section 3.3.1). For the four-jet events, used in obtaining the angular correlations, an additional cut was placed on the energy of each jet: $`E>E_{\mathrm{min}}=3`$GeV.

To obtain as much information as possible, we used all four angular variables in our analysis. The shape of the distributions of these angular correlations at LO and NLO are hardly distinguishable. Fig. 2 shows the normalised NLO theoretical predictions fitted to the data.

Although the various angular correlations are not entirely independent, the correlations (see Section 3.4) are not large as can be seen from Table 1, where the maximum values of the magnitudes of the bin-wise correlations between the observables are given. We also performed fits using a single angular variable together with one of the jet-related variables and we observed that these fits give very scattered central values. However, in these cases we found very strong correlations between the fitted parameters ($``$ 1), indicating that the distributions do not contain sufficient information to constrain all of the QCD parameters simultaneously.

## 3  Analysis Procedure and Results

### 3.1  The OPAL detector

A detailed description of the OPAL detector can be found in Ref. . This analysis relies mainly on the reconstruction of charged particle trajectories and on the measurement of energy deposited in the electromagnetic and hadronic calorimeters.

The central detector contains a silicon micro-vertex detector and three drift chamber devices: an inner vertex chamber, a large jet chamber and a surrounding $`z`$-chamber. The central detector is located inside a solenoidal magnet which provides a uniform axial222In the OPAL coordinate system the $`+x`$ axis points towards the centre of the LEP ring, the $`y`$ axis points approximately upwards and the $`z`$ axis points in the direction of the electron beam. The polar angle $`\theta `$ and the azimuthal angle $`\varphi `$ are defined with respect to $`z`$ and $`x`$, respectively, while $`r`$ is the distance from the $`z`$-axis. magnetic field of approximately 0.435 T. Tracking of the charged particles is performed with this detector. Most of the tracking information is obtained from the jet chamber which provides up to 159 measured space points per track, and close to 100 % tracking efficiency in the region $`|\mathrm{cos}\theta |<0.98`$. The average angular resolution is approximately 0.1 mrad in $`\varphi `$ and 1 mrad in $`\theta `$.

Electromagnetic energy is measured by lead glass calorimeters surrounding the magnet coil, separated into a barrel ($`|\mathrm{cos}\theta |<0.82`$) and two end-cap ($`0.81<|\mathrm{cos}\theta |<0.98`$) sections. The electromagnetic calorimeter consists of 11 704 lead glass blocks with a depth of 24.6 radiation lengths in the barrel and more than 22 radiation lengths in the end-caps. This calorimeter is surrounded by the hadronic calorimeter of the sampling type measuring the energy of hadrons emerging from the electromagnetic calorimeter.

### 3.2  Data samples

The analysis presented in this paper is based on 4.1 million hadronic events recorded within 3 GeV of the $`\mathrm{Z}`$ peak by the OPAL detector between 1991 and 1995. The descriptions of the OPAL trigger system and the offline multihadronic event selection are given in Refs.  and . In this analysis, additional criteria were applied to eliminate poorly measured tracks and obtain well contained events.

Each track was required to have (i) at least 40 measured points in the jet chamber, (ii) transverse momentum in the $`r`$$`\varphi `$ plane greater than 0.15 GeV/$`c`$, (iii) measured momentum less than 60 GeV/$`c`$, (iv) $`|\mathrm{cos}\theta |<0.97`$, (v) distance of closest approach to the origin in the $`r`$$`\varphi `$ plane of no more than 5 cm, (vi) and less than 25 cm in the $`z`$ direction. Hadronic events were required to contain at least five tracks to reduce contamination from $`e^+e^{}\tau ^+\tau ^{}`$ events.

Clusters in the electromagnetic calorimeter were accepted if they had more than 0.1 GeV energy in the barrel section or more than 0.25 GeV in the end-cap section. The corresponding clusters were required to span at least two lead glass blocks. Clusters in the hadronic calorimeter were required to have energies larger than 1 GeV.

For the calculation of the visible energy, $`E_{\mathrm{vis}}`$, the tracks are assumed to have the pion mass, and clusters are treated as photons. To determine the energy of each cluster measured in the electromagnetic and hadronic calorimeters, we used a matching algorithm . This algorithm reduces double counting of energy and gives a better resolution both in energy and angle than the calorimeters alone. We did not impose energy-momentum conservation, but checked that such a constraint does not affect our results.

The event thrust axis was determined using the accepted tracks and clusters. In order that the event be well contained, we required $`|\mathrm{cos}\theta _{\mathrm{thrust}}|<0.9`$.

The final event sample contained about 3.6 million hadronic $`\mathrm{Z}`$ decays from which about 250 000 were identified as four-jet events.

### 3.3 Data Correction

We choose as the theoretical reference the distributions given in Eqs. (22), (28) of the Appendix which are obtained from a parton level calculation. In order to compare our (detector level) data to these (parton level) theoretical predictions the measured distributions were corrected for the effects of the detector and hadronization.

First, we combined the distributions of the six variables $`\mathrm{cos}\chi _{\mathrm{BZ}},\mathrm{cos}\mathrm{\Theta }_{\mathrm{NR}},\mathrm{cos}\mathrm{\Phi }_{\mathrm{KSW}},\mathrm{cos}\alpha _{34},D_2`$ and $`R_4`$ into a one-dimensional distribution, with 128 bins, 20 for each angular correlation and 24 each for the four-jet rates (covering the range of $`0.001<y_{cut}<1.0`$) and the differential two-jet rates (in the range of $`1.0<\mathrm{ln}y_{23}<5.8`$). Next, we applied bin-by-bin corrections for detector distortions and hadronization, as described below. The correction factors were obtained from the same distributions computed from Monte Carlo events.

#### 3.3.1 Monte Carlo simulation

The correction factors were estimated using two different modes of the JETSET Monte Carlo program : the parton shower mode and the matrix element mode. In the parton shower approach, the evolution of the parton system is treated as a branching process based on the leading logarithmic approximation starting from the original $`\mathrm{q}\overline{\mathrm{q}}`$ pair. In the matrix element mode, up to four partons are generated according to the second-order ERT calculation .

It is well known that the parton shower versions of the available Monte Carlo programs provide a good description of jet rates, but cannot describe the structure of the four-jet events properly . Therefore, to determine correction factors for the four-jet angular correlations, we used Monte Carlo events with only four partons ($`\mathrm{q}\overline{\mathrm{q}}`$gg and $`\mathrm{q}\overline{\mathrm{q}}`$$`\mathrm{q}\overline{\mathrm{q}}`$) in the final state. These events were generated using the second-order matrix element (ME) calculation, subjected to hadronization, as implemented in JETSET version 7.4 assuming standard QCD colour factors and the set of parameters ERT-MC-1 of Ref. . 560 000 Monte Carlo events were generated. To avoid singularities in the four-parton generation of the LO matrix elements, an intrinsic cutoff $`y_0`$ has to be applied in the Monte Carlo simulation. In JETSET, $`y_0=0.01`$ is used, calculated according to the invariant mass resolution definition, which for four-jet events corresponds to a Durham resolution of about 0.004 (yielding the same four-jet rates). Our resolution criterion of $`y_{\mathrm{cut}}=0.008`$ is safely larger than the intrinsic cut-off. For normalised angular correlations, the distributions at leading and next-to-leading order are almost the same . Thus for the purpose of determining the corrections the ERT Monte Carlo events are sufficient.

For the correction of the differential two-jet rates and the four-jet rate, we used 2 million events generated with the parton shower (PS) version of JETSET, using parameters tuned to OPAL data at $`\sqrt{s}=`$91 GeV described in Ref. . For both models, hadronization of the parton system is based on the Lund string-fragmentation scheme . The events of both types of samples were passed through the OPAL detector simulation , processed by the same reconstruction programs, and subjected to the same event selection criteria as the data.

#### 3.3.2 Correction procedure

To correct the data for the finite acceptance and resolution of the detector, we determined bin-by-bin correction factors for every bin of the distribution from Monte Carlo samples with hadronization included, both before (hadron level, $`H_{i}^{}{}_{}{}^{\mathrm{MC}}`$) and after (detector level, $`D_{i}^{}{}_{}{}^{\mathrm{MC}}`$) the detector simulation and selection cuts. The hadron level treats all charged and neutral particles with lifetimes greater than $`3\times 10^{10}`$ s as stable and has no initial state photon radiation. The correction factor for detector effects for the $`i^{th}`$ bin of the distribution is given by:

$$C_{i}^{}{}_{}{}^{\mathrm{det}}=\frac{H_{i}^{}{}_{}{}^{\mathrm{MC}}}{D_{i}^{}{}_{}{}^{\mathrm{MC}}}.$$
(4)
We examined the effects of initial state photon radiation by switching the radiation on in the hadron level Monte Carlo simulation. The correction factors were found to be identical within statistical uncertainties with or without photon radiation over the entire phase space.

We estimated the hadronization effects by computing the distributions of Monte Carlo events at the parton ($`P_{i}^{}{}_{}{}^{\mathrm{MC}}`$) and hadron levels ($`H_{i}^{}{}_{}{}^{\mathrm{MC}}`$) (before and after hadronization) and defined the correction factor for hadronization by:

$$C_{i}^{}{}_{}{}^{\mathrm{had}}=\frac{P_{i}^{}{}_{}{}^{\mathrm{MC}}}{H_{i}^{}{}_{}{}^{\mathrm{MC}}}.$$
(5)
Using the two correction factors $`C_{i}^{}{}_{}{}^{\mathrm{det}}`$ and $`C_{i}^{}{}_{}{}^{\mathrm{had}}`$, we correct the measured distribution to the parton level according to

$$D_{i}^{}{}_{}{}^{\mathrm{corr}}=C_{i}^{}{}_{}{}^{\mathrm{det}}C_{i}^{}{}_{}{}^{\mathrm{had}}D_{i}^{}{}_{}{}^{\mathrm{meas}}.$$
(6)
We show the distributions of the total correction factors, $`C_{i}^{}{}_{}{}^{\mathrm{tot}}=C_{i}^{}{}_{}{}^{\mathrm{det}}C_{i}^{}{}_{}{}^{\mathrm{had}}`$, in Figs. 1 and 2.

### 3.4 Fit procedure

Having prepared the corrected distribution of six variables, we performed a $`\chi ^2`$ minimization to determine the most probable values of the variables $`\eta =\alpha _sC_F/(2\pi )`$, $`x=C_A/C_F`$ and $`y=T_R/C_F`$ (see Section A) with the program MINUIT . The function $`\chi ^2`$ was defined by:

$`\chi ^2=_{ij}\delta _i\sigma _{ij}^1\delta _j`$ ,

where $`\delta _i`$ is the bin-wise difference between the theoretical QCD prediction and the corrected data distribution, while $`\sigma _{ij}`$ is the covariance matrix taking into account the statistical error of the data, the finite Monte Carlo statistics, and the correlations between the bins of the distributions. The data and Monte Carlo event samples were each divided into 90 subsamples333The number of subsamples is a result of optimization for the given data and number of bins. to compute $`\sigma _{ij}`$, using the standard formula:

$$\sigma _{ij}=\frac{1}{N(N1)}\underset{n=1}{\overset{N}{}}\left(D_i^{(n)}\overline{D}_i\right)\left(D_j^{(n)}\overline{D}_j\right),$$
(7)
where $`D_i^{(n)}`$ is the value of the corrected distribution in bin $`i`$ of the $`n`$th subsample, $`\overline{D}_i`$ is the value in that bin for the whole sample and $`N`$ is the number of subsamples. We included bins in the fit only if the total correction factors were in the range $`0.9<C_i^{\mathrm{tot}}<1.1`$, except for the four-jet rate, where $`0.85<C_{i}^{}{}_{}{}^{\mathrm{tot}}<1.15`$ was used.444There are not enough bins in the smaller range to obtain a reliable fit. This bin choice is motivated by our expectation that the Monte Carlo estimate is less reliable if hadronization and/or detector distortions are large. The total number of bins used in the fit was 82, leaving 79 degrees of freedom since there are three fitted parameters ($`\eta `$, $`x`$ and $`y`$, see Section 3.5).

### 3.5 Fit results

The central values of the fit results for the parameter $`\eta `$ and colour factor ratios $`x`$ and $`y`$ are listed in the upper section of Table 2. The renormalization scale was fixed at $`x_\mu =1`$. The errors and correlations ($`\rho `$) in the table are purely statistical.

We observe a strong correlation between $`x`$ and $`y`$. We tested other variables (for instance, the three-jet rate) in an attempt to decrease these correlations, but did not observe a significant improvement. The correlations were taken into account when the fit results were converted to the standard QCD parameters $`\alpha _s(M_Z)`$ and colour factors using the definitions in Eqs. (10) and (12) and $`T_R=1/2`$. These latter results are given in the lower section of Table 2.

## 4  Systematic uncertainties

The systematic uncertainties were evaluated by considering the following effects:

* the measurement process and accuracy of the Monte Carlo detector simulation,
* dependence on the model of hadronization,
* theoretical uncertainties arising from the choice of the renormalization scale,
* variation of the matching scheme in the case of the differential two-jet distribution,
* variation of the fit range,
* the background from five parton events,
* variation of the parameter $`y_{\mathrm{cut}}`$.

In the systematic checks we always used the same covariance matrix that was optimized for the standard analysis.

### 4.1 Detector effects

The uncertainty related to the detector simulation and the measurement process was estimated by repeating the analysis using either tracks only or clusters only for the data and the detector level Monte Carlo samples. For the tracks and clusters we applied the same selection procedure as described in Section 3.2, except that all cluster related criteria were left out in the former case and all track related criteria were left out in the latter. Furthermore, the entire analysis was repeated using $`|\mathrm{cos}\theta _{\mathrm{thrust}}|<0.7`$ (rather than $`|\mathrm{cos}\theta _{\mathrm{thrust}}|<0.9`$), or $`N_{ch}8`$ (rather than $`N_{ch}5`$). The results obtained for these analyses are given in Table 3.

### 4.2 Hadronization models

Several combinations of different Monte Carlo models were considered to estimate the uncertainty associated with the hadronization correction. Since we used two different Monte Carlo samples to determine the correction factors (ME and PS), we examined the effects induced by variations of these models separately.

The angular correlations were corrected using the ME version of JETSET with the Lund string fragmentation scheme, as stated in Section 3.3.1. In the systematic studies we used the ME version of the HERWIG generator , which produces four-parton final states according to the leading order QCD matrix elements, and then a parton shower is started from each of the four partons followed by hadronization based on the cluster fragmentation.555The HERWIG-ME has not been tuned to OPAL data. We used CLSMR(1)=0.4, PSPLT(2)=0.33, DECWT=0.7, and parameters not listed here were left at their default values. We also varied the principal hadronization parameters of the JETSET ME Monte Carlo program, $`\sigma _q`$ and $`a`$. As the hadronization parameter $`b`$ is strongly correlated with $`a`$ we kept it fixed. We increased $`a`$ or decreased $`\sigma _q`$, alternately, by 10 %. These changes correspond to about one standard deviation in the parameter values allowed for the JETSET ME generator at $`Z`$ peak, as found by the LEP experiments . The fits were redone using data corrected by the factors calculated from these Monte Carlo samples at the generator level (reevaluating $`C_i^{\mathrm{had}}`$ only) and using the standard $`C_i^{\mathrm{det}}`$ correction factors for detector effects.

The effect of the hadronization process on the four-jet rates and differential two-jet rates was studied by using the PS version of HERWIG instead of the PS version of JETSET and by varying the Lund hadronization parameters $`\sigma _q`$ and $`a`$ as for the angular correlations. The parameter $`Q_0`$ defining the end of the parton shower cascade was also varied from its standard value of $`1.9`$ GeV to $`2.4`$ GeV. The theoretical predictions assume zero quark masses. To assess the effect of this assumption, we determined the hadronization correction, $`C^{\mathrm{had}}`$, using JETSET Monte Carlo events with only light primary quarks $`(u,d,s,c)`$ at the parton level but with all available flavours $`(u,d,s,c,b)`$ at the hadron level. The results provided by these variations are presented in Table 4.

### 4.3 Renormalization scale ambiguity

The arbitrariness in the choice of the renormalization scale $`\mu `$ leads to an ambiguity in the theoretical prediction. We perform our default analysis at the scale factor $`x_\mu =\mu /E_{\mathrm{cm}}=1`$.

To estimate the systematic error due to the scale ambiguity, we varied the renormalization scale within the range $`0.5x_\mu 2`$. Varying $`x_\mu `$ probes the effect of missing higher orders, and these could be very different for the NLO and NLLA calculations. We performed two variations separately: varying $`x_\mu `$ for the jet-related variables while keeping it fixed at the default value for the four-jet angular correlations, or vice versa. Table 5 presents the results of this systematic check. As expected, the normalised angular correlations were not sensitive to the scale choice, but the jet related variables showed significant sensitivity. In Fig. 3 we show the dependence of the $`\chi ^2`$ of the fit on $`x_\mu `$, for the case of jet-related variables with $`x_\mu =1`$ kept fixed for the angular correlations. The minimum of the $`\chi ^2`$ is seen to lie close to $`x_\mu =1`$.

### 4.4 Matching ambiguity

For the differential two-jet rates, we performed the fit using the $`R`$-matching scheme to combine the NLL and NLO approximations rather than the $`\mathrm{ln}R`$-matching. The result of this variation is listed in Table 5.

### 4.5 Fit range

The uncertainty related to the number of bins selected for the fit was evaluated by repeating the fit using those bins for which the values of the correction factors were in the range of $`0.8<C_{i}^{}{}_{}{}^{\mathrm{tot}}<1.2`$ instead of the default range $`0.9<C_{i}^{}{}_{}{}^{\mathrm{tot}}<1.1`$ or $`0.85<C_{i}^{}{}_{}{}^{\mathrm{tot}}<1.15`$. The results are presented in Table 5. The variation of the fit range for the jet-related variables does not influence the value of the $`\chi ^2`$ significantly. The fit-range variation for the angular variables $`\mathrm{cos}\chi _{\mathrm{BZ}}`$ and $`\mathrm{cos}\alpha _{34}`$ is the main source of the much increased value of the $`\chi ^2`$.

### 4.6 Five-parton background

Five-parton events which are reconstructed as four-jet events at the parton level are taken into account by the NLO theory. We examined the effect of those five-parton events which yield four jets at the detector level but five at the parton level. To reduce this five-parton background, we applied a cut on $`y_{45}`$, the value of the resolution variable which separates the four- and five-jet configurations of any given event. We optimized the value of this cut by studying five-parton events generated with the PENTAJET Monte Carlo program . This program generates a five-parton configuration using the same value of the intrinsic $`y_0`$ as the JETSET-ME program, then hadronizes the partons using JETSET string fragmentation with the ERT-MC-3 parameters of Ref. . The events were passed through the OPAL detector simulation program. Of 20 000 five-parton events, 2 517 were found to be four-jet events at detector level, but five-jet events at parton level, for $`y_{\mathrm{cut}}=0.008`$. The number of these events reduces to 376 when we keep events with $`y_{45}<0.003`$, eliminating 85 % of the five-parton background, while leaving 60 % of the data.

The results of the fit with this additional requirement are presented in Table 5.

### 4.7 Dependence on the parameter $`y_{\mathrm{cut}}`$

The value of the jet resolution parameter $`y_{\mathrm{cut}}`$, used in selecting four-jet events to define the angular correlations, can be chosen arbitrarily, as long as it is safely larger than the intrinsic $`y_0`$ of the JETSET-ME program. In order to have a high statistics sample, our default choice was the smallest possible value, which allows us to examine the stability of the results when $`y_{\mathrm{cut}}`$ is increased. The four-jet event sample decreases very rapidly with increasing $`y_{\mathrm{cut}}`$, which constrains the region in $`y_{\mathrm{cut}}`$ where a reliable systematic check can be performed. We repeated the entire analysis using $`y_{\mathrm{cut}}=0.01`$. The results of this fit are presented in Table 5.

### 4.8 Additional systematic checks

We tested the influence of the two- and three-parton events that give rise to four-jet events at the detector level by repeating the analysis requiring $`E_{\mathrm{min}}=5`$GeV on the minimum energy of the jets instead of the default value $`E_{\mathrm{min}}=3`$GeV. A higher value of the jet energy reduces the number of two- and three-parton events which satisfy the four-jet selection criteria in the detector. We generated 1 million events with two, three and four partons in the initial state using the ME version of JETSET with renormalization scale set to $`x_\mu ^2=0.002`$, which provides the correct proportion of two, three and four jet events for $`y_{\mathrm{cut}}=0.008`$. We then counted the proportion of the two- and three-parton events which passed the four-jet selection at detector level. With $`E>E_{\mathrm{min}}=5`$GeV, there are no background events from the two- and three-parton events. With $`E>E_{\mathrm{min}}=3`$GeV, a background of 1 % appears (519 two- or three-parton events in the sample of 53 302 four-jet events selected at the detector level). We found negligible effects on the final results of the fit, therefore, we did not include the effect of this check in the systematic errors.

As stated in Section 2, the fits to individual angular correlations give rather ambiguous results. Therefore, we choose to use all four angular correlation variables in our standard analysis. As a consistency check, we performed the fit using only five variables instead of six, removing one of the six variables at a time. The results of these fits are presented in Table 6. The errors are statistical only. We found that leaving out any one of the observables yields consistent results for the measured parameters within the statistical errors of the measurements. Therefore, we did not include the effect of this check in the systematic errors.

The number of light quark flavours in a fixed order calculation is usually set to $`N_f=5`$ at LEP energies because the hard scattering scale is assumed to be the $`e^+e^{}`$ center of mass energy. However, one may argue that in the four-quark channel the $`gb\overline{b}`$ splitting is suppressed compared to the $`gq\overline{q}`$ (q=u, d, s, c) splittings, because the scale of that process is much lower than the scale at the primary vertex. To study the effect of suppressed $`gb\overline{b}`$ splitting on our results, we also used $`N_f=4`$ light flavours at the secondary vertex in the fixed order prediction for the angular correlations. The results of this check are also shown in Table 6. The colour factor ratios increase slightly, remaining well within the statistical errors. The quality of the fit is almost unchanged. We did not include the effect of this check in the systematic errors, because it is not customary to do so.

### 4.9 Total systematic error

To combine the effects of the systematic variations we employed a Bayesian method that was developed in a similar analysis by the ALEPH collaboration , where the explicit formulae we used can be found. The basic idea is that the information obtained from a fit is accepted or rejected based on the quality of the fit, i.e. the magnitude of the $`\chi ^2`$. From a Bayesian point of view we presume that all models are equally well suited to the analysis. However, a large $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.`$ indicates that the probability of a model is low and, therefore, it should contribute to the systematic errors with a small weight. Accordingly, the systematic error corresponds to an increase of the $`\chi ^2`$ by 1.

We could not use the Bayesian method in two cases. In the variation of the fit range the number of bins changed, while in checking the renormalization scale uncertainty for the angular correlations we did not find a pronounced minimum of the $`\chi ^2`$ curve. Furthermore we did not use the Bayesian method to estimate the systematic uncertainty from the five-parton background because the event statistics for this check were modified significantly, by about 40%, relative to the standard analysis. In these cases we estimated the systematic error by the deviation of the results from the standard values.

The errors obtained from the different systematic checks are summarized in Table 7. The theoretical uncertainty refers to the larger of the errors from the scale uncertainty of the jet-related variables, the matching ambiguity (both of those estimate the effect of unknown higher order perturbative contributions), and the scale uncertainty of the four-jet angular variables. The different sources were added in quadrature to define the total systematic error.666If the Bayesian approach were not adopted and the systematic errors were estimated by the largest deviation from the standard values, then we obtain the following systematic errors: $`\mathrm{\Delta }\eta =0.0027,\mathrm{\Delta }x=0.27,\mathrm{\Delta }y=0.08`$. For the basic QCD parameters these errors are: $`\mathrm{\Delta }\alpha _s=0.030,\mathrm{\Delta }C_A=0.72,\mathrm{\Delta }C_F=0.31`$.

Taking into account the systematic correlations presented in Table 7, we converted these values to the errors of the standard QCD parameters $`\alpha _s(M_Z)`$ and colour factors using the definitions in Eq. (10) and (12) and $`T_R=1/2`$. The result of the conversion is presented in Table 7. The errors are dominated by the theoretical uncertainty.

## 5 Results and Discussion

The procedure outlined above yields

$$C_A/C_F=2.25\pm 0.08(\mathrm{stat}.)\pm 0.14(\mathrm{syst}.),T_R/C_F=0.37\pm 0.04(\mathrm{stat}.)\pm 0.06(\mathrm{syst}.).$$

Repeating the entire analysis with $`\alpha _s(M_Z)`$ fixed to its world average value 0.1184  leads to similar results:

$$C_A/C_F=2.29\pm 0.06(\mathrm{stat}.)\pm 0.14(\mathrm{syst}.),T_R/C_F=0.38\pm 0.03(\mathrm{stat}.)\pm 0.06(\mathrm{syst}.).$$

To compare the results to previous measurements, we show in Fig. 4 the two dimensional 68% and 95% confidence level contour plots of the colour factor ratios based on total uncertainties obtained by our analysis.

The DELPHI Collaboration  has performed a least-squares fit of leading order predictions to the two-dimensional distribution in the variables $`\mathrm{cos}\mathrm{\Theta }_{\mathrm{NR}}`$ and $`\mathrm{cos}\alpha _{34}`$ in order to measure of the colour factor ratios. The DELPHI Collaboration  also performed this analysis with b-tagging to separate quark jets from gluon jets which increases sensitivity. The OPAL Collaboration fitted the LO prediction to a three-dimensional distribution of the $`\mathrm{cos}\mathrm{\Theta }_{\mathrm{NR}}`$ , $`\mathrm{cos}\alpha _{34}`$ and $`\mathrm{cos}\chi _{\mathrm{BZ}}`$ variables . In a second measurement the OPAL Collaboration used event shape distributions of LEP1 data for fitting $`\alpha _s`$ and one of the colour factors . Furthermore, another recent measurement of the QCD colour factors and $`\alpha _s`$ using event shape distributions at several energy points and power corrections was carried out . The ALEPH Collaboration  was the first to present a simultaneous measurement of the strong coupling and the colour factors. The distribution $`D_2`$ and all four angular distributions defined in Section 2 have been employed. The results of the ALEPH experiment were determined using leading order predictions for the four-jet observables. The systematic uncertainties include contributions from variation of the renormalization scale, the matching ambiguity, from hadronization and detector effects, and finally from the estimation of mass effects. Our results are in agreement with these previous results. The main new feature of our analysis is that we obtained our results using NLO predictions for four-jet angular correlations. Furthermore, our systematic checks contain additional items arising from the estimate of the effects of the five-parton background and the variation of the $`y_{\mathrm{cut}}`$ parameter. We also checked that assuming four rather than five light flavours at the secondary vertex does not influence our results. This effect was not examined in any of the previous studies.

With the inclusion of the higher order corrections, the colour factor ratios became better constrained, especially $`T_R/C_F`$. We found an increase of about 15% on $`T_R/C_F`$ using the NLO theoretical predictions instead of the LO one, while $`C_A/C_F`$ was only slightly affected (about 3%).

The diamond symbol in Fig. 4 indicates the result one would expect if a leading order perturbative theory with light gluinos was the correct theory. This point falls well outside our 95% C.L. contour, indicating that the data do not favour the existence of light gluinos. A complete analysis to test for the existence of gluinos should use hadronization models and NLO predictions that include the effect of the gluinos . The latter is possible using the DEBRECEN partonic Monte Carlo program we used for generating the NLO results (see the Appendix), but presently none of the existing hadronization models include gluinos.

Finally, we converted our measured parameters to the standard QCD parameters (using $`T_R=1/2`$), which leads us to our main results:

$`\alpha _s=0.120\pm 0.011(\mathrm{stat}.)\pm 0.020(\mathrm{syst}.)`$
$`C_A=3.02\pm 0.25(\mathrm{stat}.)\pm 0.49(\mathrm{syst}.)`$
$`C_F=1.34\pm 0.13(\mathrm{stat}.)\pm 0.22(\mathrm{syst}.)`$
Other techniques  provide precise determinations of $`\alpha _s`$ . In this study the measurement of $`\alpha _s`$ serves as a crosscheck and a consistent value was found. The larger uncertainty we obtain for $`\alpha _s`$ compared to other techniques can be traced to the functional dependence on the other two parameters.

## 6 Summary

In this paper, a test of perturbative QCD at LEP has been presented. With only the assumption of non-abelian gauge symmetry and standard hadronization models, we measured the eigenvalues of quadratic Casimir operators — the colour factors — of the underlying gauge group, together with the coupling of QCD. The measurement was based on the comparison of next-to-leading order perturbative predictions for four-jet angular variables, and the resummation improved next-to-leading order perturbative predictions for the multi-jet rates, to OPAL data taken between 1991 and 1995 at the $`\mathrm{Z}`$ peak. The results for the colour factors are as follows:

$$C_A=3.02\pm 0.25(\mathrm{stat}.)\pm 0.49(\mathrm{syst}.),C_F=1.34\pm 0.13(\mathrm{stat}.)\pm 0.22(\mathrm{syst}.),$$

from a global fit to the angular correlations, four-jet rates, and differential two-jet rates. These values are in agreement with SU(3) values of $`C_A=3`$ and $`C_F=4/3`$ as well as with measured values obtained previously at LEP. The corresponding strong coupling value,

$$\alpha _s=0.120\pm 0.011(\mathrm{stat}.)\pm 0.020(\mathrm{syst}.),$$

is in agreement with the world average.

## Appendix A Appendix: Theoretical Predictions

The NLO perturbative expansion of the differential two-jet rates has the following expression:

$$D_2(y_{23})\frac{1}{\sigma _{\mathrm{tot}}}\frac{\mathrm{d}\sigma }{\mathrm{d}y_{23}}=\eta (\mu )A_{D_2}^{}(y_{23})+\eta (\mu )^2\left[B_{D_2}^{}(y_{23})+\beta _0\mathrm{ln}(x_\mu )A_{D_2}^{}(y_{23})\right],$$
(8)
where $`y_{23}`$ is the $`y_{\mathrm{cut}}`$ value which separates the two- and three-jet configurations of an event. The corresponding formula in the case of the four-jet rates takes a similar form:

$$R_4(y_{\mathrm{cut}})\frac{\sigma _{4\mathrm{jet}}}{\sigma _{\mathrm{tot}}}=\eta (\mu )^2B_{R_4}(y_{\mathrm{cut}})+\eta (\mu )^3\left[C_{R_4}(y_{\mathrm{cut}})+\beta _0\mathrm{ln}(x_\mu ^2)B_{R_4}(y_{\mathrm{cut}})\right].$$
(9)
In Eqs. (8) and (9) $`\sigma _{\mathrm{tot}}`$ is the total hadronic cross section at O($`\eta `$) accuracy and

$$\eta (\mu )\frac{\alpha _s(\mu )C_F}{2\pi }=\frac{\eta (M_Z)}{w(\mu ,M_Z)}\left(1\frac{\beta _1}{\beta _0}\eta (M_Z)\frac{\mathrm{ln}(w(\mu ,M_Z))}{w(\mu ,M_Z)}\right).$$
(10)
In the second equality in Eq. (10), we used the two-loop expression for the running coupling with

$$w(\mu ,M_Z)=1\beta _0\eta (\mu )\mathrm{ln}\left(\frac{M_Z}{\mu }\right),$$
(11)
and $`\mu `$ denotes the renormalization scale ($`x_\mu =\mu /\sqrt{s}`$, where $`\sqrt{s}`$ is the total c.m. energy). We define the colour factor ratios as

$$x=\frac{C_A}{C_F},y=\frac{T_R}{C_F},y_f=N_fy,$$
(12)
where $`N_f`$ is the number of active flavours. Then the coefficients in the expansion of the $`\beta `$ function can be written as

$$\beta _0=\frac{11}{3}x\frac{4}{3}y_f,\beta _1=\frac{17}{3}x^22y_f\frac{10}{3}xy_f.$$
(13)
The functions $`A_{D_2}^{}`$, $`B_{R_4}`$ and $`B_{D_2}^{}`$, $`C_{R_4}`$ in Eqs. (8) and (9) are the perturbatively calculable coefficient functions in the Born approximation and the radiative correction, respectively, which are independent of the renormalization scale. The $`A_{D_2}^{}`$ function is independent of $`x`$ and $`y`$. The $`B`$ functions are linear forms and the $`C_4`$ function is a quadratic form of these colour charge ratios :

$$B_O=B_0^{(O)}+B_x^{(O)}x+B_y^{(O)}y,O=D_2,\mathrm{or}R_4$$
(14)
$$C_{R_4}=C_0^{(4)}+C_x^{(4)}x+C_y^{(4)}y+C_z^{(4)}z+C_{xx}^{(4)}x^2+C_{xy}^{(4)}xy+C_{yy}^{(4)}y^2,$$
(15)
where

$$z=\frac{C_C}{N_cC_F^3},$$
(16)
with

$$C_C=\underset{a,b,c=1}{\overset{N_c^21}{}}\mathrm{Tr}(T^aT^bT^c)\mathrm{Tr}(T^cT^bT^a)$$
(17)
being a cubic Casimir invariant of the gauge group.777The results of the Monte Carlo integration for $`C_z^{(4)}`$ are compatible with zero , therefore, $`z`$ is not included as a free parameter in our fits. We calculated these coefficient functions of the differential two-jet rates in the $`6\mathrm{ln}(y_{23})0`$ $`y_{23}`$ range and those of the four-jet rates in the $`3\mathrm{log}_{10}(y_{\mathrm{cut}})0`$ $`y_{\mathrm{cut}}`$ range using the partonic Monte Carlo program DEBRECEN . All theoretical predictions were obtained for $`N_f=5`$ light-quark flavours and with the normalization $`T_R=\frac{1}{2}`$. The latter choice is arbitrary and influences the values of the coefficient functions in such a way that the physical prediction remains unchanged.

The resummation formulas for various event shapes are given in Ref.  at the NLL accuracy for the cumulative cross section,

$$R(y)=_0^ydy_{23}\frac{1}{\sigma }\frac{\mathrm{d}\sigma }{\mathrm{d}y_{23}},$$
(18)
in the following form:

$$R^{NLL}(L)=(1+C_1\eta +C_2\eta ^2)e^{Lg_1(\eta L)+g_2(\eta L)},$$
(19)
where $`L=\mathrm{ln}(y_{23})`$. The functions $`Lg_1(\eta L)`$ and $`g_2(\eta L)`$ represent the sums of the leading and next-to-leading logarithms, respectively, to all orders in the strong coupling. The colour decomposition of the $`g_1`$ and $`g_2`$ functions is given in Ref. . The coefficient $`C_1=5/2+\pi ^2/66\mathrm{ln}2`$ is known exactly  and $`C_2`$ can be determined from numerical integration of the $`\mathrm{O}(\alpha _\mathrm{s}^2)`$ matrix elements (see below).

At NLO (i.e. $`\mathrm{O}(\alpha _\mathrm{s}^2)`$ accuracy), the cumulative cross section for the differential two-jet rates has the form:

$$R^{\mathrm{O}(\alpha _\mathrm{s}^2)}(L)=1+A(L)\eta +B(L)\eta ^2.$$
(20)
In order to have the best possible theoretical description, we fit the experimental data to the matched NLL and NLO results, which is expected to describe the data in a wider kinematical range than any of the two separately. In the case of $`D_2`$, our default matching procedure is the $`\mathrm{ln}R`$-matching, which combines the NLL and NLO approximations according to the following formula:

$`\mathrm{ln}R(L)=\mathrm{\Sigma }(L,\eta )(G_{11}L+G_{12}L^2)\eta (G_{22}L^2+G_{23}L^3)\eta ^2`$ (21)
$`+A(L)\eta +\left(B_0(L)+xB_x(L)+yB_y(L){\displaystyle \frac{1}{2}}A(L)^2\right)\eta ^2.`$
where $`\mathrm{\Sigma }(L,\eta )=Lg_1(\eta L)+g_2(\eta L)`$ and the following two terms represent the two lowest order terms in the $`\eta `$-expansion of $`\mathrm{ln}R`$. In the actual analysis, we used the differential distribution $`\mathrm{d}\sigma /\mathrm{d}L`$ that is obtained from Eq. (21) as

$`{\displaystyle \frac{1}{\sigma }}{\displaystyle \frac{\mathrm{d}\sigma }{\mathrm{d}L}}=\mathrm{exp}(\mathrm{ln}R)[{\displaystyle \frac{\mathrm{d}\mathrm{\Sigma }(L,\eta )}{\mathrm{d}L}}(G_{11}+2G_{12}L)\eta (2G_{22}L+3G_{23}L^2)\eta ^2`$ (22)
$`+A^{}(L)\eta +(B_0^{}(L)+xB_x^{}(L)+yB_y^{}(L)A(L)A^{}(L))\eta ^2],`$
with $`A^{}(L)`$ and $`B^{}(L)`$ given in Eq. (8).

Eq. (22) contains the cumulative coefficients $`A(L)`$ and $`B_i(L)`$, the latter ones implicitly in $`R`$, that cannot be directly calculated using a partonic Monte Carlo program. We know, however, that

$$A(L)=C_1+\underset{m=1}{\overset{2}{}}R_{1m}L^m+D_1(L),$$
(23)
and

$$B_i(L)=C_2^{(i)}+\underset{m=1}{\overset{4}{}}R_{2m}^{(i)}L^m+D_2^{(i)}(L),i=0,x,y,$$
(24)
where the colour decomposition of the $`R_{2m}^{(i)}`$ logarithmic coefficients can be found in Ref. , except for $`R_{21}^{(i)}`$. This latter coefficient as well as the $`C_2^{(i)}`$ constants and the $`D_1(L)`$, $`D_2^{(i)}(L)`$ functions can be determined by differentiating Eqs. (23) and (24) with respect to $`L`$ and fitting to the coefficient functions as obtained from the partonic Monte Carlo program DEBRECEN. Following the fitting procedure of Refs.  we obtained $`R_{21}^{(0)}=16.531\pm 2.97(\mathrm{stat}.)`$, $`R_{21}^{(x)}=4.293\pm 1.03(\mathrm{stat}.)`$, $`R_{21}^{(y)}=4.54\pm 1.05(\mathrm{stat}.)`$, $`C_2^0=25.56\pm 4.04(\mathrm{stat}.)`$, $`C_2^x=1.53\pm 1.12(\mathrm{stat}.)`$, $`C_2^y=3.565\pm 1.17(\mathrm{stat}.)`$.

In order to account for the dependence on the renormalization scale $`x_\mu `$, one has to make the following substitutions in Eq. (22):

$`B(L)B(L)+A(L)\beta _0\mathrm{ln}x_\mu ,`$ (25)
$`\mathrm{\Sigma }(\lambda ,\eta )\mathrm{\Sigma }(\lambda ,\eta )+\lambda ^2g_1^{^{}}\mathrm{ln}x_\mu ,`$ (26)
$`G_{22}G_{22}+G_{12}\beta _0\mathrm{ln}x_\mu .`$ (27)
The other widely used matching scheme is $`R`$-matching, which will serve as a systematic check of our matching procedure. In this case the NLL and NLO predictions are combined according to the following formula :

$$R^{\mathrm{R}\mathrm{match}}=R^{\mathrm{NLL}}+\left[\eta ^2\left(BB^{\mathrm{NLL}}\right)+\eta ^3\left(CC^{\mathrm{NLL}}\right)\right],$$
(28)
where $`B^{\mathrm{NLL}}`$ and $`C^{\mathrm{NLL}}`$ are the coefficients in the $`\eta `$-expansion of $`R^{\mathrm{NLL}}`$. In the case of multi-jet rates, the NLL approximation does not exponentiate, therefore, the viable matching scheme for the four-jet rates is $`R`$-matching. The resummation formulas for the four-jet rates can be found in Refs. .

The general form of the NLO differential cross section for the four-jet angular correlations $`O_4`$ is given by :

$$\frac{1}{\sigma _0}\frac{d\sigma }{\mathrm{d}O_4}(O_4)=\eta (\mu )^2B_{O_4}(O_4)+\eta (\mu )^3\left[C_{O_4}(O_4)+B_{O_4}(O_4)\beta _0\mathrm{ln}(x_\mu ^2)\right],$$
(29)
where $`\sigma _0`$ denotes the Born cross section for the process $`e^+e^{}\overline{q}q`$. In Eq. (29) the Born and correction functions $`B_{O_4}`$ and $`C_{O_4}`$ have an analogous decomposition to Eqs. (14) and (15). To obtain the distributions normalised to unity in perturbation theory, we calculated the coefficient functions $`B_i^{(O_4)}`$ and $`C_i^{(O_4)}`$ in the linear and quadratic forms using the partonic Monte Carlo program DEBRECEN, and normalised those with the total O($`\eta ^3`$) four-jet cross section of the same $`y_{\mathrm{cut}}`$.


format: text

text:
# The Three Neutrino Scenario Talk presented at Europhysics Neutrino Oscillation Workshop (NOW2000), Otranto, Italy, September 9-16, 2000.

## 1 INTRODUCTION

Despite the current trend that many people jumped into the four neutrino scheme (see  for complehensive references), the three-flavor mixing scheme of leptons, together with the three-flavor mixing scheme of quarks, still constitutes the most promising standard model for the structure of the (known to date) most fundamental matter in nature. It is worth to note that the two of the evidences for the neutrino oscillation, one compelling  and the other strongly indicative , can be nicely fit into the three-flavor mixing scheme of neutrinos. I want to call the scheme the standard 3 $`\nu `$ mixing scheme in my talk.

I would like to address in this talk some aspects of the three-flavor mixing scheme of neutrinos which remain to be explored untill now. While I have started with a brief remark on robustness of the standard 3 $`\nu `$ mixing scheme, I do not repeat it here because I have described it elsewhere . In this manuscript I discuss mainly two key issues, (1) the sign of $`\mathrm{\Delta }m_{13}^2`$, and (2) how to measure leptonic CP violation. I will use, throughout this manuscript, the MNS matrix  in the convention by Particle Data Group.

Let me start by raising the following question;

”Suppose that the standard 3 $`\nu `$ mixing scheme is what nature exploits and the atmospheric and the solar neutrino anomalies are the hints kindly provided by her to lead us to the scheme. Then, what is left toward our understanding of its full structure?”

It is conceivable that the future atmospheric and solar neutrino observations as well as currently planned long baseline experiments will determine four parameters, $`\mathrm{\Delta }m_{23}^2`$, $`\mathrm{\Delta }m_{12}^2`$, $`\theta _{23}`$, and $`\theta _{12}`$, to certain accuracies. I then cite four things as in below which will probably be unexplored in the near future:

(i) $`\theta _{13}`$

(ii) the sign of $`\mathrm{\Delta }m_{13}^2m_3^2m_1^2`$, the normal versus inverted mass hierarchies

(iii) the CP violating leptonic Kobayashi-Maskawa phase $`\delta `$

(iv) the absolute masses of neutrinos

I make brief comments on (i) and (iv) one by one before focusing on (ii) and (iii):

Measuring $`\theta _{13}`$ is one of the goals of the currently planned long baseline experiments , and therefore I do not discuss it further. Among them the expected sensitivity in JHF, the recently approved experiment in Japan, is one of the best examined case .

Measuring absolute masses of neutrinos is certainly the ultimate challenge for neutrino experiments, but it is not clear at this moment how one can do it. Presumably, neutrinoless double $`\beta `$ decay experiments are the most promising. We, however, do not discuss it further, but just recommend the interested readers to look into a report at this conference .

## 2 SIGN OF $`\mathrm{\Delta }m_{13}^2`$

Nunokawa and I recently discussed  that the features of neutrino flavor transformation in supernova (SN) is sensitive to the sign of $`\mathrm{\Delta }m_{13}^2m_3^2m_1^2`$, making contrast between the normal ($`\mathrm{\Delta }m_{13}^2>0`$) vs. inverted ($`\mathrm{\Delta }m_{13}^2<0`$) mass hierarchies. Therefore, one can obtain insight on the sign of $`\mathrm{\Delta }m_{13}^2`$ by analyzing neutrino events from supernova. With use of the unique data at hand from SN1987A , we have obtained a strong indication that the normal mass hierarchy, $`m_3m_1m_2`$, is favored over the inverted one, $`m_1m_2m_3`$.

The point is that there are always two MSW resonance points in SN for neutrinos with cosmologically interesting mass range, $`m_\nu <100`$ eV. The higher density point, which I denote the H resonance, plays a deterministic role. If the H resonance is adiabatic the feature of $`\nu `$ flavor transformation in SN is best characterized as $`\nu _e\nu _{heavy}`$ exchange, as first pointed out in Ref. . Here, $`\nu _{heavy}`$ collectively denotes $`\nu _\mu `$ and $`\nu _\tau `$, which are physically indistinguishable in SN. See also Ref.  for a recent complehensive treatment of 3 $`\nu `$ flavor conversion in SN.

Now the question is: how does the sign of $`\mathrm{\Delta }m_{13}^2`$ make difference? The answer is: if $`\mathrm{\Delta }m_{13}^2`$ is positive (negative) the neutrino (antineutrino) undergoes the resonance. Then, if the inverted mass hierarchy is the case and assuming the adiabaticity of H resonance, the $`\overline{\nu }_e`$ which to be observed in terrestrial detectors comes from original $`\overline{\nu }_{heavy}`$ in neutrinosphere. It is widely recognized that, due to their weaker interactions with surrounding matter, $`\nu _{heavy}`$ and $`\overline{\nu }_{heavy}`$ are more energetic than $`\nu _e`$ and $`\overline{\nu }_e`$. Since the $`\overline{\nu }_e`$ induced CC reaction is the dominant reaction channel in water Cherenkov detectors, the effect of such flavor transformation would be sensitively probed by them.

We draw in Fig. 1 equal likelihood contours as a function of the heavy to light $`\nu `$ temperature ratio $`\tau T_{\overline{\nu }_x}/T_{\overline{\nu }_e}=T_{\nu _x}/T_{\overline{\nu }_e}`$ on the space spanned by $`\overline{\nu }_e`$ temperature and total neutrino luminosity by giving the neutrino events from SN1987A. The data comes from Kamiokande and IMB experiments . In addition to it we introduce an extra parameter $`\eta `$ defined by $`L_{\nu _x}=L_{\overline{\nu }_x}=\eta L_{\nu _e}=\eta L_{\overline{\nu }_e}`$ which describes the departure from equipartition of energies to three neutrino species and examine the sensitivity of our conclusion against the change in $`\eta `$.

Fig.1: Contours of constant likelihood which correspond to 95.4 % confidence regions for the inverted mass hierarchy under the assumption of adiabatic H resonance. From left to right, $`\tau T_{\overline{\nu }_x}/T_{\overline{\nu }_e}=T_{\nu _x}/T_{\overline{\nu }_e}=2,1.8,1.6,1.4,1.2`$ and 1.0 where $`x=\mu ,\tau `$. Best-fit points for $`T_{\overline{\nu }_e}`$ and $`E_b`$ are also shown by the open circles. The parameter $`\eta `$ parametrizes the departure from the equipartition of energy, $`L_{\nu _x}=L_{\overline{\nu }_x}=\eta L_{\nu _e}=\eta L_{\overline{\nu }_e}(x=\mu ,\tau )`$, and the dotted lines (with best fit indicated by open squares) and the dashed lines (with best fit indicated by stars) are for the cases $`\eta =0.7`$ and 1.3, respectively. Theoretical predictions from supernova models are indicated by the shadowed box.

At $`\tau =1`$, that is at equal $`\overline{\nu _e}`$ and $`\nu _e`$ temperatures, the 95 $`\%`$ likelihood contour marginally overlaps with the theoretical expectation  represented by the shadowed box in Fig. 1. When the temperature ratio $`\tau `$ is varied from unity to 2 the likelihood contour moves to the left, indicating less and less consistency between the standard theoretical expectation and the observed feature of the neutrino events. This is simply because the observed energy spectrum of $`\overline{\nu }_e`$ must be interpreted as that of the original one of $`\overline{\nu }_{heavy}`$ in the presence of the MSW effect in $`\overline{\nu }`$ channel. It implies that the original $`\overline{\nu }_e`$ temperature must be lower by a factor of $`\tau `$ than the observed one, leading to stronger inconsistency at larger $`\tau `$.

The solid lines in Fig. 1 are for the case of equipartition of energy into three flavors, $`\eta =1`$, whereas the dotted and the dashed lines are for $`\eta =0.7`$ and 1.3, respectively. We observe that our result is very insensitive against the change in $`\eta `$.

We conclude that if the temperature ratio $`\tau `$ is in the range 1.4-2.0 as the SN simulations indicate, the inverted hierarchy of neutrino masses is disfavored by the neutrino data of SN1987A unless the H resonance is nonadiabatic, i.e., unless $`s_{13}^2<`$ a few $`\times 10^4`$ .

## 3 HOW TO MEASURE SIGN OF $`\mathrm{\Delta }m_{13}^2`$ AND CP VIOLATION IN NEUTRINO OSCILLATION EXPERIMENTS ?

The possibility that SN can tell about the sign of $`\mathrm{\Delta }m_{13}^2`$ is, I think, interesting and in fact it is the unique available hint on the question at this moment. We, the authors of Ref. , feel that our argument and the analysis done with the SN1987A data is reasonably robust. But, of course, it would be much nicer if we can have independent confirmation by terrestrial experiments. With regard to the CP violating effect mentioned in (iii) it appears, to my understanding, that the best place for its measurement is long ($`>`$$``$ 10 km) baseline neutrino oscillation experiments.

We develop an analytic method by which we can explore various regions of experimentally variable parameters to illuminate at where CP violating effects are large and how one can avoid serious matter effect contamination. Actually we formulate below a perturbative framework to have a bird-eye view of at where the sign of $`\mathrm{\Delta }m_{13}^2`$ is clearly displayed and the CP violating phase manifests itself. Some of the earlier attempts to formulate perturbative treatment to explore the verious regions may be found in .

We rewrite the Schrödinger equation by using the basis $`\stackrel{~}{\nu }`$ defined by ($`\mathrm{\Gamma }`$ is a CP phase matrix) $`\stackrel{~}{\nu _\alpha }=\left[e^{i\lambda _5\theta _{13}}\mathrm{\Gamma }_\delta ^+e^{i\lambda _7\theta _{23}}\right]_{\alpha \beta }\nu _\beta `$, into the form

$$i\frac{\mathrm{d}}{\mathrm{d}x}\stackrel{~}{\nu _\alpha }=\left(H\right)_{\alpha \beta },\stackrel{~}{\nu _\beta }$$
(1)
where Hamiltonian H contains the following three terms:

$`H=`$ $`\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& \frac{\mathrm{\Delta }_{13}}{2E}\end{array}\right]+a\left(x\right)\left[\begin{array}{ccc}c_{13}^2& 0& c_{13}s_{13}\\ 0& 0& 0\\ c_{13}s_{13}& 0& s_{13}^2\end{array}\right]`$ (8)
$`+{\displaystyle \frac{\mathrm{\Delta }_{12}}{2E}}\left[\begin{array}{ccc}s_{12}^2& c_{12}s_{12}& 0\\ c_{12}s_{12}& c_{12}^2& 0\\ 0& 0& 0\end{array}\right]`$ . (12)
We first note the order of magnitude of a relevant quantity to observe the hierarchies of various terms in the Hamiltonian:

$$\frac{\mathrm{\Delta }m^2}{E}=10^{13}\left(\frac{\mathrm{\Delta }m^2}{10^3\text{eV}^2}\right)\left(\frac{E}{10\text{GeV}}\right)^1\text{eV}.$$
(13)
It may be compared with the matter potential $`a\left(x\right)=\sqrt{2}G_FN_e\left(x\right)`$ where $`N_e`$ denotes electron number density in the earth;

$$a\left(x\right)=1.04\times 10^{13}\left(\frac{\rho }{2.7\text{g/cm}^3}\right)\left(\frac{Y_e}{0.5}\right)\text{eV},$$
(14)
where $`Y_eN_p/\left(N_p+N_n\right)`$ is the electron fraction.

In view of these results one can identify two typical cases, the high and low energy options with $`\nu `$ beam energies $``$ 10 GeV and $``$ 100 MeV, respectively, each with a hierarchy of energy scales:

(1) High energy option

$$\frac{\mathrm{\Delta }m_{13}^2}{E}a\left(x\right)\frac{\mathrm{\Delta }m_{12}^2}{E}$$
(15)
(2) Low energy option

$$\frac{\mathrm{\Delta }m_{13}^2}{E}a\left(x\right)\frac{\mathrm{\Delta }m_{12}^2}{E}$$
(16)
Now let us discuss the high and low energy options one by one. The focus will be on the sign of $`\mathrm{\Delta }m_{13}^2`$ in the former and the CP violation in the latter.

### 3.1 High energy option: matter enhanced $`\theta _{13}`$ mechanism

In the high energy option one can formulate perturbation theory by regarding the 1st and the 2nd terms in the Hamiltonian in (12) as unperturbed part and the 3rd term as perturbation; solar $`\mathrm{\Delta }m^2`$ perturbation theory. The unperturbed system is essentially the two-flavor MSW system and it is well known that it leads to the matter enhanced $`\theta _{13}`$ mechanism in neutrino (if $`\mathrm{\Delta }m_{13}^2>0`$) or antineutrino (if $`\mathrm{\Delta }m_{13}^2<0`$) channels. Therefore, the high energy option is advantagious if $`\theta _{13}`$ is extremely small.

In leading order one can easily compute the oscillation probability in matter under the adiabatic approximation. It reads

$$P\left(\nu _\mu \nu _e\right)=s_{23}^2\mathrm{sin}^22\theta _{13}^M\mathrm{sin}^2\left(\xi _{HE}\frac{\mathrm{\Delta }m_{13}^2L}{4E}\right)$$
(17)
$$\mathrm{sin}2\theta _{13}^M=\frac{\mathrm{sin}2\theta _{13}}{\xi _{HE}},$$
(18)
where

$$\xi _{HE}=\sqrt{\left(\mathrm{cos}2\theta _{13}\pm \frac{2Ea}{\mathrm{\Delta }m_{13}^2}c_{13}^2\right)^2+\mathrm{sin}^22\theta _{13}}$$
(19)
where $`\pm `$ refers to antineutrino and neutrino channels, respectively.

Let us expand the oscillation probability by the parameter $`\frac{\mathrm{\Delta }m_{13}^2L}{4E}`$. In fact, it is a small parameter in most of the practical cases;

$`{\displaystyle \frac{\mathrm{\Delta }m_{13}^2L}{4E}}`$ $`=`$ $`0.127`$ (20)
$`\times `$ $`\left({\displaystyle \frac{\mathrm{\Delta }m_{13}^2}{10^3\text{eV}^2}}\right)\left({\displaystyle \frac{L}{1000\text{km}}}\right)\left({\displaystyle \frac{E}{10\text{GeV}}}\right)^1.`$
Then, the oscillation probability reads to next to leading order as

$`P\left(\nu _\mu \nu _e\right)`$ $`=`$ $`s_{23}^2\mathrm{sin}^22\theta _{13}\left({\displaystyle \frac{\mathrm{\Delta }m_{13}^2L}{4E}}\right)^2`$ (21)
$`\times `$ $`\left[1{\displaystyle \frac{1}{3}}\xi _{HE}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{13}^2L}{4E}}\right)^2\right]`$
The first term in (21) is identical with the vacuum oscillation probability $`P_{vac}`$ under the small $`\frac{\mathrm{\Delta }m_{13}^2L}{4E}`$ approximation. It is the simplest version of the vacuum mimicking mechanism discussed in Ref.  where a much more extensive version including the CP (or T) violating piece is uncovered.111 In passing I have a few comments on the vacuum mimicking mechanism. It might be curious that it works at the MSW resonance point because the mixing angle is certainly exhanced. But it works in such a way that there is a prolongation of oscillation length which exactly cancels the exhacement of the mixing angle . But the phenonenon of vacuum mimicking is more general which occurs not only off resonance but also in nonresonant channel as far as neutrino path length is shorter than the vacuum oscillation length. This mechanism has triggered some interests quite recently .

If the measurement is done in both neutrino and antineutrino channel, one may obtain the difference

$`\mathrm{\Delta }P`$ $``$ $`P\left(\nu _\mu \nu _e\right)P\left(\overline{\nu _\mu }\overline{\nu _e}\right)`$ (22)
$``$ $`{\displaystyle \frac{8}{3}}\left({\displaystyle \frac{Ea}{\mathrm{\Delta }m_{13}^2L}}\right)\left({\displaystyle \frac{\mathrm{\Delta }m_{13}^2L}{4E}}\right)^2P_{vac}`$
If I use $`\mathrm{\Delta }m_{13}^2=3\times 10^3`$ eV2, $`\mathrm{\Delta }P0.1P_{vac}`$ for baseline $`1000`$ km and energy $`10`$ GeV since the first parencesis is of order unity in the high energy option. Thus, the sign of $`\mathrm{\Delta }m_{13}^2`$ is determined as the sign of $`\mathrm{\Delta }P`$ .

My next and the last message about the high energy option is that the CP violating effect is small. It is obvious in leading order that no CP violating effect is induced; it is a two-flavor problem and hence there is no room for CP violation even in matter. Therefore, we have to go beyond the leading order to have CP violation. Then, the CP and T violating effect is always accompanied by the suppression factor $`\frac{\mathrm{\Delta }m_{12}^2}{\mathrm{\Delta }m_{13}^2}0.10.01`$ which comes from the energy denominator. Therefore, CP-odd effect is small in the high energy option.222It appears that this statement made some of the neutrino factory workers unhappy; probably they felt it difficult to reconsile this statement with the reported enomous sensitivities that extends toward a very small value of $`\mathrm{sin}^2\theta _{13}`$ which will be achieved by neutrino factory . However, it appears that they are actually consistent because the sensitivity is in fact achieved by the CP conserving $`\mathrm{cos}\delta `$ term not by the CP violating term in the oscillation probability at least for $`L<1000`$ km .

### 3.2 Low energy option: matter enhanced $`\theta _{12}`$ mechanism

The high energy option is certainly advantagious for the determination of the sign of $`\mathrm{\Delta }m_{13}^2`$ thanks to larger matter effects by available longer baseline due to better focusing of the $`\nu `$ beam. On the other hand, I will show that the low energy option is the natural place to look for genuine CP violation.

Because of the hierarchy in the energy scale (16), the first term in the Hamiltonian in (12) is the unperturbed term and the matter and the $`\mathrm{\Delta }m_{12}^2`$ terms are small perturbations. It is important to recognize that it is a degenerate perturbation theory because of the degeneracy in the unperturbed Hamiltonian. Then, one must first diagonalize the degenerate subspace to obtain zeroth order wave function and the first order correction to the energy eigenvalues. Then, the zeroth order wave function contains the CP violating phase effect. This is the reason why the low energy option allows large CP violation unsuppressed by the hierarchical mass ratio $`\frac{\mathrm{\Delta }m_{12}^2}{\mathrm{\Delta }m_{13}^2}`$, which is to my knowledge the unique case.

In this setting one can derive the oscillation probability $`P\left(\nu _\mu \nu _e\right)`$ as follows:

$`P`$ $`\left(\nu _\mu \nu _e\right)`$ (23)
$`=`$ $`4s_{23}^2c_{13}^2s_{13}^2\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{13}^2L}{4E}}\right)`$
$`+`$ $`c_{13}^2\mathrm{sin}2\theta _{12}^M[(c_{23}^2s_{23}^2s_{13}^2)\mathrm{sin}2\theta _{12}^M`$
$`+`$ $`2c_{23}s_{23}s_{13}\mathrm{cos}\delta \mathrm{cos}2\theta _{12}^M\left]\mathrm{sin}^2\right(\xi _{LE}{\displaystyle \frac{\mathrm{\Delta }m_{12}^2L}{4E}})`$
$``$ $`2J_M(\theta _{12}^M,\delta )\mathrm{sin}\left(\xi _{LE}{\displaystyle \frac{\mathrm{\Delta }m_{12}^2L}{2E}}\right)`$
where $`\xi _{LE}=\xi _{HE}\left(\theta _{13}\theta _{12},\mathrm{\Delta }m_{13}^2\mathrm{\Delta }m_{12}^2\right)`$, and $`J_M`$ is the matter enhanced Jarlskog factor. The probability (23) represents, apart from the $`\mathrm{cos}\delta `$ term which is small due to the factor $`s_{13}`$, the vacuum mimicking mechanism in its most extensive form including the CP violating Jarlskog term. To check how well the system mimics vacuum oscillations see Ref .

The number of appearance events in water Cherenkov detector for a beam energy E = 100 MeV is estimated by assuming 10 times stronger $`\nu _\mu `$ flux at $`L=250`$ km than the K2K design flux (despite lower energy!) and 100% conversion of $`\nu _\mu `$ to $`\nu _e`$ as

$$N6300\left(\frac{L}{100\mathrm{km}}\right)^2\left(\frac{V}{1\mathrm{Mton}}\right)\left(\frac{\mathrm{F}_{250}}{10F_{K2K}}\right).$$
(24)
where $`F_{250}`$ and $`F_{K2K}`$ are the assumed flux at a detector at $`L=250`$ km and the design neutrino flux at SK in K2K experiment, respectively. The latter is approximately, $`3\times 10^6\left({\displaystyle \frac{\mathrm{POT}}{10^{20}}}\right)`$ cm-2 where POT stands for proton on target.

To estimate the optimal distance we compute the expected number of events in neutrino and anti-neutrino channels as well as their ratios as a function of distance by taking into account of neutrino beam energy spread. For definiteness, we assume that the average energy of neutrino beam $`E`$ = 100 MeV and beam energy spread of Gaussian type with width $`\sigma _E=10`$ MeV. We present our results in Fig. 2.

Fig.2: Expected number of events for (a) neutrinos, $`N\left(\nu _\mu \nu _e\right)`$, (b) anti-neutrinos, $`N\left(\overline{\nu }_\mu \overline{\nu }_e\right)`$, and (c) their ratio $`RN\left(\nu _\mu \nu _e\right)/N\left(\overline{\nu }_\mu \overline{\nu }_e\right)`$ with a Gaussian type neutrino energy beam with $`E_\nu =`$ 100 MeV with $`\sigma `$ = 10 MeV are plotted as a function of distance from the source. Neutrino fluxes are assumed to vary as $`1/L^2`$ in all the distance range we consider. $`\mathrm{sin}^22\theta _{13}`$ is taken as 0.1, a ”maximal value” allowed by the CHOOZ limit. The remaining mixing parameters used are of the LMA MSW solution; see . The error bars are only statistical.

While this particular proposal may have several experimental problems it is sufficiantly illuminative of the fact that the low energy option is in principle more appropriate for experimental search for CP violating effect. There is a large CP violation and the matter effect is small or controllable. The remaining question is of course how to develop a feasible experimental proposal. A possibility which employs medium energy ($`12`$ GeV) conventional $`\nu `$ beam is raised by an eminent experimentalist and triggered much interests .

There were many debates between supporters of high and low energy options in the workshop. I have concluded with my personal best three flavor scenario; we measure CP violation by low energy superbeam in Japan, and you measure $`\delta `$ by neutrino factory in Europe, and then let us compare the results!

## ACKNOWLEDGMENTS

I express deep gratitude to G. L. Fogli for cordial invitation to such a focused and the well-organized workshop, in which I was able to enjoy the stimulating atmosphere of the Europian neutrino physics community. I thank Hiroshi Nunokawa for collaboration and his kind help in dealing with the elsevier latex format. This work was supported by the Grant-in-Aid for Scientific Research in Priority Areas No. 11127213, Japan Ministry of Education, Culture, Sports, Science and Technology.


text:
# Superconductivity of metallic boron in MgB2

## Abstract

Boron in MgB2 forms layers of honeycomb lattices with magnesium as a space filler. Band structure calculations indicate that Mg is substantially ionized, and the bands at the Fermi level derive mainly from B orbitals. Strong bonding with an ionic component and considerable metallic density of states yield a sizeable electron-phonon coupling. Using the rigid atomic sphere approximation and an analogy to Al, we estimate the coupling constant $`\lambda `$ to be of order 1. Together with high phonon frequencies, which we estimate via zone-center frozen phonon calculations to be between 300 and 700 cm$`^1,`$ this produces a high critical temperature, consistent with recent experiments. Thus MgB2 can be viewed as an analog of the long sought, but still hypothetical, superconducting metallic hydrogen.

Before the discovery of high-temperature superconductors much effort was devoted to boosting the critical temperature of conventional, BCS-Eliashberg superconductors. An exotic and appealing idea going back to the early 60’s, was that of metallic hydrogen. The arguments were very simple: due to the light mass, the phonon frequencies in metallic hydrogen would be very high, of the order of several thousand K, and the prefactor in the BCS formula would be very large, so that even a moderate coupling constant would provide a sizeable $`T_c.`$ This idea can be quantified as follows: for monoatomic solids, the electron-phonon coupling (EPC) constant, $`\lambda ,`$ which enters the BCS equation, can be written in the so-called McMillan-Hopfield form , $`\lambda =N(0)I^2/M\omega ^2`$, where $`N(0)`$ is the density of states (DOS) at the Fermi level per spin per atom, $`I^2`$ is the properly averaged electron-ion matrix element squared, $`M`$ is the atomic mass and $`\omega ^2`$ is (again, properly averaged) the phonon frequency. The product $`M\omega ^2`$ does not depend on the mass, but on the force constants only, while $`\eta =N(0)I^2,`$ also known as the Hopfield factor, is a purely electronic property. Correspondingly, light elements, everything else being the same, are beneficial for superconductivity.

Lacking metallic hydrogen, attention was focused upon compounds with light elements: carbides, nitrides (arguably, the superconductivity in fullerenes was a discovery along this road). Indeed, many of them were ‘high-$`T_c`$ superconductors’ on the contemporary scale: 10-15 K. It was pointed out that $`I^2`$ is rather large in these materials due to the relatively high ionicity (although not as high as in MgB2), but $`N(0)`$ is rather small. This led to the suggestion of cubic MoN, which would have a larger $`N(0)`$ than existing nitrides and carbides, as a hypothetical superconductor with $`T_c>30`$ K.

The recently discovered medium-$`T_c`$ superconductor MgB2  with $`T_c`$39 K is clearly a continuation of the same idea. The main component, B, is even lighter than C and N. Furthermore, electronic structure calculations show that the compound is not only quite ionic with a reasonable DOS, but also has strong covalent B-B bonding (the bonding-antibonding splitting due to in-plane B-B hopping is about 6 eV) and thus exhibits strong electron-phonon interactions. Interestingly, unlike carbides and nitrides, and similar to metallic hydrogen, electrons at the Fermi level (and below) are predominantly B-like. Mg $`s`$-states are pushed up by the B $`p_z`$ orbitals and fully donate their electrons to the boron-derived conduction bands. In the following we will describe the physics of such ‘metallic’ boron in detail, present an estimate of the EPC constant, and propose some routes for optimizing $`T_c`$ in this kind of compound.

MgB2 occurs in the so-called AlB2 structure. Borons form a primitive honeycomb lattice, consisting of graphite-type sheets stacked with no displacement. The borons form hexagonal prisms with the base diameter of 3.5Å nearly equal to the height. This creates large, nearly spherical pores for Mg. As in graphite, the intraplanar B-B bonds are much shorter than the distance between the planes, and hence the B-B bonding is strongly anisotropic. However, the interplane bonds are only twice as long as the intraplane ones, as compared to the ratio of 2.4 in graphite, allowing for a significant interplane hopping.

We have calculated the electronic structure of MgB2 using a general potential LAPW code. For the Rigid Atomic Spheres calculations we used the Stuttgart LMTO-TB code. For the exchange-correlation potential, the Generalized Gradient Approximation of Ref. was employed. Despite the rather simple crystal structure, very few electronic structure calculations for MgB2 have been reported (a model TB calculation of Burdett and Miller and a recent full-potential LMTO study), and these have concentrated mainly on chemical bonding, paying hardly any attention to transport and electronic properties. The results of our LAPW calculations are shown in Figs. 1 and 2. We note first that there is almost no valence charge inside the Mg MT sphere (less than 0.2 $`e`$). About half of the total valence charge resides inside the B spheres, and about the same amount in the interstitials. This is partially due to the fact that the chosen LAPW setup employs rather small MT spheres for Mg. For the LMTO calculations we used an atomic sphere of nearly the size of the free Mg atom (up to 3.13$`a_B),`$ and obtained, as expected, a larger charge of 2.8 electrons. However less than 25% of the charge has $`s`$-character. The remaining charge of $`p`$-, $`d`$-, and $`f`$\- character arises not from Mg electrons but rather from the tails of the B $`p`$-orbitals and contributions from the interstitials. In fact, one can say that Mg is fully ionized in this compound, however the electrons donated to the system are not localized on the anion, but rather are distributed over the whole crystal.

The resulting band structure can be easily understood in terms of the boron sublattice. We plotted in Fig. 1 the character of the bands (we show only the B $`p`$ character, since contributions from other orbitals near the Fermi level are very small). We observe two B band systems: two bands are derived from B $`p_z`$ states and four from B $`p_{x,y}`$. All these bands are highly dispersive (light), the former being quite isotropic and the latter more two-dimensional. Both $`p_z`$ bands cross the Fermi level (in different parts of the Brillouin zone), but only two bonding px,y bands do so, and only near the $`\mathrm{\Gamma }`$ point. They form two small cylindrical Fermi surfaces around the $`\mathrm{\Gamma }`$-$`A`$ line (Fig. 3). However, due to their 2D character, they contribute more than 30% to the total $`N(0)`$.

In contrast, the $`p_z`$ bands have 3D character, since the smaller intraplane distance compensates for a smaller (pp$`\pi `$ $`vs.`$ pp$`\sigma )`$ hopping. In the nearest neighbor tight binding (TB) model their dispersion is $`\epsilon _𝐤=\epsilon _0+2t_{pp\sigma }\mathrm{cos}ck_z\pm t_{pp\pi }\sqrt{3+2\mathrm{cos}𝐚_1𝐤+2\mathrm{cos}𝐚_\mathrm{𝟐}𝐤+2\mathrm{cos}𝐚_3𝐤}`$ where a1,2,3 are the smallest in-plane lattice vectors. The on-site parameter $`\epsilon _0`$ can be found from the eigenvalue at the $`K`$ point and is $`1.5`$ eV above the Fermi energy. We estimated $`t_{pp\sigma }`$ and $`t_{pp\pi }`$ from the LMTO calculations as $`2.5`$ eV and $``$1.5 eV, respectively. This model gives a very good description of the band structure near and below the Fermi, although the antibonding band acquires some additional dispersion by hybridizing with the Mg $`p`$ band. The role of Mg in forming this band structure can be elucidated by removing the Mg atoms from the lattice entirely and repeating the calculations in this hypothetical structure. The in-plane dispersion of the both sets of bands at and below the Fermi level changes very little (for the $`pp\pi `$ bands there is hardly any change at all, while the $`pp\sigma `$ in-plane dispersion changes by $``$ 10%). The $`k_z`$ dispersion of the $`p_z`$ bands is increased in MgB2 as compared with the hypothetical empty B2 lattice by about 30%, and these bands shift down with respect to the $`p_{x,y}`$ bands by approximately 1 eV. This shift, as well as the additional dispersion, comes mainly from the hybridization with the empty Mg $`s`$ band, which is correspondingly pushed further up, increasing the effective ionicity. Substantial $`k_z`$ dispersion of the $`p_z`$ bands produces the Fermi surface which is is approximately mirror-reflected with respect to a plane between the $`k_z=0`$ and $`k_z=\pi /c`$ planes, one pocket (electron-like) coming from the antibonding and the other (hole-like) from the bonding $`p_z`$ band. The two surfaces touch at one point on the K-H line and form a honeycomb tubular network, replicating in reciprocal space the boron lattice in real space.

The resulting bands are fairly 3D: the average Fermi velocities are $`v_{x,y}=4.9010^7`$ cm/s and $`v_z=4.7610^7`$ cm/s. The plasma frequencies are $`\omega _{px,y}=7.1`$ eV and $`\omega _{pz}=6.9`$ eV. Correspondingly, we predict fairly isotropic electrical resistivity, with the linear slope $`d\rho /dT0.08\lambda _{tr}`$ $`\mu \mathrm{\Omega }`$ cm$`/K`$. The total DOS at the Fermi level is $`N(0)=0.36`$ states/spin$``$f.u.

The above described band structure is typical for an $`sp`$metal. What is not typical is that this particular $`sp`$metal is held together by covalent bonding with a substantial ionic component, which inevitably leads to a strong electron-phonon interaction. In the following we present a semiquantitative estimate of the corresponding coupling constant, and argue that the fortunate combination of strong bonding, reasonable $`N(0)`$, and high phonon frequency is responsible for the high transition temperature in this compound.

In its most rigorous formulation the McMillan-Hopfield formula reads:

$$\lambda =n^1\underset{ij}{}NI_{i\alpha }I_{j\beta }\left(\mathrm{\Phi }^1\right)_{i\alpha ,j\beta },$$
(1)
where the indices $`i,j`$ run over all $`n`$ ions in the crystal; $`\alpha ,\beta `$ are Cartesian indices; $`NI_{i\alpha }I_{j\beta }`$ is the electron-ion matrix element averaged over the Fermi surface; and $`\mathrm{\Phi }_{i\alpha ,j\beta }`$ is the standard force matrix. While this expression is exact, a number of simplifications are needed to make it more practical. First, because of the large size of the Fermi surface, one can neglect the nondiagonal terms and reduce Eq. 1 to one unit cell:

$$\lambda \stackrel{~}{\underset{i}{}}NI^2_i\mathrm{\Phi }_{ii}^1=\stackrel{~}{\underset{i}{}}\eta _i\left(\mathrm{\Phi }_{ii}\right)^1.$$
(2)
The quantity $`\mathrm{\Phi }_{ii}=^2E_{tot}/R_i^2`$ is a local quantity, which can be calculated from the total energy differences of frozen zone-center phonons. The Hopfield factor can be calculated in the Rigid Muffin Tin (or Rigid Atomic Sphere) Approximation (RMTA), assuming that the change of the crystal potential due to an ion’s displacement can be described by shifting the electronic charge distribution rigidly inside the corresponding atomic sphere. Although expressions for $`\eta `$ have been derived for arbitrary site symmetry, these are complicated and we will use here simplified formulas formally correct for cubic site symmetry. With this simplification, $`\eta `$ can be readily calculated from the LMTO potential parameters and partial DOS’s. We show the results in Table I.

In order to get an estimate of the phonon spectrum of the system, we calculated all zone-center modes, using the full-potential LAPW method. There are four distinct modes : one silent mode, $`B_{1g}`$ (two borons displaced along $`z`$ in opposite directions), one doubly degenerate Raman mode, $`E_{2g}`$ (in-plane displacements of borons), and two infrared-active modes, which do not involve changes of in-plane bonds: $`A_{2u}`$ (B and Mg planes moving against each other), and a doubly degenerate $`E_{1u}`$ mode (B and Mg planes sliding along $`x,y).`$ Their frequencies are, respectively, 690 cm$`^1,`$ 470 cm-1, 390 cm-1, and 320 cm$`^1,`$ and their force constants are 390, 180, 70, and 44 mRy/$`a_B^2`$/atom. All calculations were carried out at the experimental lattice constant and $`c/a`$ ratio .

Given the physics of the electronic structure described above, it seems likely that the lowest mode couples little with the electrons and that its softness is derived from its in-plane acoustic character. Thus, we excluded it from the calculations of $`\lambda `$ below. Correspondingly, the average inverse force matrix for boron is $`\left(\mathrm{\Phi }_B\right)^1=M_B^1\omega ^2=7.1a_B^2`$/Ry, corresponding to $`\omega ^2^{1/2}400`$ cm-1. Together with the above value for $`\eta `$ this gives $`\lambda 0.7`$ and the logarithmically averaged frequency of the same three modes is $`(690470^2390)^{0.25}=\omega _{\mathrm{log}}500`$ cm$`{}_{}{}^{1}700`$ K. We can now estimate the critical temperature according to the McMillan formula, $`T_c=\frac{\omega _{\mathrm{log}}}{1.2}\mathrm{exp}[1.02(1+\lambda )/(\lambda \mu ^{}\mu ^{}\lambda )].`$ Using for the Coulomb pseudopotential $`\mu ^{}`$ the commonly accepted value 0.1, we obtain $`T_c`$22 K.

It is instructive to compare the calculation above with a typical $`sp`$ superconductor, Al (Table II). RMTA notoriously underestimates the electron-ion scattering in $`sp`$ metals due to the weaker screening than in $`d`$-metals. Also in Al, $`\lambda `$ is underestimated by a factor of 2. The screening properties of MgB2 are similar to those of Al (they have similar Thomas-Fermi screening lengths), and the electronic properties are similar as well. It is therefore tempting to scale the calculated $`\lambda `$ for MgB2 by the same factor. The scaled results are also shown in Tab. II and the corresponding $`T_c`$ is approximately 70 K. Although this table looks impressively quantitative, these numbers should be regarded only as rough estimates. The RMTA in MgB2 is clearly a worse approximation than in Al. The boron site symmetry is far from cubic and, moreover, large differences in the atomic sphere radii lead to artificial potential jumps at the sphere boundary, a problem for which there is no remedy in RMTA. In other words, with regard to EPC, our calculations should be considered as a qualitative indication of a strong electron-phonon interaction. We can nevertheless be confident of our main qualitative conclusions.

Our main conclusion is that MgB2, being essentially a metallic boron held together by covalent B-B and ionic B-Mg bonding, is electronically a typical $`sp`$ metal with a typical DOS. Strong bonding induces strong electron-ion scattering and hence strong electron-phonon coupling. An additional benefit is the high frequency of the boron vibrations (while the force constants remain reasonably soft). Superconductivity is mainly due to boron. The light mass and correspondingly large zero-point vibrations ($`>0.1a_B`$) suggest a possibility of anharmonic and/or nonlinear EPC, and possible deviation of the isotope effect from 100%. Isovalent doping may be beneficial if it increases the density of states $`N(0)`$. Lattice expansion due to Ca doping should lead to an overall increase of the density of states, and may provide the additional benefit of reduced $`p_z`$-$`s`$-$`p_z`$ hopping. Another interesting dopant is Na, which should not only expand the lattice, but also decrease the Fermi level, exposing more of the $`p_{x,y}`$ bands, which may provide an additional contribution to $`\lambda `$.

Finally, let us outline the directions for further theoretical investigation. First, EPC calculations beyond the RMTA (e.g., in the linear response formalism) are highly desirable and computationally feasible. Second, calculations with full structure optimization for (hypothetical) CaB2 and BeB$`_2,`$ and virtual crystal calculations for Na doping should elucidate the effect of isovalent and hole doping, giving some hints toward further optimizing $`T_c`$. Work along these lines is currently in progress.

We are thankful to J.E. Pask, C.S. Hellberg and D.J. Singh for critical reading of the manuscript. This research was supported in part by ONR and by DOE under Contract No. W-7405-82.