text:
# 1 Introduction

## 1 Introduction

Intensive experimental studies of $`B\pi \pi `$ and other charmless two-body hadronic $`B`$ decays are currently under way at $`B`$ factories . The major goal is to complement the already observed CP-violation in the $`BJ/\psi K`$ channel by measuring other CP-violating effects and quantifying them in terms of the angles $`\alpha `$ and $`\gamma `$ of the CKM unitarity triangle (for a recent comprehensive review, see, e.g. Ref. ). The extraction of these fundamental parameters from experimental data on $`B\pi \pi ,K\pi ,K\overline{K}`$ demands reliable theoretical estimates of hadronic matrix elements of the effective weak Hamiltonian between the initial $`B`$ and the final two-light-meson states. These matrix elements are not yet accessible in lattice QCD, and thus other QCD methods of their calculation are actively being developed.

The factorization ansatz in hadronic two-body $`B`$ decays has recently been put on a more solid ground. It has been shown , to $`O(\alpha _s)`$, that in the limit $`m_b\mathrm{}`$ the amplitude for $`B\pi \pi `$ factorizes into a product of the $`B\pi `$ form factor and the pion decay constant. The nonfactorizable effects due to hard gluon exchanges are directly calculated combining perturbative QCD with certain nonperturbative inputs. The latter include the $`B\pi `$ form factor and the light-cone distribution amplitudes (DA) of pions and $`B`$ meson. According to QCD factorization, various soft nonfactorizable contributions to $`B\pi \pi `$, such as exchanges of low virtuality gluons between the “emission” pion and the remaining $`B\pi `$ part of the process, vanish in the $`m_b\mathrm{}`$ limit. An open phenomenological problem for QCD factorization is to obtain quantitative estimates of soft nonfactorizable effects at a finite, physical $`b`$ quark mass. Among these effects, the contributions of current-current operators in the penguin topology (e.g., the “charming penguins”) could be important, as was argued in Ref. .

Recently, a new method to calculate the $`B\pi \pi `$ hadronic matrix elements from QCD light-cone sum rules (LCSR)  has been suggested by one of us . The main advantage of LCSR is that the hadronic matrix elements of heavy-to-light transitions can be calculated including simultaneously soft and hard effects. Thus, the $`B\pi `$ form factor including both soft (end-point) contributions and hard-gluon exchanges was obtained from LCSR  providing the main input for the factorizable $`B\pi \pi `$ amplitude. Furthermore, in Ref.  it was shown that with LCSR one achieves a quantitative control over $`1/m_b`$ effects of soft nonfactorizable gluons in $`B\pi \pi `$. The soft-gluon effect in the matrix element of the current-current operator in the emission topology turns out to be of the same size as the hard-gluon nonfactorizable contributions obtained from QCD factorization. It is important therefore to investigate other nonfactorizable effects, related to the penguin and dipole operators and/or to non-emission topologies for the current-current operators. Especially interesting is to compare the size of soft- and hard-gluon contributions.

To continue a systematic study in this direction we apply in this paper the LCSR method to the $`B\pi \pi `$ hadronic matrix element of the chromomagnetic dipole operator $`O_{8g}`$ (gluonic penguin). Apart from being phenomenologically interesting by itself, the $`\pi \pi |O_{8g}|B`$ matrix element provides a very useful study case for the LCSR approach. Indeed, as we shall see, the hard-gluon contribution is given by one-loop diagrams and is therefore relatively simple. The soft-gluon effect corresponding to a tree-level diagram is also easily calculable. One is therefore able to estimate both hard and soft gluonic-penguin contributions using one and the same method and input.

The paper is organized as follows. In section 2 we introduce the relevant three-point correlation function and calculate it using operator-product expansion (OPE) near the light-cone and taking into account both hard- ($`O(\alpha _s)`$) and soft-gluon contributions. In section 3, following the procedure described in Ref.  we match the result of this calculation to the dispersion relations in $`\pi `$ and $`B`$ channels combined with quark-hadron duality, and obtain the sum rule for the $`\pi \pi |O_{8g}|B`$ matrix element. In section 4 we investigate the $`m_b\mathrm{}`$ limit of the sum rule. Section 5 contains the numerical analysis. Furthermore, in section 6 our prediction is discussed from the phenomenological point of view and compared, in section 7, with the corresponding result of QCD factorization. We conclude in section 8. The appendix contains some useful formulae: the decompositions of the vacuum-to-pion matrix elements in terms of the light-cone DA of various twist.

## 2 Correlation function

Following Ref.  we start with defining a vacuum-pion correlation function:

$$F_\alpha (p,q,k)=d^4xe^{i(pq)x}d^4ye^{i(pk)y}\mathrm{\hspace{0.17em}0}T\{j_{\alpha 5}^{(\pi )}(y)O_{8g}(0)j_5^{(B)}(x)\}\pi ^{}(q),$$
(1)
where the quark currents $`j_{\alpha 5}^{(\pi )}=\overline{u}\gamma _\alpha \gamma _5d`$ and $`j_5^{(B)}=m_b\overline{b}i\gamma _5d`$ interpolate $`\pi `$ and $`B`$ mesons, respectively, and the quark-gluon chromomagnetic dipole operator $`O_{8g}`$ relevant for $`B\pi \pi `$ has the following standard expression:

$$O_{8g}=\frac{m_b}{8\pi ^2}\overline{d}\sigma ^{\mu \nu }(1+\gamma _5)\frac{\lambda ^a}{2}g_sG_{\mu \nu }^ab.$$
(2)
This operator enters the $`\mathrm{\Delta }B=1`$ effective weak Hamiltonian

$$H_W=\frac{G_F}{\sqrt{2}}\left\{\lambda _u\left[\left(c_1(\mu )+\frac{c_2(\mu )}{3}\right)O_1(\mu )+2c_2(\mu )\stackrel{~}{O}_1(\mu )\right]+\mathrm{}+\lambda _tc_{8g}(\mu )O_{8g}(\mu )\right\},$$
(3)
where $`\lambda _q=V_{qb}V_{qd}^{}`$, $`q=u,t`$ are the combinations of CKM parameters, and we have singled out the relevant current-current operators:

$$O_1=(\overline{d}\mathrm{\Gamma }_\mu u)(\overline{u}\mathrm{\Gamma }^\mu b),\stackrel{~}{O}_1=(\overline{d}\mathrm{\Gamma }_\mu \frac{\lambda ^a}{2}u)(\overline{u}\mathrm{\Gamma }^\mu \frac{\lambda ^a}{2}b),$$
(4)
denoting by ellipses all remaining operators which we do not consider in this paper. In Eqs. (2)-(4), $`\text{Tr}(\lambda ^a\lambda ^b)=2\delta ^{ab}`$, $`c_{1,2,8g}`$ are the Wilson coefficients, $`\mu m_b`$ is the normalization scale and $`\mathrm{\Gamma }_\mu =\gamma _\mu (1\gamma _5)`$. Furthermore, the operator $`O_2=(\overline{u}\mathrm{\Gamma }_\mu u)(\overline{d}\mathrm{\Gamma }^\mu b)`$ has been Fierz transformed: $`O_2=O_1/3+2\stackrel{~}{O}_1`$.

The correlator (1) is a function of three independent momenta, which are chosen to be $`q`$, $`pk`$, $`k`$, and, for simplicity, $`p^2=k^2=0`$. Similar to Ref. , we introduce the unphysical external momentum $`k`$ flowing into the $`O_{8g}`$ vertex. At $`k0`$, the momentum $`pq`$ in the $`B`$ channel is independent of the total momentum $`P=pqk`$ of the light-quark state formed after the $`b`$ quark decay. That will allow us later to avoid contributions of the “parasitic” light-quark states in the dispersion relation in the variable $`(pq)^2`$. In the final sum rule $`k`$ vanishes in the ground-state $`B\pi \pi `$ contribution. In Eq.(1) the pion is on shell, $`q^2=m_\pi ^2`$. We will work in the chiral limit and set $`m_\pi =0`$ everywhere, except in the enhanced combination $`\mu _\pi =m_\pi ^2/(m_u+m_d)`$. The Lorentz-decomposition of the correlation function (1) contains four invariant amplitudes:

$$F_\alpha =(pk)_\alpha F+q_\alpha \stackrel{~}{F}_1+k_\alpha \stackrel{~}{F}_2+ϵ_{\alpha \beta \lambda \rho }q^\beta p^\lambda k^\rho \stackrel{~}{F}_3,$$
(5)
depending on three kinematical invariants $`(pk)^2`$, $`(pq)^2`$ and $`P^2=(pqk)^2`$. In what follows only the amplitude $`F`$ is needed.

In the region of large virtualities: $`(pk)^2,(pq)^2,P^2<0`$, $`|(pk)^2|,|(pq)^2|,|P^2|\mathrm{\Lambda }_{QCD}^2`$, the correlation function (1) is calculated using light-cone OPE, in a form of the sum of short-distance coefficient functions (hard amplitudes) convoluted with the pion DA of growing twist and multiplicity. Only a few first terms of this expansion are usually retained, the components with higher twist/multiplicity are neglected being suppressed by inverse powers of large virtualities.

The gluon emitted from the vertex $`O_{8g}`$ in the correlation function contributes in two different ways. First, a highly virtual (hard) gluon is absorbed by one of the quark lines at a short light-cone separation from the emission point. The corresponding one-loop diagrams (see Fig. 1) are then calculated perturbatively, as an $`O(\alpha _s)`$ part of the hard amplitude. The second possibility is that the emitted gluon together with the quark-antiquark pair forms the quark-antiquark-gluon DA of the pion (see Fig. 2a). We will call these low-virtuality gluons soft. Thus, in the correlation function the hard- and soft-gluon contributions belong to different terms of the light-cone OPE and can be calculated separately. In what follows, we take into account all $`O(\alpha _s)`$ hard-gluon effects of twist 2 and 3, and, in addition the effects of quark-antiquark-gluon DA of twist 3 and 4 in zeroth order in $`\alpha _s`$. This approximation corresponds to the accuracy of the light-cone expansion adopted for the LCSR calculation of the $`B\pi `$ form factor.

The higher twist terms corresponding to the soft-gluon contributions are subleading in the light-cone OPE, being suppressed by inverse powers of large virtualities. However, in the correlation function (1) the lowest twist contribution is of $`O(\alpha _s)`$, therefore both hard and soft contributions can be equally important. We may again refer to the LCSR for the $`B\pi `$ form factor where the $`O(\alpha _s)`$ twist 2 and the zeroth-order in $`\alpha _s`$ twist 4 contributions are at the same level numerically.

In the following three subsections the calculation of the various contributions to the correlation function (1) is presented.

### 2.1 Hard gluon contribution

We begin with calculating the leading $`O(\alpha _s)`$ contribution to the correlator (1) generated by hard gluons. To obtain the relevant diagrams one has to contract the $`b`$-quark fields and, in addition, the $`d`$\- and $`\overline{d}`$quark fields from the operators $`j_{\alpha 5}^{(\pi )}`$ and $`O_{8g}`$, respectively, inserting the free-quark propagators for both contractions.

The remaining on-shell $`\overline{u}`$ and $`d`$ fields form the quark-antiquark DA of the pion. The decomposition of the relevant vacuum-pion matrix elements in terms of DA of the lowest twists 2 and 3 is presented in the appendix. Twist 4 components are neglected, because the $`O(\alpha _s)`$ twist-4 effects are beyond the accuracy of our calculation. The gluon emitted by $`O_{8g}`$ is absorbed by one of the four quark lines leading to the four diagrams shown in Fig. 1. Note that we disregard the possibility to contract the $`\overline{d}`$ quark from $`O_{8g}`$ with the $`d`$ quark from $`j_5^{(B)}`$. That type of contraction leads to “penguin-annihilation” diagrams where the heavy-light loop is connected by a single gluon with the pion part of the correlation function. An additional gluon has to be added in this case to obey color conservation, making this contribution either $`O(\alpha _s^2)`$, with both gluons hard, or $`O(\alpha _s)`$, with one hard gluon and one soft gluon, the latter entering the pion DA. Both effects are beyond the adopted approximation for the correlation function and presumably very small.

The one-loop diagrams are calculated employing dimensional regularization. Among other tools we used the FORM program . The invariant amplitude $`F`$ is obtained by taking the coefficient at $`(pk)_\alpha `$. The result for $`F`$ is ultraviolet divergent; however, this divergence does not play any role in the resulting sum rule because the $`1/ϵ`$ terms vanish after Borel transformations. On the other hand, we find that the amplitude $`F`$ does not contain infrared-collinear divergences. This is in accordance with the fact that the perturbative expansion for the correlation function starts, in twists 2 and 3, with $`O(\alpha _s)`$. If one takes into account an additional perturbative gluon correction in the correlation function, the resulting two-loop diagrams will yield infrared-collinear divergences which should be absorbed by the evolution of the pion DA, as in the case of the $`O(\alpha _s)`$ correction to the LCSR for the $`B\pi `$ form factor .

According to the procedure of the sum rule derivation explained in Ref.  we have to retain in the final answer for the amplitude $`F`$ only those terms which, in the limit of large $`|P^2|m_B^2`$, have a nonvanishing double imaginary part in the variables $`(pk)^2`$ and $`(pq)^2`$, taken in the duality region $`\mathrm{\hspace{0.33em}0}<(pk)^2<s_0^\pi `$, $`m_b^2<(pq)^2<s_0^B`$. We find that only two diagrams, in Fig. 1a,b, contain nonvanishing terms. Retaining only these terms we obtain the following analytical expression for the hard-gluon contribution to the amplitude $`F`$:

$$F^{hard}(s_1,s_2,P^2)=F^{tw2}(s_1,s_2,P^2)+F^{tw3}(s_1,s_2,P^2),$$
(6)
where the twist 2 and 3 parts are, respectively

$`F^{tw2}(s_1,s_2,P^2)={\displaystyle \frac{\alpha _sC_Fm_b^2f_\pi }{16\pi ^3}}{\displaystyle \underset{0}{\overset{1}{}}}𝑑u\phi _\pi (u)`$
$`\times \left\{{\displaystyle \frac{s_2}{m_b^2\overline{u}s_2}}\text{ln}\left({\displaystyle \frac{s_1}{\mu ^2}}\right)+{\displaystyle \frac{(s_1+s_2P^2)(m_b^2s_2)^2}{u(\overline{u}P^2us_1)s_2^2}}\text{ln}\left({\displaystyle \frac{m_b^2s_2}{m_b^2}}\right)\right\},`$ (7)
and

