text:
# I INTRODUCTION

## I INTRODUCTION

Exclusive nonleptonic $`B`$ meson decays provide information of important dynamics such as penguin contributions, CP asymmetries, and weak and strong phases. However, it is difficult to analyze these processes because of their nonperturbative origin. To simplify the analysis, the factorization approximation (FA) has been applied to decay amplitudes, under which nonfactorizable and annihilation contributions are neglected, and factorizable contributions are expressed as products of Wilson coefficients, meson decay constants, and hadronic transition form factors. In this paper we shall compare three approaches to nonleptonic charmless $`B`$ meson decays, which are more or less related to FA: the FA approach , the perturbative QCD (PQCD) approach , and the Beneke-Buchalla-Neubert-Sachrajda (BBNS) approach . We take the $`B\pi \pi `$ decays as an example, and explain how factorizable, nonfactorizable, and annihilation contributions, infrared-cutoff and scale dependences, and strong phases are treated in these approaches.

The effective Hamiltonian for the $`B\pi \pi `$ decays is given by

$$H_{\mathrm{eff}}=\frac{G_F}{\sqrt{2}}\underset{q=u,c}{}V_q\left[C_1(\mu )O_1^{(q)}(\mu )+C_2(\mu )O_2^{(q)}(\mu )+\underset{i=3}{\overset{10}{}}C_i(\mu )O_i(\mu )\right],$$
(1)
with the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements $`V_q=V_{qd}^{}V_{qb}`$ and the four-fermion operators

$`O_1^{(q)}=(\overline{d}_iq_j)_{VA}(\overline{q}_jb_i)_{VA},O_2^{(q)}=(\overline{d}_iq_i)_{VA}(\overline{q}_jb_j)_{VA},`$ (2)
$`O_3=(\overline{d}_ib_i)_{VA}{\displaystyle \underset{q}{}}(\overline{q}_jq_j)_{VA},O_4=(\overline{d}_ib_j)_{VA}{\displaystyle \underset{q}{}}(\overline{q}_jq_i)_{VA},`$ (3)
$`O_5=(\overline{d}_ib_i)_{VA}{\displaystyle \underset{q}{}}(\overline{q}_jq_j)_{V+A},O_6=(\overline{d}_ib_j)_{VA}{\displaystyle \underset{q}{}}(\overline{q}_jq_i)_{V+A},`$ (4)
$`O_7={\displaystyle \frac{3}{2}}(\overline{d}_ib_i)_{VA}{\displaystyle \underset{q}{}}e_q(\overline{q}_jq_j)_{V+A},O_8={\displaystyle \frac{3}{2}}(\overline{d}_ib_j)_{VA}{\displaystyle \underset{q}{}}e_q(\overline{q}_jq_i)_{V+A},`$ (5)
$`O_9={\displaystyle \frac{3}{2}}(\overline{d}_ib_i)_{VA}{\displaystyle \underset{q}{}}e_q(\overline{q}_jq_j)_{VA},O_{10}={\displaystyle \frac{3}{2}}(\overline{d}_ib_j)_{VA}{\displaystyle \underset{q}{}}e_q(\overline{q}_jq_i)_{VA},`$ (6)
where $`i,j`$ are color indices and $`\mu `$ a renormalization scale.

Decay amplitudes are written as

$``$ $``$ $`\pi (P_2)\pi (P_3)|H_{\mathrm{eff}}|B(P_1),`$ (7)
$``$ $`{\displaystyle \underset{i}{}}C_i(\mu )\pi (P_2)\pi (P_3)|O_i(\mu )|B(P_1),`$ (8)
with the meson momenta $`P_i`$. For convenience, we assume $`P_1=(P_1^+,P_1^{},𝐏_{1T})=(M_B/\sqrt{2})(1,1,0_T)`$, $`P_2=(M_B/\sqrt{2})(1,0,0_T)`$, and $`P_3=(M_B/\sqrt{2})(0,1,0_T)`$, $`M_B`$ being the $`B`$ meson mass. The dependences on the renormalization scale $`\mu `$ cancel between the Wilson coefficients $`C_i(\mu )`$ and the hadronic matrix elements $`\pi \pi |O_i(\mu )|B`$, such that the decay amplitudes are scale-independent (and also scheme-independent). In the practical derivation of the effective Hamiltonian, an infrared cutoff for loop integrals needs to be introduced, since four-fermion amplitudes may contain soft divergences, which are related to nonperturbative dynamics involved in the matrix elements. If the infrared cutoff appears in the form that external fermions are off-shell, the derivation will even become gauge dependent . These potential infrared-cutoff and gauge dependences are usually hidden in the matrix elements.

In the FA approach the matrix element of, say, $`O_2`$ is expressed as

$`\pi (P_2)\pi (P_3)|O_2(\mu )|B(P_1)`$ $``$ $`\pi (P_2)|(\overline{d}_iq_i)_{VA}|0\pi (P_3)|(\overline{q}_jb_j)_{VA}|B(P_1),`$ (9)
$``$ $`f_\pi F^{B\pi }(q^2=M_\pi ^2),`$ (10)
where $`f_\pi `$ is the pion decay constant and $`F^{B\pi }(M_\pi ^2)`$ the $`B\pi `$ transition form factor evaluated at $`q^2=M_\pi ^2`$, $`M_\pi `$ being the pion mass. The form factor itself is parametrized by models. Though analyses are simple under FA, a serious problem arises. When decay amplitudes are expressed as the products of $`\mu `$-dependent Wilson coeficients and $`\mu `$-independent hadronic form factors, they depend on the scale $`\mu `$. A cure of this scale dependence has been proposed in , which will be discussed in Sec. II.

The PQCD approach to exclusive $`B`$ meson decays has been developed some time ago . PQCD is a method to separate hard components from a QCD process, which are treated by perturbation theory. Nonperturbative components are organized in the form of hadron wave functions, which can be extracted from experimental data. This formalism, so-called factorization theorem, has been successfully applied to various semileptonic decays  and nonleptonic (charmed and charmless) decays . For the $`B\pi \pi `$ decays, the hard amplitude involves six external on-shell quarks, four of which correspond to the four-fermion operators and two of which are the spectator quarks in the $`B`$ and $`\pi `$ mesons. Since nonperturbative dynamics has been factored out, one can evaluate all possible Feynman diagrams for the six-quark amplitude straightforwardly, which include both factorizable and nonfactorizable contributions. Factorizable and nonfactorizable annihilation diagrams are also included. That is, FA is not necessary in the PQCD approach. Without FA, the scale independence can be achieved easily . It has been shown that nonperturbative meson wave functions, being universal, are the same for various topologies of diagrams in all decay modes, which contain the $`B`$ and $`\pi `$ mesons.

Recently, Beneke et al. proposed an alternative formalism for two-body charmless $`B`$ meson decays . They claimed that factorizable contributions (the form factor $`F^{B\pi }`$ in this case) are not calculable in PQCD, but nonfactorzable contributions are in the heavy quark limit. Hence, the former are treated in the same way as in the FA approach, and expressed as products of Wilson coefficients and $`F^{B\pi }`$. The latter, calculated as in the PQCD approach, are written as the convolutions of hard amplitudes with three $`(B,\pi ,\pi )`$ meson wave functions. Annihilation diagrams are neglected as in FA. Therefore, this formalism can be regarded as a mixture of the FA and PQCD approaches. Values of form factors at maximal recoil $`q^2=M_\pi ^2`$ and nonperturbative meson wave functions are all treated as free parameters. The BBNS approach then involves more parameters compared to FA and PQCD. This approach has been applied to the $`B\pi \pi `$ and $`K\pi `$ decays .

In this work we shall make a detailed comparision of the above three approaches to nonleponic charmless $`B`$ meson decays. We discuss their treatments of the infrared-cutoff and scale dependences in Sec. II, factorizable contributions in Sec. III, nonfactorizable contributions in Sec. IV, and annihilation contributions and strong phases in Sec. V. It will be shown that the CP asymmetry of the $`B_d^0\pi ^\pm \pi ^{}`$ decays predicted in PQCD is much larger than in the FA and BBNS approaches. Hence, the measurement of CP asymmetries will provide an appropriate test of these predictions. Section VI is the summary. Note that final-state-interaction (FSI) effects are always assumed to be absent in the above three approaches. The neglect of FSI effects has been argued in .

## II SCALE DEPENDENCE

As stated in the Introduction, a serious problem of FA is the dependence of physical predictions on the renormalization scale $`\mu `$. A plausible solution to the aforementioned problem has been proposed in . The idea is to extract a $`\mu `$-dependent evolution factor $`g(\mu )`$ from the matrix element $`\pi \pi |O(\mu )|B`$, and combine it with the $`\mu `$-dependent Wilson coefficient $`C(\mu )`$:

$$C(\mu )\pi \pi |O(\mu )|B=C(\mu )g(\mu )\pi \pi |O|B_{\mathrm{tree}}C^{\mathrm{eff}}\pi \pi |O|B_{\mathrm{tree}},$$
(11)
where the effective Wilson coefficient $`C^{\mathrm{eff}}`$ is $`\mu `$-independent. The extraction of $`g(\mu )`$ involves the mixing of effective operators below the scale $`\mu `$ and above the lower scale $`m_b`$, $`m_b`$ being the $`b`$ quark mass. After making a physical amplitude explicitly scale-independent, FA is applied to the matrix element at the tree level,

$$\pi \pi |O|B_{\mathrm{tree}}f_\pi F^{B\pi }.$$
(12)
This strategy has also been adopted in the BBNS approach.

The above issue is in fact much subtler. One-loop corrections to on-shell $`b`$ quark decay amplitudes are usually infrared divergent . To remove these infrared divergences, the conventional treatment is to assume that external quarks are off-shell by $`p^2`$. This assumption, however, results in a gauge dependence of one-loop corrections. The evolution factor $`g(\mu )`$ then depends on an infrared cutoff $`p^2`$ and on a gauge parameter, which are originally implicit in the matrix element. Therefore, the modified FA, though removes the scale dependence of a decay amplitude, often introduces the infrared-cutoff and gauge dependences.

The controversy of the scale, infrared-cutoff and gauge dependences can be resolved in the PQCD approach . In this formalism partons, i.e., external quarks, are assumed to be on shell, and both ultraviolet and infrared divergences in radiative corrections are isolated using the dimensional regularization. With external quarks being on shell, gauge invariance of a decay amplitude is guaranteed under radiative corrections. The ultraviolet poles are subtracted in a renormalization scheme. The infrared poles, which correspond to nonperturbative dynamics of decay processes, are absorbed into meson wave functions. How much finite piece of radiative corrections is absorbed into wave functions along with the infrared poles depends on a factorization scheme, or equivalently, on a factorization scale $`\mu _f`$. In the PQCD approach $`\mu _f`$ has been chosen as the inverse of the meson transverse extent, $`1/b`$. In other words, dynamics below (above) the factorization scale is absorbed into a nonperturbative wave function (perturbative hard amplitude), and a decay amplitude is factorized into convolution of these two subprocesses.

For nonleptonic $`B`$ meson decays, contributions above the factorization scale involve two scales: $`M_W`$, at which the matching conditions of the effective weak Hamiltonian to the full Hamiltonian are defined, and the typical scale $`t`$ of a hard amplitude, which reflects specific dynamics of a decay mode. The scale $`t`$, corresponding to the lower scale $`m_b`$ in the modified FA, is chosen as the maximal virtuality of internal particles in a hard amplitude. The analysis of next-to-leading-order corrections to the pion form factor  has indicated that this choice minimizes higher-order corrections.

