text:
# An interpretation of saturation phenomena as Glauber-Gribov multiple parton scatterings

## 1 pQCD, Glauber-Gribov parton rescatterings <br>and the colour dipole

The model of Ref. , hereafter labeled “pQCD+Glauber”, assumes minijet production in $`pA`$ collisions to be dominated by semihard parton multiple scatterings. By “semihard scattering” we mean a process with a minimum transverse momentum exchange $`p_012`$ GeV, described by leading order (LO) pQCD parton-parton cross section. The model assumes the S-matrix for a collision of one parton on $`n`$ partons from the target to be factorizable in terms of S-matrices for parton-parton elastic-scattering, and assumes generalized pQCD factorization. Considering only gluons for simplicity, the minijet transverse spectrum is then given by

$`{\displaystyle \frac{d\sigma _{pA}^{mj}}{d\eta d^2p_T}}=G(x,p_T^2){\displaystyle \frac{d\sigma _{hard}^A}{d^2p_T}}+AG(x^{},p_T^2){\displaystyle \frac{d\sigma _{hard}^p}{d^2p_T}},`$ (1)
where $`G(x,Q^2)`$ is the proton distribution function of a gluon with fractional momentum $`x`$, and $`x(x^{})(p_T/\sqrt{s})\mathrm{exp}(\pm \eta )`$. In the second term of Eq. (1), the A-nucleus partons are assumed to undergo a single scattering on the proton with cross section

$`{\displaystyle \frac{d\sigma _{hard}^p}{d^2p_T}}={\displaystyle _{x_{min}^{}}^1}𝑑x^{}G(x^{},p_T^2){\displaystyle \frac{d\widehat{\sigma }}{d^2p_T}}.`$ (2)
The LO pQCD cross section for gluon-gluon scattering is $`\frac{d\widehat{\sigma }}{d^2p_T}\frac{9}{2}\alpha _s^2\frac{1}{(p_T^2+p_0^2)^2}`$ and the limits of integration on $`x_i^{}`$ are given by parton-parton kinematic constraints. The first term of Eq. (1) accounts for multiple semihard scatterings of the proton partons on the nucleus. The proton is considered pointlike at an impact parameter $`b`$. Nuclear effects are assumed to be due only to multiple semihard scatterings, and are included in $`d\sigma _A^{hard}/d^2p_T`$, the transverse momentum distribution of a proton parton who suffered at least one semihard scattering. This is written as :

$`{\displaystyle \frac{d\sigma _{hard}^A}{d^2p_t}}={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}{\displaystyle d^2b\frac{d\sigma _{hard}^p}{d^2k_1}T_A(b)\times \mathrm{}\times \frac{d\sigma _{hard}^p}{d^2k_n}}`$ $`T_A(b)e^{\sigma _{hard}^p(p_0)T_A(b)}`$
$`\times \delta ^{(2)}\left({\displaystyle 𝐤_i}𝐩_𝐭\right),`$ (3)
where $`\sigma _p^{hard}(p_0)=d^2k\frac{d\sigma _{hard}^p}{d^2k}`$ is the integrated gluon-nucleon cross section, which depends explicitly on the infrared regulator $`p_0`$. $`T_A(b)`$ is the target nucleus thickness function. The exponential factor in Eq. (3) represents the probability that the parton suffered no semihard scatterings after the $`n`$-th one. In such a way, unitarity is explicitly implemented at the nuclear level, as discussed in Ref. . The sum over $`n`$ may be explicitly performed in Fourier space. The result reads:

$`{\displaystyle \frac{d\sigma _{hard}^A}{d^2p_T}}={\displaystyle \frac{d^2r_T}{4\pi ^2}e^{i\stackrel{}{k}_T\stackrel{}{r}_T}S_{hard}^A(r_T;p_0)},`$ (4)
where

$`S_{hard}^A(r_T;p_0)={\displaystyle d^2b\left[e^{\stackrel{~}{\sigma }_{hard}^p(r_T,p_0)T_A(b)}e^{\sigma _{hard}^p(p_0)T_A(b)}\right]}`$ (5)
and

$`\stackrel{~}{\sigma }_{hard}^p(r_T;p_0)={\displaystyle d^2k\left[1e^{i\stackrel{}{k}\stackrel{}{r}_T}\right]\frac{d\sigma _{hard}^p}{d^2k}}.`$ (6)
Note that $`\sigma _{hard}^p(r_T)r_T^2`$ as $`r_T0`$ and $`\sigma _{hard}^p(r_T)\sigma _{hard}^p`$ as $`r_T\mathrm{}`$. This suggests the interpretation of $`\sigma _{hard}^p(r)`$ as a dipole-nucleon “hard” cross section. This dipole is of mathematical origin, and comes from the square of the scattering amplitude written in the Fourier variable $`r_T`$, which represents the transverse size of the dipole. Then, we can interpret $`S_{hard}^A`$ as the dipole-nucleus “hard” cross section. Eq. (5) clearly incorporates Glauber-Gribov multiple scatterings of the colour dipole. Note that no nuclear effects on PDF’s are included, but shadowing is partly taken into account by the dipole multiple scatterings.

## 2 CGC and the colour dipole

The Colour Glass Condensate is an effective theory for the nucleus gluon field at small-$`x`$, which describes the high-density regime where gluon saturation comes into play to modify the free-nucleon PDF’s .

In this framework, minijet production in $`pA`$ collisions may be related to the cross section for the scattering of a colour dipole on a nucleus . The basic assumption is that the density of projectile partons is low enough for parton correlations inside the proton to be due only to DGLAP evoultion. In this case one is allowed to treat the proton in the collinear factorization limit, exactly as in the pQCD+Glauber model of the previus section, see Eqs. (1) and (4) compared to Eqs. (36) and (32) of Ref. .

However, the dipole-nucleus interaction is computed entirely in the CGC model for the target nucleus. By using the so-called “Gaussian approximation” , in which parton correlation are assumed only to be Gaussian, one may write the dipole-nucleon cross section appearing in Eq. (4) as follows:

$`\stackrel{~}{\sigma }_{CGC}^p(r_T)=4\pi \alpha _sN_c{\displaystyle \frac{d^2k_T}{(2\pi )^2}\frac{\mu _\tau (k_T)}{k_T^4}\left[1e^{ik_Tr_T}\right]},`$ (7)
where $`N_c`$ is the number of colours, and $`\tau =\mathrm{ln}(1/x)`$ is called “rapidity”. Next,, $`\mu _\tau `$ is interpreted as the unintegrated gluon distribution and has two limits :

$`\mu _\tau =\{\begin{array}{cc}{\displaystyle \frac{4\pi ^2}{N_c^21}}{\displaystyle \frac{dxG(x,k_T^2)}{d\mathrm{ln}(k_T^2)}}\hfill & k_T^2Q_s^2\hfill \\ {\displaystyle \frac{\delta }{\alpha _s}}k_T^2\mathrm{ln}{\displaystyle \frac{Q_s^2}{k_T^2}}\hfill & k_T^2Q_s^2\hfill \end{array}`$ (10)
These high- and low-$`k_T`$ limits are defined relatively to the “saturation momentum” $`Q_s^2(\tau )=Q_0^2e^{\mathrm{\hspace{0.17em}4.84}\overline{\alpha }_s(\tau \tau _0)}`$, where $`\tau _0`$ is the minimum rapidity at which gluon saturation occurs and $`Q_0`$ is the corresponding saturation scale. Unfortunately, the constants $`\delta `$, $`\tau _0`$ and $`Q_0`$ cannot be computed in this approximation.

## 3 CGC is pQCD+Glauber

The similarity of the CGC model at high-$`k_T`$ with the pQCD+Glauber model becomes evident by comparing Eqs. (7) and (10) with Eqs. (5), (6) and (2). To make this relationship more precise we proceed in four steps.

Step 1. In the CGC, we separate “soft” and “hard” interactions:

$`\mu _\tau =\mu _\tau ^S+\mu _\tau ^H\mathrm{with}\{\begin{array}{c}\mu _\tau ^H={\displaystyle \frac{4\pi ^2}{N_c^21}}{\displaystyle \frac{dxG(x,k_T^2)}{d\mathrm{ln}(k_T^2)}}\times {\displaystyle \frac{k_T^4}{(k_T^2+Q_s^2)^2}}\hfill \\ \mu _\tau ^S=\mu _\tau \mu _\tau ^H\hfill \end{array}`$ (13)
so that $`\mu _\tau ^H`$ vanishes quickly for $`k_T<Q_s`$. Accordingly, we have a soft and hard dipole-nucleon cross section: $`\stackrel{~}{\sigma }_{CGC}^p=\stackrel{~}{\sigma }_{\mathrm{eff}}^S+\stackrel{~}{\sigma }_{\mathrm{eff}}^H.`$

Step 2. In the dipole-nucleus cross section, we isolate contributions from processes with purely soft rescatterings, and processes with at leat one hard scattering. These define, respectively, the soft and hard dipole-nucleus cross section:

$`S_{CGC}^A=`$ $`S_{\mathrm{eff}}^S+S_{\mathrm{eff}}^H={\displaystyle d^2be^{\stackrel{~}{\sigma }_{\mathrm{eff}}^S(r_T)T_A(b)}e^{\sigma _{\mathrm{eff}}^HT_A(b)}}`$
$`+{\displaystyle d^2be^{\stackrel{~}{\sigma }_{\mathrm{eff}}^S(r_T)T_A(b)}\left[e^{\stackrel{~}{\sigma }_{\mathrm{eff}}^H(r_T)T_A(b)}e^{\sigma _{\mathrm{eff}}^HT_A(b)}\right]}.`$
The first term includes any number of soft scatterings but no hard ones, and the second term processes with at least one hard scattering.

Step 3. Study observables for which it is possible to neglect the soft part, e.g., hadron spectra at large $`p_T`$ (Cronin effect), or integrated charged multiplicities at large $`\sqrt{s}`$, where contribution from hard processes should become dominant.

Step 4. Following Ref. , we use the DGLAP equation to approximate the unintegrated gluon distribution in Eq. (13) with the integrated one:

$$\frac{dxG(x,k_T^2)}{d\mathrm{ln}k_T^2}\frac{\alpha _sN_c}{\pi }_x^1𝑑x^{}G(x^{},k_T^2)$$

This approximation is valid at large $`Q^2Q_s^2`$.

Result. Collecting all the steps, we obtain

$`S_{CGC}^A`$ $`{\displaystyle d^2b\left[e^{\stackrel{~}{\sigma }_{\mathrm{eff}}^H(r_T)T_A(b)}e^{\sigma _{\mathrm{eff}}^HT_A(b)}\right]}`$ (14)
$`\stackrel{~}{\sigma }_{\mathrm{eff}}^H(r_T)`$ $`={\displaystyle d^2k_T_x^1𝑑x^{}G(x^{},k_T^2)\frac{(9/2)\alpha _s^2}{(k_T^2+Q_s^2)^2}\left[1e^{ik_Tr_T}\right]}.`$ (15)
Eqs. (14) and (15) are equal to the pQCD+Glauber Eqs. (5) and (6) with $`p_0=Q_s`$ and $`x=x_{min}^{}`$.

## 4 What do RHIC data have to say?

Under the assumptions stated in the previous section, it is possible to quantitatively compute some observables of interest. While a detailed applicaton of Eq. (4) to hadron $`p_T`$-spectra and the Cronin effect is under investigation , we may use it to study RHIC Au-Au data on the centrality dependence of charged particle pseudorapidity densities . The model was derived in the case of $`pA`$ collisions, but is generalizable to AB collisions under the assumption that in the regime of interest both nuclei are composed of a diluted enough system of partons.

As a first step, we need the average minijet multiplicity at fixed impact parameter and pseudorapidity. This is obtained by integrating Eq. (3) over $`p_T`$, and introducing the thickness function $`T_B`$ of projectile nucleus:

$`{\displaystyle \frac{dN^{mj}}{d\eta }}(b;p_0)={\displaystyle d^2\beta G(x,p_0^2)T_B(b\beta )\left[1e^{K\sigma _{hard}^p(p_0)T_A(\beta )}\right]}+AB,`$ (16)
where the K-factor that simulates higher-order corrections to the pQCD cross section is explicitly shown.

The integrand in Eq. (16) may be interpreted as the average density of projectile partons (at a given $`x`$) times the probability of having at least one semi-hard scattering against the target. At low values of $`p_0`$ the semihard cross section is large and the target becomes mor and more black to the projectile partons: the probability of scattering at least once becomes so high that nearly every projectile parton scatters and is extracted from the incoming nuclear wave-function. In this regime even if we use a lower cutoff $`p_0`$ no more partons are there to be extracted. For this reason the minijet multiplicity tends to saturate . We call saturation cutoff $`p_{sat}`$ the largest value of $`p_0`$ at which this happens, and evaluate minijet multiplicities from Eq. (16) with $`p_0=p_{sat}`$. We refer to  for more details.

To apply this minijet-level computation to integrated charged multiplicities, we write the charged particle multiplicity per unit rapidity as the sum of a soft and a semi-hard part . For the soft part we use the wounded-nucleon model, which postulates a scaling of $`N^{ch}`$ with the number of participants. The semi-hard part is assumed to be completely computable from the saturation criterion for minijet production described above. To convert the minijet multiplicity to charged particle multiplicity, we further assume isentropic expansion of the initially produced minijet plasma, and parton-hadron duality. Our final formula reads