$`F^{tw3}(s_1,s_2,P^2)={\displaystyle \frac{\alpha _sC_Fm_b^3f_\pi \mu _\pi }{32\pi ^3}}[{\displaystyle \underset{0}{\overset{1}{}}}du\phi _p(u)\{{\displaystyle \frac{3}{m_b^2\overline{u}s_2}}\text{ln}\left({\displaystyle \frac{s_1}{\mu ^2}}\right)`$
$`+{\displaystyle \frac{(s_13s_2P^2)(m_b^2s_2)}{(\overline{u}P^2us_1)s_2^2}}\text{ln}\left({\displaystyle \frac{m_b^2s_2}{m_b^2}}\right)\}`$
$`+{\displaystyle \frac{1}{6}}{\displaystyle \underset{0}{\overset{1}{}}}du\phi _\sigma (u)\{({\displaystyle \frac{2}{\overline{u}(m_b^2\overline{u}s_2)}}+{\displaystyle \frac{s_2}{(m_b^2\overline{u}s_2)^2}})\text{ln}\left({\displaystyle \frac{s_1}{\mu ^2}}\right)`$
$`{\displaystyle \frac{(m_b^2s_2)}{s_2^2}}({\displaystyle \frac{2(s_1+s_2P^2)}{u(\overline{u}P^2us_1)}}+{\displaystyle \frac{(s_1P^2)(s_1+5s_2P^2)}{(\overline{u}P^2us_1)^2}})\text{ln}\left({\displaystyle \frac{m_b^2s_2}{m_b^2}}\right)\}].`$ (8)
In the above, we denote $`s_1=(pk)^2`$, $`s_2=(pq)^2`$ and $`\overline{u}=1u`$ for brevity. The terms proportional to $`\text{ln}(s_1/\mu ^2)`$ and $`\text{ln}((m_b^2s_2)/m_b^2)`$ originate from the diagrams in Fig. 1a and 1b, respectively. Note that the logarithmic dependence on the dimensional regularization scale $`\mu `$ is unimportant because it vanishes after Borel transformation in the variable $`s_1`$.

### 2.2 Soft gluon contribution

In order to reach the same accuracy in the light-cone expansion of the correlation function (1) as in the calculation of the $`B\pi `$ form factor one has to take into account the contributions of quark-antiquark-gluon three-particle DA. Being a regular part of the light-cone expansion, these contributions correspond to a nonperturbative effect of zeroth order in $`\alpha _s`$ which is qualitatively different from hard-gluon exchanges. The on-shell gluon field emitted by the penguin operator is directly absorbed in the three-particle pion DA yielding a single tree-level diagram shown in Fig. 2a. The corresponding expression is obtained by contracting $`b`$ and $`d`$ quark fields and using the decomposition of the $`0|\overline{u}(y)G_{\mu \nu }(0)d(x)|\pi `$ matrix element in terms of twist 3 and 4 quark-antiquark-gluon DA given in the appendix. We find that the contribution of twist 3 DA vanishes in the amplitude $`F`$ and the answer contains only twist 4 DA :

$`F^{soft}(s_1,s_2,P^2)={\displaystyle \frac{m_b^2f_\pi }{8\pi ^2}}{\displaystyle \frac{𝒟\alpha _i}{(m_b^2(1\alpha _1)s_2)(\alpha _2P^2(1\alpha _2)s_1)}}`$
$`\times \{2(s_1s_2P^2)[\phi _{}(\alpha _i)+\stackrel{~}{\phi }_{}(\alpha _i))]+(s_1+s_2P^2)[\phi _{}(\alpha _i)+\stackrel{~}{\phi }_{}(\alpha _i)]\},`$ (9)
where $`\alpha _i\{\alpha _1,\alpha _2,\alpha _3\}`$ and $`𝒟\alpha _i=d\alpha _1d\alpha _2d\alpha _3\delta \left(1\alpha _1\alpha _2\alpha _3\right)`$.

### 2.3 Factorizable 4-quark contribution

Continuing the light-cone expansion of the correlation function (1) one encounters also four-particle Fock states. In particular, four-quark ($`\overline{q}q\overline{q}q`$) states emerge when $`b`$ quark and hard gluon propagate at short light-cone separations and four light-quark fields are on-shell, yielding $`O(\alpha _s)`$ contribution. One of the corresponding diagrams is shown in Fig. 2b. The Fock states with $``$ 3 multiplicity contribute in two different ways. First, various multiparticle DA of higher twist can be formed. Usually, in LCSR applications their contributions are expected to be strongly suppressed (see, e.g. a more detailed discussion in Ref. ) and are therefore neglected, also in our analysis. Second, the configurations of multiple quark and gluon fields sandwiched between the vacuum and pion state, can in general be splitted into a product of two operators of lower multiplicity, one of them being vacuum averaged, and the other one forming a vacuum-pion matrix element. In order to avoid double counting, these factorizable contributions have to be subtracted from genuine multiparticle DA. In our case, the four-quark configurations emerging from the correlator (1) contain a factorizable piece, a product of vacuum-averaged quark-antiquark fields and a low-twist (2,3) quark-antiquark DA. We will take the four-quark factorizable contributions into account approximating the vacuum average of quark and antiquark fields with the standard quark condensate of dimension 3. Combined with the twist 2 and 3 DA, these contributions are effectively of twist 5 and 6 respectively, and are analogous to the factorizable twist 6 terms taken into account in the LCSR for the pion electromagnetic form factor derived in Ref. . Other factorizable multiparticle contributions, such as $`\overline{q}qGG`$ or $`\overline{q}q\overline{q}qG`$ which can be parameterized in terms of dimension 4 (gluon) and 5 (quark-gluon) condensates are neglected. The purpose of taking into account the quark-condensate contribution is twofold. First, it is enhanced by the numerically large parameter $`\mu _\pi `$, therefore the size of the quark-condensate contribution to LCSR serves as an upper limit for all neglected multiparticle contributions. Second, the QCD factorization answer for the gluonic penguin amplitude contains a certain counterpart of the quark-antiquark condensate term, with which our result will be compared.

Turning to the actual calculation, we find that only one diagram (Fig. 2b) provides nonvanishing contributions to the sum rule. The short-distance part of this diagram consists of $`b`$\- and $`u`$\- quark propagators. The momenta of the vacuum averaged $`\overline{d}d`$-quarks are neglected in the adopted condensate approximation and the vacuum averaging is done in a usual way: $`0|\overline{d_\xi ^i}d_\omega ^j|0=\delta ^{ij}\delta _{\omega \xi }\overline{q}q/12`$, where $`\overline{q}q`$ is the quark condensate density, $`q=u,d`$, and $`i,j`$, and $`\xi ,\omega `$ are the colour and spinor indices, respectively. The diagram in Fig. 2b yields:

$`F^{\overline{q}q}(s_1,s_2,P^2)={\displaystyle \frac{\alpha _sC_Fm_b^3\overline{q}qf_\pi }{12\pi }}{\displaystyle \underset{0}{\overset{1}{}}}{\displaystyle \frac{du}{s_1(\overline{u}P^2us_1)(m_b^2\overline{u}s_2)}}`$
$`\times \{2(s_1P^2)\phi _\pi (u)+{\displaystyle \frac{\mu _\pi }{2m_b}}[3\overline{u}(s_1s_2P^2)2s_1]\phi _p(u)`$
$`+{\displaystyle \frac{\mu _\pi }{12m_b}}[15(s_1P^2)+9s_2+{\displaystyle \frac{(s_1P^2)(4s_13\overline{u}s_2)5\overline{u}(s_1P^2)^2+6s_1s_2}{\overline{u}P^2us_1}}`$
$`+{\displaystyle \frac{5\overline{u}s_2(s_1P^2)+3\overline{u}s_2^2+2s_1s_2}{m_b^2\overline{u}s_2}}]\phi _\sigma (u)\},`$ (10)
where only the terms which contribute to the sum rule are retained. Adding together the contributions presented in Eqs. (8)-(10), we obtain our final result for the amplitude $`F`$:

$$F(s_1,s_2,P^2)=F^{hard}(s_1,s_2,P^2)+F^{soft}(s_1,s_2,P^2)+F^{\overline{q}q}(s_1,s_2,P^2)$$
(11)
which will be used to derive LCSR.

## 3 Dispersion relation and sum rule

The procedure of deriving the sum rule for the $`B\pi \pi `$ matrix element from the correlation function (1) is independent of the operator $`O_{8g}`$. Therefore we refer to the paper  where this procedure was explained in detail for a generic case. For completeness, let us summarize the main steps of this derivation:

1. The calculated correlation function (the amplitude $`F`$ in Eq. (11)) is matched to the dispersion relation in the pion channel, in the variable $`s_1=(pk)^2`$, at large $`s_1<0`$.
2. The Borel transformation, $`s_1M_1^2`$, is performed and quark-hadron duality is used to approximate the contribution of excited pseudoscalar and axial states (with an effective threshold $`s_0^\pi `$) yielding, at large spacelike $`P^2m_b^2`$, an expression for the matrix element of the operator product $`O_{8g}j_B`$ between pion states.
3. The matrix element $`\pi |O_{8g}j_B|\pi `$ is expanded in powers of the small ratio $`s_0^\pi /(P^2)`$ retaining the leading term of this expansion. In the resulting expression the analytic continuation from large spacelike $`P^2`$ to large timelike $`P^2=m_B^2`$ is performed. Importantly, in the adopted approximation for the correlation function the matrix element $`\pi \pi |O_{8g}j_B|0|_{P^2=m_B^2}`$ has no complex phase, as directly follows from the calculation of the imaginary part in $`s_1`$ of the relevant diagrams in Figs. 1,2. At $`0<s_1<s_0^\pi `$ the resulting expression remains real at $`P^2m_B^2`$. An imaginary part will appear if one adds hard gluon exchanges to the diagrams in Figs. 1,2, proceeding beyond the adopted order in $`\alpha _s`$. Physically, the imaginary part of $`\pi \pi |O_{8g}j_B|0|_{P^2=m_B^2}`$ should be identified, within the accuracy of the quark-hadron duality approximation, with the phase of the final-state rescattering of two pions. Our result indicates the smallness of this phase.
4. The expression obtained for $`\pi \pi |O_{8g}j_B|0`$ is then equated to the dispersion relation in the $`B`$-channel, in the variable $`s_2=(pq)^2`$. The auxiliary momentum $`k`$ vanishes in the ground state $`B`$-meson contribution. The latter contains the hadronic matrix element of $`O_{8g}`$ between $`B`$ and $`\pi \pi `$ states multiplied by the $`B`$-meson decay constant $`f_B`$.
5. The second Borel transformation, $`s_2M_2^2`$, is applied to the resulting relation and quark-hadron duality is employed again, with the corresponding threshold $`s_0^B`$, to approximate the contribution of the excited $`B`$ states. Finally, the desired sum rule for the $`\pi \pi |O_{8g}|B`$ matrix element is obtained.

For convenience, we present the resulting LCSR in a form of the sum of three separate contributions:

$$A^{(O_{8g})}(\overline{B}_d^0\pi ^+\pi ^{})\pi ^{}(p)\pi ^+(q)|O_{8g}|\overline{B}_d^0(pq)=A_{hard}^{(O_{8g})}+A_{soft}^{(O_{8g})}+A_{\overline{q}q}^{(O_{8g})},$$
(12)
where

$`A_{hard}^{(O_{8g})}`$ $`=`$ $`i{\displaystyle \frac{\alpha _sC_F}{2\pi }}m_b^2\left({\displaystyle \frac{1}{4\pi ^2f_\pi }}{\displaystyle \underset{0}{\overset{s_0^\pi }{}}}dse^{s/M_1^2}\right)({\displaystyle \frac{m_b^2f_\pi }{2m_B^2f_B}}{\displaystyle \underset{u_0^B}{\overset{1}{}}}{\displaystyle \frac{du}{u}}e^{m_B^2/M_2^2m_b^2/uM_2^2}`$ (13)
$`\times `$ $`[{\displaystyle \frac{\phi _\pi (u)}{u}}+{\displaystyle \frac{\mu _\pi }{2m_b}}\{3\phi _p(u)+{\displaystyle \frac{\phi _\sigma (u)}{3u}}{\displaystyle \frac{\phi _\sigma ^{}(u)}{6}}`$
$``$ $`\phi _p(1)\overline{u}(1+{\displaystyle \frac{3m_b^2}{um_B^2}})+{\displaystyle \frac{\phi _\sigma ^{}(1)\overline{u}}{6}}({\displaystyle \frac{5m_b^2}{um_B^2}}1)\}]),`$
$`A_{soft}^{(O_{8g})}`$ $`=`$ $`im_b^2\left({\displaystyle \frac{1}{4\pi ^2f_\pi }}{\displaystyle \underset{0}{\overset{s_0^\pi }{}}}dse^{s/M_1^2}\right)({\displaystyle \frac{f_\pi }{m_B^2f_B}}{\displaystyle \underset{u_0^B}{\overset{1}{}}}{\displaystyle \frac{du}{u}}e^{m_B^2/M_2^2m_b^2/uM_2^2}`$ (14)
$`\times `$ $`(1+{\displaystyle \frac{m_b^2}{um_B^2}})[\phi _{}(1u,0,u)+\stackrel{~}{\phi }_{}(1u,0,u)]),`$
$`A_{\overline{q}q}^{(O_{8g})}`$ $`=`$ $`i{\displaystyle \frac{\alpha _sC_F}{3\pi }}m_b^2\left({\displaystyle \frac{\overline{q}q}{f_\pi m_b}}\right)({\displaystyle \frac{m_b^2f_\pi }{2m_B^2f_B}}{\displaystyle \underset{u_0^B}{\overset{1}{}}}{\displaystyle \frac{du}{u}}e^{m_B^2/M_2^2m_b^2/uM_2^2}`$ (15)
$`\times `$ $`[{\displaystyle \frac{\phi _\pi (u)}{u}}+{\displaystyle \frac{\mu _\pi }{4m_b}}\{3(1+{\displaystyle \frac{m_b^2}{um_B^2}})\phi _p(u)+(5{\displaystyle \frac{3m_b^2}{um_B^2}})({\displaystyle \frac{\phi _\sigma (u)}{3u}}{\displaystyle \frac{\phi _\sigma ^{}(u)}{6}})\}]).`$
In the above, $`\phi _\sigma ^{}(u)=\phi _\sigma (u)/u`$, $`u_0^B=m_b^2/s_0^B`$ and $`f_B`$ is the $`B`$-meson decay constant defined as $`m_b0|\overline{q}i\gamma _5b|B=m_B^2f_B`$. Note that the first bracket in Eqs. (13) and (14) is approximately equal to $`f_\pi `$ if one uses the SVZ sum rule :

$$f_\pi ^2=\frac{1}{4\pi ^2}\underset{0}{\overset{s_0^\pi }{}}𝑑se^{s/M_1^2}\left(1+\frac{\alpha _s}{\pi }\right)+\frac{0|\alpha _s/\pi G_{\mu \nu }^aG^{a\mu \nu }|0}{12M^2}+\frac{176}{81M^4}\pi \alpha _s\overline{q}q^2,$$
(16)
retaining only the leading-order quark-loop term. The quark-condensate contribution to Eq. (16) of the form $`(m_u+m_d)\overline{q}q/M^2`$ is absent in the chiral limit, whereas the perturbative correction and the gluon and 4-quark condensate contributions are beyond our approximation.