Without assuming FA, the $`\mu `$ dependence of a decay amplitude is calculated based on six-fermion hard amplitudes in the PQCD approach. Radiative corrections produce various types of logarithms: $`\mathrm{ln}(M_W/\mu )`$, $`\mathrm{ln}(\mu /t)`$ and $`\mathrm{ln}(tb)`$. The renormalization-group (RG) summation of the first type of logarithms leads to the Wilson coefficients $`C(\mu )`$ in Eq. (11). The second and third types are summed into the evolution factor

$$g(\mu )g_1(\mu ,t)g_2(t,1/b).$$
(13)
The above expression implies that $`g`$ in fact contains two pieces, $`g_1`$ and $`g_2`$, which are governed by different anomalous dimensions. The anomalous dimension of $`g_1`$ is the same as of $`C(\mu )`$. Below the hard scale $`t`$, loop corrections associated with spectator quarks contribute, such that the anomalous dimension of $`g_2`$ differs from that of $`g_1`$ .

Insert Eq. (13) into Eq. (11), the effective Wilson coefficient is written as

$$C^{\mathrm{eff}}C(t)g_2(t,1/b),$$
(14)
which consists of two stages of RG evolutions. We have extended the Wilson evolution directly from $`M_W`$ down to $`t`$ without further considering the mixing of effective operators. The scale independence of a decay amplitude is thus constructed explicitly. The Wilson evolution is then followed by the second-stage evolution from $`t`$ down to $`1/b`$, which plays the role of the infrared cutoff $`p^2`$.

Together with the second-stage RG evolution, there also exists the Sudakov evolution $`\mathrm{exp}[s(P,b)]`$ from the resummation of double logarithms $`\mathrm{ln}^2(Pb)`$ , which arise from the overlap of collinear and soft divergences in radiative corrections to meson wave functions, $`P`$ denoting the dominant light-cone component of meson momentum. The effect of the Sudakov evolution will be explained in the next section. After summing all various large logarithms, the hard amplitude $`H(t)`$ can be evaluated at the characteristic scale $`t`$ perturbatively. $`H(t)`$ contains all possible factorizable and nonfactorizable diagrams.

At last, we address the treatment of the dependence on the factorization scale in the PQCD approach. Both the hard amplitudes and the meson wave functions $`\varphi (x,b)`$, $`x`$ being the momentum fraction associated with a valence quark, depend on factorization schemes. For example, $`\varphi (x,b)`$ and $`\varphi (x,b/2)`$ are different. However, a decay amplitude, as convolution of the two subprocesses, does not depend on factorization schemes. The wave functions, though not calculable, are universal (process-independent), since they absorb long-distance dynamics, which is insensitive to the specific decay of the $`b`$ quark into light quarks with large energy release. In the practical application $`\varphi (x,b/c)`$ for the factorization scale $`c/b`$, $`c`$ being a constant of order unity, can be determined from experimental data of some decay modes. These wave functions are then employed to make predictions for other modes, whose factorizations are constructed at the same scale $`c/b`$. With this prescription, the dependence on the factorization scale is removed. Usually, $`c`$ is chosen as unity.

A three-scale factorization formula for exclusive nonleptonic $`B`$ meson decays is then written as

$`C(t)H(t)\varphi (x,b)g_2(t,b)\mathrm{exp}\left[s(P,b)\right].`$ (15)
Note that Eq. (15) is a convolution relation, with internal parton kinematic variables $`x`$ and $`b`$ integrated out ($`t`$ depends on $`x`$ and $`b`$ as defined later). All the convolution factors, except for the wave functions $`\varphi (x,b)`$, are calculable in perturbation theory. Comparing Eq. (15) with Eq. (11), the tree-level matrix element is read off as

$$\pi \pi |O|B_{\mathrm{tree}}=H(t)\varphi (x,b)\mathrm{exp}\left[s(P,b)\right].$$
(16)
If choosing $`t`$ as $`m_b`$, Eq. (15) reduces to a simple product of the constant Wilson coefficient $`C(m_b)`$ and the hadronic matrix element.

## III FACTORIZABLE CONTRIBUTIONS

In the FA approach factorizable contributions are expressed as products of Wilson coefficients and transition form factors, which are then parametrized by some models. This treatment of factorizable contributions has also been adopted in the BBNS approach. It has been argued that the $`B\pi `$ transition form factor $`F^{B\pi }`$ can not be evaluated in PQCD, because the on-shellness of internal particles in the hard amplitude is not suppressed by meson wave functions . For example, the lowest-order hard amplitude in Fig. 1 consists of the internal $`b`$ quark and gluon propagators proportional to

$`{\displaystyle \frac{1}{[(P_1k_3)^2m_b^2](k_1k_3)^2}}`$ $`=`$ $`{\displaystyle \frac{1}{(2P_1k_3)(2k_1k_3)}}={\displaystyle \frac{1}{x_1x_3^2M_B^4}},`$ (17)
where $`k_1`$ ($`k_3`$) is the momentum of the spectator quark in the $`B`$ ($`\pi `$) meson, and $`x_1=k_1^+/P_1^+`$ ($`x_3=k_3^{}/P_3^{}`$) the corresponding momentum fraction as indicated in Fig. 1. To derive the above expression, we have neglected the mass difference $`\overline{\mathrm{\Lambda }}=M_Bm_b`$. Obviously, the pion wave function, if proportional to $`x_3(1x_3)`$, does not remove the power divergence at $`x_30`$ from the on-shell internal particles. The assumption that internal particles in a hard amplitude should be off-shell is then violated, and PQCD is not applicable to $`F^{B\pi }`$.

A different viewpoint has been adopted in the PQCD approach. Since the end-point divergence is not of the pinched type, which is absorbed into a wave function, we argue that it can be either killed by a pion wave function vanishing faster than $`x_3`$ as $`x_30`$, or smeared out by parton transverse momenta $`k_T`$ . Including $`k_T`$, Eq. (17) becomes

$`{\displaystyle \frac{1}{(x_3M_B^2+k_{3T}^2)[x_1x_3M_B^2+(𝐤_{1T}𝐤_{3T})^2]}},`$ (18)
which has no singularity at $`x_30`$. As argued in , partons can carry small amount of transverse momenta. When the end-point region with $`x0`$ is not important, it is proper to approximate Eq. (18) for the hard amplitude by Eq. (17), and the transverse degrees of freedom of partons in meson wave functions can be integrated out. When the end-point region is important as in the case of $`F^{B\pi }`$, one should keep $`k_T`$ for a consistent analysis. Once the transverse degrees of freedom of partons are taken into account, the factorization of a QCD process should be performed in the $`b`$ space, $`b`$ being the variable conjugate to $`k_T`$ . This is how the factorization scale $`1/b`$ is introduced in the previous seciton.

As stated before, the scale $`t`$ should be chosen as the maximal virtuality of internal particles in a hard amplitude,

$$t=\mathrm{max}(\sqrt{x_1x_3}M_B,\sqrt{x_3}M_B,1/b),$$
(19)
according to Eq. (18). The Sudakov factor $`\mathrm{exp}[s(P,b)]`$ suppresses the long-distance contributions in the large $`b`$ region, and vanishes as $`b=1/\mathrm{\Lambda }_{\mathrm{QCD}}`$, $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ being the QCD scale. This factor then guarantees that most contributions to decay amplitudes come from the region with small coupling constant $`\alpha _s(t)`$ . It has been shown that 97% of the contribution to $`F^{B\pi }`$ arises from the region with $`\alpha _s(t)/\pi <0.3`$ . It indicates that $`F^{B\pi }`$ at large recoil of the pion is well within the perturbative regime, and that dynamics from hard gluon exchanges dominates. This is the reason PQCD with Sudakov suppression is applicable to exclusive decays around the energy scale of the $`B`$ meson mass . The above conclusion is consistent with the following physical picture. The $`B\pi `$ transition occurs via the $`bu`$ decay with the outgoing $`u`$ quark being energetic. The spectator quark of the $`B`$ meson is more or less at rest. The probability that this spectator quark and the $`u`$ quark with large relative velocity form the pion is strongly suppressed by the pion wave function and by the Sudakov factor. Therefore, the process dominates, where a hard gluon is exchanged so that the two quarks are aligned to form the pion.

We need to illuminate another concept about the application of PQCD. The smallness of PQCD predictions compared to experimental data does not definitely imply that perturbative contribution is not important . PQCD predictions depend on wave funcitons one adopts. For example, a smooth $`B`$ meson wave function decreases PQCD predictions for $`F^{B\pi }`$ by a huge extent. However, this smooth $`B`$ meson wave function may not be correct, such that a conclusion drawn from it is not valid. In our analysis the $`B`$ meson and pion wave functions are determined from the data of the pion form factor and of the $`BD\pi (\rho )`$ decays. We then employ the extracted wave funcitons to predict $`F^{B\pi }`$ and examine the dominance of perturbative contributions. In this way the ambiguity in choosing wave functions and the model dependence of PQCD predictions are reduced.

Below we highlight the enhancement of penguin contributions observed in the PQCD approach, and its role in the explanation of the $`B\pi \pi `$ and $`BK\pi `$ data. For simplicity, we present the observation by means of the FA approach. Consider the ratios $`R`$ and $`R_\pi `$,

$`R={\displaystyle \frac{\mathrm{Br}(B^0K^\pm \pi ^{})}{\mathrm{Br}(B^\pm K^0\pi ^{})}}={\displaystyle \frac{a_K^2+2a_K\lambda ^2R_b\mathrm{cos}\varphi _3}{a_K^2}},`$ (20)
$`R_\pi ={\displaystyle \frac{\mathrm{Br}(B^0K^\pm \pi ^{})}{\mathrm{Br}(B^0\pi ^\pm \pi ^{})}}={\displaystyle \frac{a_K^2+2a_K\lambda ^2R_b\mathrm{cos}\varphi _3}{\lambda ^2R_b[R_b+2a_\pi (R_b\mathrm{cos}\varphi _3)]}},`$ (21)
with the Wolfenstein parameters $`\lambda `$ and $`R_b`$, and the Wilson coefficients

$`a_1`$ $`=`$ $`C_2+{\displaystyle \frac{C_1}{N_c}},`$ (22)
$`a_4`$ $`=`$ $`C_4+{\displaystyle \frac{C_3}{N_c}}+{\displaystyle \frac{3}{2}}e_q\left(C_{10}+{\displaystyle \frac{C_9}{N_c}}\right),`$ (23)
$`a_6`$ $`=`$ $`C_6+{\displaystyle \frac{C_5}{N_c}}+{\displaystyle \frac{3}{2}}e_q\left(C_8+{\displaystyle \frac{C_7}{N_c}}\right),`$ (24)
$`N_c`$ being the number of colors. The factor

$`a_{K(\pi )}={\displaystyle \frac{a_4+2r_{K(\pi )}a_6}{a_1}},r_{K(\pi )}={\displaystyle \frac{M_{K(\pi )}^2}{M_B(m_{s(u)}+m_d)}},`$ (25)
being negative, represents the ratio of the penguin contribution to the tree contribution in the $`K\pi (\pi \pi )`$ mode, where $`M_{K(\pi )}`$ is the kaon (pion) mass and $`m_q`$, $`q=s`$, $`u`$, and $`d`$, the $`q`$ quark masses.

It is obvious that the data $`R1`$  imply the unitarity angle $`\varphi _390^o`$, no matter what $`a_K`$, $`\lambda `$ and $`R_b`$ are. Therefore, $`R`$ is an appropriate quantity for the determination of $`\varphi _3`$. While to determine $`\varphi _3`$ from the data of the ratio $`R_\pi 4`$ , one must have precise information of $`a_K`$, $`a_\pi `$, $`\lambda `$ and $`R_b`$. It can be shown that the extraction of $`\varphi _3`$ from $`R_\pi `$ depends on these parameters sensitively. Hence, $`R_\pi `$ is not an appropriate quantity for the determination of $`\varphi _3`$. To explain $`R_\pi 4`$, a large $`|a_\pi |0.1`$ corresponding to $`m_d=2m_u=3`$ MeV and a large $`\varphi _3130^o`$ must be adopted in the FA approach . This is obvious from Eq. (21), since a large $`\varphi _3`$ leads to a constructive interference between the first and second terms in the numerator of $`R_\pi `$. In the modified FA approach with effective number of colors $`N_c^{\mathrm{eff}}`$, a large unitarity angle $`\varphi _3105^o`$ is also concluded .