$`{\displaystyle \frac{1}{N_{part}(b)/2}}{\displaystyle \frac{dN^{ch}}{d\eta }}(b)=n_{soft}+{\displaystyle \frac{1}{N_{part}(b)/2}}{\displaystyle \frac{5}{3}}{\displaystyle \frac{dN^{mj}}{d\eta }}(b;p_{sat}).`$ (17)
As the soft component $`n_{soft}=n_{soft}(\sqrt{s})`$ in Eq. (17) does not depend on the centrality of the collision, we may fit it to data for central collision only (we used data at $`\sqrt{s}`$ = 56-200 GeV  and extrapolated the fit down to 19 GeV, see Ref. ). The behaviour of the observable in non-central collisions is then completely determined by the model. In Fig. 1 we show the results of Eq. (17) compared to phobos data. Note that in our computation we used standard DGLAP evolved PDF’s from the GRV group , with no nuclear modification of any sort.

## 5 Conclusions

The pQCD+Glauber formulation gives a very well defined framework in which detailed numerical calculations may be performed, and which is well-tested in the case of nucleon-nucleon collisions. It is equivalent to the CGC model in a kinematic region where parton densities are low enough for gluon saturation effects to be negligible, but high enough for multiple scatterings and unitarity effects to become important. RHIC data on charged multiplicities at central rapidity are very satisfactorily described by the pQCD+Glauber model with free-nucleon PDF’s. Therefore, at least at central rapidity, RHIC data do not show any sign of gluon saturation.

However, integrated multiplicities have very limited sensitivity to the details of the production processes, and alone do not allow any definitive conclusion. To assess the presence or absence of nuclear effects beyond multiscatterings and unitarity, it is necessary to study more differential observables in a cleaner environment than the hot and dense medium produced in Au-Au collisions. Novel effects may then be observed as a deviation of experimental data from the baseline given by pQCD+Glauber. The Cronin effect on hadron production in $`dA`$ collisions at RHIC is an ideal candidate.


format: text

text:
# References

EFI 03-14

hep-ph/0304178

April 2003

I-SPIN, U-SPIN, AND PENGUIN DOMINANCE IN $`BKK\overline{K}`$ 111To be submitted to Physics Letters B.

Michael Gronau222Permanent Address: Physics Department, Technion – Israel Institute of Technology, 32000 Haifa, Israel. and Jonathan L. Rosner

Enrico Fermi Institute and Department of Physics

University of Chicago, Chicago, Illinois 60637

> Isospin and U-spin symmetries are useful approximations for studying penguin dominance in $`B`$ meson decays to three kaons, $`BKK\overline{K}`$. We point out certain subtleties in treating these decays in these approximations. Resulting uncertainties are discussed in determining the CP content of the final state in $`BK^+K^{}K_S`$, and in relating the CP asymmetry in this process to the value of $`\mathrm{sin}2\beta `$.

PACS codes: 12.15.Hh, 12.15.Ji, 13.25.Hw, 14.40.Nd

Recently the Belle collaboration reported branching ratio measurements for $`B`$ meson decays to three kaon states, $`K^+K^{}K^+,K^+K^{}K_S,K_SK_SK^+`$ and $`K_SK_SK_S`$ , and for $`B^+`$ decays to $`K^+K^{}\pi ^+`$ and $`K^+\pi ^+\pi ^{}`$ . The BaBar collaboration measured branching ratios for $`B^+`$ decays to all three body final states involving charged kaons or pions, $`K^+K^{}K^+,K^+K^{}\pi ^+,K^+\pi ^+\pi ^{}`$ and $`\pi ^+\pi ^+\pi ^{}`$ . The averaged values of the measured branching ratios are, in units of $`10^6`$:

$`(B^+K^+K^+K^{})`$ $`=`$ $`30.8\pm 2.1,(B^0K^+K^0K^{})=29.3\pm 5.3,`$ (1)
$`(B^+K^+K_SK_S)`$ $`=`$ $`13.4\pm 2.4,(B^0K_SK_SK_S)=4.3\pm 1.7,`$ (2)
$`(B^+K^+\pi ^+\pi ^{})`$ $`=`$ $`58.4\pm 4.4,(B^+\pi ^+\pi ^+\pi ^{})=10.9\pm 3.7,`$ (3)
$`(B^+K^+K^{}\pi ^+)`$ $`<`$ $`6.3(90\%\mathrm{c}.\mathrm{l}.).`$ (4)
These are averaged over the process shown and its CP-conjugate.

Assuming penguin dominance in $`BKK\overline{K}`$, an isospin analysis was attempted by the Belle collaboration  in order to isolate the CP-even and CP-odd components of the $`K^+K^{}K_S`$ final state. This information is useful for studying the time-dependent CP asymmetry in this channel . Similar isospin arguments were presented subsequently in , where a U-spin study relating three body $`B^+`$ decays involving charged kaons and pions was employed in order to estimate deviations from penguin dominance.

In the present Letter we will iterate the isospin analysis for $`BKK\overline{K}`$, pointing out a subtlety which was overlooked by the above two studies, thereby oversimplifying the analysis. It will be shown that these earlier studies made an implicit assumption which goes beyond isospin symmetry. We will argue that, nevertheless, in the pure penguin limit an equality holds between the amplitudes of $`B^0K^+K^{}K^0`$ and $`B^+K^0\overline{K}^0K^+`$, which is the basis for the CP argument. In order to study deviations from penguin dominance in these decays we will employ U-spin considerations which were oversimplified in . We will argue that these deviations, which are partly due to electroweak penguin contributions, may be larger than estimated, and introduce a sizable uncertainty in the value of $`\mathrm{sin}2\beta `$ determined from the CP asymmetry in $`B^0K^+K^{}K_S`$.

The effective Hamiltonian describing charmless decays $`BKK\overline{K}`$ consists of operators transforming as a sum of $`\mathrm{\Delta }I=0`$ and $`\mathrm{\Delta }I=1`$. The initial state is pure $`|I=\frac{1}{2}`$, while the final state $`|f|KK\overline{K}`$ is a superposition of three isospin states, $`|I(KK)=0,I_f=\frac{1}{2},|I(KK)=1,I_f=\frac{1}{2}`$ and $`|I(KK)=1,I_f=\frac{3}{2}`$. Consequently, these decay processes are described by five independent isospin amplitudes  corresponding to given $`KK\overline{K}`$ momenta. Suppressing the momentum dependence, we denote these amplitudes by $`A_{\mathrm{\Delta }I}^{I(KK),I_f}I(KK),I_f|\mathrm{\Delta }I|\frac{1}{2}`$, and list them as $`A_0^{0,\frac{1}{2}}`$$`A_0^{1,\frac{1}{2}}`$, $`A_1^{0,\frac{1}{2}}`$$`A_1^{1,\frac{1}{2}}`$ and $`A_1^{1,\frac{3}{2}}`$.

$`B`$ mesons decay into two kaons and an antikaon in four distinct flavor modes, $`B^+K^+K^+K^{}`$, $`B^0K^0K^0\overline{K}^0`$, $`B^+K^+K^0\overline{K}^0`$ and $`B^0K^+K^0K^{}`$. The first two processes involve identical kaons ($`K^+`$ and $`K^0`$, respectively), while in the other two processes the kaons ($`K^+`$ and $`K^0`$ in both decays) have different charges. The amplitudes for $`BKK\overline{K}`$ depend on the kaons’ momenta. In $`B^0K^+K^0K^{}`$ and in $`B^+K^+K^0\overline{K}^0`$ one expects different values for amplitudes when exchanging the $`K^+`$ and $`K^0`$ momenta. Using isospin, one can decompose decay amplitudes for $`BKK\overline{K}`$ into isospin amplitudes $`A_{\mathrm{\Delta }I}^{I(KK),I_f}`$ describing two kaons and an antikaon with given momenta. Thus, one obtains expressions for six decay amplitudes in terms of five isospin amplitudes,

$`A(K^+K^+K^{})_{p_1p_2p_3}`$ $`=`$ $`2A_0^{1,\frac{1}{2}}2A_1^{1,\frac{1}{2}}+A_1^{1,\frac{3}{2}},`$ (5)
$`A(K^0K^0\overline{K}^0)_{p_1p_2p_3}`$ $`=`$ $`2A_0^{1,\frac{1}{2}}2A_1^{1,\frac{1}{2}}+A_1^{1,\frac{3}{2}},`$ (6)
$`A(K^+K^0\overline{K}^0)_{p_1p_2p_3}`$ $`=`$ $`A_0^{0,\frac{1}{2}}A_0^{1,\frac{1}{2}}A_1^{0,\frac{1}{2}}+A_1^{1,\frac{1}{2}}+A_1^{1,\frac{3}{2}},`$ (7)
$`A(K^+K^0\overline{K}^0)_{p_2p_1p_3}`$ $`=`$ $`A_0^{0,\frac{1}{2}}A_0^{1,\frac{1}{2}}+A_1^{0,\frac{1}{2}}+A_1^{1,\frac{1}{2}}+A_1^{1,\frac{3}{2}},`$ (8)
$`A(K^+K^0K^{})_{p_1p_2p_3}`$ $`=`$ $`A_0^{0,\frac{1}{2}}+A_0^{1,\frac{1}{2}}+A_1^{0,\frac{1}{2}}+A_1^{1,\frac{1}{2}}+A_1^{1,\frac{3}{2}},`$ (9)
$`A(K^+K^0K^{})_{p_2p_1p_3}`$ $`=`$ $`A_0^{0,\frac{1}{2}}+A_0^{1,\frac{1}{2}}A_1^{0,\frac{1}{2}}+A_1^{1,\frac{1}{2}}+A_1^{1,\frac{3}{2}}.`$ (10)
On the left-hand side amplitudes are specified by the three outgoing particles and by their respective momenta. On the right-hand side we have absorbed Clebsch-Gordan coefficients in the definition of isospin amplitudes and have suppressed the momentum dependence of these amplitudes.

Let us comment briefly on Eqs. (5)–(10). The equal magnitudes and the relative signs of contributions of isospin amplitudes in pairs of processes can be easily understood from simple considerations. In Eqs. (5) and (6) the final states involve $`I(KK)=1`$ and are related to each other by an isospin reflection $`ud,\overline{d}\overline{u}`$. Consequently, the magnitudes of $`\mathrm{\Delta }I=0`$ and $`\mathrm{\Delta }I=1`$ contributions in these two processes are equal and occur with opposite and equal signs, respectively. In Eq. (7) and (8) one interchanges the $`K^+`$ and $`K^0`$ momenta in $`B^+K^+K^0\overline{K}^0`$. It then follows from Bose statistics that contributions to the two amplitudes from $`I(KK)=0`$ terms, which are antisymmetric in the isospins of the two kaons, and $`I(KK)=1`$ terms, which are symmetric in isospin, are equal and have opposite and equal signs, respectively. This corresponds to situations in which the orbital angular momentum of the $`KK`$ system in its center of mass frame (which equals the $`\overline{K}`$ angular momentum relative to this center of mass) is odd and even, respectively. The same argument applies to Eqs. (9) and (10), giving the amplitude for $`B^0K^+K^0K^{}`$ in two points of phase space where the two kaon momenta are interchanged.

An interesting consequence of the isospin decomposition is a sum rule between $`B^+`$ and $`B^0`$ decay amplitudes. The six amplitudes (5)–(10) obey one linear relation between the sum of three amplitudes for a charged $`B`$ meson and the sum of three amplitudes for a neutral $`B`$,

$`A(K^+K^+K^{})_{p_1p_2p_3}+A(K^+K^0\overline{K}^0)_{p_1p_2p_3}+A(K^+K^0\overline{K}^0)_{p_2p_1p_3}`$ $`=`$
$`A(K^0K^0\overline{K}^0)_{p_1p_2p_3}+A(K^+K^0K^{})_{p_1p_2p_3}+A(K^+K^0K^{})_{p_2p_1p_3}`$ $`=`$ $`3A_1^{1,\frac{3}{2}}.`$ (11)
The two sums, in each of which one sums over amplitudes at two points of phase space where $`K^+`$ and $`K^0`$ momenta are interchanged, are given by the $`I_f=\frac{3}{2}`$ amplitude. This relation is similar to an isospin relation among the four amplitudes for $`B^+`$ and $`B^0`$ decays to $`K\pi `$ .

So far our arguments were based purely on isospin symmetry. Let us now study the consequences of penguin dominance in $`BKK\overline{K}`$ decays, assuming that the dominant term in the $`\mathrm{\Delta }C=0,\mathrm{\Delta }S=1`$ effective Hamiltonian  contributing to these decays is a $`\overline{b}\overline{s}`$ QCD penguin operator. The assumption of penguin dominance implies that one keeps only $`\mathrm{\Delta }I=0`$ terms in Eqs. (5)–(10), neglecting $`A_1^{0,\frac{1}{2}}`$ $`A_1^{1,\frac{1}{2}}`$ and $`A_1^{1,\frac{3}{2}}`$. Note that this excludes electroweak penguin operators which contain a term transforming as $`\mathrm{\Delta }I=1`$. We will return to this point when discussing deviations from penguin dominance. In the latter approximation one has

$`A(K^+K^+K^{})_{p_1p_2p_3}`$ $`=`$ $`A(K^0K^0\overline{K}^0)_{p_1p_2p_3}=2A_0^{1,\frac{1}{2}},`$ (12)
$`A(K^+K^0K^{})_{p_1p_2p_3}`$ $`=`$ $`A(K^0K^+\overline{K}^0)_{p_1p_2p_3}=A_0^{0,\frac{1}{2}}+A_0^{1,\frac{1}{2}},`$ (13)
where $`A_0^{0,\frac{1}{2}}`$ and $`A_0^{1,\frac{1}{2}}`$ are antisymmetric and symmetric under interchange of the two kaon momenta. The equality of the two pairs of $`B^+`$ and $`B^0`$ decay amplitudes follows simply from an isospin reflection $`ud`$ in initial and final states. In this limit the two amplitudes involving $`K^+K^{}`$ in the final state are, however, different.

