text:
# Measurement of the Electroweak Penguin Process 𝑩→𝑿_𝒔ℓ⁺ℓ⁻

on leave from \]Nova Gorica Polytechnic, Nova Gorica on leave from \]University of Toronto, Toronto ON The Belle Collaboration

(August 18, 2002)

## Abstract

We report the first measurement of the branching fraction for the inclusive decay $`BX_s\mathrm{}^+\mathrm{}^{}`$, where $`\mathrm{}`$ is either an electron or a muon and $`X_s`$ is a hadronic recoil system that contains an $`s`$-quark. We analyzed a data sample of $`65.4\times 10^6`$ $`B`$ meson pairs collected with the Belle detector at the KEKB $`e^+e^{}`$ asymmetric-energy collider. We find $`(BX_s\mathrm{}^+\mathrm{}^{})=(6.1\pm 1.4(\mathrm{stat})_{1.1}^{+1.4}(\mathrm{syst}))\times 10^6`$ for dilepton masses greater than $`0.2\mathrm{GeV}/c^2`$.

preprint:   Belle Preprint 2002-28 KEK Preprint 2002-83 TIT-HPE-02-01

The first observation of the flavor-changing-neutral-current (FCNC) weak decay process $`BK\mathrm{}^+\mathrm{}^{}`$, recently reported by Belle bib:belle-kll , opens a new window for searches for physics beyond the Standard Model (SM) bib:beyond-sm . There are no SM first-order weak decay processes that can produce such decays; the dominant SM processes are second-order $`b`$-quark to $`s`$-quark processes with a virtual $`t`$-quark in a loop (electroweak penguin) or a box diagram. Non-SM processes could produce sizable modifications to the decay amplitude due to contributions from virtual non-SM particles (such as charged Higgs or SUSY particles) in the loop.

The decay amplitude is often described by an effective Hamiltonian in which non-SM contributions can modify the Wilson coefficients $`C_7`$, $`C_9`$ and $`C_{10}`$. The magnitude of allowable modifications of $`C_7`$ is constrained by the measured rate for inclusive $`BX_s\gamma `$ decays bib:belle-btosgamma ; bib:cleo-btosgamma ; bib:aleph-btosgamma . The Belle measurement of the exclusive process $`BK\mathrm{}^+\mathrm{}^{}`$ has been used to determine the best limits on non-SM contributions to $`C_9`$ and $`C_{10}`$ bib:ali-2001 . The usefulness of exclusive measurements is limited by the large theoretical uncertainties associated with the $`s`$-quark to $`K`$ meson hadronization process. Since these uncertainties are not as severe for inclusive processes, measurements of the branching fraction and the dilepton and hadronic mass spectra for $`BX_s\mathrm{}^+\mathrm{}^{}`$ will provide more stringent and less model-dependent probes for new physics. This process has not been measured; CLEO has reported a 90% confidence upper limit of $`(BX_s\mathrm{}^+\mathrm{}^{})<4.2\times 10^5`$ bib:cleo-sll .

In this Letter, we present the results of a measurement of the branching fraction for the inclusive decay $`BX_s\mathrm{}^+\mathrm{}^{}`$, where $`B`$ is either $`B^0`$ or $`B^+`$, $`\mathrm{}`$ is either an electron or a muon, and $`X_s`$ is a hadronic recoil system that contains an $`s`$-quark. Here and throughout the paper, the inclusion of charge conjugated modes is implied. We use a data sample collected at the $`\mathrm{{\rm Y}}(4S)`$ resonance with the Belle detector at the KEKB $`e^+e^{}`$ asymmetric-energy collider (3.5 GeV on 8 GeV) bib:kekb-unpublished . This sample contains $`(65.4\pm 0.5)\times 10^6`$ $`B`$ meson pairs, corresponding to an integrated luminosity of $`60\mathrm{fb}^1`$.

The Belle detector is a large-solid-angle magnetic spectrometer that consists of a three-layer silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Čerenkov counters (ACC), time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter comprised of CsI(Tl) crystals (ECL) located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-return located outside of the coil is instrumented with resistive plate counters to identify muons (KLM). The detector is described in detail elsewhere bib:belle-nim .

We reconstruct charged particle trajectories with the CDC and SVD. Electron identification is based on the position and shower shape of the cluster in the ECL, ratio of the cluster energy to the track momentum ($`E/p`$), specific energy-loss measurement ($`dE/dx`$) with the CDC and the response from the ACC. Muon identification is based on the hit positions and the depth of penetration into the ECL and KLM. We require the electron’s (muon’s) laboratory momentum to be greater than $`0.5\mathrm{GeV}/c`$ ($`1\mathrm{GeV}/c`$). We find an electron (muon) selection efficiency of $`92.5\%`$ ($`91.3\%`$) with a $`(0.2\pm 0.06)\%`$ ($`(1.4\pm 0.04)\%`$) pion to electron (muon) mis-identification probability. Charged kaon candidates are selected by using information from the ACC, TOF and $`dE/dx`$ in the CDC. The kaon selection efficiency is $`90\%`$ with a pion to kaon mis-identification probability of $`6\%`$. The remaining charged tracks are assumed to be pions. We select $`K_S^0\pi ^+\pi ^{}`$ candidates with invariant mass within $`15\mathrm{MeV}/c^2`$ of the $`K^0`$ mass. We impose additional $`K_S^0`$ selection criteria based on the distance and the direction of the $`K_S^0`$ vertex and the impact parameters of daughter tracks. We require the charged tracks other than those used in the $`K_S^0`$ reconstruction to have impact parameters with respect to the nominal interaction point of less than $`0.5\mathrm{cm}`$ in the radial direction and $`3\mathrm{cm}`$ along the beam direction.

We reconstruct photons from ECL energy clusters that have no associated charged tracks. The shower shape is required to be consistent with an electromagnetic cluster and the energy to be greater than $`50\mathrm{MeV}`$. We reconstruct $`\pi ^0\gamma \gamma `$ candidates from photon pairs with invariant mass within $`10\mathrm{MeV}/c^2`$ of the nominal $`\pi ^0`$ mass.

The $`X_s`$ system is reconstructed in 18 different combinations of either a $`K^+`$ or $`K_S^0`$ combined with 0 to 4 pions, of which up to one $`\pi ^0`$ is allowed. We combine the $`X_s`$ with two oppositely charged leptons to form a $`B`$ candidate. We identify the $`BX_s\mathrm{}^+\mathrm{}^{}`$ signal with the beam-energy constrained mass, $`M_{\mathrm{bc}}=\sqrt{(E_{\mathrm{beam}}^{\mathrm{CM}}/c^2)^2(p_B^{\mathrm{CM}}/c)^2}`$, where $`E_{\mathrm{beam}}^{\mathrm{CM}}`$ and $`p_B^{\mathrm{CM}}`$ are the beam energy and the $`B`$ candidate momentum calculated in the center-of-mass (CM) system, respectively. We find the average $`M_{\mathrm{bc}}`$ resolution is $`\sigma _{\mathrm{bc}}=2.8\mathrm{MeV}/c^2`$. The energy difference $`\mathrm{\Delta }E=E_B^{\mathrm{CM}}E_{\mathrm{beam}}^{\mathrm{CM}}`$, where $`E_B^{\mathrm{CM}}`$ is the CM energy of the $`B`$ candidate, is combined with other variables to provide background suppression.

To optimize the selection criteria and to determine their signal efficiencies, we use the following signal model. The SM calculation of Ref. bib:ali-2001  is used for the dilepton mass ($`M_\mathrm{}^+\mathrm{}^{}`$) spectrum. We require $`M_\mathrm{}^+\mathrm{}^{}`$ to be greater than $`0.2\mathrm{GeV}/c^2`$ to remove the virtual photon contribution, $`bs\gamma ^{}se^+e^{}`$ bib:ali-priv . We model the recoil mass ($`M_{X_s}`$) spectrum using a non-resonant Fermi motion model bib:ali-1997  and the JETSET hadronization program bib:jetset . For $`M_{X_s}<1.1\mathrm{GeV}/c^2`$, however, the spectrum is modeled by the sum of the exclusive $`BK^{()}\mathrm{}^+\mathrm{}^{}`$ components as predicted by the SM bib:exclusive-prediction . Using this model we find that our reconstructed final states account for $`81\%`$ of the total $`X_s`$ states. Figures 1(a) and (b) show the $`M_\mathrm{}^+\mathrm{}^{}`$ and $`M_{X_s}`$ spectra from this model.

There are two background sources that can peak in $`M_{\mathrm{bc}}`$ and $`\mathrm{\Delta }E`$. The first is hadronic $`B`$ decays into one kaon plus multiple pions ($`BX_s\pi ^+\pi ^{}`$) from abundant decays such as $`BD^{()}n\pi `$ $`(n1)`$. We estimate this background contribution by reconstructing $`BX_s\pi ^+\pi ^{}`$ events without the lepton identification criteria and multiplying by the measured momentum dependent pion mis-identification probability. We find contributions of $`2.6\pm 0.2`$ events to the $`BX_s\mu ^+\mu ^{}`$ signal and $`0.1\pm 0.05`$ events to $`BX_se^+e^{}`$. We refer to this as the fake background, and subtract it from the signal yield. The second is from $`BJ/\psi X_s`$ and $`B\psi ^{}X_s`$ where $`J/\psi `$ and $`\psi ^{}`$ decay into dileptons. These decay modes have the same final states and, in principle, interfere with the signal. For this study, these charmonium decays are explicitly vetoed. The veto windows are $`0.6`$ to $`+0.2\mathrm{GeV}/c^2`$ ($`0.35`$ to $`+0.2\mathrm{GeV}/c^2`$) around the $`J/\psi `$ mass for the $`e^+e^{}`$ ($`\mu ^+\mu ^{}`$) channel, and $`0.3`$ to $`+0.15\mathrm{GeV}/c^2`$ around the $`\psi ^{}`$ mass for both channels. We find the background from this source is negligibly small using a $`BJ/\psi X_s`$ Monte Carlo (MC) sample in which the $`M_\mathrm{}^+\mathrm{}^{}`$ spectrum reproduces that of data.

The largest background sources are random combinations of dileptons with a kaon and pions that originate from continuum $`q\overline{q}`$ ($`q=u,d,s,c`$) production or from semileptonic $`B`$ decays. We reject 83% of the $`q\overline{q}`$ background with a signal efficiency of 90% by using a Fisher discriminant bib:fisher  ($`_{\mathrm{cont}}`$) based on a modified set of Fox-Wolfram moments bib:fox-wolfram  that differentiate the event topology. In the semileptonic $`B`$ decay background, both $`B`$ mesons decay into leptons or two leptons are produced from the $`bcs,d`$ decay chain. We combine the missing mass ($`M_{\mathrm{miss}}`$) and the total visible energy ($`E_{\mathrm{vis}}`$) bib:mmiss  into another Fisher discriminant ($`_{\mathrm{sl}}`$) to reject 85% of events with two neutrinos with a signal efficiency of 91%.

We further reduce the backgrounds using $`\mathrm{\Delta }E`$ and the cosine of the $`B`$ flight direction ($`\mathrm{cos}\theta _B`$) with respect to the $`e^{}`$ beam direction in the CM frame. First, we select events with $`|\mathrm{\Delta }E|<40\mathrm{MeV}`$ $`(3\sigma _{\mathrm{\Delta }E})`$. We then calculate likelihoods $`_{S,B}=p_{S,B}^{\mathrm{\Delta }E}\times p_{S,B}^{\mathrm{cosB}}`$, where $`p_{S,B}^{\mathrm{\Delta }E}`$ and $`p_{S,B}^{\mathrm{cosB}}`$ are the uncorrelated probability density functions (PDF) for $`\mathrm{\Delta }E`$ and $`\mathrm{cos}\theta _B`$ for the signal ($`S`$) and the background ($`B`$), respectively, and form a likelihood ratio $`=_S/(_S+_B)`$. We model $`p_S^{\mathrm{\Delta }E}`$ with a Gaussian primarily determined from a signal MC sample and calibrated using $`BJ/\psi X_s`$ data; $`p_B^{\mathrm{\Delta }E}`$ is modeled with a linear function with a slope determined from a MC sample that contains $`bc`$ decays and $`q\overline{q}`$ events. The signal follows a $`1\mathrm{cos}^2\theta _B`$ distribution, while the background distributes uniformly in $`\mathrm{cos}\theta _B`$. By requiring $`>0.8`$, we retain 75% of the signal while removing 80% of backgrounds.