It is interesting to ask whether one can explain $`R_\pi 4`$ using reasonable quark masses $`m_d=2m_u10`$ MeV and a smaller $`\varphi _390^o`$, if the extraction from $`R`$ is given a higher weight. The answer is positive in the PQCD approach. An essential difference between the FA and BBNS approaches and the PQCD approach is that transition form factors can be calculated in the latter as argued above. Therefore, we do not assume that the form factors associated with the Wilson coefficients $`a_{1,4}`$ and with $`a_6`$ are the same. An explicit PQCD calculation indicates that the ratios $`a_{K,\pi }`$ are more complicate, which reduce to Eq. (25) only in some extreme kinematic conditions of partons . Moreover, the choice of the hard scale $`t`$, at which Wilson coefficients are evaluated, is also different. It is chosen arbitrarily as $`m_b`$ or $`m_b/2`$ in the FA and BBNS approaches, but as virtuality of internal particles of a hard amplitude in PQCD, such that evolution effects vary among the above form factors.

It has been noticed that $`a_1(\mu )`$, $`a_4(\mu )`$ and $`a_6(\mu )`$ have dramatically different $`\mu `$ dependences for $`\mu <m_b/2`$: $`|a_4|`$ and $`|a_6|`$ increase with the decrease of $`\mu `$ faster than $`a_1`$. The choice of the scale in PQCD makes possible that $`t`$ runs to below $`m_b/2`$, and that penguin contributions are enhanced over tree contributions. It has been observed that the evolution effect enhances the ratios $`|a_{K,\pi }|`$ by about 50% . Hence, even adopting $`m_d=2m_u10`$ MeV, we can account for $`R_\pi 4`$ using $`\varphi _390^o`$. That is, the data of $`R_\pi `$ do not demand large $`\varphi _3`$. We emphasize that the hard scale $`t`$ with an average 1.4 GeV, corresponding to $`\alpha _s(t)/\pi 0.13`$, is still large enough to guarantee the applicability of PQCD .

## IV NONFACTORIZABLE CONTRIBUTIONS

Nonfactorizable contributions are neglected in the FA approach. However, it has been known that this approximation is not always appropriate, especially for the decay modes, whose factorizable amplitudes are proportional to the Wilson coefficient $`a_2=C_1+C_2/N_c`$. Since $`a_2`$ is small, factorizable contributions and nonfactorizable contributions are of the same order, and the latter become important. This is the reason it is difficult to resolve the controversy in the branching ratios of the $`BJ/\psi K^{()}`$ decays in FA . The conventional way to introduce nonfactorizable contributions is to vary number of colors in the Wilson coeficients $`a_i`$. The resultant effective number of colors, $`N_c^{\mathrm{eff}}`$, being process-dependent parameters, render the FA approach less predictive.

To go beyond FA, it has been proposed to evaluate nonfactorizable amplitudes for two-body charmless $`B`$ meson decays using PQCD factorization theorem in the BBNS approach . The motivation is that the power divergence appearing in $`F^{B\pi }`$ does not exist in nonfactorizable amplitudes, since soft divergences cancel between the two diagrams in Fig. 2 . Hence, these diagrams contribute to the hard parts of the $`B\pi \pi `$ decays, and are calculable in PQCD. We emphasize that the above soft cancellation is not sufficient to justify the applicability of PQCD. No matter what value $`M_B`$ is, one can always obtain the soft cancellation between the two nonfactorizable diagrams. However, it is obvious that PQCD is not applicable to $`B`$ meson decays if $`M_B`$ is as low as 1 GeV. The criterion for the applicability of PQCD is stronger: PQCD is applicable only when higher-order corrections to decay amplitudes are under control.

We shall show explicitly, following the reasoning in , that the simple factorization formula in  for a nonfactorizable amplitude suffers large higher-order corrections. The expression is

$$f=C\alpha _s\frac{dx_1}{x_1}\varphi _B(x_1)\left[\frac{dx}{x}\varphi _\pi (x)\right]^2,$$
(26)
where the constant $`C`$ is related to a color factor, the denominator $`x_1x`$ comes from the hard gluon propagator, and another denominator $`x`$ comes from the internal quark propagator. If the hard part is characterized by $`M_B`$, the coupling constant $`\alpha _s`$ evaluated at the scale $`\mu M_B`$ is small and higher-order corrections are negligible. However, it is not the case. When including higher-order corrections that cause the running of $`\alpha _s`$, the logarithm $`\alpha _s(\mu )\mathrm{ln}(x_1xM_B^2/\mu ^2)`$ occurs, with $`x_1xM_B^2`$ being the invariant mass of the hard gluon. In order to diminish this logarithmic correction, the appropriate choice of the scale is $`\mu =\sqrt{x_1x}M_B`$. The coupling constant $`\alpha _s`$ in Eq. (26) becomes running and should move into the convolution relation, leading to

$$f=C\frac{dx_1}{x_1}\varphi _B(x_1)\frac{dx_2}{x_2}\varphi _\pi (x_2)\frac{dx_3}{x_3}\varphi _\pi (x_3)\alpha _s(\sqrt{x_1x_3}M_B).$$
(27)
Assume the following models for the meson wave functions:

$`\varphi _B(x)`$ $`=`$ $`N_B\sqrt{x(1x)}\mathrm{exp}\left({\displaystyle \frac{M_B^2x^2}{2\omega _B^2}}\right),`$ (28)
$`\varphi _B^{}(x)`$ $`=`$ $`N_B^{}x^2(1x)^2\mathrm{exp}\left({\displaystyle \frac{M_B^2x^2}{2\omega _B^2}}\right),`$ (29)
$`\varphi _\pi (x)`$ $`=`$ $`{\displaystyle \frac{\sqrt{6}}{2}}f_\pi x(1x).`$ (30)
The $`B`$ meson wave functions $`\varphi _B`$  and $`\varphi _B^{}`$  possess peaks at small $`x`$ that depend on the shape parameters $`\omega _B`$ and $`\omega _B^{}`$, respectively. According to , the values of the shape parameters can be fixed via the parameter $`\lambda _B=0.3`$ GeV:

$$\frac{_0^1(dx/x)\varphi _B(x)}{_0^1𝑑x\varphi _B(x)}=\frac{M_B}{\lambda _B}=17.6,$$
(31)
which gives $`\omega _B=0.65`$ GeV and $`\omega _B^{}=0.25`$ GeV. The normalization constant $`N_B`$ is related to the $`B`$ meson decay constant $`f_B`$ through $`_0^1\varphi _B(x)𝑑x=f_B/(2\sqrt{2N_c})`$. The asymptotic pion wave function in Eq. (30) has been adopted in .

Using the wave function $`\varphi _B`$($`\varphi _B^{}`$), we find that only about 20% (30%) of the full contribution to the nonfactorizable amplitude comes from the perturbative region with $`\alpha _s/\pi <0.3`$. To derive the above percentages, we have chosen $`\mathrm{\Lambda }_{\mathrm{QCD}}=250`$ MeV for the running $`\alpha _s`$, and frozen the scale $`\sqrt{x_1x_3}M_B`$ to 0.4 GeV, as it goes to below 0.4 GeV. We conclude that Eq. (27) is dominated by soft contribution from the end-point region with $`x_1,x0`$, i.e., with $`\alpha _s>1`$, where the perturbative expansion is not justified. Therefore, nonfactorizable contributions can not be evaluated reliably in the BBNS approach. This conclusion is consistent with that drawn in : PQCD without Sudakov suppression is not applicable to exclusive processes for an energy scale below 10 GeV. The formalism in  is self-consistent as $`M_B\mathrm{}`$, but not in the case with $`M_B5`$ GeV.

In the end-point region where soft contributions are important, nonfactorizable amplitudes should be treated in a more careful way. That is, the dependence on the parton transverse momentum $`k_T`$ is essential, and should be included as argued in the previous section. The argument of $`\alpha _s`$ is also set to $`t`$, the maximal virtuality of internal particles of a hard amplitude, which is similar to that in Eq. (27), but depends on the additional scale $`1/b`$. $`\alpha _s(t)`$ is then always small even in the end-point region, since $`1/b`$ is large under Sudakov suppression . It has been shown that almost 100% of full contributions to nonfactorizable amplitudes arise from the region with $`\alpha _s(t)/\pi <0.3`$. That is, the evaluation of nonfactorizable contributions is reliable in the PQCD approach.

There is another difference between the BBNS and PQCD approaches in the treatment of nonfactorizable contributions. The momentum of the light spectator quark in the $`B`$ meson has been ignored in the former , such that quark propagators in the hard part always remain time-like:

$$\frac{1}{(k_2+k_3k_1)^2+iϵ}\frac{1}{(k_2+k_3)^2}=\frac{1}{x_2x_3M_B^2}.$$
(32)
This approximation drops the strong phase from the kinematic region where the virtual quark becomes on-shell, $`(k_2+k_3k_1)^20`$. In the PQCD approach nonfactorizable contributions are complex:

$$\frac{1}{(k_2+k_3k_1)^2+iϵ}=P\left[\frac{1}{(k_2+k_3k_1)^2}\right]i\pi \delta \left[(k_2+k_3k_1)^2\right],$$
(33)
where $`P`$ denotes the principle-value prescription. As shown in , the above imaginary part is not the main source of strong phases for the $`B\pi \pi `$ decays, because nonfactorizable contributions are only few percents of factorizable ones. However, it will become important in decay modes whose factorizable amplitudes are proportional to the small Wilson coefficient $`a_2`$, or absent, such as the $`B_d^0K^\pm K^{}`$ decays .

## V STRONG PHASES

In the FA and BBNS approaches annihilation contributions have been neglected. Though annihilation amplitudes from the operator $`O_4`$ for the $`B\pi \pi `$ decays indeed vanish because of helicity suppression, those from $`O_6`$ do not. As shown in Eqs. (B1) and (B6) of , there is a large enhancing factor

$$\frac{M_B^2}{(m_b+m_d)(m_u+m_d)},$$
(34)
associated with annihilation contributions from $`O_6`$, which is much larger than the one

$$\frac{M_\pi ^2}{(m_b+m_d)(m_u+m_d)},$$
(35)
associated with usual penguin contributions. Though annihilation contributions are suppressed by a factor $`f_B/M_B`$ , they can be comparible with penguin contributions. It has been observed that the annihilation diagrams Figs. 3(e) and 3(f) contribute large imaginary parts (strong phases) in the PQCD approach .

In the approaches where annihilation contributions are neglected, strong phases must be introduced via the Bander-Silverman-Soni (BSS) mechanism shown in Fig. 4(a) . However, this mechanism, contributing through a charm quark loop, is suppressed by the charm mass threshold :

$$q^2=(x_2P_2+x_3P_3)^2=x_2x_3M_B^2>4m_c^2,$$
(36)
$`m_c`$ being the $`c`$ quark mass. Obviously, the large $`x_2`$ and $`x_3`$ regions are not favored by the pion wave funtion in Eq. (30). . Moreover, the diagram with a charm quark loop is of higher order in the PQCD formalism. An exact numerical analysis has indicated that the imaginary contribution from the BSS mechanism is smaller than that from the annihilation diagrams by a factor 10 .