When squaring the amplitudes and integrating over phase space one includes a factor $`\frac{1}{2}`$ for identical particles in the first pair of processes. The interference between the two isospin amplitudes $`A_0^{0,\frac{1}{2}}`$ and $`A_0^{1,\frac{1}{2}}`$ in the second pair of processes, corresponding to even and odd angular momenta of the $`KK`$ system, vanishes. Thus, one finds

$`\mathrm{\Gamma }(B^+K^+K^+K^{})`$ $`=`$ $`\mathrm{\Gamma }(B^0K^0K^0\overline{K}^0)=2\mathrm{\Gamma }_0^{1,\frac{1}{2}},`$ (14)
$`\mathrm{\Gamma }(B^0K^+K^0K^{})`$ $`=`$ $`\mathrm{\Gamma }(B^+K^+K^0\overline{K}^0)=\mathrm{\Gamma }_0^{0,\frac{1}{2}}+\mathrm{\Gamma }_0^{1,\frac{1}{2}},`$ (15)
where $`\mathrm{\Gamma }_0^{0,\frac{1}{2}}`$ and $`\mathrm{\Gamma }_0^{1,\frac{1}{2}}`$ are rates corresponding to the two isospin amplitudes. We conclude that, while the two rate equalities (14) and (15) follow from penguin dominance, not all four rates are equal in this approximation. In particular, the two rates involving $`K^+K^{}`$ in the final state may be different in general, contrary to arguments made in . They become equal when $`\mathrm{\Gamma }_0^{0,\frac{1}{2}}=\mathrm{\Gamma }_0^{1,\frac{1}{2}}`$ which goes beyond isospin symmetry. Experimentally, one has $`\mathrm{\Gamma }(B^+K^+K^+K^{})/\mathrm{\Gamma }(B^0K^+K^0K^{})=0.98\pm 0.19`$ which implies $`\mathrm{\Gamma }_0^{0,\frac{1}{2}}/\mathrm{\Gamma }_0^{1,\frac{1}{2}}=1.05_{0.33}^{+0.49}`$.

Eq. (13) provides the basis for attempting a separation between the CP-even and CP-odd components in the final state of the measured process $`B^0K^+K^{}K_S`$ . These components correspond to even and odd angular momentum $`K^+K^{}`$ states. Equal amplitudes in $`B^0K^+K^0K^{}`$ and $`B^+K^+K^0\overline{K}^0`$ imply that the final states in the two processes have the same angular momentum decomposition in terms of $`K^+K^{}`$ in one process and $`K^0\overline{K}^0`$ in the other. The probability for a $`K^0\overline{K}^0`$ being in an even angular momentum state, where it decays as $`K_SK_S+K_LK_L`$, is given by $`2\mathrm{\Gamma }(B^+K^+K_SK_S)/\mathrm{\Gamma }(B^+K^+K^0\overline{K}^0)`$. Using Eq. (15), this probability is given by a ratio of two measured rates $`2\mathrm{\Gamma }(B^+K^+K_SK_S)/\mathrm{\Gamma }(B^0K^+K^0K^{})=1.04\pm 0.20`$, excluding the $`\varphi K_S`$ contribution in the denominator . This is also the probability for a CP-even state in $`B^0K^+K^{}K_S`$ excluding $`\varphi K_S`$.

The above conclusion, indicating that the final state in $`B^0K^+K^{}K_S`$ (excluding $`\varphi K_S`$) is dominantly CP-even, is based on assuming penguin dominance in $`BKK\overline{K}`$. The rest of the discussion will address this issue. We will study tests for penguin dominance in $`BKK\overline{K}`$, and will evaluate deviations from this approximation in terms of measurable rates.

Penguin dominance in $`BK\pi `$, suggested in , was first tested in  by using flavor SU(3) and comparing decay rates for $`BK\pi `$ and $`B\pi \pi `$. The measured rates were also used to estimate the deviation from pure penguin dominance, given by a parameter $`|T^{}/P^{}|0.2`$. A similar analysis will be presented here in order to relate $`BKK\overline{K}`$ to $`B\pi \pi \pi `$ and $`BK\pi \pi `$ to $`BKK\pi `$. This will test penguin dominance in $`BKK\overline{K}`$. Our arguments differ in detail from those presented in . In  factorization was assumed for three body decays for which no good theoretical justification exists, while in  a subtlety in using U-spin was overlooked. We will also explain the special role of electroweak penguin contributions which were ignored in the latter study.

A useful subgroup of SU(3) permitting relations between $`B^+`$ meson decays to final states involving charged pions and kaons is U-spin , under which the pairs $`(d,s),(\overline{s},\overline{d}),(\pi ^{},K^{})`$ and $`(K^+,\pi ^+)`$ transform like doublets. The $`\mathrm{\Delta }C=0,\mathrm{\Delta }S=1`$ effective Hamiltonian transforms like a $`\overline{s}`$ component ($`\mathrm{\Delta }U_3=\frac{1}{2}`$) of a U-spin doublet, while the $`\mathrm{\Delta }C=0,\mathrm{\Delta }S=0`$ Hamiltonian transforms like a $`\overline{d}`$ component ($`\mathrm{\Delta }U_3=\frac{1}{2}`$) of another U-spin doublet. Let us consider the decays of $`B^+`$ into the final states $`K^+K^+K^{},K^+\pi ^+K^{},K^+\pi ^+\pi ^{}`$ and $`\pi ^+\pi ^+\pi ^{}`$. The initial state in these decays is pure $`U=0`$. The final states from the $`\mathrm{\Delta }U=\frac{1}{2}`$ transitions are two U-spin doublets, in which the two positively charged particles are in $`U(++)=0`$ and $`U(++)=1`$ states. Therefore, $`\mathrm{\Delta }S=1`$ and $`\mathrm{\Delta }S=0`$ amplitudes may be written separately in terms of two U-spin amplitudes corresponding to these two states . We will denote $`\mathrm{\Delta }S=1`$ amplitudes corresponding to $`U(++)=0`$ and $`U(++)=1`$ by $`A_s^0`$ and $`A_s^1`$, respectively, and analogous $`\mathrm{\Delta }S=0`$ amplitudes by $`A_d^0`$ and $`A_d^1`$.

Absorbing Clebsch-Gordan coefficients in the definition of these amplitudes and specifying the three outgoing particle momenta, one finds

$`A(K^+K^+K^{})_{p_1p_2p_3}`$ $`=`$ $`2A_s^1,`$ (16)
$`A(K^+\pi ^+\pi ^{})_{p_1p_2p_3}`$ $`=`$ $`A_s^0+A_s^1,`$ (17)
$`A(K^+\pi ^+\pi ^{})_{p_2p_1p_3}`$ $`=`$ $`A_s^0+A_s^1,`$ (18)
$`A(\pi ^+\pi ^+\pi ^{})_{p_1p_2p_3}`$ $`=`$ $`2A_d^1,`$ (19)
$`A(K^+\pi ^+K^{})_{p_1p_2p_3}`$ $`=`$ $`A_d^0+A_d^1,`$ (20)
$`A(K^+\pi ^+K^{})_{p_2p_1p_3}`$ $`=`$ $`A_d^0+A_d^1.`$ (21)
The similar forms of corresponding $`\mathrm{\Delta }S=1`$ and $`\mathrm{\Delta }S=0`$ amplitudes may be easily understood in terms of a simple U-spin reflection $`ds`$. The relative signs of terms in amplitudes in which two momenta are interchanged follow from Bose statistics.

Two amplitude relations follow from U-spin,

$`A(K^+\pi ^+\pi ^{})_{p_1p_2p_3}+A(K^+\pi ^+\pi ^{})_{p_2p_1p_3}`$ $`=`$ $`A(K^+K^+K^{})_{p_1p_2p_3},`$ (22)
$`A(K^+\pi ^+K^{})_{p_1p_2p_3}+A(K^+\pi ^+K^{})_{p_2p_1p_3}`$ $`=`$ $`A(\pi ^+\pi ^+\pi ^{})_{p_1p_2p_3}.`$ (23)
On the left-hand-side one sums over amplitudes at two points in phase space where the $`K^+`$ and $`\pi ^+`$ momenta are interchanged.

Squaring Eqs. (16)–(21) and integrating over phase space, one obtains

$`\mathrm{\Gamma }(B^+K^+K^+K^{})`$ $`=`$ $`2\mathrm{\Gamma }_s^1,\mathrm{\Gamma }(B^+K^+\pi ^+\pi ^{})=\mathrm{\Gamma }_s^0+\mathrm{\Gamma }_s^1,`$ (24)
$`\mathrm{\Gamma }(B^+\pi ^+\pi ^+\pi ^{})`$ $`=`$ $`2\mathrm{\Gamma }_d^1,\mathrm{\Gamma }(B^+K^+\pi ^+K^{})=\mathrm{\Gamma }_d^0+\mathrm{\Gamma }_d^1,`$ (25)
where $`\mathrm{\Gamma }_{s,d}^0`$ and $`\mathrm{\Gamma }_{s,d}^1`$ are rates corresponding to the two U-spin amplitudes. Contrary to arguments presented in , this does not imply an equality between $`\mathrm{\Gamma }(B^+K^+K^+K^{})`$ and $`\mathrm{\Gamma }(B^+K^+\pi ^+\pi ^{})`$, or between $`\mathrm{\Gamma }(B^+\pi ^+\pi ^+\pi ^{})`$ and $`\mathrm{\Gamma }(B^+K^+\pi ^+K^{})`$. We conclude that U-spin predictions cannot be tested in simple rate equalities. Instead, as shown in , U-spin predicts equal CP rate differences between all pairs of U-spin related decays. For instance, the CP rate differences in $`B^+K^+K^+K^{}`$ and $`B^+\pi ^+\pi ^+\pi ^{}`$ are equal in the U-spin symmetry limit. Experimental tests of such predictions are quite challenging.

In order to relate $`\mathrm{\Delta }S=1`$ and $`\mathrm{\Delta }S=0`$ processes to each other, we decompose the corresponding effective Hamiltonians into terms multiplying given CKM factors ,

$`_{\mathrm{eff}}^{\overline{b}\overline{s}}`$ $`=`$ $`V_{ub}^{}V_{us}O_u^s+V_{cb}^{}V_{cs}O_c^s,`$ (26)
$`_{\mathrm{eff}}^{\overline{b}\overline{d}}`$ $`=`$ $`V_{ub}^{}V_{ud}O_u^d+V_{cb}^{}V_{cd}O_c^d,`$ (27)
where $`O_u^{d,s}`$ and $`O_c^{d,s}`$ are two operators transforming like U-spin doublets. This implies

$`A_s^{0,1}`$ $`=`$ $`V_{ub}^{}V_{us}A_u^{0,1}+V_{cb}^{}V_{cs}A_c^{0,1},`$ (28)
$`A_d^{0,1}`$ $`=`$ $`V_{ub}^{}V_{ud}A_u^{0,1}+V_{cb}^{}V_{cd}A_c^{0,1},`$ (29)
where $`A_u^{0,1}`$ and $`A_c^{0,1}`$ are matrix elements of the U-spin doublet operators $`O_u`$ and $`O_c`$, respectively, for U-spin final states with $`U(++)=0,1`$.

Eqs. (28) and (29) may be used to test an assumption of penguin and electroweak penguin dominance in strangeness changing decays, namely an assumption that the second term in Eq. (28) dominates the amplitude. We note, however, that in contrast to the isospin analysis which assumed QCD penguin dominance, the second term includes also contributions from electroweak penguin operators.

Applying Eqs. (28) and (29) to $`B^+K^+K^+K^{}`$ and $`B^+\pi ^+\pi ^+\pi ^{}`$, respectively, which involve the amplitudes $`A_u^1`$ and $`A_c^1`$, and defining momentum-dependent amplitudes $`|P|e^{i\delta }V_{cb}^{}V_{cs}A_c^1,|T|e^{i\gamma }V_{ub}^{}V_{us}A_u^1`$, one has

$`A(B^+K^+K^{}K^+)`$ $`=`$ $`|T|e^{i\gamma }+|P|e^{i\delta },`$ (30)
$`A(B^+\pi ^+\pi ^{}\pi ^+)`$ $`=`$ $`|T|\overline{\lambda }^1e^{i\gamma }|P|\overline{\lambda }e^{i\delta },`$ (31)
where

$$\frac{V_{us}}{V_{ud}}=\frac{V_{cd}}{V_{cs}}=\overline{\lambda }\frac{\lambda }{1(\lambda ^2/2)}=0.226,(\lambda =0.22).$$
(32)
In the amplitudes for $`B^{}`$ decays to charge conjugate final states the weak phase $`\gamma `$ appears with opposite signs. Thus, one obtains for the two charge averaged partial widths:

$`\overline{\mathrm{\Gamma }}(B^\pm K^\pm K^\pm K^{})`$ $`=`$ $`{\displaystyle [|T|^2+|P|^2+2|T||P|\mathrm{cos}\gamma \mathrm{cos}\delta ]},`$ (33)
$`\overline{\mathrm{\Gamma }}(B^\pm \pi ^\pm \pi ^\pm \pi ^{})`$ $`=`$ $`{\displaystyle [|T|^2\overline{\lambda }^2+|P|^2\overline{\lambda }^22|T||P|\mathrm{cos}\gamma \mathrm{cos}\delta ]},`$ (34)
where $``$ stands for an integral over three-body phase space. We seek an upper bound on