For events with multiple candidates that pass the selection criteria, we choose the combination with the largest value of $``$. We find that the correct candidate is reconstructed in 73% of the events. We reject candidates with $`X_s`$ invariant mass greater than $`2.1\mathrm{GeV}/c^2`$. This condition removes a large fraction of combinatorial background while retaining ($`93\pm 5`$)% of the signal.

Signal MC samples are used to determine reconstruction efficiencies of $`(3.7\pm 0.4\pm 0.5)\%`$, where the first error is systematic and the second error is due to the model uncertainty. Here, an equal production rate is assumed for $`B^0\overline{B}^0`$ and $`B^+B^{}`$. Systematic errors include uncertainties from the tracking efficiency ($`2.0\%`$ per track), $`K_S^0`$ reconstruction efficiency ($`8.7\%`$ per $`K_S^0`$), $`\pi ^0`$ reconstruction efficiency ($`6.8\%`$ per $`\pi ^0`$), electron identification ($`1.8\%`$ per electron), muon identification ($`2.2\%`$ per muon), $`K^+`$ identification ($`2.5\%`$ per $`K^+`$), $`\pi ^+`$ identification ($`0.8\%`$ per $`\pi ^+`$), and the background suppression criteria ($`3.0\%`$ for $`_{\mathrm{cont}}`$, $`_{\mathrm{sl}}`$ and $``$, combined, estimated using $`BJ/\psi X_s`$ data).

The model uncertainty includes the following sources. The largest error is due to uncertainties in the SM branching fraction of the exclusive modes (11%). The uncertainty in the $`M_{X_s}`$ spectrum (4%) is due to the Fermi momentum ($`p_F`$) and spectator quark mass ($`m_q`$) parameters. We take $`p_F=0.4\mathrm{GeV}/c`$ and $`m_q=0\mathrm{GeV}/c^2`$, and vary the parameters over a range allowed by the measured heavy quark effective theory (HQET) parameters $`\lambda _1`$ and $`\overline{\mathrm{\Lambda }}`$ bib:cleo-btosgamma ; bib:cleo-hqet . The uncertainty in the fraction of the unmeasured modes (2.1%) and the uncertainty due to the fractions of $`\pi ^0`$ and $`K_S^0`$ contained in $`X_s`$ (4.6%) are estimated by comparing the inclusive hadron ($`\pi ^0`$, $`K_S^0`$, $`\eta `$, $`\varphi `$, etc.) production rates between continuum data and JETSET.

We determine the signal yield from an unbinned maximum likelihood fit to the $`M_{\mathrm{bc}}`$ distribution as shown in Fig. 2. We model the signal with a Gaussian function and the background with a threshold function bib:argus-function . The width of the signal Gaussian is obtained from $`BJ/\psi X_s`$ data. The background shape is obtained from the background MC sample. We verify that the threshold function obtained from MC reproduces the shape of data for $`BX_se\mu `$ combinations (Fig. 2(d)). The fit results are given in Table 1; we find $`60.1\pm 13.9(\mathrm{stat})_{5.4}^{+8.6}(\mathrm{syst})`$ events for the combined $`BX_s\mathrm{}^+\mathrm{}^{}`$ modes (Fig. 2(c)) with a statistical significance of $`5.4\sigma `$. Here, the systematic error is obtained from the maximum deviation when the background and signal shapes are varied by one standard deviation of the statistical error in their shape parameters. The significance is defined as $`\sqrt{2\mathrm{ln}(_0/_{\mathrm{max}})}`$, where $`_{\mathrm{max}}`$ is the maximum likelihood and $`_0`$ is the maximum likelihood when the signal yield is constrained to be zero. The fake background and the uncertainty in the background shape are included in the significance calculation. We calculate the branching fraction for $`M_\mathrm{}^+\mathrm{}^{}>0.2\mathrm{GeV}/c^2`$ to be

$$(BX_s\mathrm{}^+\mathrm{}^{})=(6.1\pm 1.4(\mathrm{stat})_{1.1}^{+1.4}(\mathrm{syst}))\times 10^6,$$

where the systematic errors in the yield, efficiency, the number of $`B`$ meson pairs and the model errors are added in quadrature to give the total systematic error. This can be compared to the SM expectation $`(BX_s\mathrm{}^+\mathrm{}^{})=(4.2\pm 0.7)\times 10^6`$ bib:ali-priv . Table 1 summarizes the branching fractions for $`BX_se^+e^{}`$ and $`BX_s\mu ^+\mu ^{}`$ separately, together with the number of candidates, signal yields, fake background estimations, efficiencies and the statistical significances.

The dilepton mass spectrum is measured by dividing the data into $`M_\mathrm{}^+\mathrm{}^{}`$ bins. For each bin, the signal yield is extracted from a fit to the $`M_{\mathrm{bc}}`$ distribution and the fake background contribution is subtracted. The result is shown in Fig. 1(c). Similarly, the recoil mass spectrum is obtained by dividing the data into $`M_{X_s}`$ bins and extracting the signal yield for each bin. The result is shown in Fig. 1(d). With the current statistics, the $`M_\mathrm{}^+\mathrm{}^{}`$ and $`M_{X_s}`$ spectra are in agreement with SM expectations. Branching fractions for each bin are given in Table 2.

In summary, we present the first measurement of the inclusive branching fraction for the electroweak penguin decay $`BX_s\mathrm{}^+\mathrm{}^{}`$. The results are in agreement with the SM expectations and can be used to constrain extensions of the SM.

We wish to thank the KEKB accelerator group for the excellent operation of the KEKB accelerator. We thank A. Ali, E. Lunghi and G. Hiller for providing us many helpful suggestions and calculations. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology of Japan and the Japan Society for the Promotion of Science; the Australian Research Council and the Australian Department of Industry, Science and Resources; the National Science Foundation of China under contract No. 10175071; the Department of Science and Technology of India; the BK21 program of the Ministry of Education of Korea and the CHEP SRC program of the Korea Science and Engineering Foundation; the Polish State Committee for Scientific Research under contract No. 2P03B 17017; the Ministry of Science and Technology of the Russian Federation; the Ministry of Education, Science and Sport of the Republic of Slovenia; the National Science Council and the Ministry of Education of Taiwan; and the U.S. Department of Energy.


format: text

text:
# Perturbative calculation of 𝑶(𝒂) improvement coefficients

## I Introduction

During the past few years the Symanzik effective field theory has been an important focus of research in lattice gauge theory. Symanzik’s idea is to describe cutoff effects in lattice field theory by a continuum effective field theory Symanzik:1979ph . One writes Symanzik:1979ph ; Symanzik:1983dc

$$_{\text{lat}}_{\text{Sym}},$$
(1)
where the symbol $``$ means that the lattice and Symanzik field theories have the same on-shell matrix elements. For lattice QCD with Wilson fermions Wilson:1975hf  the Symanzik local effective Lagrangian (LE$``$) is given by Sheikholeslami:1985ij ; Kronfeld:2002pi

$$_{\text{Sym}}=\frac{1}{2g^2}\text{tr}[F^{\mu \nu }F_{\mu \nu }]\overline{q}(/D+m)q+aK_{\sigma F}\overline{q}i\sigma _{\mu \nu }F^{\mu \nu }q+O(a^2),$$
(2)
where $`g^2`$ is a renormalized coupling, $`m`$ is a renormalized quark mass, and $`aK_{\sigma F}`$ is a short-distance coefficient. The effective field theory is useful when the scale of QCD in lattice units is small, $`\mathrm{\Lambda }a1`$, and, as used in this paper, when $`ma1`$ also. With the description in hand, the lattice field theory can be adjusted so that it approaches its continuum limit more quickly. The effective theory shows that if $`K_{\sigma F}`$ is reduced for any given on-shell matrix element, then the $`O(a)`$ term in Eq. (2) makes commensurately smaller contributions to all other on-shell matrix elements. This application of the Symanzik effective field theory is called the Symanzik improvement program Symanzik:1983dc .

A similar correspondence is set up for the vector and axial vector currents (see below), introducing further short-distance coefficients. In the last several years methods have been devised to study all of them non-perturbatively Jansen:1996ck ; Luscher:1996sc ; Martinelli:1997zc ; deDivitiis:1997ka ; Bhattacharya:1999uq . The $`O(a)`$ discretization effects violate chiral symmetry, so the key idea is to ensure that violations of chiral symmetry are at least $`O(a^2)`$. On the other hand, because of asymptotic freedom and the success of perturbative QCD, even at GeV energies Davier:1997kw , one expects perturbation theory to yield accurate estimates of the short-distance coefficients. In this paper, we compare a perturbative calculation of the currents’ short-distant coefficients to the non-perturbative results.

There are two issues that should be kept in mind when making such a comparison. First, the non-perturbative technique suffers from power corrections. Asymptotically, as $`\mathrm{\Lambda }a0`$ these are formally smaller than any error made from truncating the perturbation series. In practice, however, these effects can be significant.

Second, no two-loop results are available for the improvement coefficients considered here. Tests of perturbation theory are, therefore, not unambiguous, because different choices for the expansion parameter $`g^2`$ yield quantitatively different results. The bare coupling $`g_0^2`$ (for the Wilson gauge action) is an especially bad expansion parameter Lepage1991:zt . The obvious remedy is to rearrange the perturbative series, eliminating $`g_0^2`$ in favor of a renormalized (running) coupling, evaluated at a scale characteristic of the problem at hand. One is then faced, however, with many choices of renormalization scheme, and the question of how to determine the “characteristic scale.” In this paper we choose the Brodsky-Lepage-Mackenzie (BLM) prescription Brodsky:1983gc ; Lepage:1993xa . Once this choice is made, little subjectivity remains, so one can ask quantitatively whether one-loop BLM perturbation theory agrees with the non-perturbative method.

In the BLM method, the characteristic scale is computed from Feynman diagrams. The new information presented in this paper consists of the calculations needed to determine the BLM scales of the normalization and improvement coefficients of the vector and axial vector currents for Wilson fermions with Sheikholeslami-Wohlert action. These calculations are a by-product of our recent work on the normalization and improvement of lattice currents with heavy quarks Harada:2001fi . Details of the calculational method may be found there.

This paper is organized as follows. In Sec. II we define the lattice currents and review their description in the Symanzik effective field theory. Section III recalls the BLM prescription, focusing on points that are sometimes overlooked. Our new results for the BLM scales are given in Sec. IV. This paves the way for a systematic comparison with non-perturbative calculations of the same quantities in Sec. V. Section VI contains a few concluding remarks.

## II Lattice Currents

In this section we review the description of lattice currents with the Symanzik effective field theory. For quarks we take the Sheikholeslami-Wohlert Lagrangian Sheikholeslami:1985ij , which has an improvement coupling $`c_{\text{SW}}`$. At the tree level $`K_{\sigma F}^{[0]}=\frac{1}{4}(1c_{\text{SW}})`$, so the improvement condition $`K_{\sigma F}=0`$ requires $`c_{\text{SW}}=1+O(g^2)`$. For one-loop calculations, it is sufficient to specify $`c_{\text{SW}}`$ at the tree level. For the non-perturbative calculations cited below, $`c_{\text{SW}}1`$ is determined non-perturbatively by the methods of Ref. Luscher:1996sc .

We denote the lattice fermion field with $`\psi `$. The lattice vector and axial vector currents take the form

$`V_{\text{lat}}^\mu `$ $`=`$ $`\overline{\psi }i\gamma ^\mu \psi ac_V_{\nu }^{}{}_{\text{lat}}{}^{}\overline{\psi }\sigma ^{\mu \nu }\psi ,`$ (3)
$`A_{\text{lat}}^\mu `$ $`=`$ $`\overline{\psi }i\gamma ^\mu \gamma _5\psi +ac_A_{\text{lat}}^\mu \overline{\psi }i\gamma _5\psi .`$ (4)
The improvement couplings $`c_V`$ and $`c_A`$ should be chosen to reduce lattice artifacts, as discussed below.111The lattice currents in Eqs. (3) and (4) are useful for light quarks. For heavy quarks the “small” improvement terms become large, introducing unnecessary violations of heavy-quark symmetry. Better currents for heavy quarks are given in Refs. Harada:2001fi ; Harada:2001fj . In Symanzik’s theory of cutoff effects, the lattice currents are described by operators in a continuum effective field theory Symanzik:1979ph ; Symanzik:1983dc ; Kronfeld:2002pi ; Luscher:1996sc