Another source of strong phases in the FA and BBNS approaches is the diagrams in Fig. 4(b), whose imaginary contributions appear in the effective Wilson coefficients in Eqs.(4)-(8) of , for instance,

$`a_1^u`$ $`=`$ $`C_1+{\displaystyle \frac{C_2}{N}}+{\displaystyle \frac{\alpha _s}{4\pi }}{\displaystyle \frac{C_F}{N}}C_2F,`$ (37)
$`F`$ $`=`$ $`12\mathrm{ln}{\displaystyle \frac{\mu }{m_b}}18+f_\pi ^I+f_\pi ^{II},`$ (38)
with

$$f_\pi ^I=0.53\pi i,$$
(39)
and $`f_\pi ^{II}`$ being real. We emphasize that Eq. (39) is derived from a configuration, in which the outgoing $`\overline{u}`$ quark carries the full pion momentum. Again, this configuration is strongly suppressed by the pion wave function and by the Sudakov factor. That is, the imaginary contribution from the above source has been overestimated.

Below we make a quantative comparision among the FA, BBNS and PQCD approaches. Figure 3 collects all the Feynman diagrams for the hard amplitude of the $`B_d^0\pi ^+\pi ^{}`$ decay, including factorizable and nonfactorizable tree, penguin and annihilation diagrams, whose contributions are listed in Table I. The results obtained in the BBNS approach are also shown. The FA results are the same as the BBNS ones, but with the nonfactorizable amplitudes neglected. We observe

1. The tree contributions are almost equal in the three approaches, but the penguin contributions are larger in PQCD, corresponding to a larger ratio $`|a_\pi |`$ defined by Eq. (25) (about twice of those in FA and BBNS), i.e., corresponding to the dynamical penguin enhancement.
2. The tree and penguin contributions are real in PQCD, but complex in FA and BBNS due to the complex effective Wilson coefficients, which are the main source of strong phases in FA and BBNS.
3. The annihilation contributions are neglected in FA and BBNS, but of the same order as the penguin contributions in PQCD, which are the main source of strong phases in PQCD.
4. Nonfactorizable contributions are complex in PQCD but real in BBNS, as explained in the previous section, and neglected in FA. Note that different $`B`$ and $`\pi `$ meson wave functions have been employed as deriving the results in Table I . Nonfactorizable contributions are destructive relative to the penguin contributions, if concering only the real parts, and negligible for the $`B\pi \pi `$ branching ratios in both approaches.
5. The dominant nonfactorizable contribution $`M^P`$ in BBNS is due to the sharp hard parts in the end-point regions of momentum fractions $`x`$ as stated in Sec. III.
6. The total tree and penguin contributions, including both factorizable, nonfactorizable and annihilation ones, are

$`T_{\mathrm{PQCD}}=(71.87+3.20i)\times 10^3,P_{\mathrm{PQCD}}=(7.47+4.64i)\times 10^3,`$ (40)
$`T_{\mathrm{BBNS}}=(74.93+1.13i)\times 10^3,P_{\mathrm{BBNS}}=(1.131.27i)\times 10^3,`$ (41)
$`T_{\mathrm{FA}}=(75.70+1.13i)\times 10^3,P_{\mathrm{FA}}=(3.041.27i)\times 10^3.`$ (42)

The tree contributions obtained in the three approaches are very close. The signs of the penguin contributions are different, and the magnitude in PQCD is larger than those in FA and BBNS:

$$|P_{\mathrm{PQCD}}|5.2|P_{\mathrm{BBNS}}|,|P_{\mathrm{PQCD}}|2.7|P_{\mathrm{FA}}|.$$
(43)
The difference in the total penguin contributions is attributed to the treatment of nonfactorizable and annihilation duagrams. We emphasize that the BBNS results depend on the parameters for factorizable contributions (the transition form factors) and the parameters for nonfactorizable contributions (the meson wave funcitons). Therefore, it is possible to reverse the sign and to increase the magnitude of the penguin contribution by varying these parameters. The penguin contribution in the BBNS approach is smallest because of the cancellation between $`F^P`$ and $`M^P`$. If expressing the amplitudes of the $`B_d^0\pi ^+\pi ^{}`$ decay as

$$𝒜=V_t^{}|P|e^{i\delta }V_u^{}|T|,$$
(44)
with $`P`$ and $`T`$ given in Eq. (42), the strong phases $`\delta `$ are about

$$\delta _{\mathrm{PQCD}}123^o,\delta _{\mathrm{BBNS}}50^o,\delta _{\mathrm{FA}}113^o.$$
(45)
The predictions for the CP asymmetry of the $`B_d^0\pi ^\pm \pi ^{}`$ decays are presented in Fig. 5. Because of Eq. (43), the CP asymmetry in PQCD is larger than that ($`10\%`$ at most) in FA and much larger than that ($`4\%`$ at most) in BBNS. Because of Eq. (45), the sign of the CP asymmetry in PQCD is opposite to those in FA and BBNS. If the $`B`$ meson wave function in  with $`\omega _B=0.4`$ GeV is adopted, $`M^P`$ and also $`|P_{\mathrm{BBNS}}|`$ will increase accordingly, such that the CP asymmetry in the BBNS approach becomes a bit larger. We shall leave the resolution of the above discrepancies to future measurements of CP asymmetries. That is, experimental data of CP asymmetries will provide an appropriate test of the three approaches. Note that all the three approaches give the similar branching ratio $`B(B_d^0\pi ^\pm \pi ^{})4.6\times 10^6`$, which is consistent with the CLEO data .

As argued in Sec. III, to account for the $`BK\pi `$ and $`\pi \pi `$ data simulataneously using the unitarity angle $`\varphi _3=90^o`$, the ratio of the penguin contribution to the tree contribution must be as large as 0.1. It is easy to find

$$|P_{\mathrm{PQCD}}|/|T_{\mathrm{PQCD}}|=0.12,|P_{\mathrm{BBNS}}|/|T_{\mathrm{BBNS}}|=0.02,|P_{\mathrm{FA}}|/|T_{\mathrm{FA}}|=0.04,$$
(46)
from Eq. (42). Due to the smaller ratio in the BBNS and FA approaches, a larger angle $`\varphi _3>90^o`$ must be adopted in order to explain the $`B\pi \pi `$ data . As a consequence of the penguin enhancement, the CP asymmetry is also large in the PQCD approach.

## VI SUMMARY

We summarize the comparision among the FA, BBNS, and PQCD approaches as follows:

1. The scale independence of predictions in the FA and BBNS approaches is achieved by the generalized FA method: a $`\mu `$-dependent evolution factor is extracted from hadronic matrix elements by considering the mixing of the effective operators at a scale below $`\mu `$. This evolution factor is then combined with the $`\mu `$-dependent Wilson coefficient to form a scale-independent effective Wilson coefficient. However, the infrared-cutoff and gauge dependences appear. In the PQCD approach the scale independence is made explicit by extending the Wilson evolution down to the hard scale $`t`$ directly, which is then followed by the second-stage RG evolution evaluated based on six-quark amplitudes. Since external quarks are on-shell, the gauge invariance of PQCD predictions is guaranteed. The infrared-cutoff dependence, treated as a factorization scheme dependence, is removed by the universality of wave functons.
2. Factorizable contributions (hadronic transition form factors) are not claculable in the FA and BBNS approaches, but calculable in the PQCD approach. The inclusion of parton transverse momenta smear the end-point singularity, and Sudakov suppression guarantees that almost all contributions come from the region with $`\alpha _s/\pi <0.3`$.
3. Dynamical penguin enhancement exists in the PQCD approach, but not in the FA and BBNS approaches. This enhancement, arising from different evolution effects in the form factors associated with the Wilson coefficient $`a_1`$, $`a_4`$ and $`a_6`$, is cruical for the simultaneous explanation of the $`BK\pi `$ and $`\pi \pi `$ data using a smaller unitarity angle $`\varphi _390^o`$ .
4. Nonfactorizable amplitudes are neglected in the FA approach, but calculable in the BBNS and PQCD approaches. However, without Sudakov suppression, higher-order corrections are not under control in the BBNS approach. The end-point singularity needs to be smeared by taking into account parton transverse momenta and Sudakov suppression. On the other hand, nonfactorizable amplitudes are real in the BBNS approach, but complex in the PQCD approach. The conclusion is the same: nonfactorizable contributions are destructive relative to the penguin contributions and negligible for the evaluation of branching ratios.
5. Annihilation contributions are neglected in the FA and BBNS approaches, but calculable in the PQCD appraoch. The helicity suppression applies only to amplitudes associated with $`O_4`$, but not to those associated with $`O_6`$. It is found that the $`O_6`$ contribution is of the same order as the penguin contribution, and gives a large imaginary part.
6. Strong phases come from the BSS mechanism and the extraction of the $`\mu `$ dependence from hadronic matrix elements in the FA and BBNS approaches, but from annihilation diagrams in the PQCD approach. It has been argued that the two sources of strong phases in the FA and BBNS approaches are in fact strongly suppressed by the charm mass threshold and by the end-point behavior of meson wave functions.
7. Because of the large penguin and annihilation contributions in the PQCD approach, the predicted CP asymmetry in the $`B_d^0\pi ^\pm \pi ^{}`$ decays dominates over those in the FA and BBNS approaches as shown in Fig. 5. Future measurements of CP asymmetries can distinguish the FA, BBNS, and PQCD approaches.

Acknowledgment

We thank M. Bando, Y. Sumino and M. Yamaguchi for organizing the Summer Institute 99 at Yamanashi, Japan, whose active atmosphere stimulated this work, H.Y. Cheng, T. Mannel and A.I. Sanda for helpful discussions, and M. Kobayashi for his encouragement. The works of HNL and YYK were supported in part by the National Science Council of R.O.C. under the Grant Nos. NSC-89-2112-M-006-004 and NSC-89-2811-M-001-0053, respectively, and in part by the Monbusho visiting fellowship from the Ministry of Education, Science and Culture of Japan.

Figure Captions

1. Fig. 1: One of the lowest-order factorizable diagrams for the $`B_d^0\pi ^+\pi ^{}`$ decay.
2. Fig. 2: Lowest-order nonfactorizable diagrams for the $`B_d^0\pi ^+\pi ^{}`$ decay.
3. Fig. 3: Lowest-order hard amplitudes for the $`B_d^0\pi ^+\pi ^{}`$ decay.
4. Fig. 4: Diagrams that give strong phases in the FA and BBNS approaches.
5. Fig. 5: Direct CP asymmetries in the $`B_d^0\pi ^\pm \pi ^{}`$ decays.


format: text

text:
# A pot of gold at the end of the cosmic “raynbow”?

## Abstract

We critically review the common belief that ultrahigh energy cosmic rays are protons or atomic nuclei with masses not exceeding that of iron. We find that heavier nuclei are indeed possible, and discuss possible sources and acceleration mechanisms for such primaries. We also show detailed simulations of extensive air showers produced by “superheavy” nuclei, and discuss prospects for their detection in future experiments.

The unambiguous detection of cosmic rays (CRs) with energies above 1020 eV (see  for a survey and bibliography on the subject) is a fact of outstanding astrophysical interest. As shown in the pioneering works of Greisen, Zatsepin, and Kuzmin , the possible sources and the accelerating mechanisms are constrained by the observed particle spectra due to the interaction with the universal radiation and magnetic fields on the way to the observer. The low flux of particles at the end of the spectrum (the typical rate of CRs above 1020 eV is one event/km2/century) puts strong demands on the collection power of the experiments, such as can only be achieved by extended air shower detection arrays at ground level. This indirect method of detection bears a number of serious difficulties in determining the energy, mass and/or arrival direction of the primary particles.