$$z\left[|T|^2/|P|^2\right]^{1/2},$$
(35)
the tree-to-penguin amplitude ratio averaged over phase space.

(1) Let us assume that $`B^\pm \pi ^\pm \pi ^\pm \pi ^{}`$ is dominated by the tree amplitude and $`B^\pm K^\pm K^\pm K^{}`$ by the penguin. Using the observed value of the ratio

$$\frac{\overline{\mathrm{\Gamma }}(B^\pm \pi ^\pm \pi ^\pm \pi ^{})}{\overline{\mathrm{\Gamma }}(B^\pm K^\pm K^\pm K^{})}=0.35\pm 0.12,$$
(36)
we then find $`z=0.13\pm 0.02<0.16(90\%\mathrm{c}.\mathrm{l}.)`$. However, the possibility of tree-penguin interference weakens this bound somewhat as we show now.

(2) The Schwarz inequality implies that in the definition

$$|T||P|\mathrm{cos}\delta \xi \left[|T|^2\right]^{1/2}\left[|P|^2\right]^{1/2},$$
(37)
the magnitude of the parameter $`\xi `$ cannot exceed 1, and equals 1 when the penguin and tree amplitudes are proportional to each other and are relatively real over the entire phase space. One can then show that $`z`$ is a monotonically increasing function of $`\xi \mathrm{cos}\gamma `$ for any realistic value of $``$. An upper bound $`\xi \mathrm{cos}\gamma 0.74`$ is based on assuming $`\gamma 42^{}`$ . (The bound on $`z`$ is not very sensitive to this assumption.) Since we do not know $`\delta `$ , we shall regard $`\mathrm{cos}\delta `$ as unrestricted. One obtains the maximum value of $`z`$ for $`\xi \mathrm{cos}\gamma =0.74`$: $`z=0.19\pm 0.03`$ or $`z<0.23(90\%\mathrm{c}.\mathrm{l}.)`$.

Let us make a few remarks about deviations from penguin dominance which were neglected in Eq. (13):

1. The bound on the tree-to-penguin amplitude $`z`$ is based on applying U-spin in order to relate $`B^+K^+K^+K^{}`$ and $`B^+\pi ^+\pi ^+\pi ^{}`$. Large U-spin breaking effects are expected to affect a relation between these two processes. This may result in a value of $`z`$ as large as 0.3.
2. Corrections to penguin and electroweak penguin dominance affecting Eq. (13) may differ, in both their flavor and momentum dependence, from those estimated above for $`B^+K^+K^+K^{}`$. Therefore, the bound on $`z`$ can only be used indirectly to set an upper limit of this order on corrections to Eq. (13) from tree amplitudes. We conclude that the tree amplitude contribution to $`A(K^+K^0K^{})+A(K^0K^+\overline{K}^0)`$ could be as large as 0.3 of the equal (but opposite in sign) penguin amplitudes contributing to these two processes.
3. Electroweak penguin contributions affect the amplitude equality (13). Such terms, which were included in the denominator of Eq. (35) but not in its numerator, involve the same weak phase as the dominant penguin amplitude. A rough estimate of electroweak penguin corrections to the amplitude equality, based on Wilson coefficients or on model calculations , is about $`1020\%`$.

The combined correction to Eq. (13) from $`\mathrm{\Delta }I=1`$ tree and electroweak penguin amplitudes is hard to calculate, and depends on the interference between tree and electroweak penguin amplitudes which may be constructive except in special cases . Constructive interference could imply an overall correction of $`4050\%`$ in the most pessimistic case.

The effects of electroweak penguin and tree amplitudes on the time dependent CP asymmetry in $`B^0K^+K^{}K_S`$ are of two kinds. Both amplitudes affect the CP structure of the final state $`f=K^+K^{}K_S`$, thereby multiplying by a dilution factor the coefficient $`S_f`$ of the $`\mathrm{sin}\mathrm{\Delta }mt`$ term, where $`S_f=\mathrm{sin}2\beta `$ for a purely CP-even state excluding $`\varphi K_S`$. One way of determining the dilution factor is through a partial wave analysis in the angular momentum of the $`K^+K^{}`$ system. This may separate even and odd angular momenta corresponding to even and odd CP states.

The tree amplitude ($`T_f`$) in $`B^0K^+K^{}K^0`$, which has a weak phase different from that of the penguin and electroweak penguin amplitudes (combined to $`P_f`$), has another effect. \[We use the same convention for $`T_f`$ and $`P_f`$ as in (30)\]. It modifies $`S_f`$ from a value $`\mathrm{sin}2\beta `$ for a pure CP-even state (excluding $`\varphi K_S`$) to $`\mathrm{sin}2\beta _{\mathrm{eff}},\beta _{\mathrm{eff}}=\beta +\delta \beta `$, and introduces in the asymmetry a $`\mathrm{cos}\mathrm{\Delta }mt`$ term with coefficient $`C_f`$. These corrections depend on the tree-to-penguin ratio of amplitudes and on the relative strong phase between these amplitudes ($`\delta _f`$), both varying over phase space. We define

$$z_f\left[|T_f|^2/|P_f|^2\right]^{1/2},$$
(38)
$$|T_f||P_f|\mathrm{cos}\delta _f\xi _c\left[|T_f|^2\right]^{1/2}\left[|P_f|^2\right]^{1/2},$$
(39)
$$|T_f||P_f|\mathrm{sin}\delta _f\xi _s\left[|T_f|^2\right]^{1/2}\left[|P_f|^2\right]^{1/2}.$$
(40)
Keeping only linear terms in $`z_f`$, one obtains for $`C_f`$ and $`\delta \beta `$ expressions which generalize those obtained for two body decays

$`C_f`$ $``$ $`2z_f\xi _s\mathrm{sin}\gamma ,`$ (41)
$`\delta \beta `$ $``$ $`z_f\xi _c\mathrm{sin}\gamma .`$ (42)
We estimated that $`z_f`$ may be as large as 0.3. The largest possible effect on $`S_f`$ occurs when $`|\xi _c|=1`$, corresponding to penguin and tree amplitudes which are proportional to each other and are relaively real over the entire phase space. In this case $`|\delta \beta |`$ could be as large as $`17^{}`$, and $`S_f`$ could lie in the range $`0.21`$. The Schwarz inequality, $`\xi _c^2+\xi _s^21`$, implies $`\xi _s=0`$ when $`|\xi _c|=1`$, hence $`C_f=0`$. The other extreme but unlikely case, of a maximal $`|C_f|=2z_f\mathrm{sin}\gamma `$ and a minimal $`\delta \beta =0`$, occurs when $`\xi _s=\pm 1,\xi _c=0`$. This corresponds to penguin and tree amplitudes which are proportional to each other and are out of phase by $`\pi /2`$ over the entire phase space.

M. G. wishes to thank the CERN Theory Division and the Enrico Fermi Institute at the University of Chicago for their kind hospitality. This work was supported in part by the United States Department of Energy through Grant No. DE FG02 90ER40560.


text:
# Multilevel Monte Carlo method for simulations of fluids

## 1 INTRODUCTION

The modern theory of classical liquids is based on the statement that the macroscopic characteristics of a many-particle system can be obtained by averaging over microscopic configurations, with the probability proportional to the Gibbs distribution function. Every microscopic configuration is defined by assigning specific locations to the particles, and the main computational problem of statistical physics consists in the difficulty of averaging over the enormous space of possible configurations. In order to estimate the value of this average the Monte Carlo technique for the canonical ensemble was proposed .

However, the straightforward application of this technique is restricted to a small volume of the system under consideration, because one can consider only a relatively small number of particles in any numerical calculation. In order to minimize the surface effect on the bulk values which are to be calculated, artificial periodic boundary conditions are supposed  similar to the case of an infinite crystal .

This means that at scales comparable with the periodicity cell size the fluctuations of the particle number are cut off. In order to avoid this difficulty several approaches have been suggested using the grand canonical ensemble ,,. The main trait of these algorithms is the generation of density fluctuations in the basic cell by adding and deleting particles in accordance with the value of the chemical potential.

The Monte Carlo simulations both in canonical and grand canonical ensembles are very local, moving, e.g., one particle at a time. The operation of adding/deleting a particle is also local and needs to be followed by a local equilibration. This leads to very slow changes of large-scale features, such as averages and various types of clusters (regions of (anti)aligned dipoles, crystallized segments, etc.).The larger the scale the slower the change and longer (per particle!) is the Monte Carlo process required to produce new independent features. Since many independent features are needed for calculating accurate averages, and since very-large-scale features need to be sampled, especially in the vicinity of phase transitions, the computations often become extremely expensive, sometimes even losing practical ergodicity. In recent years, a number of novel Monte Carlo algorithms ensuring ergodicity were proposed , , . In order to avoid the slowness of the ordinary Monte Carlo simulation, steps of more collective nature are used. Nevertheless, the simulation is performed in the periodicity cell and the estimation of bulk values may only be done by extrapolation .

An approach which allows to simultaneously overcome slowness and finite size effects of the conventional Monte Carlo method consists of a multilevel view of the system, realized by multilevel algorithms , . The multilevel algorithms construct a sequence of descriptions of the system under consideration at increasingly coarser levels and transfer information back and forth between the levels in order to obtain a selfconsistent result. The efficiency of multilevel methods in solving problems of statistical physics has been shown on examples with sufficiently simple systems . It follows from these results that the effect of slowing down can be eliminated. Moreover, due to the possibility to simulate large volumes at coarse levels, the volume factor (the proportionality of the computational cost to the volume being simulated) can be suppressed as well. This means that the scale of treating a system is not restricted by the size of the basic cell. The elimination of both the slowing down and the volume factor allows one to investigate in the framework of the multilevel method long range phenomena such as phase transition. This possibility has been shown on two-dimensional examples of variable-coefficient Gaussian models  and the Ising model .

The successful application of the multilevel methods to lattice systems excites interest in adapting them to more complicated cases. The present paper treats multilevel algorithms to the investigation of fluids. Numerical results of simulations for one dimensional models are presented. These models reproduce common properties of fluids and allow comparisons to exact solutions.

## 2 MULTILEVEL MONTE CARLO APPROACH

### 2.1 Conventional Monte Carlo Method

The Monte Carlo method in the statistical theory of liquids is used to evaluate numerically the average $`\overline{A}`$ of any functional $`A`$, defined by:

$$\overline{A}=\underset{\mathrm{\Omega }}{}A(\underset{¯}{X})w(\underset{¯}{X})𝑑\underset{¯}{X}\frac{1}{m}\underset{i=1}{\overset{m}{}}A(\underset{¯}{X_i})$$
(1)
where $`w(\underset{¯}{X})`$ is the probability density of the state $`\underset{¯}{X}`$ in the configuration space $`\mathrm{\Omega }`$, and the nodes $`\underset{¯}{X_i}`$ are generated by a random walk in $`\mathrm{\Omega }`$ that satisfies detailed balance.

The simplest definition of the probability to walk from $`\underset{¯}{X}`$ to $`\underset{¯}{X^{}}`$ in detailed balance is given by :

$$\omega (\underset{¯}{X}\underset{¯}{X^{}})=\mathrm{min}[1,\frac{w(\underset{¯}{X^{}})}{w(\underset{¯}{X})}],$$
(2)
provided that the probability $`P(\underset{¯}{X^{}}\mathrm{?}|\underset{¯}{X})`$ for having chosen the state $`\underset{¯}{X^{}}`$ as the candidate to replace the current state $`X`$ is symmetric, i.e., $`P(\underset{¯}{X^{}}\mathrm{?}|\underset{¯}{X})=P(\underset{¯}{X}\mathrm{?}|\underset{¯}{X^{}})`$.

The kind of probability density given by Gibbs in the canonical ensemble for simple liquids is:

$$w(\underset{¯}{X})=\mathrm{𝑐𝑜𝑛𝑠𝑡}\mathrm{exp}(\frac{U(\stackrel{}{r}_\mathit{1},\stackrel{}{r}_\mathit{2},\mathrm{},\stackrel{}{r}_N)}{k_BT})$$
(3)
where $`\stackrel{}{r}_i`$ is the location of the $`i`$-th particle ($`i=1,\mathrm{},N`$), $`k_B`$ is the Boltzmann constant, T is the temperature and $`U`$ is the potential energy. It is usually assumed to have the form:

$$U=\underset{i<j}{}\varphi (\stackrel{}{r}_i\stackrel{}{r}_j)$$
(4)
where $`\varphi (\stackrel{}{r}_i\stackrel{}{r}_j)`$ corresponds to the energy of a two-body interaction.

In any numerical calculation it is possible to consider only relatively small number of particles (extremely small in comparison with Avogadro number). In order to minimize the surface effect in a small simulation volume, periodic boundary conditions are supposed . This means that the space is filled by particles which are located at points:

$$\stackrel{}{R}_{i,\stackrel{}{l}}=\stackrel{}{r}_i+\stackrel{}{l}L,\text{ }1iN$$
(5)
where $`\stackrel{}{l}`$ is a vector whose components are integer, particles at points $`\stackrel{}{r}_i`$ are located inside the basic cell of the linear size $`L`$, and $`N`$ is the particle number in the periodicity cell.