The second and third lines in Eq. (13) contain the contributions of the diagrams in Fig. 1a and 1b, respectively. In the latter contribution only twist 3 DA remain, the twist 2 term proportional to the small ratio $`s_0^\pi /m_B^2`$ is neglected. For the same reason, in the soft-gluon contribution (14) only two DA, $`\phi _{}`$ and $`\stackrel{~}{\phi }_{}`$ contribute.

To minimize theoretical uncertainties we will calculate the ratio of the gluonic penguin matrix element (12) to the factorizable $`B\pi \pi `$ hadronic amplitude. The latter, as shown in Ref. , simply coincides with the matrix element of the operator $`O_1`$ in the emission topology:

$$A_E^{(O_1)}(\overline{B}_d^0\pi ^+\pi ^{})\pi ^{}(p)\pi ^+(q)|O_1|\overline{B}_d^0(pq)_E=im_B^2f_\pi f_{B\pi }^+(0),$$
(17)
if one neglects the $`O(s_0^\pi /m_B^2)`$ terms and the two-gluon nonfactorizable exchanges which are beyond the adopted approximation. In the above, $`f_{B\pi }^+(0)`$ is the LCSR result for the $`B\pi `$ form factor at zero momentum transfer. For completeness we write down this sum rule :

$`f_{B\pi }^+(0)`$ $`=`$ $`{\displaystyle \frac{f_\pi m_b^2}{2m_B^2f_B}}{\displaystyle \underset{u_0^B}{\overset{1}{}}}{\displaystyle \frac{du}{u}}e^{m_B^2/M_2^2m_b^2/uM_2^2}\left(\phi _\pi (u)+{\displaystyle \frac{\mu _\pi }{m_b}}\left[u\phi _p(u)+{\displaystyle \frac{\phi _\sigma (u)}{3}}{\displaystyle \frac{u\phi _\sigma ^{}(u)}{6}}\right]\right)`$ (18)
$`+`$ $`O(q\overline{q}\text{tw4})+O(q\overline{q}G)+O(\alpha _s),`$
where, for brevity, the subleading contributions of the twist-4 quark-antiquark DA, twist-3,4 quark-antiquark-gluon DA and the $`O(\alpha _s)`$ correction to the twist 2 term are not shown explicitly but taken into account in the numerical analysis. The complete expression for LCSR (18) can be found, e.g. in the review . The recently calculated  $`O(\alpha _s)`$ correction to the twist-3 part of the sum rule (18) is very small and we neglect it here.

## 4 Heavy quark mass limit

The sum rule (12) is derived in full QCD, with the finite $`b`$-quark mass. To obtain the heavy-quark mass limit of the matrix element $`A^{(O_{8g})}(\overline{B}_d^0\pi ^+\pi ^{})`$, we substitute in LCSR the expansions of all $`m_b`$-dependent quantities:

$$m_B=m_b+\overline{\mathrm{\Lambda }},s_0^B=m_b^2+2m_b\omega _0,M_2^2=2m_b\tau ,f_B=\widehat{f}_B/\sqrt{m_b},$$
(19)
where $`\overline{\mathrm{\Lambda }}`$, $`\omega _0`$, $`\tau `$, and $`\widehat{f}_B`$ are the parameters independent of the heavy-mass scale. One also has to take into account that, at $`m_b\mathrm{}`$, $`u_B^01\omega _0/m_b`$ and only the end-point regions of the momentum fraction $`u`$ contribute in the integrals in Eqs. (13)-(15). Using a convenient integration variable $`\rho =\overline{u}m_b`$, and taking into account the end-point behavior of DA, the following substitutions have to be done in these integrals, to the leading order in $`O(1/m_b)`$ :

$`{\displaystyle \underset{u_B^0}{\overset{1}{}}}𝑑u=(1/m_b){\displaystyle \underset{0}{\overset{2\omega _0}{}}}𝑑\rho ,\phi _{\pi ,\sigma }(u)=(\rho /m_b)\phi _{\pi ,\sigma }^{}(1),\phi _p(u)=1,`$
$`\phi _{}(1u,0,u)=\stackrel{~}{\phi }_{}(1u,0,u)(\rho /m_b)\stackrel{~}{\phi }_{}^{}(0,0,1),`$ (20)
where $`\stackrel{~}{\phi }_{}^{}(0,0,1)=[\stackrel{~}{\phi }_{}(1u,0,u)/u]_{u=1}`$. With these substitutions it is easy to reproduce, e.g. the leading $`1/m_b^{3/2}`$ behavior  of the $`B\pi `$ form factor from the sum rule (18) at $`m_b\mathrm{}`$:

$$f_{B\pi }^+(0)|_{m_b\mathrm{}}=\frac{\widehat{f}_{B\pi }^+(0)}{m_b^{3/2}}+O(1/m_b^{5/2}),$$
(21)
where

$`\widehat{f}_{B\pi }^+(0)={\displaystyle \frac{f_\pi }{2\widehat{f}_B}}e^{\overline{\mathrm{\Lambda }}/\tau }{\displaystyle \underset{0}{\overset{2\omega _0}{}}}𝑑\rho e^{\frac{\rho }{2\tau }}\left[\rho \phi _\pi ^{}(1)+\mu _\pi \left(\phi _p(1){\displaystyle \frac{1}{6}}\phi _\sigma ^{}(1)\right)\right]`$ (22)
is an effective, $`m_b`$-independent form factor. Consequently, the heavy-quark mass limit for the factorizable amplitude obtained from Eq. (17) is

$$A_E^{(O_1)}(\overline{B}_d^0\pi ^+\pi ^{})|_{m_b\mathrm{}}=i\sqrt{m_b}f_\pi \widehat{f}_{B\pi }^+(0).$$
(23)
Applying Eqs. (19) and (20) to the sum rule (12) one obtains the limiting behavior for the separate contributions:

$`A_{hard}^{(O_{8g})}|_{m_b\mathrm{}}`$ $`=`$ $`i{\displaystyle \frac{\alpha _sC_F}{2\pi }}\sqrt{m_b}\left({\displaystyle \frac{1}{4\pi ^2f_\pi }}{\displaystyle \underset{0}{\overset{s_0^\pi }{}}}𝑑se^{s/M_1^2}\right)`$ (24)
$`\times `$ $`\left({\displaystyle \frac{f_\pi }{2\widehat{f}_B}}e^{\overline{\mathrm{\Lambda }}/\tau }{\displaystyle \underset{0}{\overset{2\omega _0}{}}}𝑑\rho e^{\frac{\rho }{2\tau }}\left[\rho \phi _\pi ^{}(1)+\mu _\pi \left({\displaystyle \frac{3}{2}}\phi _p(1){\displaystyle \frac{1}{12}}\phi _\sigma ^{}(1)\right)\right]\right),`$
$`A_{soft}^{(O_{8g})}|_{m_b\mathrm{}}`$ $`=`$ $`i\sqrt{m_b}\left({\displaystyle \frac{1}{4\pi ^2f_\pi }}{\displaystyle \underset{0}{\overset{s_0^\pi }{}}}𝑑se^{s/M_1^2}\right)\left({\displaystyle \frac{4f_\pi }{\widehat{f}_Bm_b^2}}e^{\overline{\mathrm{\Lambda }}/\tau }{\displaystyle \underset{0}{\overset{2\omega _0}{}}}𝑑\rho e^{\frac{\rho }{2\tau }}\rho \stackrel{~}{\phi }_{}^{}(0,0,1)\right),`$ (25)
$`A_{\overline{q}q}^{(O_{8g})}|_{m_b\mathrm{}}`$ $`=`$ $`i{\displaystyle \frac{\alpha _sC_F}{3\pi }}\sqrt{m_b}\left({\displaystyle \frac{\overline{q}q}{f_\pi m_b}}\right)`$ (26)
$`\times `$ $`\left({\displaystyle \frac{f_\pi }{2\widehat{f}_B}}e^{\overline{\mathrm{\Lambda }}/\tau }{\displaystyle \underset{0}{\overset{2\omega _0}{}}}𝑑\rho e^{\frac{\rho }{2\tau }}\left[\rho \phi _\pi ^{}(1)+\mu _\pi \left({\displaystyle \frac{3}{2}}\phi _p(1){\displaystyle \frac{1}{12}}\phi _\sigma ^{}(1)\right)\right]\right),`$
where only the leading in $`1/m_b`$ term in each of the three contributions is retained. We see that the hard-gluon contribution (24) has only $`O(\alpha _s)`$ suppression with respect to the factorizable amplitude (23). Note that in Eq. (24) only the diagram in Fig. 1a contributes, because the twist 3 contribution of the Fig. 1b diagram in Eq. (13) has an additional $`1/m_b`$ suppression. Importantly, the soft gluon contribution (25), where the pion DA is proportional to the dimensionful parameter $`\delta _\pi ^2`$, is suppressed by two powers of $`1/m_b`$ . Finally, the quark-condensate contribution (26) is both $`O(\alpha _s)`$ and $`O(1/m_b)`$ suppressed. Thus, in the heavy-quark mass limit the sum rule (12) behaves similar to the LCSR for the $`B\pi `$ form factor: the twist 2 and 3 DA survive in the asymptotic limit, whereas the higher-twist contributions are suppressed by powers of $`1/m_b`$.

Furthermore, we use the relation between the twist 3 DA: $`\phi _\sigma ^{}(u)=6(12u)\phi _p(u)`$, which follows  from QCD equations of motion in the approximation where the quark-antiquark-gluon DA are neglected. Hence, at $`u1`$, $`\phi _p(1)=\phi _\sigma ^{}(1)/6`$, and the expression in the second line in Eq. (24) (and the same in Eq. (26)) is simply equal to $`\widehat{f}_{B\pi }^+(0)`$ given by Eq. (22). As a result, at $`m_b\mathrm{}`$, the gluonic penguin amplitude is approximately equal to the heavy-quark limit (23) of the factorizable amplitude multiplied by an $`O(\alpha _s)`$ factor:

$`A^{(O_{8g})}(\overline{B}_d^0\pi ^+\pi ^{})|_{m_b\mathrm{}}=A_{hard}^{(O_{8g})}|_{m_b\mathrm{}}=i{\displaystyle \frac{\alpha _sC_F}{2\pi }}\sqrt{m_b}\left({\displaystyle \frac{1}{4\pi ^2f_\pi }}{\displaystyle \underset{0}{\overset{s_0^\pi }{}}}𝑑se^{s/M_1^2}\right)\widehat{f}_{B\pi }^+(0)`$
$`{\displaystyle \frac{\alpha _sC_F}{2\pi }}A_E^{(O_1)}(\overline{B}_d^0\pi ^+\pi ^{})|_{m_b\mathrm{}},`$ (27)
where the approximate relation (16) is used. Note that the above factorization will be violated by $`O(1/m_b)`$ corrections to Eq. (24). The latter corrections can be easily obtained expanding Eq. (13) to the next-to-leading order in the inverse heavy quark mass. Interestingly however, the quark-condensate contribution (26), being $`1/m_b`$ suppressed with respect to Eq. (24), nevertheless reveals the same factorization property:

$$A_{\overline{q}q}^{(O_{8g})}|_{m_b\mathrm{}}=i\frac{\alpha _sC_F}{6\pi }\left(\frac{\mu _\pi }{m_b}\right)\sqrt{m_b}f_\pi \widehat{f}_{B\pi }^+(0)=\frac{\alpha _sC_F}{6\pi }\left(\frac{\mu _\pi }{m_b}\right)A_E^{(O_1)}(\overline{B}_d^0\pi ^+\pi ^{})|_{m_b\mathrm{}},$$
(28)
where the PCAC relation $`\mu _\pi =2\overline{q}q/f_\pi ^2`$ is used to replace the quark-condensate density by $`\mu _\pi `$. Naturally, the soft-gluon contribution (25) also violates the factorization relation (27). Although this contribution is of order $`1/m_b^2`$ with respect to $`A_{hard}^{(O_{8g})}|_{m_b\mathrm{}}`$ , it has no $`\alpha _s`$ suppression and therefore should be considered separately from the $`O(\alpha _s/m_b^2)`$ corrections stemming from the heavy-mass expansion of hard-gluon and quark-condensate terms in LCSR beyond $`O(1/m_b)`$. A direct evaluation of the latter corrections from Eqs. (13) and (15) will, however, be incomplete because, in deriving LCSR, we have already neglected small (calculable) $`s_0^\pi /m_B^2`$ terms of the same order.

## 5 Numerical results

Let us return to the sum rule (12) at the finite $`b`$-quark mass and use it to obtain a numerical estimate for the ratio

$$r^{(O_{8g})}(\overline{B}_d^0\pi ^+\pi ^{})=A^{(O_{8g})}(\overline{B}_d^0\pi ^+\pi ^{})/A_E^{(O_1)}(\overline{B}_d^0\pi ^+\pi ^{}),$$
(29)
which determines (up to the known Wilson coefficient $`c_{8g}`$) the gluonic-penguin correction to the factorizable $`B\pi \pi `$ decay amplitude. The ratio (29) has an important advantage of being less sensitive to the input parameters than the individual matrix elements $`A^{(O_8)}`$ and $`A_E^{(O_1)}`$. Moreover, $`r^{(O_{8g})}`$ is independent of $`f_B`$, since the latter cancels in the ratio of the sum rules (12) and (18).

The parameters for the pion channel are $`f_\pi =132`$ MeV and $`s_0^\pi =0.7`$ GeV<sup>2</sup>, the latter determined from the two-point SVZ sum rule  for $`f_\pi `$. The corresponding Borel parameter interval is $`M^2=0.5\text{-}1.5`$ GeV<sup>2</sup>, accommodating also the range used in the LCSR for the pion form factor .

The inputs for the $`B`$ channel are taken from the earlier analysis of LCSR (18) for the $`B\pi `$ form factor in Refs. . In particular, to be consistent with the latter sum rule, we identify $`m_b`$ with the one-loop pole mass and adopt the interval $`m_b=4.7\pm 0.1`$ GeV which is in accordance with the most recent average of the $`\overline{\text{MS}}`$ mass $`\overline{m}_b(\overline{m}_b)=4.24\pm 0.11`$ GeV  obtained from various model-independent determinations (including QCD sum rules for the $`\mathrm{{\rm Y}}`$ system). Fixing the $`m_b`$ interval one determines the corresponding duality threshold $`s_0^B=352`$ GeV<sup>2</sup> from the QCD sum rule for $`f_B`$ taken with $`O(\alpha _s)`$ accuracy. Furthermore we take the range $`M_2^2=8\text{-}12`$ GeV<sup>2</sup> for the second Borel parameter and $`\mu _b=\sqrt{m_B^2m_b^2}2.4`$ GeV for the related scale at which the pion DA are normalized. The same scale is adopted for the $`\alpha _s`$ normalization. For the latter we also use the two-loop running with $`\overline{\mathrm{\Lambda }}^{(4)}=280`$ MeV. However, one should have in mind that in the new sum rule (12) two independent Borel parameters $`M_1`$ and $`M_2`$ are present, and the scale $`M_1`$ 1 GeV is independent of $`m_b`$. To investigate the sensitivity to the normalization scale, we will vary it in a rather wide interval, between $`\mu _b/2`$ and $`2\mu _b`$. For the quark condensate density we take $`\overline{q}q(1\text{GeV})=(240\pm 10\text{MeV})^3`$, or, equivalently, $`\mu _\pi (1\text{GeV})=1.59\pm 0.2`$ GeV. The normalization parameter of the twist 4 DA $`\delta _\pi ^2(1\text{GeV})=0.17\pm 0.05`$ GeV<sup>2</sup> is determined from the two-point QCD sum rules  (see also Ref. ).