Astrophysical mechanisms to accelerate particles to energies of up to $`10^{2122}`$ eV have been identified, but they require exceptional sites . Very recently, we have presented a comprehensive study of a possible nearby superheavy-nucleus-zevatron . We have shown that it is likely that nuclei heavier than iron with energies above a few PeV can escape from the dense core of a nearby starburst galaxy like M82, and eventually be re-accelerated to superhigh energies ($`E10^{20}`$ eV) at the terminal shocks of galactic superwinds generated by the starburst.111It is important to stress that M82 is positioned close to the arrival direction of the highest CR event detected on Earth. This was first pointed out in . This mechanism improves as the charge number $`Z`$ of the particle is increased. Furthermore, we have also shown that the nuclei may arrive on Earth. Strictly speaking, the energetic nucleus is seen to lose energy mainly as a result of its photodisintegration. In the universal rest frame (in which the microwave backgroud radiation is at $`3K`$), the disintegration rate $`R`$ of an extremely high energy nucleus with Lorentz factor $`\mathrm{\Gamma }`$, propagating through an isotropic soft photon background of density $`n`$ is given by ,

$$R=\frac{1}{2\mathrm{\Gamma }^2}_0^{\mathrm{}}𝑑E\frac{n(E)}{E^2}_0^{2\mathrm{\Gamma }E}𝑑E^{}E^{}\sigma (E^{}),$$
(1)
where primed quantities refer to the rest frame of the nucleus, and $`\sigma `$ stands for the total photon absortion cross section. Above $`10^{20}`$ eV, the energy losses are dominated by collisions with the relic photons. The fractional energy loss around this energy is $`R10^{15}`$ s-1. With this in mind, it is straightforward to check that superheavy nuclei may impact on Earth (for details see Fig. 2 of Ref. ). In the rest of this report we shall discuss the characteristics of the extensive air showers that these nuclei may produce after interaction with the atmosphere.

Golden Shower Simulations: In order to perform the simulations we shall adopt the superposition model. This model assumes that an average shower produced by a nucleus with energy $`E`$ and mass number $`A`$ is indistinguishable from a superposition of $`A`$ proton showers, each with energy $`E/A`$. We have generated several sets of 197Au air shower simulations by means of the aires Monte Carlo code .222It should be stressed that for $`A>140`$ the bulk solar–system abundance distribution peaks at $`A=195`$ . To make some estimates, we then refer our calculations to a gold nucleus. The sample was distributed in the energy range of $`10^{18}`$ up to $`10^{20.5}`$ eV. sibyll was used to reproduce hadronic collisions above 200 GeV . All shower particles with energies above the following thresholds were tracked: 750 keV for gammas, 900 keV for electrons and positrons, 10 MeV for muons, 60 MeV for mesons and 120 MeV for nucleons and nuclei. The particles were injected vertically at the top of the atmosphere (100 km.a.s.l), and the surface detector array was put at a depth of 1036 g/cm2, i.e., at sea level. Secondary particles of different types and all charged particles in individual showers were sorted according to their distance $`R`$ from the shower axis.

In Fig. 1 we show the lateral distributions of different groups of secondary particles (we have considered separately $`\gamma `$, $`e^+e^{}`$, and $`\mu ^+\mu ^{}`$). One can see that the number of muons from the gold nucleus shower is greater than the number of muons from the proton shower.

As the cascade develops in the atmosphere, it grows until a maximum size (number of particles) is reached. The location in the atmosphere where the cascade has developed the maximum size is denoted by $`X_{\mathrm{max}}`$, with units of g cm-2. For cascades of a given total energy, heavier nuclei have smaller $`X_{\mathrm{max}}`$ than nucleons because the shower is already subdivided into $`A`$ nucleons when it enters the atmosphere. At $`10^{20}`$ eV, the $`<X_{\mathrm{max}}>`$ of a proton (gold) shower is $`879`$ g/cm2 ($`777`$ g/cm2). A dust-grain has an even larger cross section, so it tends to interact sooner than protons and nuclei . In Fig. 2, we compare the longitudinal profile of showers initiated by a proton, a gold-nucleus and a dust-grain. It is clearly seen how the $`X_{\mathrm{max}}`$ decreases when increasing the mass. The simulated gold shower is partially consistent with the Fly’s Eye data. Furthermore, its longitudinal development better reproduces the data than protons or dust-grains. It should be remarked, however, that extensive air shower simulation depends on the hadronic interaction model . We also point out that for the simulation detector effects were not taken into account.

Even though the superheavy nucleus hypothesis is partially supported by data from the CASAMIA experiment , more data is certainly needed to verify this model. In order to significantly increase the statistics at the end of the spectrum, the Southern Auger Observatory is currently under construction . It will consist of a surface array which will record the lateral and temporal distribution of shower particles, and an optical air fluorescence detector, which will observe the air shower development in the atmosphere. These two techniques provide complementary methods of extracting the required information from the shower to test the ideas discussed in this paper.

This work was partially supported by CONICET, Fundación Antorchas and the National Science Foundation.


text:
# Probing the gravitational well: No supernova explosion in spherical symmetry with general relativistic Boltzmann neutrino transport

## I Back to the Past?

In the pioneering days of supernova simulations (Colgate and White , May and White ) the consideration of general relativity (GR) was standard. The observable event of a stellar core collapse followed by a supernova explosion was an ideal application for the newly derived Einstein equations in spherical symmetry (Misner and Sharp ). Lindquist  formulated the GR Boltzmann equation, and Wilson  carried out simulations based on GR Boltzmann neutrino transport using parametrized neutrino interactions. This early epoch laid the foundation leading to our current understanding of the supernova mechanism: A collapsing stellar iron core bounces at nuclear densities and launches a shock wave outwards through the infalling outer layers. The shock is energized by neutrinos diffusing out of the hot proto-neutron star, which deposit a fraction of their energy in the shock-dissociated matter via absorption on the free nucleons.

The early models underwent refinements in many respects. The equation of state evolved from simple polytropes to significantly more realistic models (Baron et al., Lattimer and Swesty ). The neutrino opacities were improved (Tubbs and Schramm , Schinder and Shapiro , Bruenn ), and sophisticated multigroup flux-limited diffusion neutrino transport schemes were developed based on this “standard” nuclear physics input (Arnett , Bowers and Wilson , Bruenn , Myra et al.).

However, improved approximations in the implementation of neutrino physics seemed to decrease the likelihood of successful explosions in spherical symmetry. This was especially true with the inclusion of the detailed neutrino energy spectrum via multigroup simulations and for simulations that took the inelastic neutrino–electron scattering into account (Bruenn ). The idea that the amount of dissociation energy expended by the shock ploughing through the outer iron core might be small enough to allow a prompt explosion had to be abandoned. The energy that is deposited behind the shock by the outstreaming neutrinos gained increased attention in light of a possible delayed shock revival (Wilson , Bethe and Wilson ). The amount of energy transfered, however, does not alone depend on the neutrino spectrum and luminosity of the neutrinos, but also on their angular distribution. The neutrino propagation angle determines the path length in a given radial shell; therefore, neutrinos with more tangential directions have a higher chance of absorption. A known problem with the multigroup flux-limited diffusion approximation to the Boltzmann transport equation is that, in the important semi-transparent region between the neutrinosphere and the heating region, the angular distribution of the neutrinos is not self-consistently determined by the transport equation. The consequence is an underestimation of the isotropy of the outstreaming neutrinos (Janka ).

The focus therefore switched on the one hand to the development of complete three-flavor Boltzmann neutrino transport (Mezzacappa and Bruenn ). In stationary comparisons between full transport and multigroup flux-limited diffusion approximations (Messer et al. , Yamada et al. ) more efficient heating was found with accurate transport. On the other hand, multidimensional phenomena were explored in the hope of finding more efficient energy transfer by neutrinos. The inclusion of neutron finger convection (see also Ref. ) produced explosions (Wilson and Mayle ) and established the delayed explosion scenario. Two-dimensional investigations of convection behind the shock and in the proto-neutron star followed with mixed results (Herant et al., Miller et al., Herant et al., Burrows et al., Janka and Müller , Keil et al., Mezzacappa et al.).

Semi-analytical investigations (Burrows and Goshy , and recently Janka ) illuminate the basic mechanisms, but are not able to self-consistently decide for or against explosions nor able to predict detailed data in either of these cases.

In the search for a robust supernova mechanism, attention was directed toward simulations in the nonrelativistic (NR) limit, because explosions in GR seemed to be less likely in the selective picture that the deeper shock formation would produce larger dissociation losses and that the neutrino luminosities would suffer gravitational redshift. Multi-dimensional simulations that approximate general relativistic effects were performed by Fryer  and Fryer and Heger  who found reduced convection with rotating progenitors.

However, there are also beneficial GR effects: Baron et al.  reported prompt explosions in general relativity with a very soft equation of state and a leakage scheme, although these explosions were not reproduced in a general relativistic multigroup flux-limited diffusion simulation by Myra et al.  when neutrino–electron scattering was included. Swesty et al.  also did not find an explosion in a systematic investigation of realistic parameters in the equation of state. Late-time beneficial effects from the decrease in the gravitational potential by neutrino energy radiation were suggested by Goldman and Nussinov . De Nisco et al.  in a quasi-static investigation pointed out that, among many detrimental effects, the hotter core in general relativistic hydrodynamics produces higher neutrino luminosities with harder spectra, resulting in a potentially beneficial impact on the heating rate. Recent dynamical simulations with GR multigroup flux-limited diffusion confirm this effect, although it does not appear to be strong enough to outweigh the negative GR contributions, i.e. the smaller heating region and greater infall velocities (Bruenn et al. ).

Recently, postbounce evolution was reexamined with Boltzmann neutrino transport without invoking multidimensional effects. As the result of an undervalued nucleon isoenergetic scattering opacity, first simulations led to an explosion of a $`13`$ M progenitor (Liebendörfer ). Simulations in Newtonian gravity with $`O(v/c)`$ Boltzmann transport and standard nuclear physics for the same progenitor do lead to an enhanced shock radius, but not to an explosion (Mezzacappa et al.). Even with the omission of the energy loss due to the escape of $`\mu `$\- and $`\tau `$ neutrinos, independent simulations of the postbounce evolution of a $`15`$ M progenitor do not produce an explosion (Rampp and Janka ). The latter simulations are based on a tangent-ray method for the $`O(v/c)`$ Boltzmann neutrino transport (see also Burrows et al. ).

It is the purpose of this paper to report on the completion of a general relativistic radiation hydrodynamics code for spherically symmetric flows that solves the detailed Boltzmann transport equation, and to present the self-consistent simulation of stellar core collapse, bounce, and postbounce evolution for a $`13`$ M star, with all neutrino flavors and “standard” nuclear physics included.

## II The Recipe

Our simulation is initiated from the stellar precollapse model of Nomoto and Hashimoto . We use the Lattimer-Swesty equation of state . Simplified silicon burning is included as material from the silicon layer surrounding the iron core is instantaneously burned under energy conservation to nuclear statistical equilibrium as soon as the temperature exceeds $`0.44`$ MeV (Mezzacappa et al. ). The simulations were carried out with a new general relativistic neutrino radiation hydrodynamics code, AGILE-BOLTZTRAN, based on a conservative formulation of GR radiation hydrodynamics in spherical symmetry and comoving coordinates (Liebendörfer et al. ).