$`V_{\text{lat}}^\mu `$ $``$ $`\overline{Z}_V^1\overline{q}i\gamma ^\mu qaK_V_\nu \overline{q}\sigma ^{\mu \nu }q+\mathrm{},`$ (5)
$`A_{\text{lat}}^\mu `$ $``$ $`\overline{Z}_A^1\overline{q}i\gamma ^\mu \gamma _5q+aK_A^\mu \overline{q}i\gamma _5q+\mathrm{},`$ (6)
where, as in Eq. (2), $`q`$ is a continuum fermion field whose dynamics is defined by $`_{\text{QCD}}`$. The ellipsis indicates operators of dimension five and higher. Further dimension-four operators are omitted from Eqs. (5) and (6), because they are linear combinations of those listed and others that vanish by the equations of motion. The short-distance coefficients in the effective Lagrangian—$`\overline{Z}_J`$ and $`K_J`$ ($`J=V`$, $`A`$)—are functions of $`g^2`$ and $`ma`$, and the improvement couplings $`c_{\text{SW}}`$ and $`c_J`$.

Symanzik improvement is achieved by adjusting $`c_J`$ so that $`K_J=0`$. Then $`\overline{Z}_VV_{\text{lat}}^\mu `$ and $`\overline{Z}_AA_{\text{lat}}^\mu `$ have the same matrix elements as $`\overline{q}i\gamma ^\mu q`$ and $`\overline{q}i\gamma ^\mu \gamma _5q`$, apart from lattice artifacts of order $`a^2`$. For light quarks one may expand $`\overline{Z}_J`$ in $`ma`$,

$$\overline{Z}_J=Z_J\left(1+mab_J\right),$$
(7)
and identify $`K_J`$ with the zeroth order of a small $`ma`$ expansion. At the tree level the coefficients of the normalization factor are $`Z_J^{[0]}=1`$, $`b_J^{[0]}=1`$. In addition, the coefficient of the lattice artifact is

$$K_J^{[0]}=c_J^{[0]}.$$
(8)
The improvement condition $`K_J=0`$ says that one should set $`c_J^{[0]}=0`$. Consequently, one-loop calculations are based solely on the first terms in Eqs. (3) and (4).

## III Brodsky-Lepage-Mackenzie Prescription

In this section we review the BLM prescription, following the argumentation from Ref. Lepage:1993xa . This material should be familiar, but some of the literature on non-perturbative improvement blurs the difference between BLM perturbation theory and other topics, such as “tadpole improvement” and mean-field estimates of the renormalized coupling, which are also discussed in Ref. Lepage:1993xa .

The problem is to find a reasonably accurate one-loop estimate of a quantity $`\zeta `$, here $`\overline{Z}_J`$ or $`K_J`$. In these cases, one gluon with momentum $`k`$ and propagator $`D(k)`$ appears. The contribution from the Feynman diagrams can be written

$$g_R^2\zeta ^{[1]}(p)=g_0^2\frac{d^4k}{(2\pi )^4}D(k)f(k,p)+\mathrm{},$$
(9)
where $`p`$ denotes $`k`$-independent parameters, such as external momenta. The ellipsis indicates higher-order terms that we would like to absorb into the renormalized coupling $`g_R^2`$. An important class of higher-order terms consists of the renormalization parts that dress the exchanged gluon. In the Fourier transform of the heavy-quark potential, for example, they turn $`g_0^2D(k)`$ into $`g_V^2(k)D(k)`$, where the potential $`V(q)=C_Fg_V^2(q)/q^2`$. Thus,

$$g_R^2\zeta ^{[1]}(p)=\frac{d^4k}{(2\pi )^4}g_V^2(k)D(k)f(k,p)+\mathrm{}$$
(10)
sums the renormalization parts. Other ways of dressing the gluon would lead to other physical running couplings, but they all are the same at order $`\beta _0g^4`$ Brodsky:1983gc , where $`\beta _0=11\frac{2}{3}n_f`$ is the one-loop coefficient of the $`\beta `$ function for $`n_f`$ light quarks.

If there is a characteristic scale $`q^{}`$, one can approximate

$`g_V^2(k)`$ $`=`$ $`{\displaystyle \frac{g_V^2(q^{})}{1+(\beta _0/16\pi ^2)g_V^2(q^{})\mathrm{ln}(k/q^{})^2}}`$ (11)
$`=`$ $`g_V^2(q^{})+{\displaystyle \frac{\beta _0}{16\pi ^2}}g_V^4(q^{})\mathrm{ln}(q^{}/k)^2+\mathrm{}.`$ (12)
The aim is to choose $`q^{}`$ so that higher-order terms are small, particularly those of order $`\beta _0g_V^4`$, which could be enhanced by a foolish choice of $`q^{}`$. Inserting Eq. (12) into Eq. (10) and setting the coefficient of $`\beta _0g_V^4`$ to zero yields

$$\mathrm{ln}q^{}a=\frac{{}_{}{}^{}\zeta _{}^{[1]}}{2\zeta ^{[1]}},$$
(13)
where $`a`$ is a reference short-distance scale (namely, the lattice spacing), and

$${}_{}{}^{}\zeta _{}^{[1]}(p)=\frac{d^4k}{(2\pi )^4}\mathrm{ln}(ka)^2D(k)f(k,p).$$
(14)
Thus, the BLM prescription is to set $`g_R^2=g_V^2(q^{})`$ in the one-loop approximation.

If one prefers a different renormalized coupling, one must change the scale in the appropriate way. The coupling in scheme “$`S`$” is related to the $`V`$ scheme by

$$\frac{1}{g_S^2(q)}=\frac{1}{g_V^2(q)}+\frac{\beta _0b_S^{(1)}+b_S^{(0)}}{16\pi ^2}+O(g^2),$$
(15)
where $`b_S^{(0)}`$ and $`b_S^{(1)}`$ are constants independent of $`n_f`$. The BLM scale $`q_S^{}`$ for this scheme is given by

$$\mathrm{ln}q_S^{}=\mathrm{ln}q^{}\frac{1}{2}b_S^{(1)}.$$
(16)
For example, for the modified minimal subtraction ($`\overline{\mathrm{MS}}`$) scheme, $`b_{\overline{\mathrm{MS}}}^{(0)}=8`$ and $`b_{\overline{\mathrm{MS}}}^{(1)}=5/3`$, $`q_{\overline{\mathrm{MS}}}^{}=e^{5/6}q^{}=0.435q^{}`$. With Eq. (16) one recovers the summary statement of Ref. Brodsky:1983gc , namely to absorb into $`q_S^{}`$ the $`n_f`$ dependence of the two-loop term, which enters through $`\beta _0`$.

The BLM prescription has several features that make it a natural choice in matching calculations, such as those considered in this paper. The effective field theory framework suggests using a renormalized coupling, in particular one that has a (quasi-)physical definition in both the underlying theory (here lattice gauge theory) and in the effective theory (here the Symanzik effective field theory). For quantitative purposes it is more interesting to note that

$$\frac{1}{g_S^2(q_S^{})}=\frac{1}{g_V^2(q^{})}+\frac{b_S^{(0)}}{16\pi ^2}+O(g^2),$$
(17)
so the numerical difference in the BLM expansion parameters is small, as long as $`g^2b_S^{(0)}/16\pi ^2`$ is small.

## IV Perturbative Results

In Ref. Harada:2001fi  we found for gauge group SU(3) and $`c_{\text{SW}}=1`$

$`Z_V^{[1]}`$ $`=`$ $`0.129423(6),`$ (18)
$`Z_A^{[1]}`$ $`=`$ $`0.116450(5),`$ (19)
in excellent agreement with previous work Gabrielli:1990us ; Capitani:2001xi . (Reference Capitani:2001xi  gives precise results as a polynomial in $`c_{\text{SW}}`$.) We also found (with $`c_J^{[0]}=0`$)

$`b_V^{[1]}`$ $`=`$ $`0.153239(14),`$ (20)
$`b_A^{[1]}`$ $`=`$ $`0.152189(14),`$ (21)
$`K_V^{[1]}`$ $`=`$ $`c_V^{[1]}+0.016332(7),`$ (22)
$`K_A^{[1]}`$ $`=`$ $`c_A^{[1]}+0.0075741(15),`$ (23)
which agree perfectly with Ref. Sint:1997jx . Solving the improvement condition $`K_J=0`$ at this order gives

$`c_V^{[1]}`$ $`=`$ $`0.016332(7),`$ (24)
$`c_A^{[1]}`$ $`=`$ $`0.0075741(15).`$ (25)
We also directly obtained

$$b_V^{[1]}b_A^{[1]}=0.0010444(16),$$
(26)
which is more accurate than the difference of the two numbers quoted above. In taking the difference, large contributions from the self energy cancel, but, even so, the near equality of $`b_V^{[1]}`$ and $`b_A^{[1]}`$ is a bit astonishing. The mass dependence of $`\overline{Z}_J`$ shows that $`b_V^{[1]}b_A^{[1]}`$ is not so small for the Wilson action Harada:2001fi .

In our method for computing the improvement coefficients it is easy to weight the integrands with $`\mathrm{ln}(ka)^2`$ and, thus, obtain the BLM scales. We find

$`{}_{}{}^{}Z_{V}^{[1]}`$ $`=`$ $`0.270691(19),`$ (27)
$`{}_{}{}^{}Z_{A}^{[1]}`$ $`=`$ $`0.243086(09),`$ (28)
$`{}_{}{}^{}b_{V}^{[1]}`$ $`=`$ $`0.321556(35),`$ (29)
$`{}_{}{}^{}b_{A}^{[1]}`$ $`=`$ $`0.318108(21),`$ (30)
$`{}_{}{}^{}b_{V}^{[1]}{}_{}{}^{}b_{A}^{[1]}`$ $`=`$ $`0.0034247(51),`$ (31)
$`{}_{}{}^{}c_{V}^{[1]}`$ $`=`$ $`0.0222383(15),`$ (32)
$`{}_{}{}^{}c_{A}^{[1]}`$ $`=`$ $`0.0147825(62),`$ (33)
and hence

$`q_{Z_V}^{}a`$ $`=`$ $`2.846,`$ (34)
$`q_{Z_A}^{}a`$ $`=`$ $`2.840,`$ (35)
$`q_{Z_A/Z_V}^{}a`$ $`=`$ $`2.898,`$ (36)
$`q_{b_V}^{}a`$ $`=`$ $`2.855,`$ (37)
$`q_{b_A}^{}a`$ $`=`$ $`2.844,`$ (38)
$`q_{b_Vb_A}^{}a`$ $`=`$ $`5.153,`$ (39)
$`q_{c_V}^{}a`$ $`=`$ $`1.975,`$ (40)
$`q_{c_A}^{}a`$ $`=`$ $`2.653.`$ (41)
The scales are in the expected range. The higher scale for $`b_Vb_A`$ means simply that the difference between these renormalization constants arises from very short distances. These numerical results are new; they have been obtained from two independent computer programs. As a further check, we have reproduced the values of $`q_{Z_V}^{}a`$ and $`q_{Z_A}^{}a`$ for the Wilson action ($`c_{\text{SW}}=0`$), given in Ref. Bernard:1998sx .

The dominant contributor to the “large” one-loop normalization constants, Eqs. (18)–(21), is the tadpole diagram (in Feynman gauge) of the self energy. One might expect perturbation theory to work better for quantities in which the effects of tadpole diagrams largely cancel (albeit in a gauge-invariant way). For example, $`Z_A/Z_V`$ and $`b_Vb_A`$ are tadpole free and have smaller one-loop coefficients.

Another way to remove the tadpoles is to write

$`Z_J=u_0\stackrel{~}{Z}_J,`$ (42)
$`b_J=\stackrel{~}{b}_J/u_0,`$ (43)
where $`u_0`$ is any convenient tadpole-dominated quantity. Then one can take $`u_0`$ from a non-perturbative Monte Carlo calculation and use perturbation theory for $`\stackrel{~}{Z}_J`$ and $`\stackrel{~}{b}_J`$. The corresponding one-loop coefficients are

$`\stackrel{~}{Z}_J^{[1]}=Z_J^{[1]}u_0^{[1]},`$ (44)
$`\stackrel{~}{b}_J^{[1]}=b_J^{[1]}+u_0^{[1]}.`$ (45)
Similarly, to get the BLM scale

$`{}_{}{}^{}\stackrel{~}{Z}_{J}^{[1]}={}_{}{}^{}Z_{J}^{[1]}{}_{}{}^{}u_{0}^{[1]},`$ (46)
$`{}_{}{}^{}\stackrel{~}{b}_{J}^{[1]}={}_{}{}^{}b_{J}^{[1]}+{}_{}{}^{}u_{0}^{[1]},`$ (47)
where $`{}_{}{}^{}u_{0}^{[1]}`$ is the BLM numerator \[cf. Eq. (13)\] for $`u_0`$. Below we take $`u_0^4`$ to be the average value of the plaquette, with $`u_0^{[1]}=1/12=0.08\overline{3}`$ and $`{}_{}{}^{}u_{0}^{[1]}=0.204049(1)`$. A glance at Eqs. (44)–(47) shows immediately that tadpole improvement reduces the one-loop coefficients. With tadpole improvement the BLM scales become