It follows from the periodicity condition that the real system is replaced by a super-lattice with the same value of the particle number density $`\rho =N/L^D`$ (where $`D`$ being the dimension of the space) in each cell.

The motion of particles is continuous and in any configuration the particles have arbitrary locations. The transition between states is made by the shift of one particle at a time  by a small amount via the following displacement:

$$\stackrel{}{r}_i^{new}=\stackrel{}{r}_i^{old}+\delta _{max}(\stackrel{}{1}2\stackrel{}{\xi })$$
(6)
where $`\stackrel{}{r}_i^{old}`$ and $`\stackrel{}{r}_i^{new}`$ are old and new locations of the $`i`$-th particle, $`\delta _{max}`$ is the maximum possible displacement of the particle along a coordinate axes, $`\stackrel{}{1}`$ is the vector all of whose components are $`1`$, and $`\stackrel{}{\xi }`$ is a vector whose components are random numbers distributed uniformly on the interval $`[0,1]`$.

For shifting the particle with the number $`i`$, say, one can see from (2) that it is enough to use, instead of the Gibbs function (3) with the energy (4), the conditional probability defined by:

$$P(\stackrel{}{r}_i𝐑_𝐢)=\mathrm{𝑐𝑜𝑛𝑠𝑡}\mathrm{exp}(u_i(𝐑_𝐢)/k_BT)$$
(7)
The last relation gives us the probability to find the particle with the number $`i`$ at the position $`\stackrel{}{r}_i`$ when the locations of all other particles, defined by the set $`𝐑_𝐢=\{\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_{i1},\stackrel{}{r}_{i+1},\mathrm{},\stackrel{}{r}_N\}`$, are fixed. The one-particle energy $`u_i(𝐑_𝐢)`$ in (7) is defined by:

$$u_i(𝐑_𝐢)=\underset{j,(ji)}{}\varphi (\stackrel{}{r}_i\stackrel{}{r}_j)$$
(8)
It follows from (6) that the conventional Monte Carlo simulation is the local process. Therefore one can expect the high speed convergence only in the case of systems with short-range correlations.

### 2.2 Coarse variables

A slowing-down is inherent not only in the conventional Monte Carlo algorithm, it is a common problem for all local processes (e.g. Gauss- Seidel relaxation for discretizied partial differential equations). The solution to this problem lies in introducing system changes of more collective nature. In the case of partial differential equations fast convergence of solutions had been attained by multigrid algorithms . These algorithms are looking for solutions on a sequence of lattices with increasingly larger meshsize (coarser scales) by combining local processing at each scale with various inter-scale (inter-lattice) interactions. A similar technique can be applied to the simulation of liquids.

A possible way to introduce a coarse description of liquid consists in the discretization of space. The periodicity cell is divided into $`M`$ disjoint parts (e.g. cubes) $`V_i^1`$ of equal volume with linear size $`h_1`$, $`1iM`$ (each $`V_i^1`$ being associated with a gridpoint $`i`$ of the first coarse-level lattice). Configurations of the finest (particle) level are mapped to the first coarse level by the operation of coarsening, this operation creates the coarse-level variable set. For example, at any instant the corresponding coarse-level variable can be defined by coarsening the particle number:

$$n_i^1=\text{ Number of particles in }V_i^1$$
(9)
with $`\underset{i=1}{\overset{M}{}}n_i=N`$, where $`N`$ is the total number of particles in the periodicity cell. The set $`\{n_i^1\}`$ defines the current configuration on the first coarse-level: instead of particle locations the occupation numbers at gridpoints are used.

Generally, the aim of coarsening is the creation of configurations represented by collective variables which describe collective particle motions at different scales. The variable at each lattice point at each coarse level is defined as a local spatial average (an average or actually a sum over a certain neighborhood of the lattice point) of similar variables at the next finer level. The total value of each such variable is well defined for each configuration of particles given at the finest level.

The extension of the coarsening operation (9) to coarser levels leads to the following definition of the coarse-variable at the level $`k`$:

$$n_j^k=\underset{V_i^{k1}V_j^k}{}n_i^{k1},\text{ }k>1$$
(10)
for each volume element $`V_j^k`$ of level $`k`$, assuming it to be a union of volume elements of the level $`k1`$. The coarsening can be repeated till the coarsest level, whose choice depends on the scale of the phenomena one wants to compute.

The coarser level simulations are performed in periodicity cells of larger sizes. In this way finite size effects are suppressed and it is possible to reach the macroscopic description of the system. No slowdown should occur provided that the coarsening ratio (the ratio between a coarse meshsize and the next finer meshsize), as well as the average number of of original particles per mesh volume of the finest lattice, are suitably low. The typical meshsize ratio is 2, typical number of particles per finest lattice mesh is between 2 and 10 (being usually larger at a higher dimension). More aggressive coarsening ratios would require much longer simulations to produce information for coarser levels.

There are many possible ways to choose the set of coarse variables. A general criterion for the quality of this set is the speed of equilibration of a compatible Monte Carlo (CMC). By this we mean a Monte Carlo process on the fine level which is restricted to the subset of fine-level configurations whose local spatial averages coincide with a fixed coarse-level configuration. A fast CMC equilibration implies that up to local processing all equilibrium configurations are fully determined by their coarse-level representations (their local spatial averages), which is the main desired property of coarsening.

The compatible Monte Carlo equilibration speed can be tested at the fine level, in the following way: after local thermalization of the ordinary Monte Carlo process, for a given configuration in the equilibrium the corresponding coarse variables are calculated by (9) or (10). Then an initial configuration in the original variables is created in accordance with the coarse variables set. A following ensemble of compatible Monte Carlo processes yields an accurate estimate of the rate of approach to the equilibrium. The interpolation from the first coarse level configuration to the finest (particle) level needs a special consideration owing to the mapping of variables defined on a grid into the continuous space of particle coordinates.

### 2.3 CP Tables and Multilevel Algorithm

The main idea of the multilevel approach is to equilibrate on each level only modes with short (comparable with the meshsize) wave lengths. Long wave modes with slow convergence at a given level are equilibrated at coarser levels where their wave lengths are comparable with the meshsize. As a result, the multilevel process leads to fast equilibration of all modes.

In the framework of this multilevel Monte Carlo algorithm, only a local process is performed at each level, defined in terms of the corresponding variables. In order to calculate transition probabilities (2) conditional probabilities similar to (7) should be derived for each coarse level. Such conditional probabilities are expressed in the form of a Conditional Probability (CP) table, which in principle tabulates numerically the probability distribution of any pair of neighboring coarse-level variables given the values of all others. Of course, not all other variables should in practice be taken into account: only the immediate neighborhood of gridpoint pair under test counts, due to the near locality property of the conditional probability (cf. also discussion of near locality in ).

For example, in terms of variables (9), (10) defined at gridpoints a conditional probability tables $`P_k(n_i^k,n_j^ks_i^1,\mathrm{},s_i^l)`$ can be constructed from a given sequence of configurations in equilibrium on the next finer level $`k1`$. These tables give us the dependence of the probability to find $`n_i^k`$ and $`n_j^k`$ particles at the two neighboring gridpoints $`i`$ and $`j`$ on $`l`$ values in their neighborhood by:

$$s_{ij}^m=\underset{qneighborhoodofiandj}{}\alpha _q^mn_q^k,1ml$$
(11)
where $`\alpha _q^m`$ are preassigned, suitably chosen coefficients, with $`\alpha _i^m=\alpha _j^m=0`$ for all $`m`$.

For example, in one dimension $`j=i+1`$ and a possible choice, for $`l=1`$, $`\alpha _{i1}^1=\alpha _{i+2}^1=1`$, otherwise $`\alpha _q^m=0`$.

In the coarse level Monte Carlo run, each trial move consists of particle exchange between two neighboring gridpoints, i.e. $`n_i^kn_i^k=n_i^k+\mathrm{\Delta }n,n_j^kn_j^k=n_j^k\mathrm{\Delta }n`$.In accordance with (2), the acceptance probability for this move is:

$$\omega (\underset{¯}{X}\underset{¯}{X}^{})=\mathrm{min}[1,\frac{P_k(n_i^k,n_j^ks_i^1,\mathrm{},s_i^l)}{P_k(n_i^k,n_j^ks_i^1,\mathrm{},s_i^l)}]$$
(12)
The CP tables for any coarse level $`k`$ are calculated by gathering appropriate statistics during Monte Carlo simulations at the next finer level $`k1`$. Because of the near-locality property, no global equilibration is needed; local equilibration is enough to provide the correct CP values for any neighborhood for which enough cases have appeared in the simulation. Thus, the fine-level simulation can be done in a relatively small periodicity cell.

However, since the fine-level canonical ensemble simulations use only a small periodicity cell, many types of neighborhoods that would be typical at some parts of a large volume (e.g., typical at parts with average densities different than that used in the periodicity cell) will not show up or will be too rare to have sufficiently accurate statistics. Hence, simulations at some coarse level may run into a situation in which the CP table being used has flags indicating that values one wants to extract from it start to have poor accuracy. In such a situation, a temporary local return to finer levels should be made, to accumulate more statistics, relevant for the new local conditions.

To return from a coarse level to the next finer level one needs first to interpolate, i.e., to produce the fine level configurations represented by the current coarse level configuration, with correct relative probabilities. The interpolation is performed by CMC sweeps at the fine level (few sweeps are enough, due to the fast CMC equilibration). This fast equilibration also implies that the interpolation can be done just over a restricted subdomain, serving as a window: In the window interior good equilibrium is reached. Additional passes can then be made of ordinary (not compatible) MC, to accumulate in the interior of the window the desired additional CP statistics, while keeping the window boundary frozen (i.e., compatible). The window can then be coarsened (by the local spatial averaging) and returned to the coarse level, where simulations can now resume with the improved CP table.

Iterating back and forth between increasingly coarser levels and window processing at finer levels whenever missing CP statistics is encountered, one can quickly converge the required CP tables at all levels of the system, with only relatively small computational domains employed at each level. The size of those domains needs only be several times larger than the size of the neighborhoods being used (with a truncation error that then decreases exponentially with that size). However, larger domains are better, since they provide sampling of a richer set of neighborhoods (diminishing the need for returning later to accumulate more statistics), and since the total amount of work at each level depends anyway only on the desired amount of statistics, not on the size of the computational domain.

## 3 NUMERICAL TESTS

Thermodynamic properties of a one-dimensional fluid can be described exactly under the assumption of the nearest neighbor finite range interparticle interaction ,. One can cite as an example of the short-range potential the truncated Lennard-Jones potential:

$$\varphi _{LJ}(r)=\{\begin{array}{cc}4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right],\hfill & \text{for }rr_{cut}\hfill \\ 0,\hfill & \text{for }r>r_{cut}\hfill \end{array}$$
(13)
where $`\sigma `$ is related to the diameter of atoms, $`ϵ`$ is the measure of interparticle interaction and $`r_{cut}`$ is the cut-off distance.

At high temperatures the attractive part of the potential can be neglected and a good approximation for the interparticle interaction is the hard rods model:

$$e^{\frac{\varphi _H(r)}{k_BT}}=\{\begin{array}{cc}0\hfill & \text{}rd\hfill \\ 1\hfill & \text{}r>d\hfill \end{array}$$
(14)
where $`d`$ is the diameter of every particle.

In order to test the multilevel algorithm it was applied to the simulation of one-dimensional fluids. First of all, the speed of the compatible Monte Carlo equilibration was tested in the case of the interparticle interaction (13) with $`{}_{cut}{}^{}=2.5\sigma `$. After locally equilibrating an ordinary Monte Carlo process, a sequence of equilibrium configurations was picked out. For each configuration the coarse-level variable set is defined by (9) and an initial configuration for compatible Monte Carlo simulation is constructed in accordance with the frozen equilibrium coarse level structure.

For the Lennard-Jones fluid a convenient quantity measuring the approach to thermal equilibrium is the energy of the system. However, this quantity is subject to fluctuations during a single run of the equilibration. In order to smooth out simulation results a number of compatible Monte Carlo simulations are performed using the same initial configuration but different sequences of random numbers. The relaxation function is defined by the averaging over these CMC runs:

$$R_{U_0}(\tau )=\frac{<U(\tau )><U(\tau _f)>}{U_0<U(\tau _f)>}$$
(15)
where $`U(\tau )`$ is the energy of the system after sweep $`\tau `$, subscript $`f`$ means final configurations, and $`U_0=U(0)`$.

The ensemble average relaxation function (15) is used to estimate the relaxation time:

$$\tau _{eq}(U_0)=\underset{0}{\overset{\tau _f}{}}R_{U_0}(\tau )𝑑\tau $$
(16)
Here the relaxation time depends on the initial configuration of CMC runs, therefore additional averaging of the relaxation time (16) is done over all the equilibrium configurations which were chosen during the ordinary Monte Carlo process as the initial configuration set for CMC:

$$\tau _{eq}=<\tau _{eq}(U_0)>_{U_0}$$
(17)
Definitions (16) and (17) lead to the exact result in the case of a single classical relaxation behavior. Nevertheless, even if the relaxation behavior is polydispersive, the integral (16) still gives an unambiguous measure of the relaxation rate .

The temperature dependence of the relaxation time estimated in accordance with (17) is shown in Fig.1. In contrast to the ordinary Monte Carlo process the relaxation of CMC is not sensitive to the size of the periodicity domain and a fast equilibration implies that coarsening the particle number satisfies (in the case of the interaction potential (13) and, hence, in the high temperature limit (14) the requirements for the quality of the coarse variables set.

In the special case of the one-dimensional fluid, the main contribution to the distribution of the particle number $`n_i^k`$ at any lattice point is expected to come from its nearest neighbors. Therefore a simple possible form of the conditional probability table is given by:

$$P_k(n_i^k,n_j^ks_{ij}^1,\mathrm{},s_{ij}^l)=P_k(n_i^k,n_{i+1}^kn_{i1}^k,n_{i+2}^k)$$
(18)
Quantities of $`n_i^k`$ were calculated at each lattice point $`i`$ for each produced configuration at the next finer level, and from them the probability distributions (18) were accumulated, forming the CP table. Then, multilevel Monte Carlo simulation was performed using this CP tables.

A suitable quantity for comparing simulation results at different levels is the fluctuation $`\nu _k`$ of the particle number in the subdomain of size $`V^k`$ (corresponding to the gridpoint of the lattice with meshsize $`h_k=h_12^{k1}`$):

$$\nu _k=\frac{<n^{k^2}><n^k>^2}{<n^k>}$$
(19)
where the mean values are defined by probabilities $`p_k(n)`$ to display the particle number $`n^k`$ in a volume of size $`V^k`$: $`<n^k>=\underset{n^k}{}n^kp_k(n^k)`$, $`<n^{k^2}>=\underset{n^k}{}n^{k^2}p_k(n^k)`$.

The fluctuation of the particle number does not depend on the volume if it is sufficiently large (as compared with the correlation length) and is related with the isothermal compressibility $`\kappa =(\frac{1}{V}\frac{V}{P})_{T=\mathrm{𝑐𝑜𝑛𝑠𝑡}}`$ . For 1D hard core system the equation of state is found analytically (it is known as the Tonk’s equation ) and the isothermal compressibility is given by:

$$\nu _{\mathrm{}}=\rho k_BT\kappa =(1\rho d)^2$$
(20)
Results of particle number fluctuation multilevel measurements in volumes of different sizes for hard rods systems are presented in Fig.2. The dependence of this quantity on the subdomain size is explained by the size effects and coincides with the expression of the grand canonical ensemble :

$$\nu _k=\nu _{\mathrm{}}+\frac{C}{V^k}$$
(21)
where $`\nu _{\mathrm{}}`$ is the bulk value of the particle number fluctuation (20). The coefficient $`C`$ as well as $`\nu _{\mathrm{}}`$ depend on the particle number density. One can see from Fig.2 that at high temperatures (the interaction model (14) five levels are enough in order to achieve the bulk behavior in subdomains connected with gridpoints of the coarsest level.

In the case of the Lennard-Jones fluid the particle number fluctuations depends also on the temperature. The dependence on the meshsize (the subdomain volume) for system with the interaction potential (13) is shown in Fig.3 for different temperatures. One can see that at high temperatures the properties of the Lennard-Jones fluid are similar to the hard rods system (see Fig.2). At low temperatures the behavior of fluctuations drastically changes. The results obtained indicate that in contrast to the hard rods system the Lennard-Jones system loss homogeneity at low temperatures on fine scales.

## 4 CONCLUSION

It follows from the present results that in the framework of the multilevel Monte Carlo method a suitable choice of coarsening allows to attain the macroscopic behavior. The particle number fluctuation in subdomains of intermediate size tends to the exact value on increasingly coarser levels. The accuracy of averages (including the particle number fluctuation) depends only on the CP tables adequacy, which is achieved by the correct selection of the variable set.

The advantage of the multilevel Monte Carlo method consists in fast convergence of measured mean values due to the selfconsistent equilibration on different levels modes of wave lengths comparable with the simulation domain. The computational work on each level is proportional to the number of gridpoints and is independent of the particle number associated with the gridpoint. It leads to the high speed of the method as compared with conventional algorithms.

The equilibrium on fine levels (and in particular the frequency of the appearance of a given particle number in the simulation domain) is defined by the canonical ensemble configurations on coarsest level. The particle number in fine level simulation domains is variable and its distribution imitates the result of the grand canonical ensemble simulation . The multilevel Monte Carlo method is not concerned with the equilibrium in accordance with the value of a chemical potential. It opens new way for the development of the one-phase approach to the phase-transition problem .

## 5 ACKNOWLEDGMENTS

The research has been supported by Israel Absorption Ministry, project No. 6682, by the U.S. Air Force Office of Scientific Research, contract No. F33615-97-D-5405, by the European Office of Aerospace Research and Development (EOARD) of the U.S. Air FORCE, contract No. F61775-00-WE067, by Israel Science Foundation Grant No. 696/97 and by the Carl F.Gauss Minerva Center for Scientific Computation at the Weizmann Institute of Science.


text:
# A New Symmetric Expression of Weyl Ordering

## 1 Introduction

The problem for ordering of operators in Quantum Physics has been investigated then and now, see, for example, , , ,  and their references.

Quantization

$$\{f,g\}_{\text{Pb}}\frac{i}{\mathrm{}}[\widehat{f},\widehat{g}]$$

is a Lie algebra homomorphism, but not a ring homomorphism, where $`\{,\}_{\text{Pb}}`$ denotes a Poisson bracket and $`[,]`$ a commutator of operators on a Fock space. That is

$$fg=gf\text{but}\widehat{f}\widehat{g}\widehat{g}\widehat{f}.$$

This describes the noncommutativity in Quantum Physics. Therefore there exist various operator orderings and we need to select a suitable one depending on a problem in question.

Let $`a^{}`$ and $`a`$ be the creation operator and the annihilation operator of a harmonic oscillator respectively. The relation is

$$[a,a^{}]=1.$$

We put the number operator $`Na^{}a`$.

There are various useful orderings in Quantum Physics . In this paper, we deal with normal ordering, anti-normal ordering, and Weyl ordering. We denote them as follows:

1. Normal Ordering $`:(a^{})^na^m:_N(a^{})^na^m`$,
2. Anti-Normal Ordering $`:(a^{})^na^m:_{AN}a^m(a^{})^n`$,
3. Weyl Ordering
$`:(a^{})^na^m:_W\left(\begin{array}{c}n+m\\ n\end{array}\right)^1(\text{sum of all symmetric products of }na^{}\text{ and }ma)`$,

For example, when $`m=n`$,

$`:a^{}a:_W`$ $`=`$ $`{\displaystyle \frac{1}{2}}(a^{}a+aa^{}),`$ (1.1)
$`:(a^{}a)^2:_W`$ $`=`$ $`{\displaystyle \frac{1}{6}}(a^{}a^{}aa+a^{}aa^{}a+a^{}aaa^{}+aa^{}a^{}a+aa^{}aa^{}+aaa^{}a^{}),`$ (1.2)
$`:(a^{}a)^3:_W`$ $`=`$ $`{\displaystyle \frac{1}{20}}(a^{}a^{}a^{}aaa+a^{}a^{}aa^{}aa+a^{}a^{}aaa^{}a+a^{}a^{}aaaa^{}+`$ (1.3)
$`a^{}aa^{}a^{}aa+a^{}aa^{}aa^{}a+a^{}aa^{}aaa^{}+a^{}aaa^{}a^{}a+`$
$`a^{}aaa^{}aa^{}+a^{}aaaa^{}a^{}+aa^{}a^{}a^{}aa+aa^{}a^{}aa^{}a+`$
$`aa^{}a^{}aaa^{}+aa^{}aa^{}a^{}a+aa^{}aa^{}aa^{}+aa^{}aaa^{}a^{}+`$
$`aaa^{}a^{}a^{}a+aaa^{}a^{}aa^{}+aaa^{}aa^{}a^{}+aaaa^{}a^{}a^{}),`$
here we write $`:(a^{})^na^n:_{}`$ as $`:(a^{}a)^n:_{}`$ for simplicity ($``$=N, AN, W).

We often make proper use of these orderings. Therefore it is important to investigate what are the relations among them. The intertwining formulas of these orderings have been known by ,  as a formula of “$`s`$-ordered power-series expansions”;

$$\{(a^{})^na^m\}_s=\underset{k=0}{\overset{\mathrm{min}\{n,m\}}{}}k!\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}m\\ k\end{array}\right)\left(\frac{ts}{2}\right)^k\{(a^{})^{nk}a^{mk}\}_t,$$
(1.4)
where

$$\{(a^{})^na^m\}_s\frac{^{n+m}}{\alpha ^n(\alpha ^{})^m}\text{e}^{\alpha a^{}\alpha ^{}a}\text{e}^{s|\alpha |^2/2}|_{\alpha =0},$$

and the ordering specified $`s,t=+1,0,1`$ are, respectively, normal, Weyl, and anti-normal. On the other hand, the intertwining formulas among three orderings are used in the Weyl calculus. Recently, Omori et al.  have investigated strange phenomena of $``$-exponential functions of quadratic forms in the Weyl algebra. They have written the intertwining formulas as a form of linear operators by using the product formulas of the $``$-exponential functions;

$$\text{e}^{\frac{1}{2}_a^{}_a}:f(a^{},a):_N=:f(a^{},a):_W,$$
(1.5)
$$\text{e}^{\frac{1}{2}_a^{}_a}:f(a^{},a):_{AN}=:f(a^{},a):_W,$$
(1.6)
where $`f(a^{},a)`$ is a suitable 2-variable function $`f(\overline{z},z)`$ substituting $`a^{},a`$ for $`\overline{z},z`$ respectively.

Under these circumstances, we investigate expressions of Weyl ordering with the intertwining formulas of orderings. The motivation is as follows; The definition of Weyl ordering is clear, but it is hard to give an explicit expression of it, like (1.1), (1.2), (1.3), especially when $`n`$ is large. Therefore it is important to search a simple expression of it. If we study by using the formula (1.4) or (1.5), (1.6), though (1.4) give some explicit expressions of Weyl ordering, it is not “symmetric” that Weyl ordering would be. Therefore, we study this problem in another point of view, namely in the combinatorial theory and in the difference theory. Then we find a symmetric expression of Weyl ordering by using the number operator $`N`$, which is new and would be useful in Quantum Physics. Since $`a^{}a=N`$ and $`aa^{}=N+1`$, we shall show Weyl ordering $`:(a^{}a)^n:_W`$ as a symmetric polynomial with respect to $`N`$ and $`N+1`$.

For example, (1.1), (1.2), (1.3) are

$`:a^{}a:_W`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{N+(N+1)\},`$ (1.7)
$`:(a^{}a)^2:_W`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{N^2+(N+1)^2\},`$ (1.8)
$`:(a^{}a)^3:_W`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{N^3+(N+1)^3\}+{\displaystyle \frac{1}{4}}\{N+(N+1)\}.`$ (1.9)
In fact, for all $`n=1,2,\mathrm{},`$ we can show $`:(a^{}a)^n:_W`$ as a symmetric form like these examples in section 6. It is a new and very clear-cut expression of Weyl ordering.

The contents of this paper are as follows. In section 2, we prepare the difference theory and the Stirling number of the first kind. In section 3, first we prepare some facts about normal and anti-normal ordering. Next, we give a simple proof for the intertwining formulas of orderings as above for a help of the following sections. In section 4, we interpret these expressions by the difference theory. Then we find that the noncommutativity of harmonic oscillators is nothing but the the Newton expansion w.r.t. variable $`N`$. In section 5, we describe nontrivial relations including the Stirling number of the first kind. These relations are derived from the noncommutativity of harmonic oscillators and are needed in the following section. In section 6, we describe the main result, that is a new symmetric expression of Weyl ordering. In section 7, we summarize our results.

## 2 Mathematical Preliminaries

In this section, we prepare the difference theory and the Stirling number of the first kind.

First we define factorial of degree $`n`$ with an increment $`ϵ`$

$$x^{\underset{¯}{n}}x(xϵ)\mathrm{}(x(n1)ϵ)$$
(2.1)
and

$$x^{\overline{n}}x(x+ϵ)\mathrm{}(x+(n1)ϵ).$$
(2.2)
We remark

$$(x+ϵ)^{\overline{n}}=(x+ϵ)\mathrm{}(x+nϵ)=(x+nϵ)^{\underset{¯}{n}}.$$
(2.3)
Let $`p_n(x)`$ be a polynomial of degree $`n`$ w.r.t. $`x`$, and $`\mathrm{\Delta }_{+x}`$ the forward difference operator defined by $`\mathrm{\Delta }_{+x}f(x)ϵ^1\{f(x+ϵ)f(x)\}.`$ Then we have the Newton expansion

$$p_n(x)=\underset{k=0}{\overset{n}{}}\frac{x^{\underset{¯}{k}}}{k!}\mathrm{\Delta }_{+x}^kp_n(0)=\underset{k=0}{\overset{n}{}}\frac{x^{\underset{¯}{nk}}}{(nk)!}\mathrm{\Delta }_{+x}^{nk}p_n(0).$$
(2.4)
We also remark that

$$\mathrm{\Delta }_{+x}x^{\underset{¯}{n}}=nx^{\underset{¯}{n1}}.$$
(2.5)
Next we define the Stirling number of the first kind and calculate some of them explicitly.

###### Definition 2.1.

The Stirling number of the first kind $`s(n,i)(i=1,\mathrm{},n)`$ is defined by the following equation ;

$$x^{\underset{¯}{n}}=x(xϵ)\mathrm{}(x(n1)ϵ)=\underset{i=1}{\overset{n}{}}s(n,i)ϵ^{ni}x^i$$

and $`s(n,0)0(n1)`$, $`s(j,i)0(j<i)`$.

We remark

$$x^{\overline{n}}=x(x+ϵ)\mathrm{}(x+(n1)ϵ)=\underset{i=1}{\overset{n}{}}(1)^{ni}s(n,i)ϵ^{ni}x^i.$$
(2.6)
###### Remark 2.2.

We can obtain the Stirling number of the first kind explicitly by the following recursion formula and initial conditions.

$`s(n+1,i)`$ $`=`$ $`s(n,i1)ns(n,i)(ni1),`$ (2.7)
$`s(n,n)`$ $`=`$ $`1(n1),`$ (2.8)
$`s(n,0)`$ $`=`$ $`0.`$ (2.9)
For example, if we put $`i=n`$ in (2.7), then

$`s(n+1,n)`$ $`=`$ $`s(n,n1)ns(n,n)`$
$`=`$ $`s(n,n1)n.`$
Therefore

$$s(n,n1)=\underset{k=1}{\overset{n1}{}}k=\frac{n(n1)}{2}.$$
(2.10)
Next, since

$`s(n+1,n1)`$ $`=`$ $`s(n,n2)ns(n,n1)`$
$`=`$ $`s(n,n2)+{\displaystyle \frac{1}{2}}n^2(n1),`$
we have

$$s(n,n2)=\underset{k=1}{\overset{n1}{}}\frac{1}{2}k^2(k1)=\frac{1}{24}n(n1)(n2)(3n1).$$
(2.11)
## 3 A Simple Proof for the Intertwining Formulas of Orderings

First we prepare some facts about normal and anti-normal ordering. For a while, we put $`[a,a^{}]=ϵ`$ for more general situation.

###### Lemma 3.1.

We have following formulas using the number operator $`N=a^{}a`$;

$$:(a^{}a)^n:_N=N(Nϵ)\mathrm{}(N(n1)ϵ)=N^{\underset{¯}{n}},$$
(3.1)
$$:(a^{}a)^n:_{AN}=(N+ϵ)(N+2ϵ)\mathrm{}(N+nϵ)=(N+ϵ)^{\overline{n}}.$$
(3.2)
###### Proof.

By using simple relations

$$aN=(N+ϵ)a,a^{}N=(Nϵ)a^{}\text{or}Na=a(Nϵ),Na^{}=a^{}(N+ϵ),$$

we have an induction foumula

$`(a^{})^na^n`$ $`=`$ $`a^{}\mathrm{}a^{}a^{}a^{}aaa\mathrm{}a`$
$`=`$ $`a^{}\mathrm{}a^{}a^{}Naa\mathrm{}a`$
$`=`$ $`a^{}\mathrm{}a^{}(Nϵ)a^{}aa\mathrm{}a`$
$`=`$ $`\mathrm{}`$
$`=`$ $`(N(n1)ϵ)(a^{})^{n1}a^{n1}.`$
Therefore we obtain (3.1) by the mathematical induction. We also obtain (3.2) in a similar way. ∎

As we mentioned in section $`1`$, since Omori et al have proved (1.5), (1.6) by using the product formulas of the $``$-exponential functions, we give an another simple proof of them. We have only to prove the next proposition. This proof helps account for our theory to be described in the following sections.

###### Proposition 3.2.

$$\text{e}^{\frac{1}{2}ϵ_a^{}_a}:(a^{})^na^m:_N=:(a^{})^na^m:_W,$$
(3.3)
$$\text{e}^{\frac{1}{2}ϵ_a^{}_a}:(a^{})^na^m:_{AN}=:(a^{})^na^m:_W.$$
(3.4)
###### Proof.

By using the generating function of the Weyl ordering,

$$\text{e}^{\alpha a^{}+\beta a}=\underset{n,m=0}{\overset{\mathrm{}}{}}\frac{\alpha ^n\beta ^m}{n!m!}:(a^{})^na^m:_W$$

and the fundamental Baker-Campbell-Hausdorff formula

$$\text{e}^{A+B}=\text{e}^{\frac{1}{2}[A,B]}\text{e}^A\text{e}^B\text{whenever }([A,[A,B]]=[B,[A,B]]=0\text{), }$$
(3.5)
we have

$`\text{e}^{\alpha a^{}+\beta a}`$ $`=`$ $`\text{e}^{\frac{1}{2}ϵ\alpha \beta }\text{e}^{\alpha a^{}}\text{e}^{\beta a}(\text{normal order})`$ (3.6)
$`=`$ $`\text{e}^{\frac{1}{2}ϵ_a^{}_a}\text{e}^{\alpha a^{}}\text{e}^{\beta a}.`$
If we expand both sides of (3.6) and compare the coefficients of $`\alpha ^n\beta ^m`$, we obtain (3.3). Similarly, we have

$`\text{e}^{\alpha a^{}+\beta a}`$ $`=`$ $`\text{e}^{\frac{1}{2}ϵ\alpha \beta }\text{e}^{\beta a}\text{e}^{\alpha a^{}}(\text{anti-normal order})`$ (3.7)
$`=`$ $`\text{e}^{\frac{1}{2}ϵ_a^{}_a}\text{e}^{\alpha a^{}}\text{e}^{\beta a}.`$
If we expand both sides of (3.7) and compare the coefficients of $`\alpha ^n\beta ^m`$, we obtain (3.4). ∎

By the proposition 3.2 and the lemma 3.1, we obtain the following corollary.

###### Corollary 3.3.

We can show Weyl ordering $`:(a^{})^na^m:_W`$ in terms of normal ordering.

$$:(a^{})^na^m:_W=\underset{k=0}{\overset{\mathrm{min}\{n,m\}}{}}\frac{ϵ^kk!}{2^k}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}m\\ k\end{array}\right)(a^{})^{nk}a^{mk}.$$
(3.8)
Especially, in case of $`m=n`$, if we use (3.1), we obtain an expression of $`:(a^{}a)^n:_W`$ in terms of the number operator $`N`$,

$$:(a^{}a)^n:_W=\underset{k=0}{\overset{n}{}}\frac{ϵ^kk!}{2^k}\left(\begin{array}{c}n\\ k\end{array}\right)^2N^{\underset{¯}{nk}}.$$
(3.9)
By a similar calculation, we obtain the next expression using anti-normal ordering.

###### Corollary 3.4.

$$:(a^{})^na^m:_W=\underset{k=0}{\overset{\mathrm{min}\{n,m\}}{}}\frac{(ϵ)^kk!}{2^k}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}m\\ k\end{array}\right)a^{nk}(a^{})^{mk}.$$
(3.10)
Especially, in case of $`m=n`$, if we use (3.2), we obtain an another expression of $`:(a^{}a)^n:_W`$ in terms of the number operator $`N`$,

$$:(a^{}a)^n:_W=\underset{k=0}{\overset{n}{}}\frac{(ϵ)^kk!}{2^k}\left(\begin{array}{c}n\\ k\end{array}\right)^2(N+ϵ)^{\overline{nk}}.$$
(3.11)
###### Remark 3.5.

In , , K.E.Cahill and R.J.Glauber introduced “$`s`$-ordered power-series expansions” as follows;
Define the $`s`$-ordered displacement operator $`D(\alpha ,s)`$ by

$$D(\alpha ,s)=\text{e}^{\alpha a^{}\alpha ^{}a}\text{e}^{s|\alpha |^2/2},\text{where }\alpha 𝐂\text{ and }\alpha ^{}\text{ is its complex conjugate.}$$
(3.12)
By the fundamental Baker-Campbell-Hausdorff formula (3.5), for the three discrete values of $`s=+1,0,1`$, the operator $`D(\alpha ,s)`$ can be written as an exponential which is, respectively, normal order

$$D(\alpha ,1)=\text{e}^{\alpha a^{}}\text{e}^{\alpha ^{}a},$$

Weyl order

$$D(\alpha ,0)=\text{e}^{\alpha a^{}\alpha ^{}a},$$

and anti-normal order

$$D(\alpha ,1)=\text{e}^{\alpha ^{}a}\text{e}^{\alpha a^{}}.$$

They defined the $`s`$-ordered product $`\{(a^{})^na^m\}_s`$ by means of the Taylor series

$$D(\alpha ,s)=\underset{n,m=0}{\overset{\mathrm{}}{}}\{(a^{})^na^m\}_s\frac{\alpha ^n(\alpha ^{})^m}{n!m!}$$

or, equivalently,

$$\{(a^{})^na^m\}_s\frac{^{n+m}D(\alpha ,s)}{\alpha ^n(\alpha ^{})^m}|_{\alpha =0}.$$

The intertwining formula between $`s`$-ordered product and $`t`$-ordered product is

$$D(\alpha ,s)=\text{e}^{(st)|\alpha |^2/2}D(\alpha ,t)$$

or

$$\{(a^{})^na^m\}_s=\underset{k=0}{\overset{\mathrm{min}\{n,m\}}{}}k!\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}m\\ k\end{array}\right)\left(\frac{ts}{2}\right)^k\{(a^{})^{nk}a^{mk}\}_t.$$
(3.13)
If we put $`(s,t)=(0,1),(0,1)`$, we recover (3.8), (3.10) in case of $`ϵ=1`$ respectively.

We remark that since they did not write in the form of (3.3) and (3.4) in ,  explicitly, though the idea of the proof was almost the same, we gave the proof of them in this section.

## 4 A Relation Between the Noncommutativity of Harmonic Oscillators and the Difference Theory

In this section, we describe a relation between the noncommutativity of harmonic oscillators and the difference theory. The noncommutativity of harmonic oscillators is as follows;

$$\text{e}^{t\alpha \beta }\text{e}^{\alpha a^{}}\text{e}^{\beta a}=\text{e}^{s\alpha \beta }\text{e}^{\beta a}\text{e}^{\alpha a^{}},$$
(4.1)
where $`ts=ϵ`$. If we expand both sides of (4.1), then we have

$$\underset{k,l,m}{}\frac{t^k\alpha ^k\beta ^k}{k!}\frac{\alpha ^l(a^{})^l}{l!}\frac{\beta ^ma^m}{m!}=\underset{k,l,m}{}\frac{s^k\alpha ^k\beta ^k}{k!}\frac{\beta ^la^l}{l!}\frac{\alpha ^m(a^{})^m}{m!}.$$
(4.2)
We choose coefficients of $`\alpha ^n\beta ^n`$ of (4.2), then the left hand side is

$`{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{t^k}{k!}}{\displaystyle \frac{(a^{})^{nk}}{(nk)!}}{\displaystyle \frac{a^{nk}}{(nk)!}}`$ $`=`$ $`{\displaystyle \frac{1}{(n!)^2}}{\displaystyle \underset{k=0}{\overset{n}{}}}k!t^k\left(\begin{array}{c}n\\ k\end{array}\right)^2(a^{})^{nk}a^{nk}`$
$`=`$ $`{\displaystyle \frac{1}{(n!)^2}}{\displaystyle \underset{k=0}{\overset{n}{}}}k!t^k\left(\begin{array}{c}n\\ k\end{array}\right)^2N^{\underset{¯}{nk}}.`$
By a similar calculation of the right hand side, we obtain

$$\underset{k=0}{\overset{n}{}}k!t^k\left(\begin{array}{c}n\\ k\end{array}\right)^2N^{\underset{¯}{nk}}=\underset{k=0}{\overset{n}{}}k!s^k\left(\begin{array}{c}n\\ k\end{array}\right)^2(N+ϵ)^{\overline{nk}}.$$
(4.3)
###### Remark 4.1.

In the case of $`t={\displaystyle \frac{ϵ}{2}}`$ and $`s={\displaystyle \frac{ϵ}{2}}`$, (4.3) is the equation that (3.9) and (3.11) are equal. Moreover, since $`(a^{})^na^n=N^{\underset{¯}{n}}`$ and

$$_a^{}(a^{})^na^n=n(a^{})^{n1}a^n=nN^{\underset{¯}{n1}}a,\mathrm{\Delta }_+N^{\underset{¯}{n}}=nN^{\underset{¯}{n1}},(\mathrm{\Delta }_+=\mathrm{\Delta }_{+N})$$
(4.4)
relations using various orderings of $`(a^{})^{n1}a^n={\displaystyle \frac{_a^{}}{n}}(a^{})^na^n`$ are equivalent to the relation (4.3) acted by $`\mathrm{\Delta }_+`$ with some $`t`$. Similarly, for any $`n,m`$, relations using various orderings of $`(a^{})^na^m`$ are equivalent to the relation (4.3) acted by some power of $`\mathrm{\Delta }_+`$ and some $`t`$.

###### Remark 4.2.

From (4.3),

$$\frac{d}{dt}\underset{k=0}{\overset{n}{}}k!t^k\left(\begin{array}{c}n\\ k\end{array}\right)^2N^{\underset{¯}{nk}}=\frac{d}{dt}\underset{k=0}{\overset{n}{}}k!(tϵ)^k\left(\begin{array}{c}n\\ k\end{array}\right)^2(N+ϵ)^{\overline{nk}}$$

$$\underset{k=0}{\overset{n}{}}k!kt^{k1}\left(\begin{array}{c}n\\ k\end{array}\right)^2N^{\underset{¯}{nk}}=\underset{k=0}{\overset{n}{}}k!k(tϵ)^{k1}\left(\begin{array}{c}n\\ k\end{array}\right)^2(N+ϵ)^{\overline{nk}}.$$

Therefore, if we put $`t={\displaystyle \frac{1}{2}}`$ and $`ϵ=1`$, we obtain the following equation

$$\underset{k=0}{\overset{n}{}}\frac{k!k}{2^{k1}}\left(\begin{array}{c}n\\ k\end{array}\right)^2\left\{N^{\underset{¯}{nk}}+(1)^k(N+1)^{\overline{nk}}\right\}=0.$$
(4.5)
Next we show the equation of noncommutativity of harmonic oscillators (4.3) by using the difference theory. We put

$$p_n(N)=\underset{k=0}{\overset{n}{}}k!s^k\left(\begin{array}{c}n\\ k\end{array}\right)^2(N+ϵ)^{\overline{nk}}=\underset{k=0}{\overset{n}{}}k!s^k\left(\begin{array}{c}n\\ k\end{array}\right)^2(N+(nk)ϵ)^{\underset{¯}{nk}},$$
(4.6)
and we expand (4.6) using the forward difference operator $`\mathrm{\Delta }_+`$. Since

$`\mathrm{\Delta }_+^{nl}p_n(0)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{n}{}}}k!s^k\left(\begin{array}{c}n\\ k\end{array}\right)^2(nk)(nk1)\mathrm{}(lk+1)(0+(nk)ϵ)^{\underset{¯}{lk}}`$
$`=`$ $`{\displaystyle \underset{k=0}{\overset{n}{}}}k!s^k\left(\begin{array}{c}n\\ k\end{array}\right){\displaystyle \frac{n!}{k!(nk)!}}(nk)(nk1)\mathrm{}(lk+1){\displaystyle \frac{(nk)!}{(nl)!}}ϵ^{lk}`$
$`=`$ $`l!\left(\begin{array}{c}n\\ l\end{array}\right){\displaystyle \underset{k=0}{\overset{n}{}}}s^k\left(\begin{array}{c}n\\ k\end{array}\right)\left({\displaystyle \frac{d}{dϵ}}\right)^{nl}ϵ^{nk}`$
$`=`$ $`l!\left(\begin{array}{c}n\\ l\end{array}\right)\left({\displaystyle \frac{d}{dϵ}}\right)^{nl}(s+ϵ)^n`$
$`=`$ $`l!\left(\begin{array}{c}n\\ l\end{array}\right){\displaystyle \frac{n!}{l!}}t^l,`$
by (2.4), we have

$`p_n(N)`$ $`=`$ $`{\displaystyle \underset{l=0}{\overset{n}{}}}\mathrm{\Delta }_+^{nl}p_n(0){\displaystyle \frac{N^{\underset{¯}{nl}}}{(nl)!}}`$
$`=`$ $`{\displaystyle \underset{l=0}{\overset{n}{}}}l!\left(\begin{array}{c}n\\ l\end{array}\right){\displaystyle \frac{n!}{l!}}t^l{\displaystyle \frac{N^{\underset{¯}{nl}}}{(nl)!}}`$
$`=`$ $`{\displaystyle \underset{l=0}{\overset{n}{}}}l!t^l\left(\begin{array}{c}n\\ l\end{array}\right)^2N^{\underset{¯}{nl}}.`$
Therefore, the noncommutativity of harmonic oscillators is nothing but the the Newton expansion w.r.t. variable $`N`$. More precisely, the noncommutative parameter $`ϵ`$ in $`[a,a^{}]=ϵ`$ corresponds to the increment of the difference operator $`\mathrm{\Delta }_+=\mathrm{\Delta }_{+N}`$.

## 5 Relations Including the Stirling Number of the First Kind

We shall describe nontrivial relations including the Stirling number of the first kind. These relations are derived from the noncommutativity of harmonic oscillators. Hereafter, we put $`ϵ=1`$ again.

###### Theorem 5.1.

We have

$$\underset{k=0}{\overset{2j+1}{}}\frac{k!}{2^k}\left(\begin{array}{c}n1\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)s(nk,n2j1)=0(j=0,1,\mathrm{}).$$
(5.1)
For example,

$`{\displaystyle \underset{k=0}{\overset{1}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n1\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)s(nk,n1)`$ $`=`$ $`s(n,n1)+{\displaystyle \frac{(n1)n}{2}}s(n1,n1)`$
$`=`$ $`{\displaystyle \frac{n(n1)}{2}}+{\displaystyle \frac{(n1)n}{2}}`$
$`=`$ $`0.`$
This theorem is used in a new symmetric expression of $`:(a^{}a)^n:_W`$ with respect to $`a^{}a=N`$ and $`aa^{}=N+1`$.

###### Proof.

From (4.3),

$$\underset{k=0}{\overset{n}{}}k!t^k\left(\begin{array}{c}n\\ k\end{array}\right)^2\mathrm{\Delta }_+N^{\underset{¯}{nk}}=\underset{k=0}{\overset{n}{}}k!(t1)^k\left(\begin{array}{c}n\\ k\end{array}\right)^2\mathrm{\Delta }_+(N+nk)^{\underset{¯}{nk}}$$

$$\underset{k=0}{\overset{n}{}}k!t^k\left(\begin{array}{c}n\\ k\end{array}\right)^2(nk)N^{\underset{¯}{nk1}}=\underset{k=0}{\overset{n}{}}k!(t1)^k\left(\begin{array}{c}n\\ k\end{array}\right)^2(nk)(N+nk)^{\underset{¯}{nk1}}.$$

Therefore, if we put $`t={\displaystyle \frac{1}{2}}`$ , we obtain the following equation.

$$\underset{k=0}{\overset{n}{}}\frac{k!(nk)}{2^k}\left(\begin{array}{c}n\\ k\end{array}\right)^2\left\{N^{\underset{¯}{nk1}}(1)^k(N+nk)^{\underset{¯}{nk1}}\right\}=0.$$
(5.2)
If we multiply (5.2) by $`(N+1)`$, then $`(N+nk)^{\underset{¯}{nk1}}(N+1)=(N+1)^{\overline{nk}}`$ and the left hand side is

$`{\displaystyle \underset{k0}{}}{\displaystyle \frac{k!}{2^k}}n\left(\begin{array}{c}n1\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)\left\{(N+1)^{\underset{¯}{nk}}+(1)^{k+1}(N+1)^{\overline{nk}}\right\}`$
$`=`$ $`n{\displaystyle \underset{k0}{}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n1\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)\left\{{\displaystyle \underset{i=1}{\overset{nk}{}}}s(nk,i)(N+1)^i+(1)^{k+1}{\displaystyle \underset{i=1}{\overset{nk}{}}}(1)^{nki}s(nk,i)(N+1)^i\right\}`$
$`=`$ $`n{\displaystyle \underset{i1}{}}{\displaystyle \underset{k0}{}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n1\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)s(nk,i)\left\{1+(1)^{n1i}\right\}(N+1)^i`$
$`=`$ $`2n{\displaystyle \underset{j0}{}}{\displaystyle \underset{k0}{}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n1\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)s(nk,n2j1)(N+1)^{n2j1}.`$
Therefore we obtain the theorem 5.1. ∎

###### Remark 5.2.

More generally, by the equation

$$\text{e}^{\frac{1}{2}_a^{}_a}:(a^{})^{n+m}a^n:_N=:(a^{})^{n+m}a^n:_W=\text{e}^{\frac{1}{2}_a^{}_a}:(a^{})^{n+m}a^n:_{AN},$$

we obtain the next theorem.

###### Theorem 5.3.

For $`n>0`$ and $`m=n+1,n+2,\mathrm{},1,0`$, we have

$`\{(1)^{n+mi}+1\}{\displaystyle \underset{k=0}{\overset{n+m+1i}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n+m\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)s(n+m+1k,i)`$
$`+`$ $`{\displaystyle \underset{l=i+1}{\overset{n+m+1}{}}}(m+1)^{li}\left(\begin{array}{c}l1\\ i1\end{array}\right){\displaystyle \underset{k=0}{\overset{n+m+1l}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n+m\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)s(n+m+1k,l)=0`$
$`(\text{for}i=1,2,\mathrm{},n+m+1),`$
and for $`n>0`$ and $`m=0,1,2,\mathrm{}`$, we have

$`\{(1)^{ni}+1\}{\displaystyle \underset{k=0}{\overset{n+1i}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n+m\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)s(n+1k,i)`$
$`+`$ $`{\displaystyle \underset{l=i+1}{\overset{n+1}{}}}(1m)^{li}\left(\begin{array}{c}l1\\ i1\end{array}\right){\displaystyle \underset{k=0}{\overset{n+1l}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n+m\\ k\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right)s(n+1k,l)=0`$
$`(\text{for}i=1,2,\mathrm{},n+1).`$
## 6 A Symmetric Expression of Weyl Ordering

We have expressed $`:(a^{}a)^n:_W`$ as a polynomial of the number operator $`N`$ in two ways, namely normal ordering expression (3.9) and anti-normal one (3.11). Moreover, we have a new symmetric expression of $`:(a^{}a)^n:_W`$ with respect to $`a^{}a=N`$ and $`aa^{}=N+1`$.

###### Theorem 6.1.

$$:(a^{}a)^n:_W=\underset{j=0}{\overset{[(n1)/2]}{}}\alpha (n,2j)\{N^{n2j}+(N+1)^{n2j}\},$$
(6.1)
where the constants $`\alpha (n,i)`$ are defined by

$$\alpha (n,i)\frac{1}{2}\underset{k=0}{\overset{i}{}}\frac{k!}{2^k}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}n1\\ k\end{array}\right)s(nk,ni)(i=1,2,\mathrm{}),$$

and we have $`\alpha (n,2j+1)=0(j=0,1,\mathrm{})`$.

For example,

$$\alpha (n,0)=\frac{1}{2},$$
(6.2)
and by using (2.8), (2.10) and (2.11), we have easily

$`\alpha (n,2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=0}{\overset{2}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}n1\\ k\end{array}\right)s(nk,n2)`$ (6.3)
$`=`$ $`{\displaystyle \frac{1}{2}}\left\{s(n,n2)+{\displaystyle \frac{1}{2}}n(n1)s(n1,n2)+{\displaystyle \frac{1}{8}}n(n1)^2(n2)s(n2,n2)\right\}`$
$`=`$ $`{\displaystyle \frac{1}{24}}n(n1)(n2).`$
We remark that if we put $`n=1,2,3`$ in (6.2), (6.3), we recover (1.7), (1.8) and (1.9).

By a similar calculation, we have

$`\alpha (n,4)={\displaystyle \frac{1}{2880}}n(n1)(n2)(n3)(n4)(5n7),`$
$`\alpha (n,6)={\displaystyle \frac{1}{725760}}n(n1)(n2)(n3)(n4)(n5)(n6)(35n^2147n+124),`$
where we used Mathematica.

###### Proof.

First, by (4.5), we remark that we have

$`{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}n1\\ k1\end{array}\right)\left\{N^{\underset{¯}{nk}}+(1)^k(N+1)^{\overline{nk}}\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2n}}{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{k!k}{2^{k1}}}\left(\begin{array}{c}n\\ k\end{array}\right)^2\left\{N^{\underset{¯}{nk}}+(1)^k(N+1)^{\overline{nk}}\right\}`$
$`=`$ $`0.`$
By the proposition 3.2 and the corollary 3.3, 3.4,

$`:(a^{}a)^n:_W`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\text{e}^{\frac{1}{2}_a^{}_a}:(a^{}a)^n:_N+\text{e}^{\frac{1}{2}_a^{}_a}:(a^{}a)^n:_{AN})`$
$`=`$ $`{\displaystyle \frac{1}{2}}\left\{{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)^2N^{\underset{¯}{nk}}+{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{(1)^kk!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)^2(N+1)^{\overline{nk}}\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)\left\{\left(\begin{array}{c}n1\\ k\end{array}\right)+\left(\begin{array}{c}n1\\ k1\end{array}\right)\right\}\left\{N^{\underset{¯}{nk}}+(1)^k(N+1)^{\overline{nk}}\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=0}{\overset{n1}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}n1\\ k\end{array}\right)\left\{N^{\underset{¯}{nk}}+(1)^k(N+1)^{\overline{nk}}\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=0}{\overset{n1}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}n1\\ k\end{array}\right)\{{\displaystyle \underset{i1}{}}s(nk,i)N^i`$
$`+(1)^k{\displaystyle \underset{i1}{}}(1)^{nki}s(nk,i)(N+1)^i\}`$
$`=`$ $`{\displaystyle \underset{i1}{}}\{{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=0}{\overset{ni}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}n1\\ k\end{array}\right)s(nk,i)N^i`$
$`+(1)^{ni}{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=0}{\overset{ni}{}}}{\displaystyle \frac{k!}{2^k}}\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}n1\\ k\end{array}\right)s(nk,i)(N+1)^i\}`$
$`=`$ $`{\displaystyle \underset{i1}{}}\alpha (n,ni)\{N^i+(1)^{ni}(N+1)^i\},`$
where we used $`s(nk,i)=0(nk<i)`$.
If $`i=n2j1(j=0,1,\mathrm{})`$, $`\alpha (n,2j+1)`$ vanish because of (5.1).
And if $`i=n2j(j=0,1,\mathrm{})`$, since $`(1)^{ni}=1`$, we obtain the theorem 6.1. ∎

## 7 Discussion

We obtained a new symmetric expression of $`:(a^{}a)^n:_W`$ w.r.t. $`a^{}a=N`$ and $`aa^{}=N+1`$ by using the difference theory. We also found that the noncommutative parameter $`ϵ`$ corresponded to the increment of the difference operator w.r.t. variable $`N`$. Therefore, quantum (noncommutative) calculations of harmonic oscillators were done by classical (commutative) ones of the number operator by using the difference theory. As a by-product, nontrivial relations including the Stirling number of the first kind (5.3), (5.3) were also obtained.

## Acknowledgements

We are very grateful to Shin’ichi Nojiri for valuable discussions.