AGILE is an implicit GR hydrodynamics code using an adaptive grid technique to conservatively implement shift vectors (Liebendörfer and Thielemann ). It maximizes, with respect to the number of required grid points, the resolution in regions with large density gradients and allows a smooth propagation of the shock through the outer layers (see graph (a) in Figs. (1-5)). BOLTZTRAN is an implicit three-flavor multigroup Boltzmann neutrino transport solver (Mezzacappa and Bruenn , Mezzacappa and Messer ) that was consistently coupled to AGILE, enabled for adaptive gridding, and extended to GR flows under omission of the small gravitational neutrino backreaction (Liebendörfer ). We choose $`103`$ spatial zones, and discretize the neutrino-momentum phase space, as in Mezzacappa et al. , with $`6`$-point Gaussian quadrature and $`12`$ neutrino energy groups, ranging from $`5`$ MeV to $`300`$ MeV.

We emphasize at this point that the energy evolution in Lagrangian radiation hydrodynamics obeys the very simple conservation equation

$`{\displaystyle \frac{}{t}}\left[\mathrm{\Gamma }e+{\displaystyle \frac{2}{\mathrm{\Gamma }+1}}\left({\displaystyle \frac{1}{2}}u^2{\displaystyle \frac{m}{r}}\right)+\mathrm{\Gamma }J+uH\right]`$ (1)
$`+{\displaystyle \frac{}{a}}\left[4\pi r^2\alpha up+4\pi r^2\alpha \rho \left(uK+\mathrm{\Gamma }H\right)\right]`$ $`=`$ $`0.`$ (2)
The variables are defined as in Lindquist  with the exceptions that $`J`$, $`H`$, and $`K`$ represent the zeroth, first, and second angular moments of the specific radiation energy; $`\rho `$, $`e`$, and $`p`$ the fluid rest mass density, internal energy, and pressure; and, $`a`$, $`r`$, and $`u`$ the enclosed rest mass, the radius, and the time derivative of the radius devided by the lapse function $`\alpha `$. The integration of Eq. (2) over the total rest mass $`a_{tot}=1.55`$ M in the computational domain provides a check on the conservation of total energy, which can easily be monitored. Energy conservation presents a significant challenge. The cancellation of the gravitational binding energy with the internal energy sets the scale of the total energy, which is two orders of magnitude smaller than either contribution separately and comparable to the $`10^{51}`$ erg ($`1\%`$ of radiation energy) expected to be deposited in the heating region. Our total energy drift is no more than $`5\times 10^{49}`$ erg in the first $`200`$ ms in which the explosion outcome is decided. Lepton number is conserved to within a fraction of a percent.

## III The Failed Explosion of a $`13`$ M Star

We discuss the evolution of our simulation at five time slices: At bounce, and at $`1`$ ms, $`10`$ ms, $`100`$ ms, and $`500`$ ms after bounce. At each time slice we show eight graphs that compare the general relativistic evolution (GR, fat lines) to the nonrelativistic evolution (NR, thin lines). We have already reported on the NR simulation in Mezzacappa et al. . The two runs are synchronized at bounce. The graphs are (a) the velocity profile and the location of the adaptive grid points (for the GR profiles only); (b) the rest mass density profile; (c) the entropy per baryon profile; (d) the electron fraction profile. The latter three quantities enter the equation of state and determine the composition and thermodynamic state of the local fluid element. For the convenience of interpretation, we additionally show the mass fraction of heavy nuclei, free nucleons, and alpha particles in graph (f). The remaining three graphs are devoted to neutrino radiation quantities: (e) the luminosities and (h) the rms energies defined by

$$E_{rms}=\left(\frac{E^2n_\nu (E)𝑑E}{n_\nu (E)𝑑E}\right)^{\frac{1}{2}}.$$
(3)
The integrand $`n_\nu (E)dE`$ is the number density of neutrinos in the comoving frame energy interval $`[E,E+dE]`$. Graph (g) shows the mean flux factor defined by the quotient of the neutrino energy flux, $`H`$, and the neutrino energy density, $`J`$: $`FH/(cJ)`$. We also show the energy-averaged neutrinospheres in this plot because they separate the interior diffusion regime, with small flux factors, from the exterior regime, where large flux factors describe increasingly forward-peaked neutrino propagation through semitransparent and free-streaming regions. AGILE-BOLTZTRAN solves directly for the specific neutrino distribution function, $`F(t,a,E,\mu )`$, at each time $`t`$. The angular resolution in the neutrino propagation angle cosine $`\mu `$ is set by the order of a Gaussian quadrature. This limits the representation of a strongly forward peaked radiation field, as can be observed in the mean flux factor at large radii in the graphs (g) of Fig. (1-5). However, our choice of 6-point quadrature allows an accurate representation of the mean flux factor out to a radius of at least four times the radius of the neutrinosphere. This range includes the cooling and heating region for all times in the simulation, and any influence of this numerical restriction on the dynamics is excluded. To define an energy-averaged neutrinosphere, we write the energy dependent optical depth at enclosed rest mass, $`a`$, as

$$\tau (a,E)=_a^{a_{tot}}\underset{i}{}\frac{\sigma _i(E)}{4\pi r^2}\frac{n_i}{\rho }da.$$
(4)
The variables are defined as in Eq. (2). The index $`i`$ labels different reactions, with cross sections $`\sigma _i`$ and target number densities $`n_i`$. We draw in graph (g) the energy-averaged neutrinospheres at the location $`r(a)`$ where the condition

$$\frac{\tau (a,E)n_\nu (a,E)𝑑E}{n_\nu (a,E)𝑑E}=\frac{2}{3}$$
(5)
for the enclosed mass $`a`$ holds. The radiation quantities are given for the electron neutrino and the electron antineutrino separately. The $`\mu `$ and $`\tau `$ neutrinos do not exchange energy with the heating region. Their luminosities and rms energies are shown in later figures. In our discussion, we follow the general relativistic simulation and supply numbers for the nonrelativistic case in parentheses.

The first prominent differences between the GR and NR simulations appear at bounce in Fig. (1). At this time, the inner core passes maximum compression with a central density of $`3.87\times 10^{14}`$ g/cm3 (NR: $`3.05\times 10^{14}`$ g/cm3) as shown in graph (b). The homologous collapse of the causally connected inner core is abruptly halted when the central region reaches nuclear density. The short range nuclear forces suddenly increase the stiffness of the equation of state. A pressure wave emerges from the center and steepens into a shock wave at the border of the inner iron core, where it faces the disconnectedly infalling material from the outer iron core. We observe the shock formation in the velocity profile in graph (a) at an enclosed mass of $`0.53`$ M (NR: $`0.65`$ M). The conversion of kinetic energy into internal energy at the shock front becomes apparent in the discontinuity of the entropy profile in graph (c). The additional discontinuity at $`1.2`$ M stems from continued silicon burning by compressional heating in the infalling material. From the composition in graph (f), we know that the shock starts to dissociate the heavy nuclei as they fall in. The shock in the NR case encloses a region of undissociated heavy nuclei at high density. The extent of this region is limited by the phase transition to bulk nuclear matter in the Lattimer-Swesty equation of state. These heavy nuclei are not present in the GR case because the shock is formed at a smaller radius, almost coincident with the location of the phase transition. The electron fraction profile in graph (d) exhibits values from $`0.499`$ in the outer layers down to $`0.29`$ (NR: $`0.28`$) in the central region. The lepton fraction is set by neutrino trapping during core collapse (Sato , Mazurek ). Thereafter, smaller changes in the electron fraction are allowed according to an establishing equilibrium between electron- and neutrino-capture at constant lepton fraction (see also Mezzacappa and Bruenn , Messer ). Determined by this equilibrium are the rms energies for the electron neutrinos, as shown in graph (h). The mean flux factor in graph (g) reflects an isotropic neutrino distribution in the innermost $`0.7`$ M and outwardly directed radiation at the border of the computational domain. The peaks in the luminosity in graph (e) are caused by numerical noise at the shock front. The small flux factor in graph (g) for the electron neutrinos shows that the short range energy flux at the peak is small in the sense that it does not exceed one percent of the prevailing neutrino energy density multiplied by the velocity of light. An influence on the dynamics is therefore excluded.

One millisecond later, in Fig. (2), the shock has moved from the edge of the homologous core to densities $`10^{12}`$ g/cm3 as shown in graph (b). In graph (f) we observe a region of dissociated material behind the shock. The heavy nuclei still present in the NR simulation between the shock heated matter and bulk nuclear matter have a mass of $`0.15`$ M. They leave the shock an additional energy $`2\times 10^{51}`$ erg with respect to the GR simulation. The NR shock in graph (a) still contains kinetic energy, while the GR shock has almost completely turned into an accretion shock. In graph (g) we see that the shock approaches the neutrinospheres. Plotted are the energy-averaged neutrinospheres located in the broad region where neutrinos with different energies emanate from the diffusion regime. As the shock passes this region at about $`4`$ ms after bounce, an energetic neutrino burst will be released from the hot shocked material, rendering it “neutrino-visible” to the outside world. A rise in the electron neutrino luminosity at the neutrinospheres at this time is manifest in graph (e). A dip in the rms energy profiles of the electron antineutrino in graph (h) at a radius $`11`$ km results from the fact that there are two distinct regions producing these neutrinos: the cold compressed center of the star and the hot shocked mantle.

The time slice in Fig. (3) corresponds to $`10`$ ms after bounce, after the launch of the neutrino burst. The luminosity peak in graph (e) has already propagated to a radius of $`1800`$ km. This is consistent with the time of shock breakout, $`t=10`$ ms - $`1800`$ km/$`c`$ $`4`$ ms. The propagation of the neutrinos is also visible in their breakout rms energies. In graph (h), the neutrinos from collapse, with rms energies $`10`$ MeV, are replaced with the burst neutrinos having an rms energy of almost $`15`$ MeV. The electron neutrinos in the burst were produced by electron capture on free protons. The corresponding deleptonization behind the shock is dramatic in graph (d). The energy carried off with the neutrino burst completely drains the shock in both the NR and the GR cases. The pure accretion shock continues to propagate outwards as infalling material is dissociated and layered on the hot mantle (graph (a)). This stage is the definitive end of a “prompt,” i.e., purely hydrodynamic, explosion.

In the standard picture, the ensuing evolution is driven by electron neutrino heating. Electron flavor neutrinos are diffusing out of the cold unshocked core as well as created in the accreting and compactifying matter around the neutrinospheres in a hot shocked mantle. On their escape from the semi-transparent region, they deposit energy behind the shock by absorption on free nucleons. This situation is well established at $`100`$ ms after bounce, corresponding to Fig. (4). At this time, the accretion shock has passed the border of the original iron core. In graph (a), it has stalled at a radius of roughly $`200`$ km. The GR shock shows a smaller stall radius than the NR shock, and the infall velocity behind the shock is about two times larger than in the NR case. The incomplete dissociation of the infalling silicon layer at the shock in graph (f) illustrates the weakening of the shock. Half of the mass remains bound in alpha particles when it crosses the shock and is dissociated later during the fall through the heating region. We additionally notice that the heavy nuclei in the cold unshocked core have either been dissociated by heating or have made the phase transition to bulk nuclear matter. Thus, the same amount of mass has been dissociated in the NR and GR cases and the corresponding earlier advantage of the NR case in the energy balance behind the shock does not persist. We see in graph (g) that the neutrino radiation becomes increasingly forward peaked outside of the neutrinospheres. There is an extended region between $`50`$ and $`200`$ km where the neutrinos can deposit energy by absorption. In the interior part of this region, the net cooling region, this process is more than outweighed by neutrino emission of the infalling material via electron and positron capture, as the material is compressed, heated, and layered onto the proto-neutron star. In the outer net heating region, however, electron flavor neutrino absorption just exceeds neutrino emission, leading to comparatively small net heating. The heating can indirectly be observed in graph (c) where the entropy per baryon starts to increase behind the shock. We read from graph (g) that the neutrinospheres in the more compact and hotter GR core are at smaller radii. This leads to emission of neutrinos at higher rms energies (graph (h)). The discontinuity in the rms energies at the shock are due to the Doppler shift across it. All radiation quantities are those measured by comoving observers. This choice also leads to the luminosity step in graph (e): the infalling observer heading towards the center observes a higher luminosity than the observer almost at rest behind the shock. Although the GR luminosities in the heating region are up to $`20\%`$ higher than in the NR case, the heating is not more effective in the GR case. The shocked material settles onto the proto-neutron star with a higher velocity in the deeper gravitational potential and therewith stays for a shorter time in the net heating region.