For our initial numerical illustration all pion light-cone DA are taken in the asymptotic form, that is, the coefficients $`a_2,f_{3\pi },\omega _{3\pi },ϵ_\pi `$ are put to zero in the expressions for DA presented in the appendix. The numerical results for the ratio $`r^{(O_{8g})}`$, together with the hard-, soft-gluon and quark-condensate contributions, calculated at the central values of the input parameters and for the asymptotic form of the pion DA are plotted in Fig. 3 as a function of Borel variables. We observe a nice stability in both $`M_1,M_2`$. The most important feature of our calculation which is clearly seen in Fig. 3 is that the soft-gluon part of the gluonic penguin amplitude, while being suppressed by two powers of $`m_b`$, numerically amounts to about 50% of the hard-gluon contribution and has an opposite sign. A relatively large magnitude of the soft-gluon term (14) can be traced back to the large numerical factor of 20 originating from the contributing twist-4 DA and compensating the suppression factor $`\delta _\pi ^2/m_b^2`$. The quark-condensate contribution, on the other hand, is subleading also numerically, not larger than 30% of the hard-gluon contribution.

To take into account deviations from the asymptotic form, we have repeated our calculation with the nonasymptotic part of the pion DA taken to the next-to-leading order in the conformal spin . In this case we adopt in the pion twist-2 DA $`a_2(1\text{GeV})=0.24\pm 0.14\pm 0.08`$ obtained  by fitting the LCSR for the pion form factor to the experimental data. The relevant parameters in the twist 3,4 DA have been estimated at the scale of $`1`$ GeV using two-point QCD sum rules: $`f_{3\pi }=0.0035`$, $`\omega _{3\pi }=2.88`$  and $`ϵ_\pi =0.5`$  . These estimates have an uncertainty of $`\pm 30\%`$. The corresponding LO anomalous dimensions can be found in the appendix.

The difference between the values of $`r^{(O_{8g})}`$ obtained with nonasymptotic and asymptotic DA is inessential for the hard-gluon and quark-condensate contributions. On the contrary, in the soft-gluon term, the effect of the nonasymptotic part determined by the parameter $`ϵ_\pi `$ is substantial. The reason is that in the soft-gluon term described by the diagram in Fig. 2a the gluon carries the dominant fraction $`\alpha _3=u`$ ($`u_0^B<u<1`$) of the pion momentum, the $`d`$-quark carries $`\alpha _1=1u`$, whereas the $`\overline{u}`$-quark momentum fraction is restricted to small values: $`\alpha _2s_0^\pi /m_B^2`$. The integration over the variables $`\alpha _i`$ in this subspace of the integration region $`_i\alpha _i=1`$ is strongly influenced by the presence of the second polynomial in the conformal expansion of the twist 4 DA $`\phi _{}(\alpha _i)`$ and $`\stackrel{~}{\phi }_{}(\alpha _i)`$. As a result, the soft-gluon contribution, being $`50\%`$ of the hard-gluon part in the case of asymptotic DA, drops to $`(10\%)\text{-}(+20\%)`$ at $`ϵ_\pi =0.35\text{-}0.65`$.

Finally, we get the following numerical estimate for the ratio of the gluonic-penguin and factorizable amplitudes including nonasymptotic effects:

$$r^{(O_{8g})}(\overline{B}_d^0\pi ^+\pi ^{})=0.035\pm 0.015.$$
(30)
In obtaining the above range, the uncertainties caused by the variation of all parameters within their allowed intervals are added linearly. The uncertainty in Eq. (30) is mainly due to the sensitivity to the nonperturbative parameters $`ϵ_\pi `$, $`\overline{q}q`$ and $`\delta _\pi ^2`$, whereas the choice of the $`b`$-quark mass, Borel parameters and $`a_2`$ has a smaller impact. There is an additional uncertainty of approximately $`\pm 20\%`$ from varying the normalization scale in the limits $`\mu _b/2\text{-}\mathrm{\hspace{0.17em}2}\mu _b`$. The LCSR prediction for the $`B\pi `$ form factor with the same input and the same conservative procedure of estimating the uncertainty is $`f_{B\pi }^+=0.28\pm 0.06`$.

One can also investigate the quality of the $`m_b\mathrm{}`$ limit for the ratio $`r^{(O_8g)}`$. Dividing Eq. (23) by Eq. (27) we find that, within uncertainties of our input, $`r^{(O_{8g})}`$ varies between 0.6 and 1.1 of its heavy quark limit. Future improvements of nonperturbative parameters determining the vacuum-pion matrix elements of twist 3,4 will allow to decrease the uncertainty of our numerical estimates, in particular, of the subleading in $`1/m_b`$ terms.

## 6 Gluonic penguins and the decay amplitude

We are now in a position to assess the role of the gluonic-penguin matrix element in $`\overline{B^0}\pi ^+\pi ^{}`$. The decay amplitude can be represented in the following form:

$`𝒜(\overline{B}_d^0\pi ^+\pi ^{})`$ $``$ $`\pi ^+\pi ^{}|H_W|\overline{B_d^0}`$ (31)
$`=`$ $`i{\displaystyle \frac{G_F}{\sqrt{2}}}f_\pi f_{B\pi }^+(0)m_B^2\{\lambda _u[c_1(\mu )+{\displaystyle \frac{c_2(\mu )}{3}}+2c_2(\mu )r_E^{(\stackrel{~}{O}_1)}(\overline{B}_d^0\pi ^+\pi ^{})]`$
$`+`$ $`\mathrm{}+\lambda _tc_{8g}(\mu )r^{(O_{8g})}(\overline{B}_d^0\pi ^+\pi ^{})\},`$
where we denote by ellipses all terms that are not relevant for the present discussion: the hadronic matrix elements of $`O_{1,2}`$ in the annihilation, penguin and penguin-annihilation topologies; the contributions generated by current-current operators with $`c`$-quarks (“charming penguins”) and, finally, all effects of quark- and electroweak-penguin operators (for a model-independent classification of all contributions and topologies see Ref. ). The second line in Eq. (31) contains the factorizable amplitude and the nonfactorizable correction due to the operator $`\stackrel{~}{O}_1`$. The latter is parametrized by the ratio

$$r_E^{(\stackrel{~}{O}_1)}(\overline{B}_d^0\pi ^+\pi ^{})A_E^{(\stackrel{~}{O}_1)}(\overline{B}_d^0\pi ^+\pi ^{})/A_E^{(O_1)}(\overline{B}_d^0\pi ^+\pi ^{}),$$
(32)
where $`A_E^{(\stackrel{~}{O}_1)}`$ is the $`\overline{B}_d^0\pi ^+\pi ^{}`$ matrix element of $`\stackrel{~}{O}_1`$ in the emission topology. Finally, in the third line in Eq. (31) we added the gluonic penguin correction calculated above.

The matrix element of $`\stackrel{~}{O}_1`$ contains the hard $`O(\alpha _s)`$ and soft $`O(1/m_b)`$ contributions, so that

$$r_E^{(\stackrel{~}{O}_1)}(\overline{B}_d^0\pi ^+\pi ^{})=r_E^{(\stackrel{~}{O}_1,hard)}(\overline{B}_d^0\pi ^+\pi ^{})+\frac{\lambda _E(\overline{B}_d^0\pi ^+\pi ^{})}{m_B},$$
(33)
where the soft contribution was calculated from LCSR in Ref.  with an estimate $`\lambda _E(\overline{B}_d^0\pi ^+\pi ^{})=0.1\pm 0.05`$ GeV. The hard nonfactorizable correction $`r_E^{(\stackrel{~}{O}_1,hard)}`$ is also accessible with LCSR but demands two-loop calculation. Following Ref.  we simply use the QCD factorization estimate: $`r_E^{(\stackrel{~}{O}_1,hard)}(\overline{B}_d^0\pi ^+\pi ^{})=\alpha _sC_FF/(8\pi N_c)`$, where the expression for $`F`$ is given in Ref. . Since we are only interested in an order-of-magnitude estimate, we neglect the twist-3 part of $`r_E^{(\stackrel{~}{O}_1,hard)}`$ which in QCD factorization  has to be regularized introducing some additional parameter. For consistency we use the same normalization scale $`\mu _b`$ and twist-2 pion DA as in LCSR. We obtain:

$$r_E^{(\stackrel{~}{O}_1)}(\overline{B}_d^0\pi ^+\pi ^{})=[(0.002\text{-}0.016)+0.045i]+[0.01\text{-}0.03],$$
(34)
where the first and second brackets contain the hard-gluon (QCD factorization) and soft-gluon (LCSR) contributions. The relatively large uncertainty of the hard-gluon part is due to the variation of $`a_2(1\text{GeV})`$ in the adopted range 0.02-0.46. Substituting the estimates (30) and (34) in Eq. (31) and using also the LCSR estimate for the $`B\pi `$ form factor we obtain:

$`𝒜(\overline{B}_d^0\pi ^+\pi ^{})=i{\displaystyle \frac{G_F}{\sqrt{2}}}f_\pi m_B^2(0.28\pm 0.05)\{\lambda _u(1.03`$
$`+[0.001\text{-}0.009+0.025i][0.005\text{-}0.015])+\mathrm{}\lambda _t[0.003\text{-}0.008]\},`$ (35)
where, for consistency, the Wilson coefficients  are taken at the scale $`\mu _b`$. We use $`c_1(\mu _b)=1.124`$ and $`c_2(\mu _b)=0.272`$ in NLO and in the NDR scheme and, correspondingly, $`c_{8g}(\mu _b)=0.166`$ in LO. In Eq. (35) the first (second) bracket in the term proportional to $`\lambda _u`$ is due to the hard (soft) parts of $`r_E^{(\stackrel{~}{O}_1)}`$. Note that some of the uncertainties in the $`B\pi `$ form factor and in the nonfactorizable contributions are correlated.

Both nonfactorizable corrections in Eq. (35) turn out to be very small, not exceeding a one-percent level. However, before one can claim that $`\overline{B}_d^0\pi ^+\pi ^{}`$ has a predominantly factorizable amplitude, it is necessary to calculate all remaining contributions, having in mind that the sum of all nonfactorizable contributions may accumulate into a noticeable effect. As both analyses of Ref.  and this paper show, in this calculation the soft-gluon effects are indispensable.

Note that the situation with the colour-suppressed $`\overline{B}_d^0\pi ^0\pi ^0`$ channel is quite different. Here the factorizable part is proportional to the small coefficient $`c_2(\mu _b)+c_1(\mu _b)/30.1`$, therefore the contribution of the $`\stackrel{~}{O}_1`$ matrix element multiplied by $`2c_1`$ is equally important. Also the gluonic-penguin correction in this channel can reach a few % of the total amplitude. One may expect that in some other penguin-dominated channels of charmless $`B`$ decays, such as $`BK\pi `$, the effects of gluonic penguins may even be more important.

## 7 Comparison with QCD factorization

It is interesting to compare our estimate of the hadronic matrix element $`\pi \pi |O_{8g}|B`$ expressed via the ratio $`r^{(O_{8g})}`$ with the same ratio derived from QCD factorization. We can only compare the $`O(\alpha _s)`$ contributions in our calculation because, in the QCD factorization scheme, the soft-gluon effects have not yet been investigated quantitatively.

To simplify the comparison, we set to zero all Wilson coefficients in $`H_W`$ except $`c_{8g}`$. In this case, the hard-scattering kernel in $`O(\alpha _s)`$ calculated within QCD factorization approach contains a single diagram with a hard gluon emitted by the $`O_{8g}`$ vertex and converted into a quark-antiquark pair. This diagram (the sixth one in Fig. 2 of the second paper in Ref. ) is simply reproduced from the one in Fig. 1a of this paper if one puts on shell the quarks emitted by the $`B`$-meson and pion currents. With the same procedure of putting on shell the quark lines applied to the Fig. 1b diagram, one gets a diagram which describes a hard-gluon exchange between the gluonic-penguin vertex and the spectator quark in the $`B`$-meson, with a subsequent hadronization of the light-quark pair into a pion pair. The long-distance part of this process can be described by a two-pion DA, an object different from the one-pion DA used before. In the framework of QCD factorization this particular diagram has not been taken into account. A convincing reason for that is that the corresponding contribution is $`1/m_b`$ suppressed, which is also confirmed by the heavy-mass expansion of LCSR discussed above. We conclude, that, apart from soft-gluon effects, at finite $`m_b`$ there are certain hard-gluon effects in LCSR, such as the contribution of the Fig. 1b diagram, which do not appear in the QCD factorization.

After this qualitative discussion we turn to the quantitative comparison with the gluonic penguin contribution in the QCD factorization. In our toy scenario with only $`c_{8g}0`$ it is easy to extract this contribution from the expression for the effective coefficients $`a_i`$ obtained in Refs.  for $`B\pi \pi `$. The only effective coefficients where $`c_{8g}`$ contributes are $`a_4`$ and $`a_6`$, so that the ratio $`r^{(O_{8g})}`$ in the QCD factorization scheme is :

$$r^{(O_{8g})}(\overline{B}_d^0\pi ^+\pi ^{})_{QCDfact.}=\frac{\alpha _sC_F}{2\pi N_c}\left(_0^1𝑑u\frac{\phi _\pi (u)}{1u}+\frac{2\mu _\pi }{m_b}\right),$$
(36)
together with the $`O(\mu _\pi /m_b)`$ correction which is the only $`1/m_b`$ effect retained in QCD factorization  because of the large numerical value of $`\mu _\pi `$. At $`m_b\mathrm{}`$ the above ratio, for the asymptotic DA $`\phi _\pi =6u(1u)`$, coincides with the heavy-quark mass limit (27) obtained from LCSR, if in the latter the approximation (16) is used. Thus, the LCSR prediction for the gluonic-penguin matrix element reproduces the QCD factorization result for the same matrix element in the limit of infinitely heavy $`b`$ quark.

The two methods, however, substantially differ in subleading terms. This difference manifests itself already at the $`O(\alpha _s/m_b)`$ level. In addition to the penguin-spectator exchange effect corresponding to the diagram in Fig. 2b and discussed above, there is a factor of 2 difference between the second, subleading term in Eq. (36) and the heavy quark mass limit (28) of the quark-condensate contribution in LCSR. Furthermore, the fact that the soft gluon corrections in LCSR are suppressed by two powers of $`1/m_b`$ is in a general agreement with the expectations of QCD factorization. However, as we already noted, this effect which is absent in the QCD factorization prediction, turns out to be essential at finite $`m_b`$.