$`q_{\stackrel{~}{Z}_V}^{}a`$ $`=`$ $`2.061,`$ (48)
$`q_{\stackrel{~}{Z}_A}^{}a`$ $`=`$ $`1.803,`$ (49)
$`q_{\stackrel{~}{b}_V}^{}a`$ $`=`$ $`2.317,`$ (50)
$`q_{\stackrel{~}{b}_A}^{}a`$ $`=`$ $`2.289.`$ (51)
The scales are lower than without tadpole improvement, but still ultraviolet.

It is perhaps worthwhile emphasizing the difference between tadpole improvement and the BLM prescription. The aim of tadpole improvement is to re-sum large contributions appearing at order $`g^2`$ and higher, replacing the sum with a non-perturbative estimate ($`u_0`$, for example). The aim of the BLM prescription is to re-sum potentially large renormalization parts into the renormalized coupling. Although the aims are similar, they are not identical. They are not mutually exclusive, and neither is a substitute for the other.

## V Comparison to Non-perturbative Calculations

With the BLM scales in hand we can compare the prediction of one-loop BLM-improved perturbation theory with non-perturbative determinations of the improvement coefficients. We shall make the comparison in two ways. First we compare the numerical values directly, at two values of the bare coupling. Here there are two methods in the literature, one based on finite-size techniques and the Schrödinger functional Luscher:1996ck ; Luscher:1996ax ; Guagnelli:1997db , and another based on large volumes with hadronic matrix elements Bhattacharya:2001pn . The difference between these two illustrates how large power corrections to the improvement coefficients are. We also compare our results graphically, as a function of coupling, to Padé approximants given in Refs. Sint:1997jx ; Luscher:1996ck ; Luscher:1996ax . These graphs are helpful for seeing whether discrepancies in the one-loop and non-perturbative estimates arise from two-loop or power corrections.

We obtain $`\alpha _V(q^{})`$ as follows. First we compute

$$\alpha _{1\times 1}=\frac{3}{4\pi }\mathrm{ln}\mathrm{},$$
(52)
where $`\mathrm{}`$ is the ensemble average of the plaquette. Then we follow Ref. Lepage:1993xa  and take $`\alpha _V`$ to be

$$\alpha _V(3.402/a)\frac{2\alpha _{1\times 1}}{1+\sqrt{14.741\alpha _{1\times 1}}},$$
(53)
which agrees with the standard definition of $`\alpha _V`$ with an accuracy of order $`\alpha _s^3`$. The scale $`3.402/a`$ is the BLM scale for $`\mathrm{}`$. We then run from $`3.402/a`$ to $`q^{}`$ with the two-loop evolution equation. Of course, once two-loop perturbation theory is available, one would have to extend the accuracy of Eq. (53) and of the evolution.

Table 1 gives results from our perturbative calculation with non-perturbative results from the Alpha Collaboration Luscher:1996ck ; Luscher:1996ax ; Guagnelli:1997db  and from Bhattacharya *et al.* Bhattacharya:2001pn , at $`\beta =6.2`$. Table 2 gives the same comparison at $`\beta =6.0`$. Above (below) the horizontal line, we have applied the BLM prescription without (with) tadpole improvement. The error bars on the entries from Refs. Luscher:1996ck ; Luscher:1996ax ; Guagnelli:1997db ; Bhattacharya:2001pn  are statistical, and compiled in Ref. Bhattacharya:2001pn .

For the normalization factors $`Z_V`$ and $`Z_A`$, BLM perturbation theory and the non-perturbative methods agree well, within 3–4%. The difference between the two non-perturbative values of $`Z_V`$ exceeds the reported errors, but is easily explained by power correction of order $`(\mathrm{\Lambda }a)^2`$. For the tadpole-free ratio $`Z_A/Z_V`$ and for the tadpole-improved quantities $`u_0\stackrel{~}{Z}_J`$, BLM perturbation theory lies very close to the non-perturbative range. These impressions are strengthened by Fig. 1, which shows $`Z_V`$ and $`Z_A`$ as functions of $`g_0^2`$.

Circles show BLM perturbation theory, and the thin solid (dashed) lines indicate how two-loop contributions of $`\pm \alpha _V^2`$ ($`\pm 2\alpha _V^2`$) could modify the result. We show the result with and without tadpole improvement in Figs. 1(b,d) and (a,c), respectively. For the non-perturbative method, a heavy (blue) line shows the Padé approximants Luscher:1996ax

$`Z_V`$ $`=`$ $`{\displaystyle \frac{10.7663g_0^2+0.0488g_0^4}{10.6369g_0^2}},`$ (54)
$`Z_A`$ $`=`$ $`{\displaystyle \frac{10.8496g_0^2+0.0610g_0^4}{10.7332g_0^2}},`$ (55)
which deviate from the underlying calculations negligibly for $`g_0^21`$. The shaded bands behind the Padé curves show a power-correction of $`\pm (\mathrm{\Lambda }a)^2`$, with $`\mathrm{\Lambda }500\text{GeV}`$. The finite-volume result also suffers from power corrections of order $`(a/L)^2`$. They are estimated to be small by comparing calculations on lattices with $`a/L=1/8`$ and $`1/12`$ Luscher:1996ax . Also, they are parametrically smaller, because Ref. Luscher:1996ax  holds $`L\mathrm{\Lambda }2`$ for all $`g_0^2`$.

Next let us turn to the $`O(ma)`$ corrections to the normalization factors, $`b_V`$ and $`b_A`$. There is only one calculation of $`b_A`$ Bhattacharya:2001pn , so let us concentrate first on $`b_V`$. The two non-perturbative results for $`b_V`$ agree perfectly with each other (see the Tables), but they deviate significantly from one-loop BLM perturbation theory. Some insight can be gleaned from Fig. 2, which shows $`b_V`$ as a function of $`g_0^2`$.

The non-perturbative method is represented with the Padé approximant Sint:1997jx

$$b_V=\frac{10.7613g_0^2+0.0012g_0^40.1136g_0^6}{10.9145g_0^2},$$
(56)
with light (blue) shading for a power correction $`\pm \mathrm{\Lambda }a`$. In finite volume there is also a power correction to $`b_V`$ of order $`a/L`$; by construction it applies to $`b_V1`$ Luscher:1996ax , but now $`L`$ with $`a`$ varies such that $`a/L=1/8`$ for all $`g_0^2`$. We model this effect as $`(b_V1)(1\pm \frac{1}{8})`$, shown in the darker (pink) shading in Fig. 2. Judging from its size and shape, the deviation seen in Fig. 2 looks less like a two-loop effect than a combination of power corrections of order $`a/L`$ and $`\mathrm{\Lambda }a`$. (Similar conclusions are reached in Ref. Bhattacharya:2001pn .) There is almost no difference whether one applies tadpole improvement to $`b_V`$ or not, once the BLM prescription is applied. These two approximations truncate higher orders of the perturbation series differently, substantiating the idea that the discrepancy is a power correction.

The non-perturbative calculation of $`b_A`$ agrees with one-loop BLM perturbation theory. Note, however, that Ref. Bhattacharya:2001pn  obtains $`b_V`$ and $`b_Vb_A`$ directly, and then $`b_A=b_V(b_Vb_A)`$. The agreement between BLM perturbation theory and the non-perturbative results for $`b_V`$ and $`b_Vb_A`$ is not good, so the agreement for $`b_A`$ may be an accident. Since the coefficient $`b_V^{[1]}b_A^{[1]}`$ in Eq. (26) is remarkably small, the two-loop contribution could be as large as the one-loop term. Furthermore, inspection of Fig. 14 in Ref. Bhattacharya:2001pn  suggests that a fit to the three smallest masses would yield a smaller value of $`b_Vb_A`$. We consider the comparison of $`b_A`$ and $`b_Vb_A`$ to be unsettled pending a two-loop calculation and a more robust non-perturbative calculation.

In any case, the mild disagreement on $`b_V`$ and $`b_Vb_A`$ is not of much practical importance. For the sake of argument, suppose $`ma<0.1`$, which holds for the light quarks for which the currents were designed. Then power corrections in $`b_J`$, at fixed $`a`$, lead to an uncertainty in a decay constant or a form factor of only a few per cent. After a continuum limit extrapolation, these uncertainties will not be important.

Now let us turn to the coefficients $`c_J`$ of the improvement terms in Eq. (3) and (4). At the tabulated values of $`\beta `$, the non-perturbative and BLM calculations of $`c_A`$ do not agree at all. At $`\beta =6.0`$ (Table 2) the two non-perturbative calculations also do not agree with each other. Figure 3(a) shows $`c_A`$ as a function of $`g_0^2`$,

using the Padé approximant Luscher:1996ck

$$c_A=0.00756g_0^2\frac{10.748g_0^2}{10.977g_0^2}$$
(57)
to represent the non-perturbative calculations. The disagreement between BLM perturbation theory and Eq. (57) sets in for $`g_0^2>0.9`$. There are two reasons to suspect that the discrepancy stems from a power correction of order $`\mathrm{\Lambda }a`$ to the results of Ref. Luscher:1996ck . First, Fig. 3(a) shows that it has the shape and size of such a power correction. Second, the extracted value of $`c_A`$ depends on the lattice derivative used to define the current Collins:2001mm . Note Bhattacharya:2001pn  that errors in $`c_A`$ propagate to $`c_V`$, because in the Ward identities $`c_A`$ is multiplied by large hadronic matrix elements such as $`am_K^2/m_sa\times 2.5\text{GeV}`$. This enhancement also explains why Eq. (57) leads to worse scaling in $`f_\pi `$ Collins:2001mm . Figure 3(a) also includes the non-perturbative results of Refs. Bhattacharya:2001pn ; Collins:2001mm . The difference between those points and BLM perturbation theory could be a modest two-loop effect or a small power correction.

For $`c_V`$, the two non-perturbative results agree neither with each other, nor with BLM perturbation theory. The Alpha Collaboration has only a preliminary calculation Guagnelli:1997db . We have taken the liberty of extracting results from Fig. 3 of Ref. Guagnelli:1997db  and fitting them to a Padé formula. The leading behavior is fixed to Eq. (24), and we obtain

$$c_V=0.01633g_0^2\frac{10.257g_0^2}{10.963g_0^2}.$$
(58)
Figure 3(b) plots Eq. (58), the underlying points Guagnelli:1997db , the non-perturbative results from hadronic correlation functions Bhattacharya:2001pn , and BLM perturbation theory. As usual we show possible power corrections to Eq. (58) of order $`\pm \mathrm{\Lambda }a`$, as well as the size of typical two-loop effects. At small $`g_0^2`$, there is good agreement with (BLM) perturbation theory, but once $`g_0^2>0.9`$, there is a sharp turnover. It is probably a power correction, possibly exacerbated by power corrections to $`c_A`$ as modeled by Eq. (57). With hadronic correlation functions Bhattacharya:2001pn  the non-perturbative value of $`c_V`$ is half or a third as large. It is not clear at present whether the discrepancy between Ref. Bhattacharya:2001pn  and BLM perturbation theory is a power correction to the former or a sizable two-loop correction to the latter.

We should also mention that BLM perturbation theory works better than several forms of mean-field perturbation theory (let alone bare perturbation theory). In Table 3 we list several choices for $`\alpha _s`$:

$`\alpha _0=g_0^2/4\pi ,`$ (59)
$`\stackrel{~}{\alpha }_0=\alpha _0/u_0,`$ (60)
as well as $`\alpha _{1\times 1}`$ \[Eq. (52)\] and $`\alpha _{\overline{\mathrm{MS}}}(q_{\overline{\mathrm{MS}}}^{})`$ \[Eq. (17)\]. With only one-loop expansions available, the mean-field choices $`\stackrel{~}{\alpha }_0`$ and $`\alpha _{1\times 1}`$ give smaller corrections, and one-loop perturbation theory falls short even when power corrections are negligible. The consistency of BLM-$`V`$ perturbation theory for $`Z_V`$, $`Z_A`$, and $`Z_A/Z_V`$ indicates that the BLM prescription does indeed re-sum an important class of higher-order contributions. On the other hand, the coupling $`\alpha _{\overline{\mathrm{MS}}}(q_{\overline{\mathrm{MS}}}^{})`$ seems, empirically, to work less well. In continuum perturbative QCD, it usually does not matter whether one adopts $`\alpha _V(q_V^{})`$, $`\alpha _{\overline{\mathrm{MS}}}(q_{\overline{\mathrm{MS}}}^{})`$ or some other renormalized coupling (at the BLM scale), once two-loop effects are included. It would not be surprising for the same to hold for short-distance quantities in lattice gauge theory, such as improvement coefficients.