Finally, we show the situation at $`500`$ ms after bounce in Fig. (5). The shock in graph (a) has receded to $`60`$ km in the GR case and to $`90`$ km in the NR case. The accretion shock, sitting deep in the gravitational potential, experiences a high ram pressure from the infalling material. All nuclei are completely dissociated when they pass the shock (graph (f)). The infall velocities behind the shock are as large as $`5000`$ km/s in the NR simulation and almost double that in the GR simulation. This produces high accretion luminosities in the GR case, as shown in graph (e). Higher rms energies in the GR case are evident in graph (h), as expected from the deeper position of the neutrinospheres. The net heating is not very efficient on the rapidly infalling material, and the entropy per baryon increases only very slowly at $`1.7k_B`$/($`100`$ ms). This rise is due to the infall of higher entropy material and the retraction of the shock radius to deeper points in the gravitational potential such that a larger kinetic energy from the infalling material is dissipated. The entropy per baryon reaches $`20k_B`$ in graph (c). Overall, most features in the GR simulation are enhanced relative to the NR simulation, and the GR case exhibits a more compact structure in the proto-neutron star and its surroundings. The central density reached is $`5.15\times 10^{14}`$ g/cm3 in graph (b).

The trajectories of mass shells containing $`0.01`$ M are shown for the relativistic simulation in Fig. (6). Note that the high velocities (steep gradients of mass traces) in the heating region can hardly be distinguished from the ten times larger infall velocities outside of the shock. As discussed above, mass elements are almost stopped at the shock front and do not directly fall onto the proto-neutron star surface, as might be inferred from Fig. (6). The shock position is drawn for both the GR and the NR simulation. The GR shock barely reaches $`200`$ km, the NR shock reaches $`230`$ km.

Graphs (a)-(d) of Fig. (7) show the neutrino heating rates (dashed lines with positive values), neutrino cooling rates (dashed lines with negative values), and their superpositions (solid lines). These energy transfer rates are shown separately for the electron neutrinos (a,b) and electron antineutrinos (c,d). We see that the net energy transfer rate is the difference between much larger cooling and heating rates, except very close to the shock, where heating clearly dominates the cooling. The location where the net heating changes sign defines the gain radius that separates the net cooling region (neutrinosphere to gain radius) from the net heating region (gain radius to shock radius). The peak net heating rate occurs rather close to the gain radius. Note the sharp cutoff at the shock in the heating by electron neutrino absorption with respect to the smoother decrease in electron antineutrino absorption. This is due to the steep drop in the free neutron mass fraction. Free neutrons are less abundant in the presence of neutron rich heavy nuclei outside the shock than the targets for electron antineutrino absorption, the free protons.

At $`100`$ ms after bounce we observe a net heating rate $`10^{20}`$ erg/g/s outside of the gain radius at $`100`$ km. Note that the net heating by electron antineutrino absorption is higher than the net heating by electron neutrino absorption. This is due to the fact that the increasing electron chemical potential for the infalling material, at comparable heating rates, favors cooling by electron capture compared to cooling by positron capture. This example illustrates that the net energy transfer or the net change in the electron fraction by weak interactions cannot be estimated simply based on luminosities and absorption cross sections without a detailed consideration of the compensating inverse reaction rates. We found increasing electron fractions that exceed a value of $`0.5`$ only in our earlier exploding models where the decrease in the chemical potential in the hot expanding layers allowed for comparable electron and positron capture at a slightly favored electron neutrino absorption rate (e.g. in Ref. ).

Graph (e) shows the rest mass flux $`4\pi r^2u\rho /\mathrm{\Gamma }`$ through spherical surfaces measured at constant radii $`r`$ at $`100`$ ms after bounce. The time slices at $`70`$ ms, (g), and at $`130`$ ms, (h), are added to sketch the time dependence of the mass flux when the shock starts to recede. The largest mass flux occurs at small radii where the high density proto-neutron star slowly adjusts its size to the newly accreted material. Most interesting for shock revival is the region between the neutrinospheres (shown as vertical bars) and the shock radius, which enclose the cooling and heating region visible in graph (a) and (b). The discontinuity in graph (g) in the mass flux at the shock front indicates that the shock is still moving outwards. The material at the gain radius and at the shock radius show a comparable mass flux at $`70`$ ms after bounce. However, at $`100`$ ms after bounce in graph (e), and even more so at $`130`$ ms after bounce in graph (h), a slope in the mass flux develops. The mass flux through the neutrinospheres is larger than the mass flux through the gain radius, and both are larger than the mass flux through the shock. The latter is determined by the infalling material, i.e., the density profile of the progenitor and the gravitational potential. The material in the heating region is drained from below (Janka ), and the conditions for heating become inefficient because of the shortened time the infalling matter spends in the heating region. In the general relativistic simulation we observe this process earlier than in the nonrelativistic simulation. The mass flux and the energy transfer rates are about a factor of two larger.

At $`500`$ ms after bounce, the neutrinospheres and the shock have both retracted to smaller radii. The mass flux in graph (f) is almost constant down to the neutrinospheres. It is mainly set by the close to free infall of the layers around the iron core of the progenitor. Although the infall velocities have reached very large values, the actual mass flux has decreased because of the decreasing density in the outer layers. The gain radius has adjusted to a much smaller location $`35`$ km (NR $`60`$ km). In this compact configuration, the energy transfer rates have increased by more than an order of magnitude, as shown in graphs (b) and (d). At this late-time, we observe a peak in the heating rate across the shock front. Part of it is due to the increased heating rate caused by the Doppler shift in the rms energies and luminosities that applies to the material in the shock front when it is not yet completely decelerated but already dissociated. Another part is caused by numerical artifacts in the shock front: e.g., the unphysical width of the feature is set by the artificial viscosity. However, if we convolute the heating rate per mass with the density decline in the shock, we find a smooth decrease in the heating rate per volume across the shock and exclude relevant perturbations in the results owing to this feature.

We close the description of our simulations with two figures showing the temporal evolution of the GR radiation quantities. In Fig. (8) we show the neutrino luminosities at $`500`$ km radius as a function of time. The electron luminosity is slowly rising during collapse and decreases as the core reaches maximal density. It remains suppressed for the $`4`$ ms the shock needs to propagate to the electron neutrinosphere. The most prominent feature is the electron neutrino burst, reaching $`3.5\times 10^{53}`$ erg/s at shock breakout, and declining afterwards. Graph (b) shows the luminosity spectra at $`100`$ ms after bounce, just before the shock starts to recede. Figure (9) shows the long-term evolution with better resolution, sampled at a radius of $`500`$ km. The maximum luminosities in the GR simulation are $`3.9\times 10^{52}`$ erg/s, $`3.6\times 10^{52}`$ erg/s, and $`2.3\times 10^{52}`$ erg/s for the electron neutrino, electron antineutrino, and $`\mu `$ and $`\tau `$ neutrinos respectively. This is roughly $`10\%`$ more than the maximum luminosities reached (at different times) in the NR simulation. The GR results also show higher rms energies, but no qualitatively different features.

## IV Conclusion and Outlook

We have completed the construction of a general relativistic radiation hydrodynamics code, AGILE-BOLTZTRAN, and simulated the spherically symmetric, general relativistic core collapse, bounce, and postbounce evolution of a $`13`$ M progenitor. We investigate the confluence of (i) matter and radiation in a deeper effective gravitational potential, (ii) a GR core hydrodynamic structure that acts as a more intense neutrino source, and (iii) an increased heating efficiency obtained from accurate three-flavor Boltzmann neutrino transport. However, we find that the combination of these ingredients does not result in a supernova explosion. Our model shares this outcome with recent simulations that investigated a subset of these issues (Rampp and Janka , Mezzacappa et al. , Bruenn et al. ).

(i) During core collapse, the homologous core is smaller in the general relativistic gravitational potential than in Newtonian gravity. This leads to a smaller enclosed mass at shock formation and higher initial dissociation losses when the shock moves to larger radii. The neutrinos in curved spacetime propagate on trajectories with nearly constant $`\epsilon =(\mathrm{\Gamma }+u\mu )E`$ and $`b=r\sqrt{1\mu ^2}/(\mathrm{\Gamma }+u\mu )`$, where $`\mu `$ and $`E`$ are the neutrino propagation-angle cosine and energy measured by comoving observers. This latter effect is of little importance for our low-mass progenitor. The GR corrections for redshift and curvature between the neutrinospheres and the heating region do not exceed $`3\%`$ until shock recession, and increase only later, when the neutrinospheres have receded to smaller radii. (ii) Nonlinear effects in GR enhance the self-gravitation in the high-density domain of the proto-neutron star. The latter becomes more compact and exhibits higher internal energies. This in turn leads to higher core and accretion luminosities with harder spectra in all neutrino flavors (Bruenn et al.). (iii) The solution of the Boltzmann equation for the neutrino transport reproduces accurately the angular distribution of the neutrino radiation field behind the shock. An increased angular spread in the semi-transparent regime keeps the outstreaming neutrinos longer in the heating region and, therefore, increases the absorption efficiency with respect to the multigroup flux-limited diffusion approximation (Messer et al., Yamada et al. ).

The main difference between the GR and NR simulations stems from the difference in the size of the proto-neutron star, which is caused by the nonlinear GR effects at very high densities in the center of the star. At radii of order $`100`$ km and larger, the gravitational potential becomes comparable in the GR and NR cases. However, large differences arise if the steep gravitational well is probed at different positions in the GR and NR simulations, as happens with accretion down to the neutrinospheres, which represent the outer “boundary” of the proto-neutron star. The neutrinospheres in the more compact GR case receive material that has traversed a larger potential difference and, setteling deeper in the well, produces more energetic accretion luminosities. After the shock has stalled, the surrounding layers (cooling-heating region, shock radius) adjust to a smaller radius and settle to a stationary equilibrium in the spirit of Burrows and Goshy . In the GR simulation, this occurs at a smaller radius, deeper in the gravitational well, with higher infall velocities, higher accretion luminosities, and higher heating rates, which sustain equilibrium between the pressure behind the shock and the increased ram pressure ahead of it. The tight feedback between infall velocity, accretion luminosity and heating leads to higher luminosities in more compact configurations where the shock stalls at a smaller radius. Moreover, this self regulation has observable consequences: The gravitational well is probed when the proto-neutron star radius is set, tuning the directly observable neutrino luminosities and spectra.

From the comparison between the NR and GR shock trajectories in Fig. (6), we might be tempted to conclude that the chances for a vigorous explosion in the more realistic GR case are more pessimistic than in the NR case. We however feel that such a conclusion would be premature if only based on a postbounce evolution that fails to reproduce an explosion. An actual boost in the explosion energy due to general relativistic effects has been seen by Liebendörfer  in a simulation that led to an explosion because of incorrect nucleon isoenergetic scattering opacities. As in the simulation with the correct opacities, a more compact proto-neutron star was created and higher luminosities with harder spectra developed in the GR case. However, in the exploding case, a turning point occurred where the infall velocities behind the shock almost vanished. This happened $`100`$ ms after bounce, before the heating region had the time to step deeper into the gravitational well. The consequence was a more efficient heating with the GR enhanced luminosities and spectra, which lead to a more energetic supernova than in the corresponding NR simulation. Our conclusion for low-mass progenitors with small redshifts is: GR effects do not cause explosions in our current simulations. An essential physical ingredient is missing. However, GR effects may eventually drive an explosion more efficiently when it has been launched.