To complete our comparison with other methods let us parenthetically note that LCSR predictions do not support the PQCD approach  to charmless $`B`$ decays. In particular, according to LCSR and contrary to PQCD, the soft contributions dominate in the $`B\pi `$ form factor.

## 8 Conclusion

In this paper, we have presented the first calculation of the hadronic $`B\pi \pi `$ matrix element of the gluonic penguin operator $`O_{8g}`$ using QCD LCSR. This matrix element is a very suitable study object for the sum rule approach because both hard- and soft-gluon contributions are calculable within one procedure using the same input. The results obtained for $`\pi ^+\pi ^{}|O_{8g}|\overline{B}_d^0`$ from LCSR clearly indicate that at finite $`m_b`$ the soft-gluon and other $`1/m_b`$ suppressed effects are important motivating further investigation in this direction, in particular, calculating the penguin matrix elements of current-current operators. Importantly, in the heavy-quark mass limit the LCSR prediction is consistent with the QCD factorization. Similar to the nonfactorizable emission, the gluonic-penguin contributions to $`B\pi \pi `$ are suppressed either by $`\alpha _s`$ or by powers of $`1/m_b`$. On the phenomenological side, our result clearly indicates that gluonic penguins play an insignificant role in the $`\overline{B}_d^0\pi ^+\pi ^{}`$ decay amplitude but become noticeable in $`\overline{B}_d^0\pi ^0\pi ^0`$. Our prediction can also be used in the future analyses of penguin-dominated $`BK\pi `$ modes.

### Acknowledgments

This work is supported by the DFG Forschergruppe “Quantenfeldtheorie, Computeralgebra und Monte Carlo Simulationen”, by the German Ministry for Education and Research (BMBF), and by the KBN grant 5P03B09320. P.U. acknowledges the support of the Humboldt Foundation.

## Appendix: Vacuum-pion matrix elements and distribution amplitudes

The following convenient decomposition of the vacuum-pion matrix element near the light-cone is used in the calculation of the correlation function (1):

$`0|\overline{u}_\xi ^i(x_2)d_\omega ^j(x_1)|\pi ^{}(q)={\displaystyle \frac{i\delta ^{ij}}{12}}f_\pi {\displaystyle \underset{0}{\overset{1}{}}}due^{iuqx_1i\overline{u}qx_2}[(\overline{)}q\gamma _5)_{\omega \xi }\phi _\pi (u)`$
$`+(\gamma _5)_{\omega \xi }\mu _\pi \phi _p(u){\displaystyle \frac{1}{6}}(\sigma _{\beta \tau }\gamma _5)_{\omega \xi }q_\beta (x_1x_2)_\tau \mu _\pi \phi _\sigma (u)]+\text{twist }\text{ 4},`$ (A1)
where the quark fields are taken near the light-cone ($`x_i=u_ix`$, $`x^2=0`$), $`i,j=1,2,3`$ and $`\xi ,\omega =1\text{-}4`$ are the quark color and spinor indices, respectively. Only the twist-2 DA $`\phi _\pi `$ and the twist 3 DA $`\phi _p`$ and $`\phi _\sigma `$ are retained. The familiar definitions of these DA are obtained by multiplying both parts of this equation by the corresponding combinations of $`\gamma `$ matrices and taking Dirac and color traces. For the twist-2 pion DA, in the next-to-leading accuracy in the conformal spin  adopted in this paper, only the first nonasymptotic term is retained in the expansion in Gegenbauer polynomials $`C_{2n}^{3/2}`$:

$$\phi _\pi (u,\mu )=6u\overline{u}\left[1+a_2(\mu )C_2^{3/2}(u\overline{u})\right].$$
(A2)
The LO scale-dependence of the nonasymptotic part is determined by

$$a_2(\mu _2)=\left[L(\mu _2,\mu _1)\right]^{\gamma _2^{(0)}/\beta _0}a_2(\mu _1),$$
(A3)
where $`L(\mu _2,\mu _1)=\alpha _s(\mu _2)/\alpha _s(\mu _1)`$, $`\beta _0=11\frac{2}{3}N_F`$, and $`\gamma _2^{(0)}=50/9`$.

The analogous decomposition for the quark-antiquark-gluon matrix element in terms of the pion three-particle DA reads :

$`0|\overline{u}_\xi ^i(x_2)g_sG_{\mu \nu }^a(x_3)d_\omega ^j(x_1)|\pi ^{}(q)={\displaystyle \frac{\lambda _{ji}^a}{32}}{\displaystyle 𝒟\alpha _ie^{iq(\alpha _1x_1+\alpha _2x_2+\alpha _3x_3)}}`$
$`\times [if_{3\pi }(\sigma _{\lambda \rho }\gamma _5)_{\omega \xi }(q_\mu q_\lambda g_{\nu \rho }q_\nu q_\lambda g_{\mu \rho })\phi _{3\pi }(\alpha _i)`$
$`f_\pi (\gamma _\lambda \gamma _5)_{\omega \xi }\left\{(q_\nu g_{\mu \lambda }q_\mu g_{\nu \lambda })\phi _{}(\alpha _i)+{\displaystyle \frac{q_\lambda (q_\mu x_\nu q_\nu x_\mu )}{(qx)}}\left(\phi _{}(\alpha _i)+\phi _{}(\alpha _i)\right)\right\}`$
$`+{\displaystyle \frac{if_\pi }{2}}ϵ_{\mu \nu \delta \rho }(\gamma _\lambda )_{\omega \xi }\{(q^\rho g^{\delta \lambda }q^\delta g^{\rho \lambda })\stackrel{~}{\phi }_{}(\alpha _i)+{\displaystyle \frac{q_\lambda (q^\delta x^\rho q^\rho x^\delta )}{(qx)}}(\stackrel{~}{\phi }_{}(\alpha _i)+\stackrel{~}{\phi }_{}(\alpha _i))\}].`$ (A4)
In the above

$$\phi _{3\pi }=360\alpha _1\alpha _2\alpha _3^2\left[1+\frac{\omega _{3\pi }}{2}(7\alpha _33)\right]$$
(A5)
is the twist-3 quark-antiquark-gluon DA, with nonperturbative parameters $`f_{3\pi }`$ and $`ϵ_\pi `$ defined via matrix elements of the following local operators:

$`0|\overline{u}\sigma _{\mu \nu }\gamma _5G_{\alpha \beta }d|\pi ^{}(q)=if_{3\pi }\left[(q_\alpha q_\mu g_{\beta \nu }q_\beta q_\mu g_{\alpha \nu })(q_\alpha q_\nu g_{\beta \mu }q_\beta q_\nu g_{\alpha \mu })\right],`$ (A6)
$$0|\overline{u}\sigma _{\mu \lambda }\gamma _5[D_\beta ,G_{\alpha \lambda }]d\frac{3}{7}_\beta \overline{u}\sigma _{\mu \lambda }\gamma _5G_{\alpha \lambda }d|\pi ^{}(q)=\frac{3}{14}f_{3\pi }\omega _{3\pi }q_\alpha q_\beta q_\mu .$$
(A7)
The scale dependence of the twist 3 parameters is given by:

$`\mu _\pi (\mu _2)=\left[L(\mu _2,\mu _1)\right]^{\frac{4}{\beta _0}}\mu _\pi (\mu _1),f_{3\pi }(\mu _2)=\left[L(\mu _2,\mu _1)\right]^{\frac{1}{\beta _0}\left(\frac{7C_F}{3}+3\right)}f_{3\pi }(\mu _1),`$ (A8)
$`(f_{3\pi }\omega _{3\pi })(\mu _2)=\left[L(\mu _2,\mu _1)\right]^{\frac{1}{\beta _0}\left(\frac{7C_F}{6}+10\right)}(f_{3\pi }\omega _{3\pi })(\mu _1),`$ (A9)
The corresponding expressions for the twist-3 quark-antiquark DA are:

$`\phi _p(u)=1+30{\displaystyle \frac{f_{3\pi }}{\mu _\pi f_\pi }}C_2^{1/2}(u\overline{u})3{\displaystyle \frac{f_{3\pi }\omega _{3\pi }}{\mu _\pi f_\pi }}C_4^{1/2}(u\overline{u}),`$
$`\phi _\sigma (u)=6u(1u)\left(1+5{\displaystyle \frac{f_{3\pi }}{\mu _\pi f_\pi }}\left(1{\displaystyle \frac{\omega _{3\pi }}{10}}\right)C_2^{3/2}(u\overline{u})\right).`$ (A10)
Finally, the expressions for the four twist-4 DA entering the decomposition (A4) are:

$`\phi _{}(\alpha _i)=120\delta _\pi ^2ϵ_\pi (\alpha _1\alpha _2)\alpha _1\alpha _2\alpha _3,`$
$`\phi _{}(\alpha _i)=30\delta _\pi ^2(\alpha _1\alpha _2)\alpha _3^2\left[{\displaystyle \frac{1}{3}}+2ϵ_\pi (12\alpha _3)\right],`$
$`\stackrel{~}{\phi }_{}(\alpha _i)=120\delta _\pi ^2\alpha _1\alpha _2\alpha _3\left[{\displaystyle \frac{1}{3}}+ϵ_\pi (13\alpha _3)\right],`$
$`\stackrel{~}{\phi }_{}(\alpha _i)=30\delta _\pi ^2\alpha _3^2(1\alpha _3)\left[{\displaystyle \frac{1}{3}}+2ϵ_\pi (12\alpha _3)\right],`$ (A11)
where the nonperturbative parameters $`\delta _\pi ^2`$ and $`ϵ_\pi `$ are defined as

$$0|\overline{u}\stackrel{~}{G}_{\alpha \mu }\gamma ^\alpha d|\pi ^{}(q)=i\delta _\pi ^2f_\pi q_\mu ,$$
(A12)
and (up to twist 5 corrections):

$$0|\overline{u}[D_\mu ,\stackrel{~}{G}_{\nu \xi }]\gamma ^\xi d\frac{4}{9}_\mu \overline{u}\stackrel{~}{G}_{\nu \xi }\gamma ^\xi d|\pi ^{}(q)=\frac{8}{21}f_\pi \delta _\pi ^2ϵ_\pi q_\mu q_\nu ,$$
(A13)
with the scale-dependence:

$$\delta _\pi ^2(\mu _2)=\left[L(\mu _2,\mu _1)\right]^{\frac{8C_F}{3\beta _0}}\delta _\pi ^2(\mu _1),(\delta _\pi ^2ϵ_\pi )(\mu _2)=\left[L(\mu _2,\mu _1)\right]^{\frac{10}{\beta _0}}(\delta _\pi ^2ϵ_\pi )(\mu _1).$$
(A14)


format: text

text:
# Numerical and Monte Carlo Bethe ansatz method: 1D Heisenberg model

## I Introduction

The study of exactly solvable models is a very important field in condensed matter physics, which began with Bethe’s solution of the isotropic Heisenberg chain HABethe31 . In general, the Bethe ansatz (BA) solution of a model has several drawbacks: it has a complex mathematical structure; the excitations are not immediately available; and most important, it does not give explicit results even for the thermodynamic quantities of the system. It was only when Yang and Yang CNYang69  presented a strategy to study the thermodynamics of BA solvable systems that the temperature dependence of quantities such as the specific heat and the magnetic susceptibility become available. The method is now designated as thermodynamic Bethe ansatz and has known many developments in the last thirty years MTakahashib . Additionally, correlation functions, such as the conductivity, can not be obtained from the BA equations alone, and combination of BA results with other methods is required for their calculation prl .

The BA ansatz method has been applied to Bose, Fermi EHLieb63 ; CNYang67 ; EHLieb68 , and spin systems HABethe31 ; ROrbach58 . It is a general feature of the BA solution, first proved by Yang and YangCNYang69  for the Bose case, that a given eigenstate of the model is characterized by a unique set of quantum numbers $`\{I_j\}`$. Further, it also can be shown that all configurations of these quantum numbers $`I_j`$ exhaust the Hilbert space of a given model. Since the energy eigenvalues are functions of the above quantum numbers, instead of using TBA and quantum Monte Carlo approaches, we can study BA solvable models in quantum number space by classical Monte Carlo. Furthermore, for a small system (but for larger systems than those available to exact diagonalization methods), it is possible to solve the BA equations for all eigenvalues. Therefore, the expectation value of an Hermitian operator in thermal equilibrium can be computed.

In this paper, we shall introduce two numerical approaches for computing thermodynamic quantities of Bethe ansatz solvable models. The methods are ilustrated with the 1D isotropic Heisenberg model, since this model is well studied in the literature. Furthermore, the study of the Heisenberg model is it self relevant, since this system predicts many properties of quasi-one-dimensional materials NMotoyama96 ; STakagi96 ; KRThurber01 . This model has been investigated by many kinds of methods. For example, the low temperature behaviors are quite well understood by a combination of Bethe ansatz RJBaxterb  and conformal field theory AABelavin84 ; JLCardyb . A strong logarithm singularity in the susceptibility at low temperature was first found by the Bethe ansatz calculation of quantum transfer matrix SEggert94  and then verified experimentally NMotoyama96 ; STakagi96 . The thermodynamics of the model has been studied by TBA MTakahashib ; MTakahashi71 ; MGaudin71 ; MTakahashi85 ; PSchlottmann85 ; MTakahashibb ; MShiroishi02  as well as by QTM MSuzuki87 ; TKoma89 ; AKlumper93 ; AKlumper98 .

The paper is organized as follows. In Sec. II, we first briefly review the BA solution of the isotropic Heisenberg model. In Sec. III and Sec. IV, we introduce the basic idea of numerical Bethe ansatz (MBA) and Monte Carlo Bethe ansatz (MCBA). In Sec.V we check the effectiveness of these two methods computing the specific heat and the magnetic susceptibility in the absence of an external magnetic field and compare our results with those obtained from the TBA. We than use our methods to study the two quantities above in the presence of an external magnetic field. A brief summary is given in Sec. VI.

## II Isotropic Heisenberg Model

Now let us first review the Bethe ansatz solution of the 1D Heisenberg chain, which can be found in the book of Takahashi MTakahashib . The Hamiltonian of the isotropic Heisenberg model is

$`=J{\displaystyle \underset{l=1}{\overset{N}{}}}\left(S_l^xS_{l+1}^x+S_l^yS_{l+1}^y+S_l^zS_{l+1}^z\right),`$ (1)
where $`N`$ is the number of sites, $`S_l^x,S_l^y,S_l^z`$ are spin $`1/2`$ operators at site $`l`$ and $`J=1,1`$ represent anti-ferromagnetic and ferromagnetic cases, respectively. The solution with periodic boundary condition $`\stackrel{}{S}_{N+1}=\stackrel{}{S}_1`$ using the string hypothesis takes the form

$`N\theta (x_\gamma ^n/n)=2\pi I_\gamma ^{(n)}+{\displaystyle \underset{m,\beta n,\gamma }{}}\mathrm{\Theta }_{nm}(x_\gamma ^nx_\beta ^m).`$ (2)
Here $`\theta (x)=2\mathrm{tan}^1(x)`$, and