## VI Conclusions

In this paper we have compared non-perturbative calculations of several improvement coefficients to perturbation theory with the BLM prescription. Previously this could not be done, because the “BLM numerators” in Eqs. (27)–(33) were not available. We find that, for the scale-independent quantities considered here, the integration of the $`\mathrm{log}k^2`$-weighted integrals is numerically straightforward.

BLM perturbation theory for the current normalization factors $`Z_J`$ agrees very well with non-perturbative calculations of the same quantities. Here the leading power correction is only of order $`(\mathrm{\Lambda }a)^2`$, and the small deviations can probably be removed with a two-loop calculation. Note that generalizations of the BLM method for higher-order perturbation theory have been considered in continuum perturbative QCD Brodsky:1994eh  and in lattice gauge theory Hornbostel:2000ey .

For the improvement coefficients $`b_J`$ and $`c_J`$, the leading power corrections are of order $`\mathrm{\Lambda }a`$ (and in the Schrödinger functional also of order $`a/L=1/8`$), while some of the one-loop coefficients are small. It is consequently difficult to diagnose the discrepancies. By noting the size and dependence on $`g_0^2`$ of the differences, we concur with the authors of Refs. Bhattacharya:2001pn ; Collins:2001mm , namely, that power corrections contaminate the non-perturbative results. In particular, it seems unlikely that higher orders in perturbative series could explain all discrepancies between one-loop BLM perturbation theory and the results from Refs. Luscher:1996ck ; Luscher:1996ax ; Guagnelli:1997db .


text:
# Optimization of three-dimensional micropost microcavities for cavity quantum electrodynamics

## 1 Introduction

Spontaneous emission is not an intrinsic property of an isolated atom, but is rather a property of an atom coupled to its electromagnetic vacuum environment. The spontaneous emission rate is directly proportional to the density of electromagnetic states that a spontaneously emitted photon can couple to, and can be modified with respect to its value in free space by placing the atom in a cavity . The experimental demonstrations of the inhibition and enhancement of spontaneous emission rate were carried out starting in the mid-70’s , using atoms coupled to single mirrors, planar cavities, or spherical Fabry-Perot resonators. Advances in microfabrication techniques enabled the construction of high-quality semiconductor micropost and microdisk microcavities in the late 80’s and early 90’s, and ignited interest in solid-state cavity quantum electrodynamics (QED) experiments . In 1987, photonic-crystal structures were proposed as promising candidates for strong spontaneous emission modification , but the first experimental results on photonic crystal microcavities were not reported until a decade later .

Cavity-QED phenomena in the low-$`Q`$ (weak-coupling) regime, as well as in the high-$`Q`$ (strong-coupling) regime, can be used in construction of high-efficiency light-emitting diodes, low-threshold lasers, and single photon sources. A powerful property of solid-state microcavities is that a single narrow-linewidth emitter (quantum dot) can be embedded in them during the growth process, enabling cavity-field interaction with such artificial atom . Due to imperfections in fabricated structures, un-optimized structure parameters, and the inability to precisely control position of a quantum dot, only phenomena in the low-$`Q`$ regime have been observed so far.

The first successful optical characterizations of photonic-crystal microcavities with quantum dots were performed recently . $`Q`$ factors as large as $`2800`$ were reported, together with mode volumes as small as $`0.5(\lambda /n)^3`$, where $`\lambda `$ is the optical wavelength, and $`n`$ is the refractive index of the dielectric material . The possibility of improving the quality factor while preserving such a small mode volume makes these structures good candidates for cavity QED, in particular with neutral atoms (due to a strong field intensity in the air region for the optimized cavity designs) . So far, this has not been demonstrated experimentally.

The advantages of microposts relative to other microcavities are that the light escapes in the normal direction to the sample in a single-lobed Gaussian-like pattern, and that it is relatively straightforward to isolate a single quantum dot in a post. However, in order to observe such cavity QED phenomena as strong coupling with a single dot or single-dot lasing in these structures, a number of design and fabrication issues have to be addressed. In this article, we present the optimization of micropost parameters (illustrated in Fig. 1), in order to maximize the quality factor and minimize the volume of the fundamental ($`HE_{11}`$) mode (whose field pattern is shown in Fig. 2). We show that both strong coupling cavity QED with a single quantum dot, and single-quantum-dot lasing are possible in the optimized micropost microcavity.

All analyzes presented in this article are performed by the finite-difference time-domain (FDTD) method, which enables accurate modelling of the electromagnetic properties of structures with complex geometries. The rotational symmetry of micropost microcavities allows us to use a cylindrical FDTD algorithm and reduce the order of the computer memory requirements from $`N^3`$ to $`N^2`$, where $`N`$ represents a linear dimension of the computational domain. The method used is described in detail in our earlier publication .

## 2 Motivation for maximizing the ratio of quality factor to mode volume

Let us assume that a single quantum dot is isolated in a microcavity, and that the transition frequency from the one-exciton state to the zero-exciton state is on resonance with the fundamental optical cavity mode frequency $`\omega `$. Under these conditions, the system can be modelled in the same way as a single two-level atom coupled to a single cavity mode, and described by the Jaynes-Cummings Hamiltonian . The coupling parameter $`g`$ between the exciton and the cavity field reaches its maximum value equal to the vacuum Rabi frequency $`g_0`$, when the dot is located at the point of the maximum electric field intensity, and when the excitonic dipole moment is aligned with the electric field:

$$g_0=\frac{\mu }{\mathrm{}}\sqrt{\frac{\mathrm{}\omega }{2ϵ_MV}},$$
(1)
where $`ϵ_M`$ is the dielectric constant at the location of the exciton, $`\mu `$ is the dipole moment matrix element between the one-exciton and zero-exciton states, and $`V`$ is the cavity mode volume, defined as

$$V=\frac{ϵ(\stackrel{}{r})|E(\stackrel{}{r})|^2d^3\stackrel{}{r}}{\mathrm{max}\left[ϵ(\stackrel{}{r})|E(\stackrel{}{r})|^2\right]}.$$
(2)
Depending on the ratio of the coupling parameter $`g`$ to the cavity field decay rate $`\kappa =\omega /2Q`$ and the excitonic dipole decay rate $`\gamma `$, we can distinguish two regimes of coupling between the exciton and the cavity field: strong coupling, for $`g>\kappa ,\gamma `$, and weak coupling, for $`g<\kappa ,\gamma `$. In the strong-coupling case, the exciton is coherently coupled to the cavity field, spontaneous emission is reversible, and vacuum Rabi oscillation occurs. On the other hand, in the weak-coupling case, the spontaneous emission is irreversible, and the spontaneous emission decay rate $`\mathrm{\Gamma }`$ is

$$\mathrm{\Gamma }=g^2\frac{4Q}{\omega }.$$
(3)
The spontaneous emission rate of an exciton in free space, on the other hand, is given by

$$\mathrm{\Gamma }_0=\frac{\omega ^3\mu ^2}{3\pi ϵ_0\mathrm{}c^3}.$$
(4)
The ratio of $`\mathrm{\Gamma }`$ to $`\mathrm{\Gamma }_0`$ is called the Purcell factor . For an exciton positioned at the maximum of the field intensity and aligned with the electric field, the Purcell factor is equal to

$$F_0=\frac{3Q\lambda ^3ϵ_0}{4\pi ^2Vϵ_M}.$$
(5)
We usually define the Purcell factor $`F`$ as the spontaneous emission rate enhancement relative to the bulk material. The spontaneous emission rate in the bulk material with refractive index $`n_h`$ is enhanced $`n_h`$ times with respect to its value in free space, which implies that $`F=F_0/n_h`$.

If the Purcell factor is much greater than one, the exciton will radiate much faster in the cavity than in free space. The radiative-rate enhancement is proportional to the ratio of the quality factor to the volume of the cavity mode, according to Eqn. 5. The Purcell factor increases with $`Q/V`$ only to the point where the coupling parameter $`g`$ becomes larger than the decay rates of the system ($`\kappa `$ and $`\gamma `$). At that point, the coupled exciton-cavity system enters the strong-coupling regime.

Increasing $`Q/V`$ can also lead to a reduction in laser threshold. The fraction of the light emitted by an exciton that is coupled into one particular cavity mode is known as the spontaneous emission coupling factor $`\beta `$, and is related to the Purcell factor via the following expression:

$$\beta =\frac{F}{1+F}.$$
(6)
Therefore, if the emission rate of an exciton is strongly enhanced by its interaction with a cavity mode, the fraction of spontaneous emission going into all other modes ($`1\beta `$) is reduced. The fraction of spontaneous emission going into non-lasing modes is one of the fundamental losses in a laser, and by decreasing it, one can lower the laser threshold.

Of particular interest would be a single-dot laser, which represents an ultimate microscopic limit for semiconductor lasers. The realization of such a device would allow physical investigations similar to those afforded by the single-atom laser . Lasing of such a microscopic system would occur when the mean spontaneously-emitted photon number $`n_{sp}`$ in the laser mode becomes larger than one :

$$n_{sp}=\frac{\beta \tau _{ph}N_A}{\tau _{sp}}=\frac{N_A(\mathrm{\Gamma }_0n_h)}{\omega /Q}\frac{F^2}{1+F}1,$$
(7)
where $`\tau _{sp}=1/\mathrm{\Gamma }`$, $`\tau _{ph}=Q/\omega `$, and $`N_A`$ is the average probability over time that the quantum dot contains an exciton.

One of the most interesting applications of cavity QED is the construction of efficient sources of single photons . Single-photon sources are useful for quantum cryptography , quantum computation , quantum networking , and random number generators . A single quantum dot can be used to generate single photons, and the output coupling efficiency can be enhanced by cavity QED. In other words, by changing the cavity parameters ($`Q/V`$) and the quantum dot location, we can control the probability of coupling this spontaneously emitted single photon into the mode of interest, and subsequently coupling it into the communication channel.

## 3 Micropost microcavities

Micropost microcavities consist of a high-refractive-index region (spacer) sandwiched between two dielectric mirrors, as shown in Fig. 1. Confinement of light in these structures is achieved by the combined action of distributed Bragg reflection (DBR) in the longitudinal direction (along the post axis), and total internal reflection (TIR) in the transverse direction (along the post cross-section). The microposts analyzed in this paper are rotationally symmetric around the vertical axis. The DBR mirrors can be viewed as one dimensional (1D) photonic crystals generated by stacking high- and low-refractive-index disks on top of each other. The microcavity is formed by introducing a defect into this periodic structure. The periodicity of the photonic crystal is denoted as $`a`$, the thickness of the low-refractive-index disks is $`t`$, the diameter of the disks is $`D`$, and the refractive indices of the low- and high-refractive-index regions are $`n_l`$ and $`n_h`$, respectively. The defect is formed by increasing the thickness of a single high-refractive-index disk from $`(at)`$ to $`s`$, as shown in Fig. 1. The number of photonic crystal periods above and below the defect region (i.e., the number of DBR pairs) is labelled as $`MPT`$ and $`MPB`$, respectively.

The mode of interest to us is the doubly-degenerate fundamental ($`HE_{11}`$) mode, whose field pattern is shown in Fig. 2. The parallel component of the electric field is dominant in this mode, and has an antinode in the center of the spacer. Furthermore, in this central plane, the electric field is practically linearly polarized along the vertical axis of the micropost, while there is a small deviation from the linear polarization at larger distances from this axis.

The rule of thumb generally used for designing microposts is to make mirror layers one-quarter wavelength thick, and to choose the optical thickness of the spacer equal to the target wavelength. In the case of a planar DBR cavity (with $`D\mathrm{}`$), this choice of parameters leads to the maximum reflectivities of the mirrors and the maximum $`Q`$ factor of the cavity mode: the cavity operates at the Bragg wavelength, for which the partial reflections from all high- and low-refractive-index interfaces add up exactly in phase. However, the strength of the cavity QED phenomena is proportional to the ratio of the cavity $`Q`$ factor to the mode volume $`V`$, as discussed in the previous section, and we will try to design microposts in such a way that this ratio is maximized.