There are many facets of a supernova that cannot consistently be included in a spherically symmetric model. There is no doubt that convection behind the shock will occur, and significant rotation and strong magnetic fields might be present. Observations of neutron star kicks, mixing of species, inhomogeneous ejecta, and polarization of spectra suggest the presence of asymmetries in supernova explosions (e.g. Tueller et al. , Strom et al. , Galama et al. , Leonard et al. , and references therein). Motivated by such observations, various multi-dimensional explosion mechanisms have been explored (Herant et al., Miller et al., Herant et al., Burrows et al., Janka and Müller , Mezzacappa et al., Fryer , Fryer and Heger ) and jet-based explosion scenarios received new momentum (Höflich et al. , Khokhlov et al. , MacFadyen and Woosley , Wheeler et al. ). It remains a challenge for the next decade to include all relevant physics in supernova simulations in an attempt to reproduce supernova observations in detail. For the time being, one has to single out a subset of the known physics and to investigate the role each part plays in the restricted simulation. It is natural to start with ingredients that have long been believed to be essential for the explosion and to add modifiers in a systematic way until the observables can be reproduced (Mezzacappa et al. ). Spherically symmetric supernova modeling has a long tradition and is nearing a definitive point in the sense that high-resolution hydrodynamics, general relativity, complete Boltzmann neutrino transport, and reasonable nuclear and weak interaction physics are being combined to dispel remaining uncertainties. Models that are restricted to spherical symmetry allow the detailed simulation and study of radiation hydrodynamic flows in general relativity, with all feedbacks included. Moreover, they will provide a gauge for unavoidable approximations in future, more inclusive simulations. The accurate description of the high-density region in the proto-neutron star and realistic multigroup transport are prerequisits for the evaluation of new developments in the input physics. With the former well underway, research in spherically symmetric supernova models can focus on the uncertainties in the nuclear physics and neutrino interactions. A major revision in the equation of state (Shen et al. , see also Swesty et al. ) or in the neutrino-matter interactions (Burrows and Sawyer , Reddy et al. , Thompson et al. , Langanke and Martinez-Pinedo ), or the inclusion of neutrino oscillations (Fuller et al. , Abazajian et al. ), might be the last opportunity to understand the most basic supernova observable in restricted spherical models: explosion.

## Acknowledgments

We enjoyed fruitful discussions with Thomas Janka and Markus Rampp. M.L. is supported by the National Science Foundation under contract AST-9877130 and, formerly, was supported by the Swiss National Science Foundation under contract 20-53798.98. A.M. is supported at the Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the U.S. Department of Energy under contract DE-AC05-00OR22725. F.-K.T. is supported in part by the Swiss National Science Foundation under contract 20-61822.00 and as a Visiting Distinguished Scientist at the Oak Ridge National Laboratory. O.E.B.M. is supported by funds from the Joint Institute for Heavy Ion Research and a Dept. of Energy PECASE award. W.R.H. is supported by NASA under contract NAG5-8405 and by funds from the Joint Institute for Heavy Ion Research and a Dept. of Energy PECASE award. S.W.B. is supported by the NSF under contract 96-18423 and by NASA under contract NAG5-3903. Our simulations were carried out on the ORNL Physics Division Cray J90 and the National Energy Research Supercomputer Center Cray J90.


text:
# 𝑂(𝑚_𝑑-𝑚_𝑢) Effects in CP-even and CP-odd 𝐾→𝜋𝜋 Decays

## Abstract

Strong isospin-breaking (IB) effects in CP-even and CP-odd $`K\pi \pi `$ decays are computed to next-to-leading order (NLO) in the chiral expansion. The impact of these corrections on the magnitude of the $`\mathrm{\Delta }I=1/2`$ Rule and on the size of the IB correction, $`\mathrm{\Omega }_{IB}`$, to the gluonic penguin contribution to $`ϵ^{}/ϵ`$ are discussed.

thanks: E-mail: maltman@fewbody.phys.yorku.cathanks: E-mail: wolfe@niobe.iucf.indiana.edu

In the presence of IB, the standard isospin decomposition of the $`K^+\pi ^+\pi ^0`$, $`K^0\pi ^+\pi ^{},\pi ^0\pi ^0`$ decay amplitudes, $`A_{+0}`$, $`A_+`$ and $`A_{00}`$, becomes

$`A_{00}`$ $`=`$ $`\sqrt{1/3}A_0\mathrm{e}^{i\mathrm{\Phi }_0}[\sqrt{2/3}]A_2\mathrm{e}^{i\mathrm{\Phi }_2},`$
$`A_+`$ $`=`$ $`\sqrt{1/3}A_0\mathrm{e}^{i\mathrm{\Phi }_0}+[1/\sqrt{6}]A_2\mathrm{e}^{i\mathrm{\Phi }_2},`$
$`A_{+0}`$ $`=`$ $`[\sqrt{3}/2]A_2^{}\mathrm{e}^{i\mathrm{\Phi }_2^{}}.`$ (1)
In the absence of the $`I=2`$ component of electromagnetism (EM), the $`\mathrm{\Phi }_I`$ are the $`\pi \pi `$ phases. In general, $`|A_2^{}||A_2|`$ due to EM- and strong-IB-induced $`\mathrm{\Delta }I=5/2`$ contributions. $`A_0`$, $`A_2`$ can be chosen real in the absence of CP violation.

Since $`|A_0|20|A_2|`$, IB “leakage” of the large octet amplitude into the $`\mathrm{\Delta }I=3/2`$ amplitude can be numerically significant. EM leakage contributions have been computed to NLO in Ref. ; we compute the NLO strong octet IB contributions. These enter Standard Model predictions of $`ϵ^{}/ϵ`$ where the strong cancellation between gluonic penguin ($`O_6`$) and electroweak penguin ($`O_8`$) contributions is sensitive to the degree of strong-IB-induced suppression of the $`O_6`$ contribution .

At leading chiral order (LO), the computation of the octet leakage contribution is unambiguous; the magnitude of the LO weak $`27`$-plet low-energy constant (LEC) is decreased by $`\mathrm{\Omega }_{IB}=13\%`$. The corresponding $`O_6`$ suppression in $`ϵ^{}/ϵ`$ is $`1\mathrm{\Omega }_{IB}`$. Recent analyses of $`ϵ^{}/ϵ`$ employ $`\mathrm{\Omega }_{IB}=0.25\pm 0.08`$, the difference from the LO value reflecting estimates of the effect of $`\eta ^{}`$ mixing. This effect is NLO in the chiral expansion, but does not exhaust NLO contributions. A full NLO calculation can be performed using Chiral Perturbation Theory (ChPT). The importance of such a complete NLO determination can be seen from the recent discussion of NLO $`\pi \eta `$ mixing effects : the $`\eta ^{}`$ contribution (associated with the strong LEC $`L_7^r`$) turns out to be almost completely cancelled by a contribution proportional to $`L_8^r`$ . To compute the NLO IB leakage contributions one evaluates the tree and one-loop graphs of Ref. . NLO tree contributions are either proportional to the product of the LO weak octet LEC $`c^\pm `$ and a single NLO strong LEC or proportional to one of the NLO weak LEC’s. All loop graphs involve one vertex from the LO octet effective weak Lagrangian, $`c^\pm Tr\left[\lambda ^\pm _\mu U^{}^\mu U\right]`$, where the superscripts $`\pm `$ label the CP even and odd cases, respectively, $`\lambda ^+=\lambda _6`$, $`\lambda ^{}=\lambda _7`$, and $`U=exp\left(i\lambda \pi \right)`$, is the usual matrix variable. The (scale-dependent) ratio of the sum of the loop contributions to the LO octet contribution for a given amplitude is thus completely fixed; the main uncertainty lies in a lack of knowledge of the NLO weak LEC’s, for which we are forced to use models (see Refs.  for further discussion).

The contributions to $`A_2`$ and $`A_2^{}`$ associated with EM  and octet IB  leakage are given in Table 1. The errors reflect uncertainties in the estimates of the NLO LEC’s. Denoting the ratio of LO $`27`$-plet to octet weak LEC’s obtained neglecting, or including, IB by $`r_{IC}`$, or $`r_{IB}`$, respectively, we find $`R_{IB}r_{IB}/r_{IC}=0.963\pm 0.029\pm 0.010\pm 0.034`$. The errors reflect uncertainties in the weak NLO LEC combinations, the input value of $`B_0(m_dm_u)`$, and the EM contributions, respectively. The deviation from $`1`$ is significantly smaller than at LO (where $`R_{IB}=0.870`$). The $`\mathrm{\Delta }I=5/2`$ contribution (dominantly EM in character ), leads to $`|A_2|/|A_2^{}|=1.094\pm 0.0391`$, and significantly exacerbates the phase discrepancy problem for the neutral $`K`$ decays .

For the CP-odd case, $`\mathrm{\Omega }_{IB}=\omega Im\delta A_2/ImA_0`$ ($`\omega =ReA_0/ReA_222.2`$; $`\delta A_2`$ is the octet leakage contribution). At LO, $`\mathrm{\Omega }_{IB}=0.13[\mathrm{\Omega }_{IB}]_{LO}`$. At NLO $`\mathrm{\Omega }_{IB}=[\mathrm{\Omega }_{IB}]_{LO}\left[1+\frac{\mathrm{Im}\delta A_2^{(NLO;ND)}}{\mathrm{Im}\delta A_2^{(LO)}}\frac{\mathrm{Im}A_0^{(NLO;ND)}}{\mathrm{Im}A_0^{(LO)}}\right][\mathrm{\Omega }_{IB}]_{LO}\left[1+R_2R_0\right]`$. The superscript $`(NLO;ND)`$ indicates the sum of non-dispersion NLO contributions (involving NLO weak and strong LEC’s and the non-dispersive parts of loop graphs). Neither the NLO $`I=0`$ IC nor NLO $`I=2`$ IB leakage CP-odd LEC combinations are known. The NLO dispersive contributions create phases consistent with Watson’s theorem. Although the positive $`I=0`$ phases correspond to attractive FSI, NLO weak LEC corrections may, nonetheless, make $`ImA_0`$ smaller at NLO than the LO (see comments on Ref.  in Ref.  for a related discussion). If, however, NLO effects do enhance $`ImA_0`$ (decreasing the level of $`O_6`$-$`O_8`$ cancellation and increasing $`ϵ^{}/ϵ`$) $`\mathrm{\Omega }_{IB}`$ will be simultaneously suppressed, further increasing $`ϵ^{}/ϵ`$. The known NLO contributions (loops and strong LEC terms) give contributions $`0.24(0.31)`$ to $`R_2`$ and $`0.02(+0.42)`$ to $`R_0`$, at scale $`\mu =m_\eta (m_\rho )`$. Using the weak deformation model to estimate the weak NLO LEC’s, $`1+R_2R_0=0.27`$ while, for the chiral quark model, it lies between $`0.62`$ and $`1.42`$. Averaging these results, and taking their spread as a minimal indication of theoretical error, we thus obtain, at NLO, $`\mathrm{\Omega }_{IB}=0.11\pm 0.08`$, a much smaller value than conventionally employed, though with comparable errors. It is also significantly smaller than the partial NLO estimate of Ref. , $`0.16\pm 0.03`$, based on the strong LEC contributions only. The central value above, combined with conventional central values for the $`B`$-factors, leads to a $`50\%`$ increase in the predicted value for $`ϵ^{}/ϵ`$. To be conservative, we would propose using this lower central value with an even larger error estimate, in all future calculations of Standard Model values for $`ϵ^{}/ϵ`$.