$`\mathrm{\Theta }_{nm}(x)`$ $`=`$ $`\theta \left({\displaystyle \frac{x}{|nm|}}\right)+2\theta \left({\displaystyle \frac{x}{|nm|+2}}\right)+\mathrm{}+2\theta \left({\displaystyle \frac{x}{n+m2}}\right)+\theta \left({\displaystyle \frac{x}{n+m}}\right)\mathrm{For}:nm`$ (3)
$`=`$ $`2\theta \left({\displaystyle \frac{x}{2}}\right)+2\theta \left({\displaystyle \frac{x}{4}}\right)+\mathrm{}+2\theta \left({\displaystyle \frac{x}{2n2}}\right)+\theta \left({\displaystyle \frac{x}{2n}}\right)\mathrm{For}:n=m`$
and $`x_\gamma ^n`$ is the real part of the n-string which is designated by

$$x_\gamma ^{n,j}=x_\gamma ^n+i(n+12j),j=1,\mathrm{},n$$

$`I_\gamma ^n`$ is the quantum number of $`\gamma `$th n-string (note that $`n`$ and $`\gamma `$ are indices). We denote the number of the n-string by $`\alpha _n`$, thus $`\gamma =1,\mathrm{},\alpha _n`$ and the string configuration $`\{\alpha \}`$ satisfy

$`\alpha _1+2\alpha _2+\mathrm{}+(M1)\alpha _{M1}+M\alpha _M=M,`$ (4)
where $`M`$ is the number of down spins. The quantum number of n-string $`I_\gamma ^n`$ is an integer (half-odd integer) if $`N\alpha _n`$ is odd (even) and satisfy

$`|I_\alpha ^n|\left(N1{\displaystyle \underset{m=1}{\overset{M}{}}}t_{nm}\alpha _m\right)/2,`$ (5)
where $`t_{nm}2\mathrm{min}(n,m)\delta _{nm}`$. For a given set of $`\{I_\gamma ^n\}`$, Eq. (2) can be solved numerically and the energy is given by

$`E\{I_\gamma ^n\}=NJ/4+{\displaystyle \underset{n,\gamma }{}}{\displaystyle \frac{2Jn}{(x_\gamma ^n)^2+n^2}},`$ (6)
which represents the energy of the lowest weight state in SU(2) irreducible space designated by $`S=N/2M,S_z=S,S1,\mathrm{},S`$. In the presence of an external field $`h`$ a Zeeman term is added to Eq. (6). Hence the total energy of a given quantum number configuration is given by

$`E=E\{I_\gamma ^{(n)}\}h,`$ (7)
where $`=2S_z`$ is the magnetization of the state

## III Numerical Bethe Ansatz

In statistical mechanics, the expectation value of an Hermitian operator $`Q`$ in thermal equilibrium is given by

$`Q={\displaystyle \frac{1}{Z}}{\displaystyle \underset{\mu }{}}Q_\mu e^{\beta E_\mu }.`$ (8)
where $`Z`$ is known as partition function, defined as

$`Z={\displaystyle \underset{\mu }{}}e^{\beta E_\mu },`$ (9)
$`\beta `$ is inverse temperature, and $`_\mu `$ represents sum over all possible eigenstates of the Hamiltonian. It turns out that the variation of $`Z`$ with respect to temperature or any other external parameters affecting the system can virtually tell us everything we might want to know about the macroscopic behavior of the system. For example, the internal energy is given by

$`U={\displaystyle \frac{1}{Z}}{\displaystyle \underset{\mu }{}}E_\mu e^{\beta E_\mu }`$ (10)
From Eq. (9), it is easy to see that the internal energy can also be written in terms of a derivative of the partition function:

$`U={\displaystyle \frac{1}{Z}}{\displaystyle \frac{Z}{\beta }}={\displaystyle \frac{\mathrm{ln}Z}{\beta }}.`$ (11)
The specific heat is given by the derivative of the internal energy:

$`C_v={\displaystyle \frac{U}{T}}=k_B\beta ^2{\displaystyle \frac{U}{\beta }}=k_B\beta ^2{\displaystyle \frac{^2\mathrm{ln}Z}{\beta ^2}}.`$ (12)
where $`k_B`$ is the Boltzmann constant which is set to unit hereafter.

Our aim is to combine the idea of statistical mechanics mentioned above with the numerical solution of the BA equations. The main idea of the numerical Bethe ansatz method we introduce here is, first, to compute all eigenvalues of a BA solvable model from its corresponding BA equations. Then to compute the expectation value of the Hermitian operators, representing the physical observables we are interested in, by averaging those operators over all states of the system, weighting each state with its own Boltzmann weight.

It has been shown MTakahashib  that the Hilbert space of the isotropic Heisenberg model is complete under the string classification. Here we want to show how to travel through all $`C_{N/2}^N`$ states in quantum number space and illustrate it by considering a system of 6 sites.

For the case of $`M`$ down spins, the first task is how to obtain all string configurations fulfilling the restriction (4). We adopt a time-like number “$`\alpha _M:\alpha _{M1}:\mathrm{}:\alpha _2:\alpha _1`$”, where the magnitude $`\alpha _n`$ measures from 0 to $`[M/n]`$ (here $`[x]`$ returns the truncated integer value of $`x`$), just like hour and minute in “HH:MM” measure from 0 to 23 and 0 to 59 respectively. If we increase the number “$`\alpha _M:\mathrm{}:\alpha _1`$” by adding 1 to the first digit $`\alpha _1`$ step by step, we can travel through all possible values. Among all these numerical values, only those whose digits satisfy the condition (4) are what we need. Then all string configurations can be found by this procedure. Of course, these operations are realized in a computer. In order to make the method clear, let us consider a problem of 6 sites.

M=0: It is easy to have the state with all spins up, i. e. $`M=0`$, which has energy $`E=JN/4`$.

M=1: In this case, we only have one string configuration $`\alpha _1=1`$ and one quantum number $`2I_12`$, thus there are 5 states. Each of them is represented by one quantum number in the interval $`[2,2]`$. We can get all possible quantum number configurations from the following figure,

$``$
where the dot is the occupied quantum number, and the open circles represent other possible quantum numbers. Then the BA equation is just

$`6\mathrm{tan}^1x_1^1=\pi I_1,`$ (13)
which has a simple solution $`x_1=\mathrm{tan}(\pi I_1/6)`$.

M=2: Here the string configuration is characterized by $`\{\alpha _1,\alpha _2\}`$. We construct a number “$`\alpha _2:\alpha _1`$”, in which the maximum value of $`\alpha _1`$ is 2 ($`[2/1]=2`$), and $`\alpha _2`$ 1 ($`[2/2]=1`$). Increasing $`\alpha _1`$ step by step we generate all possible configurations of the “$`\alpha _2:\alpha _1`$” number, ranging from 0:0 to 1:2. Among all the generated configurations, we are only interested in those satisfying the condition $`\alpha _1+2\alpha _2=2`$. The first case is $`\alpha _1=2,\alpha _2=0`$, in which the quantum number satisfy $`1.5I_1^1,I_2^11.5`$, the second one is $`\alpha _1=0,\alpha _2=1`$, in which the number satisfy $`1I_1^21`$. They can be characterized by

$``$
and

$``$
respectively, where $``$ denotes the occupation for a quantum number of 2-string. In Table 1, we list all quantum number configurations for $`M=2`$. The BA equations for these two cases are

$`6\mathrm{tan}^1x_1^1=`$ $`\pi I_1^1+\mathrm{tan}^1{\displaystyle \frac{x_1^1x_2^1}{2}},`$
$`6\mathrm{tan}^1x_2^1=`$ $`\pi I_2^1+\mathrm{tan}^1{\displaystyle \frac{x_2^1x_1^1}{2}},`$ (14)
and

$`6\mathrm{tan}^1(x_1^2/2)=\pi I_1^2,`$ (15)
respectively.

M=3: In this case the string configuration is characterized by $`\{\alpha _1,\alpha _2,\alpha _3\}`$. In the same way as we did above, we construct a number “$`\alpha _3:\alpha _2:\alpha _1`$”, the maximum value for each digit from left to right is 1, 1, 3 respectively. Then we have 3 string configurations with the condition $`\alpha _1+2\alpha _2+3\alpha _3=3`$, which correspond to the following sequences

$`a:`$
$`b:`$ $``$
$``$
$`c:\mathrm{\S }`$
where $`a,b,c`$ have 1, 3, 1 states respectively, $`\mathrm{\S }`$ denotes the site for 3-string. And in Table 2, we list all quantum number configurations for $`M=3`$, whose BA equations are

$`6\mathrm{tan}^1x_1^1`$ $`=`$ $`\pi I_1^1+\mathrm{tan}^1(x_1^1x_2^1)`$
$`+\mathrm{tan}^1(x_1^1x_3^1),`$
$`6\mathrm{tan}^1x_2^1`$ $`=`$ $`\pi I_2^1+\mathrm{tan}^1(x_2^1x_1^1)`$
$`+\mathrm{tan}^1(x_2^1x_3^1),`$
$`6\mathrm{tan}^1x_3^1`$ $`=`$ $`\pi I_3^1+\mathrm{tan}^1(x_3^1x_1^1)`$ (16)
$`+\mathrm{tan}^1(x_3^1x_1^1).`$
$`6\mathrm{tan}^1x_1^1`$ $`=`$ $`\pi I_1^1+\mathrm{tan}^1(x_1^1x_1^2)`$
$`+\mathrm{tan}^1((x_1^1x_1^2)/3),`$
$`6\mathrm{tan}^1x_1^2`$ $`=`$ $`\pi I_1^2+\mathrm{tan}^1(x_1^2x_1^1)`$
$`+\mathrm{tan}^1((x_1^2x_1^1)/3),`$
and

$`6\mathrm{tan}^1x_1^3=\pi I_1^3.`$ (18)
respectively.

As a result we have totally $`C_3^6=20`$ distinct configurations of quantum number whose Hilbert space is complete.

Then, we compute the eigenvalue for a given quantum number configuration $`\{I_\gamma ^n\}`$ by solving the BA equations numericaly. For the Heisenberg chain, the BA equations can be solved by iteration, for other models, such as the Hubbard model, the BA equations can be solved by a gradient method.

## IV Monte Carlo Bethe Ansatz

For a system of $`N`$ sites, there are $`C_{N/2}^N`$ quantum number configurations. This number increases exponentially with the size of the system, so it is impossible to calculate all eigenvalues for a large system, such as $`N>40`$, under present computer capacity. This restriction can be overcome by a Monte Carlo method. There are many Monte Carlo methods available, and we introduce below a new method that we call Monte Carlo Bethe ansatz. This method is a classical Monte Carlo strategy applied to a quantum problem. The basic idea behind the MCBA method is to simulate the random thermal fluctuation of the system from state to state in quantum number space of the BA solution. This method is not limited by the sign problem, that may show up in the usual quantum Monte Carlo methods.

Since the energy eigenvalues are a function of both $`M`$ and of the quantum numbers $`I_\gamma ^n`$ we can follow a classical Monte Carlo strategy, by sampling the configuration space of $`M`$ and $`\{I_\gamma ^n\}`$. We now explain how to implement the Monte Carlo calculation, which follows three steps. Let us assume the present state is $`\mu `$ with a corresponding $`M_\mu `$ – the number of down spins in state $`\mu `$. From the state $`\mu `$ any other state $`\nu `$ with $`M_\nu `$, in number of $`C_{N/2}^N1`$, can be obtained.

step one: first we choose $`M_\nu `$, knowing that the number of states with $`M_\nu `$ spins down is $`C_{M_\nu }^NC_{M_\nu 1}^N`$, thus the probability of selecting $`M_\nu `$ is $`(C_{M_\nu }^NC_{M_\nu 1}^N)/C_{N/2}^N`$.

step two: selected $`M_\nu `$, all possible string configurations given $`M_\nu `$ are determined from of Eq. (4) which satisfy MTakahashib

$`{\displaystyle \underset{\alpha _1+\mathrm{}+M\alpha _M=M}{}}D(\{\alpha _n\})=(C_M^NC_{M1}^N),`$ (19)
with $`D(\{\alpha _n\})`$ is the number of states, characterized by the set of quantum numbers $`\{I_\alpha ^n\}`$ associated with the string configuration $`\{\alpha _n\}`$, and reads

$`D(\{\alpha _n\})={\displaystyle \underset{i=1}{\overset{M}{}}}C_{\alpha _j}^{N_{j=1}^Mt_{ij}\alpha _j}.`$ (20)
So, in step two, we select a string configuration with the probability $`D(\{\alpha _n\})/(C_M^NC_{M1}^N)`$.

step three: having determined the string configuration, we then select at random a quantum number configuration, which is the state $`\nu `$ we want, for the given string configuration. From the partition function $`Z`$, the probability density for a state $`\mu `$ is

$`p_\mu =(N2M_\mu +1)e^{\beta E_\mu },`$ (21)
where the degenerancy of state $`\mu `$ was taken into account. The detailed balance condition tells us the transition probability should satisfy

$`{\displaystyle \frac{p_\nu }{p_\mu }}={\displaystyle \frac{(N2M_\nu +1)}{(N2M_\mu +1)}}e^{\beta (E_\nu E_\mu )}.`$ (22)
Hence it is possible to use the Metropolis algorithm for the acceptance ratio to accept or reject the state $`\mu `$ according to

$`A(\mu \nu )=\{\begin{array}{cc}\frac{(N2M_\nu +1)}{(N2M_\mu +1)}e^{\beta (E_\nu E_\mu )},\hfill & \frac{p_\nu }{p_\mu }<1\hfill \\ 1,\hfill & \mathrm{otherwise}.\hfill \end{array}`$ (25)
The MCBA algorithm is complete and the three basic steps are repeated a number of times. After an initial equilibration time, the expectation values can be then estimated as an arithmetic mean over the repeated Markov chain

$`Q={\displaystyle \frac{1}{N}}{\displaystyle \underset{\{\mu \}}{}}Q(\mu ).`$ (26)
## V Specific Heat and Susceptibility

In order to check the validity of our approaches, we apply these two methods to the study of the specific heat and of the magnetic susceptibility of the anti-ferromagnetic and ferromagnetic Heisenberg models.

For the present model, however, because of the degeneracy in each set of quantum number configuration, Eq. (8) should be revised according to the property of the operator. For example, the internal energy and magnetization are

$`E={\displaystyle \frac{1}{Z}}{\displaystyle \underset{\mu }{}}(N2M_\mu +1)E_\mu e^{\beta E_\mu },`$
$`M={\displaystyle \frac{1}{Z}}{\displaystyle \underset{\mu }{}}{\displaystyle \underset{M_\mu ^z=N/2+M_\mu }{\overset{N/2M_\mu }{}}}2M_\mu ^ze^{\beta E_\mu }.`$ (27)
where $`Z=_\mu (N2M_\mu +1)e^{\beta E_\mu }`$. From thermodynamics it is easy to have the expression for specific heat and magnetic susceptibility per site