In our earlier work , we analyzed the $`Q`$ factor of the $`HE_{11}`$ mode in a GaAs/AlAs micropost as the cavity diameter was tuned between 0.5 $`\mu `$m and 2 $`\mu `$m. The remaining cavity parameters were chosen according to the large-cavity rule of thumb, i.e. in such a way that the cavity would operate at the Bragg wavelength for $`D\mathrm{}`$. When the cavity diameter was decreased from 2 $`\mu `$m to 0.5 $`\mu `$m, the mode volume decreased by a factor of almost ten, from $`19.2(\lambda /n_h)^3`$ to $`2(\lambda /n_h)^3`$, while the cavity $`Q`$ dropped by only a factor of two, from $`11500`$ to $`5000`$. Thus, in order to maximize the ratio of the quality factor $`Q`$ to the mode volume $`V`$, we need to explore structures with small diameters $`D`$, and try improve their $`Q`$ factors.

The reduction in $`Q`$ with decrease in $`D`$ is due to the combination of two loss mechanisms: longitudinal loss through DBR mirrors, and transverse loss due to imperfect TIR confinement in the transverse direction. Let us address the longitudinal loss first. The decrease in the post diameter $`D`$ implies a change in the dispersion relation of the 1D photonic crystal, and the size and position of its band gap, as illustrated in Fig. 3. In this figure, it is assumed that the high- and low-refractive-index regions of the photonic crystal consist of GaAs and AlAs, with refractive indices of $`n_h=3.57`$ and $`n_l=2.94`$, and thicknesses of 70 nm and 85 nm, respectively, or that they consist of GaAs and AlxGa1-xAs, with refractive indices of $`n_h=3.57`$ and $`n_l=3.125`$, and thicknesses of 70 nm and 80 nm, respectively. When the diameter $`D`$ decreases, the frequencies of the band gap edges increase, and the size of the band gap decreases. For structure diameters larger than 2 $`\mu `$m, band gap edges can be approximated by their values at $`D\mathrm{}`$. Therefore, as $`D`$ decreases, the blue shift of the cavity mode wavelength $`\lambda `$ increases relative to the target wavelength at which the 1D cavity operates . Simultaneously, the size of the photonic band-gap decreases, implying that the cavity mode is less confined in the longitudinal direction than in the planar cavity case.

The cavity mode is strongly localized in real space, and consequently delocalized in Fourier space ($`k`$-space), meaning that it consists of a wide range of wave-vector components. Some of these components are not confined in the post by TIR; i.e., they are positioned above the light line, where they can couple to radiative modes, leading to transverse loss. A cavity mode which is strongly confined in the longitudinal direction by high-reflectivity mirrors is delocalized in Fourier space and suffers large transverse loss. Similarly, a mode that is delocalized in the longitudinal direction is more localized in Fourier space and suffers less transverse loss. Therefore, when optimizing the quality factor of three-dimensional microposts, there is a tradeoff between these two loss mechanisms.

In the middle of a large band-gap, the longitudinal confinement is strongest, but the $`Q`$ factor is limited by transverse loss. By shifting the resonant wavelength away from the mid-gap (e.g., by tuning the thickness of the cavity spacer) one can delocalize the mode in real space, localizing it more strongly in Fourier space, reducing the contribution of wave-vector components above the light line, and thereby decreasing the transverse radiation loss. Eventually, as the mode wavelength approaches the band-gap edges, the loss of longitudinal confinement starts to dominate and $`Q`$ drops. Therefore, in the microposts with high reflectivity mirrors and finite diameter, it is expected that the maximum $`Q`$ will be located away from the mid-gap position.

Moreover, since the mode wavelength can be tuned from the mid-gap towards any of the two band-gap edges, two local maxima of $`Q`$ (i.e., a double peak behavior in $`Q`$ vs. mode wavelength) are expected. Besides detuning the mode wavelength from the mid-gap, we can also suppress the transverse loss by relaxing the mode slightly in the longitudinal direction, i.e., by reducing the reflectivities of photonic crystal mirrors and decreasing the band gap size. This can be achieved by shrinking the cavity diameter, or by changing the photonic crystal parameters (e.g., by reducing the refractive-index contrast).

In this article, we study both of these approaches to $`Q`$ optimization: tuning the mode wavelength away from the mid-gap by changing the spacer thickness, and tuning the mirror reflectivities by changing photonic crystal parameters or cavity diameter. We also show that the employment of very high reflectivity mirrors cannot lead to high-$`Q`$ cavities with small diameters, as the transverse radiation loss is high, resulting from very strong mode localization in the longitudinal direction.

## 4 Maximizing the ratio of quality factor to mode volume for the fundamental mode in a micropost microcavity

### 4.1 Tuning the cavity diameter and the cavity spacer

In our earlier work , we analyzed the $`Q`$ factor of the $`HE_{11}`$ mode in a GaAs/AlAs micropost as the cavity diameter was tuned between 0.5 $`\mu `$m and 2 $`\mu `$m. The remaining cavity parameters were chosen in such a way that the cavity would operate at the Bragg wavelength for $`D\mathrm{}`$ ($`a=155`$ nm, $`t=85`$ nm, $`s=280`$ nm, $`n_h=3.57`$ and $`n_l=2.94`$). The number of DBR mirror pairs on top and bottom were $`MPT=15`$ and $`MPB=30`$, respectively.

Let us first study the $`HE_{11}`$ mode as the diameter is decreased below 0.5 $`\mu `$m, keeping all other structure parameters the same as above. In order to tune the mode frequency within the band gap, we tune the spacer thickness $`s`$. Results for $`\lambda `$, $`Q`$, $`V`$, and $`Q/V`$ are shown in Figs. 4 and 5. From Fig. 3, we see that the band gaps in these structures extend from 875 nm to 969 nm, from 850 nm to 920 nm, and from 790 nm to 850 nm, for structure diameters of 0.5 $`\mu `$m, 0.4 $`\mu `$m, and 0.3 $`\mu `$m, respectively. As we have noted previously, when $`D`$ decreases, the band gap edges shift towards lower wavelengths, and the size of the band gap decreases. The cavity mode wavelength is blue-shifted in this process, as can be seen in Figure 4.

The mode volume $`V`$ is minimized when the mode wavelength is located near the middle of the band gap. For the structures with $`D`$ equal to 0.4 $`\mu `$m and 0.3 $`\mu `$m, the maximum $`Q`$ factor also occurs close to the mid-gap. Different behavior is seen for the structure with $`D`$ equal to 0.5 $`\mu `$m, which has a local minimum of $`Q`$ at mid-gap and exhibits a double-peak behavior.

The double-peak behavior was already introduced in the previous section. In the middle of the band gap, where the longitudinal mode confinement is strongest and the mode volume is minimum, the radiation loss in the transverse direction is high, and the $`Q`$ factor is degraded. By shifting the resonant wavelength away from the mid-gap, the mode is delocalized in real space, leading to a reduction in the transverse radiation loss (e.g., at the positions of the two peaks in $`Q`$). Eventually, as the mode wavelength approaches the band-gap edges, the loss of longitudinal confinement starts to dominate, $`Q`$ drops, and the mode volume increases.

To support this explanation, we analyze the same structure, with $`D=0.5`$ $`\mu `$m, but with the number of mirror pairs on top ($`MPT`$) increased from $`15`$ to $`25`$. As expected, at mid-gap, $`Q`$ does not increase significantly with $`MPT`$. The mode there is already strongly confined in the longitudinal direction, and the addition of extra pairs does not change the longitudinal loss. The modal $`Q`$ factor is determined by the radiation loss in the transverse direction, which is independent of $`MPT`$. On the other hand, the $`Q`$’s at the two peaks increase with $`MPT`$. At these points, the mode is not confined as well in the longitudinal direction, and longitudinal loss can be reduced by adding more mirror pairs.

As an even stronger demonstration of our explanation for the double-peak behavior, we separate the radiation loss into the loss above the top micropost surface ($`L_a`$), and the loss below it ($`L_b`$). The total $`Q`$ is a combination of two newly introduced quality factors, $`Q_a`$ and $`Q_b`$, which are inversely proportional to $`L_a`$ and $`L_b`$, respectively:

$$\frac{1}{Q}=\frac{1}{Q_a}+\frac{1}{Q_b}.$$
(8)
It follows from their definition that $`Q_a`$ and $`Q_b`$ are measures of the longitudinal and transverse loss, respectively. We analyze two sets of structure parameters, corresponding to the local maximum or minimum in $`Q`$. For $`s=270`$ nm and $`D=0.5`$ $`\mu `$m (local maximum), we calculate $`Q_a14500`$ and $`Q_b13910`$, while, for $`s=290`$ nm and $`D=0.5`$ $`\mu `$m (local minimum), we calculate $`Q_a16000`$ and $`Q_b5100`$. These results show that the local minimum in $`Q`$ is due to an increase in the transverse loss, manifested as a drop in $`Q_b`$.

Let us now address the single-peak behavior of $`Q`$ as a function of cavity spacer thickness, when $`D`$ is equal to 0.4 $`\mu `$m or 0.3 $`\mu `$m. Structures with smaller diameters have smaller band gaps, as illustrated in Fig. 3, and the cavity modes are more delocalized in the longitudinal direction, relative to the structure with $`D=0.5`$ $`\mu `$m. The defect modes must therefore be more localized in Fourier space, and will thus suffer less radiation loss in the transverse direction. This implies that the $`Q`$ factors are determined mostly by the longitudinal loss. They reach their maxima at the mid-gap, where the mode volume is minimum, and the longitudinal confinement is strongest.

The maximum $`Q/V`$ ratio of almost $`6000`$ (where $`V`$ is measured in cubic wavelengths in the high-refractive index material) is achieved for the structure with $`D=0.4`$ $`\mu `$m. For this structure, the $`Q`$ factor is close to $`9500`$, and the mode volume is $`1.6(\lambda /n_h)^3`$. For $`D=0.4`$ $`\mu `$m, a variation in the thicknesses of the mirror layers allows us to achieve a small increase in the $`Q`$ factor, to $`10500`$, and in the $`Q/V`$ ratio, to $`6500`$. This result is obtained for $`a=155`$ nm, $`t=75`$ nm, and $`s=290`$ nm.

In the introduction, we mentioned that an advantage of microposts, relative to other solid state microcavities, is that the light escapes from them in a single-lobed Gaussian-like pattern, normal to the sample surface. In order to show this, we calculate the far-field radiation pattern from a micropost with $`D=0.5`$ $`\mu `$m and $`s=270`$ nm. We are unable to directly compute the far field by employing the FDTD method, as we are limited by our computer memory size. However, we can estimate the far field from the Fourier transform of the near field, using the method described in Ref. . The calculated radiation pattern is shown in Fig. 6. Its resolution is limited by the resolution that we can achieve in Fourier space, or, more precisely, by the number of pixels in the light cone. This, in turn, is dictated by the size of the computational domain. The best resolution in Fourier space that we can obtain with a reasonable size of the computational domain is seven pixels per light cone radius. Nonetheless, the computed radiation pattern demonstrates that even microposts with small diameters can emit light in a Gaussian-like pattern. The FWHM of the emission lobe shown is approximately equal to $`50^o`$.

### 4.2 Other material systems

#### 4.2.1 GaAs/AlxGa1-xAs cavities

In the previous section of this article, we stated that a potential route to maximizing $`Q`$ for small micropost diameters is the construction of a photonic crystal with a small refractive-index perturbation. As the perturbation gets smaller, the cavity mode becomes more delocalized in real space, and consequently more localized in Fourier space. This, in turns, leads to reduction in the transverse radiation loss. Furthermore, the cavity resonance can be located at lower frequencies, where the density of free-space radiation modes is smaller. In order to compensate for the increased longitudinal loss, we need to put more mirror pairs on top of these structures.

We will now analyze a micropost with the following parameters: $`a=150`$ nm, $`t=80`$ nm, $`MPT=25`$, $`MPB=30`$, $`n_h=3.57`$, and $`n_l=3.125`$. This choice of refractive indices corresponds to GaAs/AlxGa1-xAs layers. Both the cavity diameter $`D`$ and the spacer thickness $`s`$ are tuned. The positions of the band gap edges as a function of $`D`$ are illustrated in Fig. 3. By comparing to the positions of the band gap edges for the GaAs/AlAs system, we confirm that the band gap in the GaAs/AlxGa1-xAs system is shifted to lower frequencies, and that its size is decreased. This affects the $`HE_{11}`$ mode dramatically, as can be seen in Figs. 7 and 8.

By comparing Fig. 8 to Fig. 5, we can see that the mode volume increases when the refractive-index contrast is reduced, as a result of the reduction in band-gap size. Even though $`Q`$ larger than 14000 can be achieved for $`D=0.5`$ $`\mu `$m, $`V`$ also increases, and the maximum $`Q/V`$ ratio is similar to that calculated for the GaAs/AlAs system. Furthermore, this $`Q/V`$ ratio can be achieved in the GaAs/AlAs system with fewer top mirror pairs. Longitudinal loss dominates in the GaAs/AlxGa1-xAs system, and $`Q`$ vs. $`s`$ plots demonstrate a single-peak behavior.

When the number of mirror pairs on top is reduced from 25 to 20, the peak $`Q`$ factor of the GaAs/AlxGa1-xAs micropost with diameter of 0.5 $`\mu `$m drops from around 14000 to 4000, showing that the the longitudinal loss is dominant in this case, and a large number of mirror pairs is necessary to achieve large $`Q`$ factors.

#### 4.2.2 GaAs/AlOx cavities

From the results already presented in this article, it is clear that a material system with a high refractive-index contrast, such as GaAs/AlOx, is not a good choice for high $`Q`$, small mode-volume microposts. High refractive index contrast can certainly produce larger band gaps, and thereby provide a better longitudinal confinement of the cavity mode. However, if the contrast is increased, the mode suffers more radiation loss in the transverse direction, which limits its $`Q`$ factor. To confirm this, we analyzed a structure with $`n_h=3.57`$, $`n_l=1.515`$, $`a=235`$ nm, $`t=165`$ nm, $`MPT=15`$, and $`MPB=30`$, for different $`D`$. We were unable to obtain good mode localization for $`D<0.8`$ $`\mu `$m, and the calculated $`Q`$ factors were under $`250`$. For $`D=0.8`$ $`\mu `$m, the mode has $`Q=600`$ and $`\lambda =947`$ nm. If we keep increasing $`D`$ to 1.3 $`\mu `$m, $`Q`$ factors remain below $`1000`$.

## 5 Cavity quantum electrodynamics with microposts

The question that we would like to address in this section is whether such cavity QED phenomena as strong coupling or single-dot lasing can be observed in the optimized microposts. Let us revisit our best design, with $`Q10^4`$, $`V=1.6(\lambda /n_h)^3`$, $`D=0.4\mu `$m, $`\lambda =885`$ nm, and the cavity field decay rate $`\kappa =(\pi c)/(\lambda Q)=106`$ GHz.

By combining Eqns. 1 and 4, the Rabi frequency $`g_0`$ of a system on resonance can be expressed as

$$g_0=\frac{\mathrm{\Gamma }_0}{2}\sqrt{\frac{V_0}{V}},$$
(9)
where $`V_0=(3c\lambda ^2ϵ_0)/(2\pi \mathrm{\Gamma }_0ϵ_M)`$. Let us assume that a quantum dot exciton without a cavity has a typical homogenous linewidth $`\gamma _h=20`$ GHz, and a radiative lifetime of 0.5 ns, corresponding to a spontaneous emission rate of $`\mathrm{\Gamma }=2`$ GHz. The free-space spontaneous emission rate is $`\mathrm{\Gamma }_0=\mathrm{\Gamma }/n_h=0.56`$ GHz. The Rabi frequency for our optimized cavity, calculated from Eqn. 9, is equal to $`g_0=400\mathrm{\Gamma }_0=224`$ GHz. If we assume that the quantum dot is located in the center of the micropost and that its dipole is aligned with the electric field, we have $`g=g_0`$. Strong coupling is therefore possible in this case, since $`g_0>\kappa ,\pi \gamma _h`$. The minimum quality factor necessary to achieve strong coupling is approximately equal to 5000. This provides a reasonable margin for $`Q`$ degradation due to fabrication imperfections.

Is strong coupling possible with larger diameter microposts, such as $`D=2\mu `$m? The mode volume in such a structure is on the order of $`20(\lambda /n_h)^3`$, as we mentioned previously. For the same quantum dot, with $`\mathrm{\Gamma }_0=0.56`$ GHz, placed in the center of this large cavity, the Rabi frequency is $`g_0=60`$ GHz. For our experimentally observed homogenous broadening $`\gamma _h=20`$ GHz, it is impossible to reach strong coupling, since $`\pi \gamma _h>g_0`$. Even if the homogenous linewidth were reduced to 2 GHz (i.e., if the homogeneous broadening were entirely due to radiative decay), the $`Q`$ factor required to achieve strong coupling would be on the order of $`2\times 10^4`$. We therefore conclude that large-diameter microposts are not promising candidates for the observation of strong coupling with a single quantum dot.

Designs of two-dimensional photonic crystal microcavities in free-standing membranes were recently proposed that allow for very strong coupling between the cavity field and a neutral atom trapped in one of photonic crystal holes . We will now address the feasibility of strong coupling with a single quantum dot in these structures, and compare them to our micropost designs. These photonic-crystal microcavities can localize light into mode volumes equal to $`1/2(\lambda /n_h)^3`$, with $`Q`$ on the order of $`10^4`$. However, since the field intensity is strongest in or around the defect air hole (where a neutral atom would be trapped), it is almost impossible to place a quantum dot at the point where its interaction with the cavity field would be strongest. For example, if the dot is placed at the point where the field intensity is $`60\%`$ of its maximum value, the Rabi frequency remains the same as for our best micropost design ($`g=224`$ GHz), despite a three-fold decrease in the mode volume. The quality factor is in the same range as for the optimized microposts, which implies that the potential of these structures to achieve strong coupling with single quantum dots is similar to that of microposts.

What about single-dot lasing in microposts? The lasing condition for such a microscopic system is given by Eqn. 7. Clearly, in order to reach the laser threshold, it is necessary to increase the Purcell factor $`F`$ and the quality factor $`Q`$. Our analysis indicates that large spontaneous-emission enhancement is possible in microposts. As an example, let us consider an unoptimized microcavity with the following parameters: $`n_h=3.57`$, $`n_l=2.94`$, $`D=0.5`$ $`\mu `$m, $`s=280`$ nm, $`a=155`$ nm, $`t=85`$ nm, $`MPT=15`$, $`MPB=30`$, $`Q=4800`$, $`\lambda =920`$ nm and $`V=2(\lambda /n_h)^3`$. (The method used for calculation is described in Refs.  and .) The Purcell factor for an emitter with zero linewidth positioned in the center of this micropost, is equal to 147. The enhancement drops to 65 for a linewidth of 100 GHz. Such a high Purcell factor would imply that $`\beta 1`$. Let us also assume that $`N_A1`$, corresponding to fast pumping. In order to observe single-dot lasing, we then need to satisfy the condition $`\tau _{ph}>\tau _{sp}`$. For a cavity with $`Q=10^4`$ operating at $`\lambda 1`$ $`\mu `$m, we have $`\tau _{ph}5.3`$ ps. Therefore, to achieve single-dot lasing, we would need $`\tau _{sp}`$ shorter than $`5`$ ps. If we again assume that the lifetime of an exciton without a cavity is $`0.5`$ ns, corresponding to a spontaneous emission rate of 2 GHz, lifetime reduction to $`5`$ ps would require a Purcell factor equal to 100. As mentioned above, such Purcell factors are possible for sufficiently narrow homogeneous linewidths. Single-dot lasing should therefore be possible in the optimized microposts.

## 6 Conclusions

Using the FDTD method, we have analyzed the fundamental ($`HE_{11}`$) mode in ideal, three-dimensional micropost cavities, for a variety of material systems (GaAs/AlAs, GaAs/AlxGa1-xAs and GaAs/AlOx). Microcavities were treated as single defects in a 1D cylindrical photonic crystal, which allowed us to push the limits of quality factors and mode volumes $`V`$ beyond those achievable by standard micropost designs. Our motivation was to maximize the $`Q/V`$ ratio of the defect mode, in order to use cavity-QED phenomena to build novel optoelectronic devices, such as single-dot lasers and high-efficiency light-emitting diodes, or to construct hardware for quantum computers and quantum communication systems, such as single-photon sources and strongly coupled quantum dot – cavity systems.

The standard approach for designing micropost microcavities is to choose the thicknesses of mirror layers and the spacer corresponding to the Bragg wavelength of a planar microcavity. We have shown that this approach does not necessarily lead to the highest $`Q`$ factors for the small cavity diameters analyzed in this article ($`D0.5`$ $`\mu `$m). Another widespread misconception is that the $`Q`$ of the cavity mode can always be improved by increasing the refractive-index contrast of the mirror layers. We have shown that this approach fails for small post diameters. Two primary loss mechanisms in three-dimensional microposts are the loss in the longitudinal direction, through DBR mirrors, and the loss in the transverse direction, due to imperfect confinement by TIR. A cavity mode which is strongly confined in the longitudinal direction by high-reflectivity mirrors is delocalized in Fourier space, leading to increased coupling to radiation modes and increased transverse loss. Similarly, a mode that is delocalized in the longitudinal direction and suffers significant longitudinal loss is more localized in Fourier space and suffers less transverse loss. When designing three-dimensional microposts, there is a tradeoff between these two loss mechanisms.

We were able to achieve $`Q`$ as high as $`10^4`$ together with mode volume as small as $`1.6(\lambda /n_h)^3`$ by optimizing structure parameters. Even though this range of values can be achieved in both GaAs/AlAs and GaAs/AlxGa1-xAs material systems, the former is a better choice from the perspective of fabrication, since the optimized structures require fewer mirror pairs on top.

We have also demonstrated that the optimized cavities can be used to observe novel cavity-QED phenomena, such as single-dot lasing or strong coupling between a single quantum dot and the cavity field. Moreover, the potential of microposts to achieve strong coupling with quantum dots is comparable to that of the largest $`Q/V`$ planar photonic crystal microcavities that are presently known .

Acknowledgement
This work was supported in part by the ARO, under the research grant 0160-G-BC575 (MURI program entitled “Single photon turnstile devices”). Financial assistance for M.P. was provided by the Stanford Graduate Fellowships program.
Present address of M. Pelton: Laboratory of Quantum Optics and Quantum Electronics, Department of Microelectronics and Information Technology, Royal Institute of Technology (KTH), Electrum 229, SE-164 40 Kista, Sweden.
Address if A. Scherer: California Institute of Technology, Pasadena, CA 91125.


text:
# Small-world properties of the Indian Railway network

## Abstract

Structural properties of the Indian Railway network is studied in the light of recent investigations of the scaling properties of different complex networks. Stations are considered as ‘nodes’ and an arbitrary pair of stations is said to be connected by a ‘link’ when at least one train stops at both stations. Rigorous analysis of the existing data shows that the Indian Railway network displays small-world properties. We define and estimate several other quantities associated with this network.

Given a chance, how would we have possibly organized our train travel? People dislike to change trains to reach their destinations. Therefore an extreme possibility would be to run a single train passing through all stations in the country so that no change of train is needed at all! An obvious disadvantage in this strategy is that the average distance between the stations become very large and so also the time needed for travel. The other limiting situation would be, to run a train between any pair of neighbouring stations and try to travel along the minimal paths. This requires a change of train at every station, which is also clearly not economically viable. Railway networks in no country in the world follow either of the two ways, actually they go mid-way. Like any other transport system the main motivation of railways is to be fast and economic. To achieve it, railways run simultaneously many trains, covering short as well as long routes so that a traveller does not need to change more than only a few trains to reach any arbitrary destination in the country.

In this paper we analyse the structure of the Indian Railway network (IRN). This is done in the context of recent investigations of the scaling properties of several complex networks e.g., social, biological, computational networks Small  etc. Identifying the stations as nodes of the network and a train which stops at any two stations as the link between the nodes we measure the average distance between an arbitrary pair of stations and find that it depends only logarithmically on the total number of stations in the country. While from the network point of view this implies the small-world nature of the railway network, in practice a traveller has to change only few trains to reach an arbitrary destination. This implies that over years, the railway network has been evolved with the sole aim in mind to make it fast and economic, eventually its structure has become a small-world network WS .

The structure and properties of several social, biological and computational networks like the World-wide web (WWW) web , network of the Internet structure Faloutsos , neural networks neural , collaboration network collab  etc. are being studied recently with much interest. In general a network has a number of ‘nodes’ and some ‘links’ connecting different pairs of nodes. Typically the following quantities are defined to characterize a network of $`N`$ nodes: (i) the diameter is the maximum distance between an arbitrary pair of nodes (ii) the clustering coefficient $`𝒞(N)`$ is the average fraction of connected triplets (iii) the probability distribution $`P(k)`$ that an arbitrarily selected node has the degree $`k`$ i.e. this node is linked to $`k`$ other nodes.