$`C={\displaystyle \frac{\beta ^2}{N}}(E^2E^2),`$
$`\chi ={\displaystyle \frac{\beta }{N}}(M^2M^2)`$ (28)
We apply NBA to 24-site system and MCBA to 60-site system, respectively. The later has $`C_{30}^{60}`$ different quantum number configurations, hence it is impossible to calculate all the eigenvalues of the system.

In Figs. 1 and 2, we show the specific heat and the magnetic susceptibility, for a 24-site system, obtained from NBA and compare our results with those obtained from TBA. It is clear that the two results match. In Fig 3, we show the specific heat and the magnetic susceptibility, for a 60-site system obtained from MCBA together with the results from TBA. They both agree to each other except at low temperature. In Table 3, we compare, for the ferromagnetic case, the two methods we introduced here with TBA, giving the explicit numerical values. It is clear that our methods work very well for the present model. Hence our conclusion is that for small systems, such as $`N38`$, and for the Heisenberg chain, it is possible to compute all eigenvalues and to obtain all possible thermodynamic quantities of interest by using Eq. (8). For temperatures larger than the finite size energy gap our results agree with TBA results exactly. For larger systems, however, results can still be obtained by using the MCBA method.

Now we study the thermodynamics of the model in the presence of a magnetic field by NBA, which has also been studied by Klümper AKlumper98 . In Fig. 4 and Fig. 5, the results for the specific heat and the magnetic susceptibility of the anti-ferromagnetic case are shown for various magnetic fields. It is clear from these two figures that there are two different behaviors at low temperature, separated by the saturation field $`h_c=1.0`$ at the ground state. In order to understand better this behavior of the antiferromagnetic case, let us use the mapping between the Heisenberg model and the spinless fermion model. This mapping is achieved by the Jordan-Wigner transformation EFradkinb , and Hamiltonian (1) can be written as

$``$ $`=`$ $`{\displaystyle \frac{J}{2}}{\displaystyle \underset{l=1}{\overset{N}{}}}(f_l^{}f_{l+1}+f_{l+1}^{}f_l)`$ (29)
$`+J{\displaystyle \underset{l=1}{\overset{N}{}}}(n_l{\displaystyle \frac{1}{2}})(n_{l+1}{\displaystyle \frac{1}{2}}).`$
where the spinless fermion operators $`f_l^{},f_l`$ obey the usual anti-commutation relation, $`n_l`$ is the usual local number operator. When $`h<h_c`$, the system is not fully polarized, that is $`_{l=1}^Nn_l>0`$, hence we always have two Fermi points $`\pm k_F`$ at the ground state. The dispersion relation of low-lying excitations is dominated by the linear-$`k`$ dependence, hence we still have the Fermi-liquid like specific heat: $`CT`$ at low temperatures. If $`hh_c`$, however, and from the point view of spinless fermions, we have $`_{l=1}^Nn_l=0`$, and the dispersion relation becomes $`k^2`$, because of the $`\mathrm{cos}k`$ dispersion-relation for the fermions in the lattice. Hence, the specific heat manifest a $`T^{1/2}`$ behavior at sufficiently low temperature for $`h=h_c`$, which can be seen in Fig. 4 (open circles). Moreover, the magnetic susceptibility presents a strong peak for $`h=h_c`$, when $`T0`$ \[see Fig. 5, panel (c)\]. This strong magnetic response is associated with a change in the nature of the elementary excitations when the line $`h=h_c`$ is crossed at zero temperature. Indeed, at $`T=0`$ and $`h_c=1.0`$, the system manifest infinite susceptibility, as can be seen from Fig. 5, panel (b). We attribute it due to the degeneracy between the state of $`[N1,1]`$ and $`[N]`$, and a small magnetic field can fully polarize the system. The phase with $`hh_c`$ share anti-ferromagnetic-like behavior \[Fig. 6, panel (b)\], while for $`h<h_c`$, the susceptibility shows a logarithm singularity AKlumper98 .

For the ferromagnetic case the specific heat and the magnetic susceptibility are ploted in Fig. 6, for different values of the magnetic field. As is known, if $`h=0`$ the ground state of the ferromagnetic case is highly degenerate with $`S=N/2,S_z=S,S+1,\mathrm{},S`$ and a very small $`h`$ can fully polarized the system. So it is easy to understand why zero temperature susceptibility is infinite. After it is magnetized (in the presence of small $`h`$), however, the susceptibility should be zero. This behavior is seen in Fig. 6, panel (b). We also show, in Fig. 6, panel (c), the susceptibility obtained by MCBA. Both the results of the two methods agree with each other perfectly.

## VI Summary

In summary, we presented two numerical approaches to discuss the thermodynamics of Bethe ansatz solvable models. The first one is the numerical Bethe ansatz which works very well for a small system. We think it is possible to obtain all eigenvalues of a system up to size $`L=38`$, for the Heisenberg model. For a relatively larger system, we also find that the Monte Carlo simulation in quasi-momentum space works well in the moderate and high temperature regions. At low temperatures the present selection method is not excellent, and a better one is required. The discovery of such a method is a challenging and interesting research problem.

There are many physical quantities of interest at finite temperature which are still not well understood, such as spin stiffness of XXZ model, important to understand the transport properties, because of the complex form of the thermodynamic equations. Our methods provide a new route to compute all these quantities directly from the Bethe ansatz equations.

This work is supported by trans-century projects and Cheung Kong projects of China Education Ministry and by the Portuguese Program PRAXIS XXI under grant number 2/2.1/FIS/302/94. We thank M. Takahashi for sending us the TBA data. SJG and NMRP want to thank the support of the Physics Department of the University of Évora, where part of this research was done. SJG would like to thank Dr. M. B. Luo for helpful discussion on Monte Carlo methods.


text:
# Global Strichartz Estimates for nontrapping Geometries: About an Article by H. Smith and C. Sogge

## 1. Introduction

Consider $`\mathrm{\Theta }B(0,R)^d`$ a smooth ($`C^3`$) obstacle, $`\mathrm{\Omega }=\mathrm{\Theta }^c`$, $`g`$ a smooth ($`C^2`$) metric on $`\overline{\mathrm{\Omega }}`$ equal to $`g_{i,j}(x)=\delta _{i,j}`$ for $`|x|R`$, and the wave equation on $`\mathrm{\Omega }`$:

(1.1)
$$\{\begin{array}{cc}\hfill (_t^2\mathrm{\Delta }_g)u(t,x)& =F(t,x)\hfill \\ \hfill u_{t=0}=f\dot{H}_D^\gamma ,& _tu_{t=0}=g\dot{H}_D^{\gamma 1}\hfill \\ \hfill u_\mathrm{\Omega }& =0\hfill \end{array}$$

where $`\mathrm{\Delta }_g`$ is the Laplace Beltrami operator on $`\mathrm{\Omega }`$, and $`\dot{H}^\gamma `$, $`\gamma `$ are the homogeneous Sobolev spaces associated to the square root of the Laplace operator with Dirichlet boundary conditions on $`\mathrm{\Omega }`$, $`\sqrt{\mathrm{\Delta }_D}`$. In  H. Smith and C. Sogge show that in odd space dimensions, if the perturbation $`(g,\mathrm{\Theta })`$ is assumed to be non trapping, then local (in time) Strichartz estimates imply global (in time) Strichartz estimates for the solutions of the wave equation. The method of proof relies on the exponential decay for the local energy of solutions of the corresponding wave equation with compactly supported initial data, a localization argument on the wave cone and a result by M. Christ and A. Kiselev  allowing to deduce the inhomogeneous estimates from the homogeneous ones. The purpose of this note is to generalize this result to the case of even dimensions of space. In this case, the exponential decay is no more true and one has to replace it by another form of decay. In fact the decay we will use (global $`L^2`$ integrability of the local energy norm for any initial data, i.e. non necessarily compactly supported) is closely related to the exponential decay of the local energy used by H. Smith and C. Sogge, since the proof of our decay relies on a weaker form of the estimates used to prove the exponential decay in scattering theory.

More precisely, the result by H. Smith and C. Sogge reads as follows:

###### Definition 1.1 ().

We say that $`1r,s2p,q+\mathrm{}`$ and $`\gamma `$ are admissible if the following two mixed norm estimates hold

* Local Strichartz estimates. For data $`f,g,F`$ supported in $`\{|x|R\}`$, and $`u`$ solution of (1.1), one has

(1.2)
$$\begin{array}{c}u_{L_t^p([0,1];L_x^q(\mathrm{\Omega })}+(u,_tu_{L_t^{\mathrm{}}([0,1];\dot{H}_D^\gamma \times \dot{H}_D^{\gamma 1}}\hfill \\ \hfill C\left(f_{\dot{H}_D^\gamma }+g_{\dot{H}_D^{\gamma 1}}+F_{L_t^r([0,1];L_x^s(\mathrm{\Omega }))}\right)\end{array}$$
* Global Minkovski Strichartz Estimates. For data $`f,g,F`$ and $`u`$ solution of (1.1) in the case $`\mathrm{\Omega }=^d`$ and $`g_{i,j}(x)=\delta _{i,j}`$, one has

(1.3)
$$u_{L_t^p(;L_x^q(^d))}C\left(f_{\dot{H}_D^\gamma }+g_{\dot{H}_D^{\gamma 1}}+F_{L_t^r(;L_x^s(^d))}\right)$$

In , H. Smith and C. Sogge show:

###### Theorem 1 (\[10, Theorem 1.1\]).

Assume that $`d3`$ is odd, $`1r,s2p,q+\mathrm{}`$ and $`\gamma `$ are admissible, $`p>r`$ and $`\gamma \frac{d1}{2}`$. Then for data $`f,g,F`$ and $`u`$ solution of (1.1), one has

(1.4)
$$u_{L_t^p(;L_x^q(\mathrm{\Omega })}C\left(f_{\dot{H}_D^\gamma }+g_{\dot{H}_D^{\gamma 1}}+F_{L_t^r(;L_x^s(\mathrm{\Omega }))}\right)$$

In this article we are going to prove:

###### Theorem 2.

Assume $`d2`$, $`1r,s2p,q+\mathrm{}`$ and $`\gamma `$ are admissible, $`p>2`$ and $`\gamma <d/2`$. Assume also if $`d=2`$ that $`\mathrm{\Theta }\mathrm{}`$. Then for data $`f,g,F`$ and $`u`$ solution of (1.1), then

(1.5)
$$u_{L_t^p(;L_x^q(\mathrm{\Omega }))}C\left(f_{\dot{H}_D^\gamma }+g_{\dot{H}_D^{\gamma 1}}+F_{L_t^r(;L_x^s(\mathrm{\Omega }))}\right)$$

###### Remark 1.2.

The assumption $`\mathrm{\Theta }\mathrm{}`$ is used to avoid low frequency problems. It could be replaced by a “$`0`$ is not a resonance” assumption (see section 2).

As explained above, the idea in this article is to replace the exponential decay of the local energy by another (weaker) form of decay. Section 2 is hence devoted to the proof of this decay. Then in section 3 we prove Theorem 2.

###### Remark 1.3.

At the final stage of this work, I have learned that the results of this article have been obtained independantly (by a slightly different method) by Jason Metcalfe

## 2. Global $`L^2`$-integrability of the local energy

We begin first by recalling the definition of the homogeneous Sobolev spaces (in the boundary case) $`\dot{H}_D^s`$. Remark that we shall restrict ourselves to the case $`s<d/2`$, hence the multiplication by a smooth function $`\xi C_0^{\mathrm{}}(^d)`$ is continuous from $`\dot{H}^s(^d)`$ to $`H^s(^d)`$ and the two norms are equivalent on functions with fixed compact support. Take $`R`$ large enough so that $`\mathrm{\Theta }B(0,R)`$ and $`g_{i,j}(x)=\delta _{i,j}`$ for $`|x|R`$. Fix $`\beta C_0^{\mathrm{}}(^d)`$ equal to $`1`$ on $`\{|x|R\}`$, and define for $`s0`$ $`\dot{H}_D^s`$ to be the set of functions $`fH_{\text{ loc}}^s(\mathrm{\Omega })`$ such that

(2.1)
$$f_{\dot{H}_D^s}^2=\beta f_{H^s(\mathrm{\Omega })}^2+(1\beta )f_{\dot{H}^s(^d)}^2$$

is bounded and which satisfy the following compatibility conditions:

(2.2)
$$\mathrm{\Delta }_g^jf_\mathrm{\Omega }=0,j;s2j>\frac{1}{2}.$$

For $`s<0`$, define $`\dot{H}_D^s=(\dot{H}_D^s)^{}`$ Remark that on $`\dot{H}_D^s`$, the norms $`(\mathrm{\Delta }_g)^{s/2}f_{L^2}`$ and $`f_{\dot{H}_D^s}`$ are equivalent (it is clear for $`s`$ an even integer and follows in general by interpolation and duality). Denote by $`\dot{}^s=\dot{H}_D^s\times \dot{H}_D^{s1}`$ and $`^s=H_D^s\times H_D^{s1}`$.

In this part we are going to prove:

###### Theorem 3.

Consider $`\chi C_0^{\mathrm{}}(^d)`$. Assume that $`d2`$ and if $`d=2`$ that $`\mathrm{\Theta }\mathrm{}`$. Then For data $`f,g,F`$ such that $`F`$ is supported in $`B_x(0,R)\mathrm{\Omega }`$ and $`u`$ solution of (1.1), one has

(2.3)
$$(\chi u,\chi _tu)_{L_t^2(;\dot{H}_D^\gamma \times \dot{H}_D^{\gamma 1})}C\left(f_{\dot{H}_D^\gamma }+g_{\dot{H}_D^{\gamma 1}}+F_{L_t^2(;\dot{H}_D^{\gamma 1})}\right)$$

###### Remark 2.1.

Remark that in the result above only $`F`$ (and neither $`f`$ nor $`g`$) is assumed to be compactly supported.

###### Remark 2.2.

As will appear below, the proof of this result is closely related to the proof of the exponential decay of the local energy (for compactly supported data). In fact we will use some estimates on the resolvent weaker than the estimates required to prove this exponential decay. The result above can be transposed to the framework of Schrödinger equations. It gives then a (global in time) “smoothing effect” .

###### Remark 2.3.

For $`d=2,\gamma 1/2`$ and $`g_{i,j}(x)=\delta _{i,j}`$, for any $`u`$ solution of (1.1) with $`F=0`$, we have (see \[10, Lemma 2.2\])

(2.4)
$$(\chi u,\chi _tu)_{L_t^2(;\dot{H}_D^\gamma \times \dot{H}_D^{\gamma 1})}C\left(f_{\dot{H}_D^\gamma }+g_{\dot{H}_D^{\gamma 1}}\right)$$

To prove Theorem 3, we first show that by a $`TT^{}`$ argument it is enough to prove a non homogeneous estimate. Indeed denote by $`A=i\left(\begin{array}{cc}0& \text{Id}\\ \mathrm{\Delta }& 0\end{array}\right)`$ the self adjoint operator on $`=H_0^1(\mathrm{\Omega })\times L^2(\mathrm{\Omega })`$ with domain $`𝒟=\{(u,v)H^2(\mathrm{\Omega })\times H_0^1(\mathrm{\Omega })\}`$. The solution of (1.1) with $`F=0`$ is given by

(2.5)
$$\left(\begin{array}{c}u\\ _tu\end{array}\right)=e^{itA}\left(\begin{array}{c}f\\ g\end{array}\right)$$

Remark that the operator $`e^{itA}`$ is a group of isometries on $`\dot{}_D^s`$ (equipped with the norm $`A^sf_{L^2}`$).

Denote by $`T=\chi e^{itA}`$. The continuity of $`T`$ from $`\dot{}^s`$ to $`L^2(_t;^s)`$ is equivalent to the continuity of $`T^{}`$ from $`L^2(_t;^s)`$ to $`\dot{}^s`$. But

(2.6)
$$T^{}G=_se^{isA}\chi G(s,)𝑑s$$

and the continuity of $`T^{}`$ is in turn equivalent to the continuity of $`TT^{}`$ from $`L^2(_t;\dot{}^s)`$ to $`L^2(_t;^s)`$. Finally

(2.7)
$$TT^{}G(t,)=_s\chi e^{i(ts)A}\chi G(s,)𝑑s=_{s<t}+_{t<s}$$

And it sufficies to prove the continuity of any of the two operators above, say for example the continuity of $`\chi TT_1^{}\chi `$ with

(2.8)
$$TT_1^{}G(t,)=_{s<t}e^{i(ts)A}G(s,)𝑑s$$

Remark that this latter result implies in turn the nonhomogeneous part of Theorem 3. Suppose that $`G`$ is supported in $`\{s>0\}`$ (which by density and translation invariance is possible). Then $`U=TT_1^{}G`$ satisfies the equation

(2.9)
$$(i_t+A)U=\chi G(t,),U_{t<0}=0$$

The Fourier transforms with respect to $`t`$ of $`U`$ and $`G`$ are holomorphic in the half plane $`\text{Im}z<0`$ (due to the support property) and satisfy there

(2.10)
$$(Az)\widehat{U}=\chi \widehat{G}$$

Taking $`z=xi\epsilon `$ and letting $`\epsilon >0`$ tend to $`0`$, Theorem 3 is a consequence of the fact that the Fourier transform is an isometry on $`L^2(_t;)`$ if $``$ is a Hilbert space and of

###### Proposition 2.4.

Under the assumptions of Theorem 2, the resolvent $`\chi (A(xi\epsilon ))^1\chi `$ is for $`x`$ and $`0<|\epsilon |<1`$, uniformly bounded on $`_D^s`$.

We are going to deduce Proposition 2.4 from the classical result:

###### Proposition 2.5.

Suppose that the obstacle $`\mathrm{\Theta }`$ is non trapping. Then the resolvent of the operator $`\mathrm{\Delta }_D`$, $`(\mathrm{\Delta }_D\lambda )^1`$ (which is analytic in $`^+`$) satisfies:

(2.11)
$$\begin{array}{c}\chi C_0^{\mathrm{}}(^2),C>0;\lambda ,0<\epsilon <<1,\\ \chi (\mathrm{\Delta }_D(\lambda \pm i\epsilon ))^1\chi _{L^2L^2}\frac{C}{1+\sqrt{|\lambda |}}\end{array}$$

###### Remark 2.6.

Proposition 2.5 was proven for $`|\lambda |>>1`$ in greater generalities by Lax and Phillips , Melrose and Sjöstrand ), Vainberg, Tang Zworski  (see also  for a self contained proof which joined with the results in  would relax the smoothness assumption required in these papers to a $`C^2`$ assumption). The proof for $`|\lambda |<<1`$ can be found in \[2, Annexe B.2\] (see also ). Remark that in  the Poincaré inequality is used to control the local $`L^2`$-norm by the local $`L^2`$-norm of the gradient of a function (wich is why $`\mathrm{\Theta }\mathrm{}`$ is required). However, the following inequality

(2.12)
$$\beta <\frac{d}{2}1,uC_0^{\mathrm{}}(^d),\frac{u}{r(1+r^2)^{\beta /2}}_{L^2}C\frac{u}{(1+r^2)^{\beta /2}}_{L^2}$$

allows to handle the argument for $`d4`$ (and the result is standard in odd dimensions).

Since

(2.13)
$$(Az)^1=\left(\begin{array}{cc}z(\mathrm{\Delta }_g+z^2)^1& i(\mathrm{\Delta }_g+z^2)^1\\ i\mathrm{\Delta }_g(\mathrm{\Delta }_g+z^2)^1& z(\mathrm{\Delta }_g+z^2)^1\end{array}\right)$$

to prove Proposition 2.4 we have to estimate uniformly for any $`s`$

(2.14) $`\chi (1+|z|)(\mathrm{\Delta }_g+z^2)^1\chi _{H_D^sH_D^s}`$
(2.15) $`\chi (\mathrm{\Delta }_g+z^2)^1\chi _{H_D^sH_D^{s+1}}`$
(2.16) $`\chi \mathrm{\Delta }_g(\mathrm{\Delta }_g+z^2)^1\chi _{H_D^sH_D^{s1}}`$
Remark that using duality and interpolation it is enough to do this for $`s`$.

For $`s=0`$, (2.14) is Proposition 2.5 and if $`u=(\mathrm{\Delta }_g+z^2)^1\chi f`$, we have

(2.17) $`{\displaystyle _\mathrm{\Omega }}\chi f\chi \overline{u}`$ $`={\displaystyle _\mathrm{\Omega }}(\mathrm{\Delta }+z^2)u\chi \overline{u}`$
$`={\displaystyle _\mathrm{\Omega }}\chi |u^2\chi u\overline{u}+z^2\chi |u|^2`$
which by Cauchy Schwartz implies (using (2.14) for $`s=0`$ and a different $`\chi `$) easily (2.15) for $`s=0`$. By duality we get also (2.15) for $`s=1`$. Applying $`\mathrm{\Delta }_g`$ to $`\chi (\mathrm{\Delta }_g+z^2)^1\chi `$ and commuting with $`\chi `$ we obtain by induction (2.15) for $`s`$ (hence by duality and interpolation for any $`s`$). Applying again $`\mathrm{\Delta }_g`$ we deduce (2.16) for any $`s`$. Finally, remark that

(2.18)
$$\chi \mathrm{\Delta }_g(\mathrm{\Delta }_g+z^2)^1\chi =\chi \text{Id}\chi +\chi z^2(\mathrm{\Delta }_g+z^2)^1\chi $$

from which we estimate uniformly

(2.19)
$$\chi z^2(\mathrm{\Delta }_g+z^2)^1\chi _{H_D^sH_D^{s1}}.$$

Interpolating this latter inequality with $`\text{(}\text{2.15}\text{)}_s`$ we get $`\text{(}\text{2.14}\text{)}_s`$.

## 3. From Energy decay to global Strichartz estimates

In this section we are going to show how the local (in time) Strichartz estimates and the global (in time) $`L^2`$\- energy decay imply global (in time) Strichartz estimates. We consider firstly the homogeneous case ($`F=0`$) and we also assume that $`\gamma 1/2`$ if $`d=2`$.

Consider $`\chi C_0^{\mathrm{}}(^d)`$ equal to $`1`$ on $`B(0,R)`$, $`u`$ solution of (1.1) with $`F=0`$ and $`v=\chi (x)u`$, $`w=(1\chi )(x)u`$ solutions of

(3.1)
$$\{\begin{array}{cc}\hfill (_t^2\mathrm{\Delta }_g)v(t,x)& =[\mathrm{\Delta }_g,\chi ]u\hfill \\ \hfill v_{t=0}& =\chi f\dot{H}^\gamma \hfill \\ \hfill _tv_{t=0}& =\chi g\dot{H}_D^{\gamma 1}\hfill \\ \hfill v_\mathrm{\Omega }& =0\hfill \end{array}$$

(3.2)
$$\{\begin{array}{cc}\hfill (_t^2\mathrm{\Delta }_g)w(t,x)& =(_t^2\mathrm{\Delta }_0)w(t,x)=[\mathrm{\Delta }_0,\chi ]u\hfill \\ \hfill w_{t=0}& =(1\chi )f\dot{H}^\gamma ,\hfill \\ \hfill _tw_{t=0}& =(1\chi )g\dot{H}^{\gamma 1}\hfill \end{array}$$

First we show that $`w`$ satisfies global Strichartz stimates. Since $`w`$ is a solution of the free wave equation (because it is supported in the set where $`g_{i,j}(x)=\delta _{i,i}`$), the contributions of $`(1\chi )f`$ and $`(1\chi )g`$ are handled using (1.3). Consequently it sufficies to study $`\stackrel{~}{w}`$ the solution of

(3.3)
$$\{\begin{array}{cc}\hfill (_t^2\mathrm{\Delta }_0)\stackrel{~}{w}(t,x)& =[\mathrm{\Delta }_0,\chi ]u\hfill \\ \hfill \stackrel{~}{w}_{t=0}& =0\hfill \\ \hfill _t\stackrel{~}{w}_{t=0}& =0\hfill \end{array}$$

For this we use the following result by M. Christ and A. Kiselev :

###### Theorem 4 (M. Christ and A. Kiselev).

Consider a bounded operator $`T:L^p(;B_1)L^q(;B_2)`$ given by a locally integrable kernel $`K(t,s)`$ with value operators from $`B_1`$ to $`B_2`$ where $`B_{1;2}`$ are Banach spaces. Suppose that $`p<q`$. Then the operator

(3.4)
$$\stackrel{~}{T}f(t)=_{s<t}K(t,s)f(s)𝑑s$$

is bounded from $`L^p(;B_1)`$ to $`L^q(;B_2)`$ by

(3.5)
$$\stackrel{~}{T}_{L^p(;B_1)L^q(;B_2)}(12^{p^1q^1})^1T_{L^p(;B_1)L^q(;B_2)}$$

According to (3.3),

(3.6)
$$(\stackrel{~}{w},_t\stackrel{~}{w})(t,x)=_0^te^{i(ts)A_0}[\mathrm{\Delta }_g,\chi ]u(s,)ds$$

According to Theorem 3 (and Remark 2.3 since $`\gamma 1/2`$ if $`d=2`$) applied to the free wave operator, for any $`\chi C_0^{\mathrm{}}(^d)`$ the operator $`T_0=\chi e^{itA_0}`$ maps $`\dot{H}^\gamma \times \dot{H}^{\gamma 1}`$ into $`L^2(_t^+;\dot{H}^\gamma \times \dot{H}^{\gamma 1})`$ and consequently its adjoint

(3.7)
$$T_0^{}G=_0^{\mathrm{}}e^{isA_0}\chi G𝑑s$$

maps $`L^2(_t^+;\dot{H}^\gamma \times \dot{H}^{\gamma 1})`$ into $`\dot{H}^\gamma \times \dot{H}^{\gamma 1}`$ and since

(3.8)
$$(\stackrel{~}{w},_t\stackrel{~}{w})=_0^te^{i(ts)A_0}(0,[\mathrm{\Delta }_g,\chi ]u)𝑑s,$$

the global Strichartz estimate follows from Theorem 4 (here we are speaking about adjoints relative to the $`L^2`$ in time and $`\dot{H}^\gamma \times \dot{H}^{\gamma 1}`$ in space norms).

We come back to the analysis of $`v=\chi u`$. Consider $`\phi C_0^{\mathrm{}}(]0,1[)`$ equal to $`1`$ on $`[1/4,3/4]`$. Denote by $`v_n=\phi (tn/2)v`$ solution of

(3.9)
$$\{\begin{array}{cc}\hfill (_t^2\mathrm{\Delta }_g)v_n(t,x)& =\phi (tn/2)[\mathrm{\Delta }_g,\chi ]u+[_t^2,\phi (tn/2)]\chi u\hfill \\ & =u_n(t,x)\hfill \\ \hfill v_n_{t<n/2}& =0\hfill \\ \hfill _tv_n_{t<n/2}& =0\hfill \\ \hfill v_{n_\mathrm{\Omega }}& =0\hfill \end{array}$$

With, according to Theorem 3,

(3.10)
$$\underset{n}{}u_n_{L^2(_t;\dot{H}_D^{\gamma 1})}^2C(f_{\dot{H}_D^\gamma }^2+g_{\dot{H}_D^{\gamma 1}}^2)$$

According to the (local in time) Strichartz inequality (1.2) (and Minkovski inequality), we obtain

(3.11)
$$v_n_{L_t^p(_t;L_x^q(\mathrm{\Omega }))}+(v_n,_tv_n_{L_t^{\mathrm{}}([0,1];\dot{H}_D^\gamma \times \dot{H}_D^{\gamma 1})}C(u_n_{L_t^1([1+n/2,2+n/2];H_x^{\gamma 1}(\mathrm{\Omega }))})$$

Since

(3.12)
$$v_{L_t^p(;L_x^q(\mathrm{\Omega }))}^p\underset{n}{}v_n_{L_t^p(_t;L_x^q(\mathrm{\Omega })}^p$$

and since $`p2`$, we obtain, using (3.10) and (3.11)

(3.13) $`v_{L_t^p(;L_x^q(\mathrm{\Omega }))}^p`$ $`{\displaystyle \underset{n}{}}v_n_{L_t^p(_t;L_x^q(\mathrm{\Omega }))}^pC{\displaystyle \underset{n}{}}u_n_{L_t^1;H_x^{\gamma 1}(\mathrm{\Omega })}^p`$
$`C\left({\displaystyle \underset{n}{}}u_n_{L_t^2;H_x^{\gamma 1}(\mathrm{\Omega })}^2\right)^{p/2}C(f_{\dot{H}_D^\gamma }^2+g_{\dot{H}_D^{\gamma 1}}^2)^{p/2}`$
which proves Theorem 2 for $`F=0`$. As recognized by H. Smith and C. Sogge , another application of Theorem 4 gives the general case.

Finally, if $`d=2`$ we have to get rid of the assumption $`\gamma 1/2`$. Interpolating between (2.4) and the energy estimate (for any $`\gamma <1`$)

(3.14)
$$(\chi u,\chi _u)_{L^{\mathrm{}}(_t;\dot{H}_D^\gamma \times \dot{H}^{\gamma 1})}(f,g)_{L^{\mathrm{}}(_t;\dot{H}_D^\gamma \times \dot{H}^{\gamma 1})}$$

we obtain for any $`q>1/(1\gamma )`$

(3.15)
$$(\chi u,\chi _tu)_{L^q(_t;\dot{H}_D^\gamma \times \dot{H}^{\gamma 1})}(f,g)_{L^{\mathrm{}}(_t;\dot{H}_D^\gamma \times \dot{H}^{\gamma 1})}$$

using this to estimate $`\stackrel{~}{w}`$ and Christ-Kiselev Theorem we can conclude as above provided $`p>1/(1\gamma )`$. Finally we remark that this latter condition is automatically fulfilled due to the scaling condition $`\frac{1}{p}+\frac{d}{q}=\frac{d}{2}\gamma `$ which is necessary for (1.2) to hold.