Watts and Strogatz WS  proposed a model of small-world network (SWN) in the context of various social and biological networks. They argued that SWNs must have small diameters which grow as $`\mathrm{ln}N`$ like random networks but should have large values of the clustering coefficients $`𝒞(N)1`$ like regular networks. On the other hand the scale-free networks (SFN) are characterized by the power law decay of the degree distribution function: $`P(k)k^\gamma `$. It was observed later that the degree distributions of nodes for two very important networks e.g., World-wide web web  which is a network of web-pages and the hyperlinks among various pages and the Internet network Faloutsos  of routers or autonomous systems have scale-free property. Barabási and Albert (BA) proposed a model for SFN which grows from an initial set of nodes and at every time step some additional nodes are introduced which are randomly connected to the previous nodes with the linear attachment probabilities barabasi . All scale-free networks are believed to display small-world properties while a small-world network is not necessarily scale-free.

Networks defined on the Euclidean space have also generated much interests in recent times. Internet, transport systems, postal networks etc. are naturally defined on two-dimensional space. In these generalised networks the attachment probabilities depend jointly on the nodal degrees as well as the lengths of the links euclid1 ; euclid2 .

A railway network is one of the most important examples of transport systems. The very complex topological structures of railway networks have attracted the attention of researchers in many different contexts. For example the fractal nature of the structure of railway networks was studied by Benguigui Benguigui . Very recently the efficiency of Boston subway network has been studied where a new measure for such networks has been proposed boston .

Our scheme is to associate first a representative graph $`G_N`$ with the IRN of $`N`$ stations in the following way. Here the stations represent the ‘nodes’ of the graph, whereas two arbitrary stations are considered to be connected by a ‘link’ when there is at least one train which stops at both the stations. These two stations are considered to be at unit distance of separation irrespective of the geographical distance between them. Therefore the shortest distance $`\mathrm{}_{ij}`$ between an arbitrary pair of stations $`s_i`$ and $`s_j`$ is the minimum number of different trains one needs to board on to travel from $`s_i`$ to $`s_j`$. Thus $`\mathrm{}_{ij}=1`$ implies that there is at least one train which stops at both $`s_i`$ and $`s_j`$. Similarly, $`\mathrm{}_{ij}=2`$ implies that there is no train which stops at both $`s_i`$ and $`s_j`$ and one has to change the train at least once in some intermediate station to board the second train to reach $`s_j`$. With this definition, if the trains $`t_1`$, $`t_2`$, $`\mathrm{}`$ $`t_n`$ pass through a station $`s_i`$, then all the stations through which these $`n`$ trains pass are unit distance away from $`s_i`$ and are considered as first neighbours of $`s_i`$. Consequently, the number $`k_i`$ of such stations is the degree of the node $`s_i`$.

Indian Railway network is a densely populated network of more than 8000 stations where the number of trains plying in this network is of the order of 10000 IRsite . However, we collected the data of IRN on a coarse-grained level following the recent Indian Railways time table ‘Trains at a Glance’ Timetable  containing the important trains and stations in India. This table contains a total of $`L=579`$ trains covering $`N=587`$ stations in a total of 86 tables. A grand rectangular matrix $`𝒢(N,L)`$ is then constructed such that the $`ij`$-th element of this matrix is 1 if the train $`j`$ stops at the station $`i`$, otherwise this element is zero. A second matrix $`𝒯(0:N,N)`$ is also constructed where the degree $`k_i`$ of the station $`i`$ is stored at the element $`𝒯(0,i)`$ and the serial numbers of the $`k_i`$ neighbours of $`i`$ are stored at the locations $`𝒯(j,i),j=1,k_i`$, rest of the elements being zero. We define and estimate the following quantities for the IRN.

Since $`G_N`$ is a connected graph, there are $`N(N1)/2`$ distinct shortest paths among the $`N`$ stations. We calculate the probability distribution of the shortest path lengths Prob$`(\mathrm{})`$. The shortest path lengths are calculated using a burning algorithm Herrmann  and using the matrix $`𝒯`$. In this algorithm the fire starts from an arbitrary node $`i`$, and burns this node at time $`t=0`$. At time $`t=1`$ the fire burns all $`k_i`$ neighbours of $`i`$. At time $`t=2`$ all unburnt neighbours of $`k_i`$ nodes are burnt and so on. The burning time of a node is the length of the shortest path of that node from the node $`i`$. This calculation has been repeated for all $`N`$ nodes to get $`N(N1)/2`$ shortest distances. In Fig. 1 inset we plot this distribution which goes to a maximum of $`\mathrm{}=5`$ implying that one needs to change at most four trains to reach any station from any station in India on the coarse-grained level. Similarly the distribution has a peak at $`\mathrm{}=2`$ implying that one can go to the majority of stations in India by changing train only once. In the graph theory the diameter of a graph is measured by the maximum distance between the pairs of nodes. Therefore according to this definition the diameter of our network is exactly equal to 5. However the average shortest path between an arbitrarily selected pair of nodes which we call as the mean distance $`𝒟(N)`$ is also a measure of the topological size of the graph and have been used by many authors to measure the size of networks as described in barabasi . We therefore measure the mean distance $`𝒟(N)`$ of the railway network of $`N`$ stations as the average shortest distance $`\mathrm{}_{ij}`$ between an arbitrary pair of stations $`s_i`$ and $`s_j`$. We obtain $`𝒟(N)2.16`$ for this network.

It is desirable to see how $`𝒟(N)`$ varies with $`N`$ amaral . Since we have the data of a single railway network, we divide the whole IRN into 25 different subsets consisting of trains and stations of 10 different states, 7 different combinations of states, 7 different railway zones and the whole IRN. As a result we obtained 25 data points (though they are not necessarily non-overlapping samples), reflecting the nature of variation of $`𝒟(N)`$ with $`N`$. In Fig. 1 we plot this data on a semi-log scale and though there is some wild fluctuations for small values of $`N`$, for large values of $`N`$ the linear behaviour is quite apparent. The whole range is fitted with $`𝒟(N)=A+B\mathrm{log}(N)`$ where $`A1.33`$ and $`B0.13`$.

The clustering coefficient $`𝒞(N)`$ is defined in the following way. Let the subgraph $`G_i`$ consisting of the neighbours of $`s_i`$ i.e., $`(s_1,s_2,s_3,\mathrm{},s_{k_i})`$ have $`E_i`$ links among them. Then the clustering co-efficient $`𝒞_i`$ of the node $`i`$ is $`2E_i/k_i(k_i1)`$ and that of the whole network is $`𝒞_i=𝒞_i`$. A direct measure of the clustering co-efficient of the whole IRN gives: $`𝒞0.69`$ (Fig. 2). The high value of the clustering coefficient is explained in the following way. The number $`n_s`$ of stations in which a particular train stops are all at unit distance from one another on the network and therefore form an $`n_s`$-clique. Therefore if only one train stops at some station $`i`$ then $`C_i=1`$. When two trains stop at the station $`i`$ and the sets $`n_s(1)`$ and $`n_s(2)`$ of stations covered by these two trains are different, $`C_i`$ is in general smaller than 1. However there may be other trains which do not stop at $`i`$ but stop at the stations which are not in both $`n_s(1)`$ and $`n_s(2)`$. These trains enhance the value of $`C_i`$. The value of $`𝒞0.69`$ is compared with a corresponding random graph network having the same number of vertices and edges as in IRN with the edges distributed randomly. It is found that the number of edges in IRN is 19603. If these edges are distributed randomly within the maximum possible edges on a graph of $`N`$=587 nodes the the clustering coefficient should be 19603/$`[N(N1)/2]0.113`$ which is the same as Prob(1). We also compute a modified clustering coefficient $`𝒞_o`$ by counting in $`E_i`$ only those links in the subgraph $`G_i`$ which pass through the node $`i`$. We obtained a value $`𝒞_o0.55`$ for the IRN.

Recently, the study of the clustering coefficient as a function of the degree of the node of some real-world network has shown an interesting feature Ravasz  $`𝒞(k)`$, defined as the clustering coeffcient of the node with degree $`k`$, shows a decrease (apparently a power law decay) with $`k`$ in several networks like the actor, language or world-wide-web networks. However in the network of internet at the router level or power grid network of the Western US, $`𝒞(k)`$ was found to be more or less a constant. In the IRN also, we find that $`𝒞(k)`$ (Fig. 3) remains at a constant value close to unity for small $`k`$ and shows a logarithmic decay at larger values of $`k`$. In all these real-world networks where $`𝒞(k)`$ remains more or less a constant, the nodes are linked by physical connections which may be responsible for this common feature. However, in this context it should also be mentioned that the scale-free Barabási-Albert network barabasi  also predicts $`𝒞(k)k^0`$ and $`𝒞(N)N^{0.75}`$. In the IRN, although $`𝒞(N)`$ shows a decrease with $`N`$, it is apparently much slower than a power law.

The degree distribution of the network, that is, the distribution of the number of stations $`k`$ which are connected by direct trains to an arbitrary station is denoted by $`P(k)`$. We plot the cumulative degree distribution $`F(k)=_k^{\mathrm{}}P(k)𝑑k`$ using a semi-log scale in Fig. 4 for the whole IRN. We see that $`F(k)`$ approximately fits to an exponentially decaying distribution $`F(k)exp(\alpha k)`$ with $`\alpha `$= 0.0085.

We also calculated the distribution $`D(n_t)`$ of the number of trains $`n_t`$ which stop at an arbitrary station. This is plotted in Fig. 5 on a semi-log scale after scaling by the average number of trains $`n_t12.06`$ along both the abscissa and the ordinate. The data is binned as before and is fitted to an exponential form: $`D_t(n_t)n_t=a\mathrm{exp}(bx)`$ with $`x=n_t/n_t`$, $`a0.47`$ and $`b0.75`$.

The distribution $`D(n_s)`$ of the number of stations through which an arbitrary train passes is plotted in Fig. 6. The data is scaled by the average number of stations $`n_s12.37`$ along both the abscissa and the ordinate. The $`D(n_s)`$ grows very fast at the beginning, reaches a maximum and then decays to zero. A numerical fit to a functional form like $`D_s(n_s)n_s=ax^4/(x^2+b)^3`$ with $`x=n_s/n_s`$, $`a0.6`$ and $`b0.096`$ turns out to be reasonably good.

We also measure the connectivity correlation of IRN following the works of Pastor . Let $`F(k^{}|k)`$ denote the conditional probability that a node of degree $`k`$ has a neighbour of degree $`k^{}`$. Then to see how the nodes of different degrees are correlated we measure the average degree $`k_{nn}(k)=\mathrm{\Sigma }_k^{}k^{}F(k^{}|k)`$ of the subset of nodes which are all neighbours to a particular node of degree $`k`$. In general this average has a variation like $`k_{nn}(k)k^\nu `$ where a non-zero $`\nu `$ reflects a non-trivial correlation among the nodes of the network. We calculated $`k_{nn}(k)`$ for IRN and plotted it in Fig. 7 on a double logarithmic scale. Almost over a decade the $`k_{nn}(k)`$ remains same on the average and is independent of $`k`$, indicating the absence of correlations among the nodes of different degrees.

A more sensitive measure for the degree correlations was proposed in Newman . Newman has defined a degree-degree correlation function $`r`$ which measures whether a vertex of high degree at one end of a link prefers a vertex of high degree (“assortative mixing”, $`r>0`$) or low degree (“disassortative mixing” $`r<0`$) at the other end. It has been observed that social networks are assortative and technological and biological networks are disassortative. We have measured for IRN the normalized correlation function following Newman  and found its values to be $`r`$ = -0.033. This indicates that the IRN is of disassortative nature, i.e. rich vertices at one end of a link show some preference towards poor vertices at the other end, and vice versa.

To summarize, we investigated the structural properties of the Indian Railway network to see if some of the general scaling behaviour obtained for many complex networks in recent times may also be present in IRN. While nodes of the network are evidently the stations, the links are defined as the pairs of stations communicated by single trains. With such a definition of link, the mean distance of the network is a measure of how good is the connectivity of the network. Indeed, we observed that the mean distance of IRN varies logarithmically with the number of nodes with a high value of the clustering coefficient. This implies that IRN behaves like a small-world network, which we believe should be typical of the railway network of any other country, which we are unable to study at present for unavailability of data.

We like to thank I. Bose for constantly encouraging us to work on this problem and also to S. Goswami for suggesting ref. . PS acknowledges financial support from DST grant SP/S2-M11/99. GM acknowledges hospitality in the S. N. Bose National Centre.