text:
# hep-ph/0008276CERN-NUFACT note 43CERN-TH/2000-213, FISIST/8-2000/CFIF Charged-Lepton-Flavour Violation in Kaon Decays in Supersymmetric Theories

## 1 Introduction

The recent Super-Kamiokande data  have triggered an upsurge of interest in extensions of the Standard Model with massive neutrinos and/or violation of the charged-lepton numbers in processes such as $`\mu e\gamma `$, $`\mu 3e`$, $`\tau \mu \gamma `$ and $`\mu e`$ conversion on heavy nuclei . The present experimental upper bounds on the most interesting of these processes are:

$`BR(\mu e\gamma )<\mathrm{\hspace{0.17em}1.2}\times 10^{11}`$ $`:`$
$`BR(\mu ^+e^+e^+e^{})<\mathrm{\hspace{0.17em}1.0}\times 10^{12}`$ $`:`$
$`R(\mu ^{}Tie^{}Ti)<\mathrm{\hspace{0.17em}6.1}\times 10^{13}`$ $`:`$
$`BR(\tau \mu \gamma )<1.1\times \mathrm{\hspace{0.17em}10}^6`$ $`:`$  (1)
Our main interest in this paper is in rare kaon decays, for which the current experimental bounds are:

$`BR(K_L^0\mu ^\pm e^{})<4.7\times \mathrm{\hspace{0.17em}10}^{12}`$ $`:`$
$`BR(K_L^0e^+e^{})=(8.7_{4.1}^{+5.7})\times \mathrm{\hspace{0.17em}10}^{12}`$ $`:`$
$`BR(K_L^0\mu ^+\mu ^{})=(7.18\pm 0.17)\times \mathrm{\hspace{0.17em}10}^9`$ $`:`$
$`BR(K^+\pi ^+\mu ^+e^{})<2.8\times \mathrm{\hspace{0.17em}10}^{11}`$ $`:`$
$`BR(K^+\pi ^+\mu ^{}e^+)<5.2\times \mathrm{\hspace{0.17em}10}^{10}`$ $`:`$
$`BR(K_L^0\pi ^0\mu ^\pm e^{})<3.1\times \mathrm{\hspace{0.17em}10}^9`$ $`:`$
$`BR(K^+\pi ^+e^+e^{})<(2.94\pm 0.05\pm 0.14)\times \mathrm{\hspace{0.17em}10}^7`$ $`:`$
$`BR(K_L^0\pi ^0e^+e^{})<4.3\times \mathrm{\hspace{0.17em}10}^9`$ $`:`$
$`BR(K^+\pi ^+\mu ^+\mu ^{})<(7.6\pm 2.1)\times \mathrm{\hspace{0.17em}10}^8`$ $`:`$
$`BR(K_L^0\pi ^0\mu ^+\mu ^{})<5.1\times \mathrm{\hspace{0.17em}10}^9`$ $`:`$
$`BR(K^+\pi ^{}\mu ^+e^+)<5.0\times \mathrm{\hspace{0.17em}10}^{10}`$ $`:`$
$`BR(K^+\pi ^{}e^+e^+)<6.4\times \mathrm{\hspace{0.17em}10}^{10}`$ $`:`$
$`BR(K^+\pi ^{}\mu ^+\mu ^+)<3.0\times \mathrm{\hspace{0.17em}10}^9`$ $`:`$  (2)
Projects are being proposed that could be used to improve these limits significantly, e.g., some of the powerful proton sources being proposed for neutrino factories  could provide intense secondary kaon beams.

Any observable rate for one of these processes would constitute unambiguous evidence for new physics. The rates for such processes remain extremely suppressed if we simply extend the Standard Model to include right-handed neutrinos, but larger rates are possible in more ambitious extensions of the Standard Model. Supersymmetry is one example of new physics that could amplify rates for some of the rare processes (1, 2), either in the minimal supersymmetric extension of the Standard Model (MSSM) or in its modification to include violation of $`R`$ parity.

In previous works, rare charged-lepton decays and anomalous $`\mu e`$ conversions on heavy nuclei have received considerable attention , leading to a good understanding of the correlations between the predicted rates and possible violations of leptonic universality and/or exotic Yukawa couplings with $`\mathrm{\Delta }L0`$. It is natural to include rare kaon decays in this analysis, because strangeness-changing decays occur in the Standard Model. As mentioned above, there is now considerable discussion of intense proton sources to be used as muon sources for neutrino factories . If the protons have high energy above $`15`$ GeV, as in some neutrino factory designs, they would also be copious sources of kaons. These might provide a new opportunity to study rare $`K`$ decays with high statistics, and we are interested to know whether these might cast additional light on neutrino masses and mixing. Specifically, as we show in this paper, rare $`K`$ decays that violate charged-lepton number allow one, in the context of the MSSM, to correlate the (s)quark and (s)lepton mixing, whilst in $`R`$-violating supersymmetry one can probe interesting products of Yukawa couplings.

This paper is organized as follows. In Section 2, we analyze the decays $`K^0\mu ^\pm e^{}`$ in the MSSM, finding rates that are relatively small, although for certain model parameters they may be within the reach of an imaginable kaon beam at a neutrino factory. In Section 3, we analyze the same decays in $`R`$-violating supersymmetric models, finding that the rates could in principle be much larger, close to the present experimental limits. Indeed, these existing upper limits on these decays already impose interesting upper limits on products of $`LL\overline{E}`$ and $`LQ\overline{D}`$ couplings. We also point out that $`\stackrel{~}{b}_L`$$`\stackrel{~}{b}_R`$ mixing in the presence of $`R`$ violation could lead to the like-sign-lepton decays $`K^\pm \pi ^{}\mathrm{}^\pm \mathrm{}^\pm `$, although the existing experimental limits on these processes do not impose interesting upper limits on couplings. Taken together, our results indicate that rare $`K`$ decays could provide interesting signatures for supersymmetric models, in particular those with $`R`$ violation.

## 2 Rare Kaon Decays in the MSSM

In this Section, we evaluate the rates for rare kaon decays in the MSSM with massive neutrinos, using the see-saw mechanism , which we consider to be the most natural way to obtain neutrino masses in the sub-eV range. In particular, we assume Dirac neutrino masses $`m_\nu ^D`$ of the same order as the charged-lepton and quark masses, and heavy Majorana masses $`M_{\nu _R}`$, leading to a light effective neutrino mass matrix of the form:

$$m_{eff}=m_\nu ^D(M_{\nu _R})^1m_\nu ^{D^T}$$
(3)
Neutrino-flavour mixing  may then occur through either the Dirac matrix $`m_\nu ^D`$ and/or the Majorana mass matrix $`M_{\nu _R}`$, which may also feed flavour violation through to the charged leptons. In non-supersymmetric models with massive neutrinos, the amplitudes for charged-lepton-flavour violation are proportional to inverse powers of the right-handed neutrino mass scale $`M_{\nu _R}`$, and thus the rates for rare decays are extremely suppressed . On the other hand, in supersymmetric models one must also take the dynamics of sneutrinos $`\stackrel{~}{\nu }`$ into account, and these processes may only be only suppressed by inverse powers of the supersymmetry-breaking scale, which characterizes $`m_{\stackrel{~}{\nu }}`$ and is at most $`𝒪(1)`$ TeV .

The magnitudes of the predicted rates depend on the details of the masses and mixings of sparticles including the sneutrinos. If their soft supersymmetry-breaking masses are non-universal at $`M_{GUT}`$, large rates are in general predicted . However, even if the soft supersymmetry-breaking masses are universal at $`M_{GUT}`$, renormalization effects in the MSSM with right-handed neutrinos spoil this diagonal form  at lower scales. This is because the Dirac neutrino and charged-lepton Yukawa couplings and masses $`m_{\mathrm{}}`$ cannot, in general, be diagonalized simultaneously. Since both these sets of lepton Yukawa couplings appear in the renormalization-group equations, the lepton and slepton mass matrices also may not be diagonalized simultaneously at low energies.

To illustrate this point, let us consider the simplest example of a model based on Abelian flavour symmetries and symmetric mass matrices . Requiring large (2-3) mixing in this model  severely constrains the possible flavour charges and thus the forms of the charged-lepton and the neutrino mass matrices. A representative example is given by Ansatz A of :

$`m_{\mathrm{}}\left(\begin{array}{ccc}\overline{ϵ}^7& \overline{ϵ}^3& \overline{ϵ}^{7/2}\\ \overline{ϵ}^3& \overline{ϵ}& \overline{ϵ}^{1/2}\\ \overline{ϵ}^{7/2}& \overline{ϵ}^{1/2}& 1\end{array}\right),m_\nu ^D\left(\begin{array}{ccc}ϵ^7& ϵ^3& ϵ^{7/2}\\ ϵ^3& ϵ& ϵ^{1/2}\\ ϵ^{7/2}& ϵ^{1/2}& 1\end{array}\right)`$ (10)
where $`\overline{ϵ}=\sqrt{ϵ}=0.2`$. We already see that the two matrices cannot be simultaneously diagonalized and indeed,

$`V_{\mathrm{}}=\left(\begin{array}{ccc}1& \overline{ϵ}^2& \overline{ϵ}^{7/2}\\ \overline{ϵ}^2& 1& \overline{ϵ}^{1/2}\\ \overline{ϵ}^{7/2}& \overline{ϵ}^{1/2}& 1\end{array}\right),V_{\nu _D}=\left(\begin{array}{ccc}1& \overline{ϵ}^4& \overline{ϵ}^7\\ \overline{ϵ}^4& 1& \overline{ϵ}\\ \overline{ϵ}^7& \overline{ϵ}& 1\end{array}\right)`$ (17)
Within this general framework, there is ambiguity in the specification of numerical coefficients in the matrix elements, which are expected to be of order unity. We return later to this point.

For a generic texture where the charged lepton and neutrino matrices are not simultaneously diagonal, the slepton mass matrix acquires non-diagonal contributions from renormalization at scales below $`M_{GUT}`$. In the basis where $`m_{\mathrm{}}`$ is diagonal, these corrections take the form :

$`\delta m_\stackrel{~}{\mathrm{}}^2{\displaystyle \frac{1}{16\pi ^2}}(3+a^2)\mathrm{ln}{\displaystyle \frac{M_{GUT}}{M_N}}\lambda _D^{}\lambda _Dm_{3/2}^2,`$ (18)
where $`\lambda _D`$ is the Dirac neutrino Yukawa coupling, $`M_N`$ is the scale where the effective neutrino-mass operator is formed, $`a`$ is related to the trilinear mass parameter $`A_{\mathrm{}}am_{3/2}`$, and $`m_{3/2}`$ is the presumed common value $`m_0`$ of the scalar masses at the GUT scale.

In the case of non-universal soft masses, these corrections are generically negligible. However, the rates for $`\mathrm{\Delta }L0`$ processes are generally too large in such non-universal models . On the other hand, models with scalar-mass universality at $`M_{GUT}`$, such as no-scale  and gauge-mediated models , may yield acceptable rates for $`\mathrm{\Delta }L0`$ processes. In such models, the contributions (18) related to neutrino masses dominate and lead to non-negligible rates for the lepton-flavour-violating processes which are determined by the off-diagonal terms in the Yukawa matrix $`\lambda _N`$. The various different solutions of the solar neutrino deficit , with a small/large mixing angle and with eV or $``$ 0.03 eV neutrinos, predict in general different rates for charged-lepton-flavour violation: the larger the $`\mu e`$ mixing, and the larger the neutrino mass scales that are required, the larger the rates. Thus, degenerate neutrinos with bimaximal mixing may be expected to yield significantly larger effects than, for instance, hierarchical neutrinos with a small vacuum mixing angle. Note, in particular, that the just-so solutions to the solar neutrino problem with $`\delta m^210^{10}`$ eV2 predict, in the case of hierarchical neutrino masses, rates that are small, even if the first/second-generation neutrino mixing is large.

In this MSSM framework, rare kaon decays are generated by box diagrams involving chargino and neutralino exchanges. For instance, for $`K^0\mu ^\pm e^{}`$ we have the diagrams 111The contribution from the first diagram when the neutralino is a purely photino state has been discussed in . of Fig. 1 222Kaon decays in left-right symmetric models have been analysed in . In principle, using the mixing in the quark sector \- in particular between $`s`$ and $`d`$ quarks - we can also generate $`K\mu e`$ by penguin diagrams. In the cases that all the quark mixing is in the down sector, or the right-handed mixing in the down sector is much bigger than the one in the left (which is the one bounded by $`V_{CKM}`$) the rates might be of relevance. However, we do not address here this model-dependent issue.. It is clear that $`K^{\pm ,0}\pi ^{\pm ,0}\mu ^\pm e^{}`$ decays can be generated in a similar way, but in this case the experimental bounds are worse by almost two orders of magnitude, and we do not discuss them in detail.

Our procedure for evaluating these contributions is as follows:

* We first find the maximal squark mixing that is allowed by the neutral-kaon mass difference $`\mathrm{\Delta }m_K`$ .
* We next find the maximal slepton mixing allowed by $`\mu e\gamma `$ and $`\mu `$$`e`$ conversion in nuclei in a model-independent way.
* Having fixed these values, finally we calculate the rates for rare kaon decays.

As we noted previously, the $`\mu e`$ mixing is constrained by the form of the neutrino textures and thus by the recent neutrino data: in general, degenerate neutrinos with large angle MSW oscillations require smaller soft masses to be consistent with the observed rates. However, as we found in , even for the small-angle MSW solutions of the solar neutrino deficit, we can obtain large rates for values of $`m_0`$ and $`m_{1/2}`$ well below 500 GeV. The latest Super-Kamiokande data  on solar neutrinos favour large mixing angles, which might suggest larger $`\mu e`$ flavour violation. Thus, considering models with small mixing angles is conservative.

For illustration, we focus here on one such model, namely variant $`A_1`$ of the texture (10, 17), which has the numerical values of the ambiguous $`𝒪(1)`$ coefficients fixed as described in . For definiteness, we choose its ‘inverted’ option with negative off-diagonal entries in the Dirac neutrino coupling matrix . Our results are presented in Figure 2, where we summarise our predictions for $`BR(K\mu e)`$, in association with the predictions for $`BR(\mu e\gamma )`$ and $`BR(\mu e)`$ for different values of $`\mathrm{tan}\beta `$ and $`m_{1/2}`$. We parameterize the supersymmetric masses in terms of the universal GUT-scale parameters $`m_0`$ and $`m_{1/2}`$, for sfermions and gauginos respectively, and use the renormalization-group equations of the MSSM to calculate the low-energy sparticle masses. Other relevant free parameters of the MSSM are the trilinear coupling $`A`$ (for which we start with the initial condition $`A_0=m_{1/2}`$), the sign of the Higgs mixing parameter $`\mu `$, and the value of $`\mathrm{tan}\beta `$. Models with different signs of $`\mu `$ give similar results: here we assume $`\mu <0`$.

As expected, the larger the value of $`\mathrm{tan}\beta `$ and the smaller the soft supersymmetric terms, the larger the branching ratios, apart from certain cancellations. In the case $`\mathrm{tan}\beta =10`$ and $`m_{1/2}=250`$ GeV, we see that, for the range $`m_0170`$ GeV where $`BR(\mu e\gamma )`$ and $`BR(\mu e)`$ conversion are consistent with the current experimental bounds (1), $`BR(K\mu e)`$ is at most $`210^{18}`$. However, for the same value of $`m_{1/2}`$, when $`\mathrm{tan}\beta =20`$ we find a significantly larger branching ratio at small values of $`m_0170`$ GeV. Moreover, for smaller $`m_{1/2}=150`$ GeV, we gain almost two orders of magnitude when we consider $`m_0`$ in the low-mass window between 100 and 150 GeV. We recall  that these lower values of $`m_{1/2},m_0`$ are consistent with accelerator constraints and generically yield cold dark matter densities in the range preferred by cosmology .

We do not discuss here other model textures for the mass matrices. Rather, our point here has been to demonstrate that, despite the limits from $`\mu e\gamma `$, $`\mu `$$`e`$ conversion and $`\mathrm{\Delta }m_K`$, the branching ratio of $`K\mu e`$ may be within the reach of the next generation of experiments, namely in the range $`10^{16}10^{18}`$, at least if $`\mathrm{tan}\beta `$ is large and the soft supersymmetry-breaking terms are small. The sensitivities (1) to $`\mu e\gamma `$ and $`\mu `$$`e`$ conversion could be improved significantly even before a neutrino factory comes into operation, and such a machine would offer enhanced prospects for probing them. It is therefore likely that the best prospects for discovering charged-lepton flavour violation may be offered by $`\mu e`$ conversion and $`\mu e\gamma `$. However, rare kaon decays provide a complementary probe which also gives information on the squark mixing, in the context of the MSSM.

## 3 Rare Kaon Decays in $`R`$-Violating Supersymmetry

We now discuss kaon decays violating charged-lepton flavour beyond the context of the MSSM. As is well known, the gauge symmetries of the MSSM allow additional dimension-four Yukawa couplings, of the form

$$\lambda L_iL_j\overline{E}_k,\lambda ^{}L_iQ_j\overline{D}_k,\lambda ^{\prime \prime }\overline{U}_i\overline{D}_j\overline{D}_k$$

where the $`L(Q)`$ are the left-handed lepton (quark) superfields, and the $`\overline{E}`$,($`\overline{D},\overline{U}`$) are the corresponding right-handed fields. If all these couplings were present simultaneously in the low-energy Lagrangian, they would generate unacceptably fast proton decay. Therefore, extra symmetries must be invoked to forbid all ($`R`$-parity ), or subsets (baryon and lepton parities ) of these couplings. In the latter case, very interesting baryon- and lepton-number-violating processes may occur .

Imposing electroweak $`SU(2)`$ and colour $`SU(3)`$ invariance, one finds that there are just 45 $`R`$-violating couplings in total. Besides proton decay, there are many experimental constraints on these couplings, both individually and in various combinations, from the non-observation of modifications to Standard Model processes and of possible exotic processes . In order to understand the possible hierarchies of $`R`$-violating couplings, models of flavour symmetries have been invoked . For instance, it was found in previous work  that theories with symmetric fermion mass textures lead to the expectation that $`R`$-violating couplings are small: $`10^3`$-$`10^4`$ for 100 GeV sfermion masses, whilst, in models with asymmetric fermion mass textures, dominance by a single coupling may be permitted, without however excluding several products of couplings from being non-negligible.

In this paper, therefore, we allow the general possibility that several $`R`$-violating operators may be large, and discuss the limits on their combinations that are obtainable from kaon decays. In this class of models, whilst $`\mu e\gamma `$ occurs at the one-loop level, $`\mu 3e`$, $`\mu e`$ conversion, $`K^0\mu ^\pm e^{}`$ and $`K^{\pm ,0}\pi ^{\pm ,0}\mu ^\pm e^{}`$ may occur at tree level via different combinations of couplings. For instance, in the case of $`LQ\overline{D}`$ couplings, $`\mu 3e`$ gives the limit

$`(L_1Q_i\overline{D}_j)(L_2Q_i\overline{D}_j)10^4\left({\displaystyle \frac{m_{\stackrel{~}{f}}}{100\mathrm{GeV}}}\right)^2`$ (19)
whilst $`\mu `$$`e`$ conversion in Titanium gives:

$`(L_2Q_1\overline{D}_k)(L_1Q_1\overline{D}_k)10^8\left({\displaystyle \frac{m_{\stackrel{~}{f}}}{100\mathrm{GeV}}}\right)^2`$ (20)
$`(L_2Q_j\overline{D}_1)(L_1Q_j\overline{D}_1)10^8\left({\displaystyle \frac{m_{\stackrel{~}{f}}}{100\mathrm{GeV}}}\right)^2`$ (21)
and the bound

$`(L_iQ_1\overline{D}_2)(L_iQ_2\overline{D}_1)10^9\left({\displaystyle \frac{m_{\stackrel{~}{f}}}{100\mathrm{GeV}}}\right)^2`$ (22)
is obtainable from $`\mathrm{\Delta }m_K`$ .

What is the connection of these results with neutrino masses? If $`R`$ parity is violated, neutrino masses are generated via one loop diagrams involving the vertices $`\nu _i\overline{d}_k\stackrel{~}{d}_j,\overline{\nu }_i^cd_j\stackrel{~}{d}_k^{}`$ for $`LQ\overline{D}`$ operators and the vertices $`\nu _i\overline{e}_k\stackrel{~}{\mathrm{}}_j,\overline{\nu }_i^c\mathrm{}_j\stackrel{~}{e}_k^{}`$ for $`LL\overline{E}`$ operators. Focusing on $`LQ\overline{D}`$ couplings and assuming that the left-right squark soft mixing terms are diagonal in the physical basis and proportional to the associated quark mass, the induced masses are given by

$$m_{\nu _{ii^{}}}\frac{n_c\lambda _{ijk}^{}\lambda _{i^{}kj}^{}}{16\pi ^2}m_{d_j}m_{d_k}\left(\frac{f(m_{d_j}^2/m_{\stackrel{~}{d}_k}^2)}{m_{\stackrel{~}{d}_k}}+\frac{f(m_{d_k}^2/m_{\stackrel{~}{d}_j}^2)}{m_{\stackrel{~}{d}_j}}\right),$$
(23)
where $`f(x)=(x\mathrm{ln}xx+1)/(x1)^2`$, $`m_{d_i}`$ is the down quark mass of the $`i`$th generation inside the loop, $`m_{\stackrel{~}{d}_i}`$ is an average of $`\stackrel{~}{d}_{Li}`$ and $`\stackrel{~}{d}_{Ri}`$ squark masses, and $`n_c=3`$ is the colour factor.

Requiring that these contributions be consistent with the neutrino data gives bounds on the associated $`R`$-violating products. We can see from the above expression that the heavier the fermions in the loop (including the associated fermion mass arising from the soft mixing term), the stricter the bounds . For example, demanding $`m_{e\mu }<2.5`$ eV for sparticle masses of 100 GeV leads to $`\lambda _{133}^{^{}}\lambda _{233}^{^{}}3.810^7`$, whilst for $`\lambda _{122}^{^{}}\lambda _{222}^{^{}}`$ the bound drops to $`2.310^4`$ . For higher sfermion masses, larger $`R`$-violating couplings are allowed. In Super-Kamiokande-friendly solutions with hierarchical neutrinos the bounds are stricter by two orders of magnitude.

Note, however, that neutrino masses do not strictly constrain $`K\mu e`$ (and in certain cases the rest of the flavour-violating-processes), since:

$``$ Neutrino masses may only constrain products of $`LL\overline{E}`$ or $`LQ\overline{D}`$ operators, not mixed $`LL\overline{E}`$-$`LQ\overline{D}`$ products.

$``$ Even for the diagrams with products of only $`LQ\overline{D}`$ operators, rare kaon decays involve quarks of the lightest and second-lightest generations. In this case the bounds from neutrino masses are significantly weaker, and the stricter limits come from the current measurements of the rare kaon decays themselves. The same is true for $`\mu e`$ conversion and even for $`\mu e\gamma `$ via $`LL\overline{E}`$ couplings. For fermions of the first two generations, the bounds from the lepton-flavour-violating processes themselves tend to dominate.

Two-body $`K^0`$ decays to muons and electrons proceed via the diagrams shown in Fig. 3.

At the quark level, the effective Lagrangian for such processes has the form 333Note that, for energies of the order of the kaon mass, both $`s`$\- and $`t`$-channel diagrams yield contact interactions. :

$`_{d\overline{s}\mathrm{}_j^{}\mathrm{}_k^+}`$ $`=`$ $`{\displaystyle \frac{1}{m_{\stackrel{~}{\nu }_i}^2}}\left[\lambda _{ijk}^{}\lambda _{i12}^{}\left(\overline{s_R}d_L\right)\left(\overline{\mathrm{}_{j}^{}{}_{L}{}^{}}\mathrm{}_{k}^{}{}_{R}{}^{}\right)+\lambda _{ikj}\lambda _{i21}^{}\left(\overline{s_L}d_R\right)\left(\overline{\mathrm{}_{j}^{}{}_{R}{}^{}}\mathrm{}_{k}^{}{}_{L}{}^{}\right)\right]`$ (24)
$``$ $`{\displaystyle \frac{\lambda _{ji1}^{}\lambda _{ki2}^{}}{2m_{\stackrel{~}{u}_i}^2}}\left(\overline{s_R}\gamma _\mu d_R\right)\left(\overline{\mathrm{}_{j}^{}{}_{L}{}^{}}\gamma ^\mu \mathrm{}_{k}^{}{}_{L}{}^{}\right).`$
The two different contributions from $`s`$-channel and $`t`$-channel diagrams put limits on different couplings , as we discuss below.

We have derived Feynman rules for the relevant effective kaon, pion and lepton interactions. Based on these Feynman rules, we have recalculated the important kaon decay processes, and update the limits on the products of $`R`$-violating couplings using the present experimental limits (2). The diagrams of Fig. 3 lead to the the following effective Lagrangian for $`K^0\mathrm{}^+\mathrm{}^{}`$ interactions:

$`_{K^0\mathrm{}_j^{}\mathrm{}_k^+}`$ $`=`$ $`{\displaystyle \frac{F_{K^0}}{2m_{\stackrel{~}{\nu }_i}^2}}\left[\lambda _{ijk}^{}\lambda _{i12}^{}\left(\overline{\mathrm{}_{j}^{}{}_{L}{}^{}}\mathrm{}_{k}^{}{}_{R}{}^{}\right)\lambda _{ikj}\lambda _{i21}^{}\left(\overline{\mathrm{}_{j}^{}{}_{R}{}^{}}\mathrm{}_{k}^{}{}_{L}{}^{}\right)\right]K^0(p_K)`$ (25)
$``$ $`{\displaystyle \frac{f_K}{4m_{\stackrel{~}{u}_i}^2}}\lambda _{ji1}^{}\lambda _{ki2}^{}p_K^\mu \left(\overline{\mathrm{}_{j}^{}{}_{L}{}^{}}\gamma _\mu \mathrm{}_{k}^{}{}_{L}{}^{}\right)K^0(p_K),`$
where $`F_{K^0}=m_{K^0}^2f_K/(m_s+m_d)`$, $`m_s+m_d0.15`$ GeV is the sum of the current masses of the $`s`$ and $`d`$ quarks, and $`f_K=0.1598`$ GeV is the kaon decay constant. The value of $`F_{K^0}`$ is related to the pseudoscalar $`<0|\overline{s}\gamma ^5d|K^0>=F_{K^0}`$ matrix element, and is obtained from $`f_K`$ by using the Dirac equations for quarks. All QCD corrections are included in this phenomenological approach. In the following, we assume that the $`R`$-violating couplings are real and that only one of their products in (24) is non-zero.

We have implemented the Feynman rules in the CompHEP package , using the effective Lagrangian (25), and have obtained the following results:

$$\mathrm{\Gamma }_{K^0\mathrm{}_j^{}\mathrm{}_k^+}(\stackrel{~}{\nu }_i)=\frac{(\lambda _{ijk}\lambda _{i12}^{})^2F_{K^0}^2}{64\pi m_{\stackrel{~}{\nu }_i}^4m_{K^0}}\left(1\frac{m_\mathrm{}_j^2+m_\mathrm{}_k^2}{m_{K^0}^2}\right)\mathrm{\Delta }(m_{K^0},m_\mathrm{}_j,m_\mathrm{}_k)$$
(26)
$$\mathrm{\Gamma }_{K^0\mathrm{}_j^{}\mathrm{}_k^+}(\stackrel{~}{u}_i)=\frac{(\lambda _{ji1}^{}\lambda _{ki2}^{})^2f_K^2}{256\pi m_{\stackrel{~}{u}_i}^4m_{K^0}}\left(m_\mathrm{}_j^2+m_\mathrm{}_k^2\frac{(m_\mathrm{}_j^2m_\mathrm{}_k^2)^2}{m_{K^0}^2}\right)\mathrm{\Delta }(m_{K^0},m_\mathrm{}_j,m_\mathrm{}_k)$$
(27)
where $`\mathrm{\Delta }(a,b,c)=\sqrt{\left[a^2(b+c)^2\right]\left[a^2(bc)^2\right]}`$ is the triangle function. As there exist two similar contributions in the $`s`$ channel, coming from terms with different couplings $`\lambda _{ijk}\lambda _{i12}^{}`$ and $`\lambda _{ikj}\lambda _{i21}^{}`$, as seen in the first line in (25), we give only one of them in (26).

We obtain the following nominal numerical results for $`K^0\mathrm{}^+\mathrm{}^{}`$ decay via $`\stackrel{~}{\nu }`$ exchange, and the corresponding limits on the $`\lambda \lambda ^{}`$ products (for the numerical results for $`\mathrm{\Gamma }_{K^0\mathrm{}^+\mathrm{}^{}}`$ we have used the nominal values $`\lambda =\lambda ^{}=1`$ and $`m_{\stackrel{~}{\nu }(\stackrel{~}{u})}=100\text{ GeV}`$):

$`\mathrm{\Gamma }_{K^0e^+\mu ^{}}=1.57\times 10^{12}\text{ GeV},\lambda _{i21}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^26.2\times 10^9`$
$`\lambda _{i12}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^26.2\times 10^9`$
$`\mathrm{\Gamma }_{K^0e^{}\mu ^+}=1.57\times 10^{12}\text{ GeV},\lambda _{i12}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^26.2\times 10^9`$
$`\lambda _{i21}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^26.2\times 10^9`$
$`\mathrm{\Gamma }_{K^0e^+e^{}}=1.72\times 10^{12}\text{ GeV},\lambda _{i11}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^21.0\times 10^8`$
$`\lambda _{i11}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^21.0\times 10^8`$
$`\mathrm{\Gamma }_{K^0\mu ^+\mu ^{}}=1.42\times 10^{12}\text{ GeV},\lambda _{i22}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^22.6\times 10^7`$
$`\lambda _{i22}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^22.6\times 10^7`$ (28)
For $`K^0\mathrm{}^+\mathrm{}^{}`$ decay via up-squark exchange we have the following limits:

$`\mathrm{\Gamma }_{K^0e^+\mu ^{}}=1.61\times 10^{15}\text{ GeV},\lambda _{2i1}^{}\lambda _{1i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^21.9\times 10^7`$
$`\mathrm{\Gamma }_{K^0e^{}\mu ^+}=1.61\times 10^{15}\text{ GeV},\lambda _{1i1}^{}\lambda _{2i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^21.9\times 10^7`$
$`\mathrm{\Gamma }_{K^0e^+e^{}}=8.25\times 10^{20}\text{ GeV},\lambda _{1i1}^{}\lambda _{1i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^24.7\times 10^5`$
$`\mathrm{\Gamma }_{K^0\mu ^+\mu ^{}}=3.19\times 10^{15}\text{ GeV},\lambda _{2i1}^{}\lambda _{2i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^25.4\times 10^6`$ (29)
We have used in our calculations the decay width $`\mathrm{\Gamma }_{exp}(K_L^0)=1.273\times 10^{17}`$ GeV and the experimental limits on $`K^0\mathrm{}^+\mathrm{}^{}`$ decay widths shown in (2).

We now discuss the diagrams for 3-body kaon decays to pions and two charged leptons, of which there are two qualitatively different kinds:

* The kaon may decay into a pion of the same charge, in which case the leptons in the final state must have opposite signs: $`K^\pm \pi ^\pm \mathrm{}^{}\mathrm{}^\pm `$ and $`K^0\pi ^0\mathrm{}^{}\mathrm{}^\pm `$. The corresponding diagram for the first process is shown in Fig. 4.
* The kaon may decay into a pion with the opposite charge, in which case the leptons in the final state must have the same signs: $`K^\pm \pi ^{}\mathrm{}^\pm \mathrm{}^\pm `$. Representative diagrams for this process are shown in Fig. 5. This process involves two heavy virtual particles, the $`W`$ boson and a down squark. One should note that decay width of this process is directly proportional to the mixing between the left- and right-handed squark states, denoted by $`\stackrel{~}{b}_L`$ and $`\stackrel{~}{b}_R`$, respectively. If there is no mixing, this same-sign-lepton process mentioned vanishes. One can expect sizeable mixing only for squarks of the third generation (and especially in the high-$`\mathrm{tan}\beta `$ region), which is why we have used $`\stackrel{~}{b}_{L,R}`$ in the diagram.

The possibilities for $`K^\pm \pi ^\pm \mathrm{}^{}\mathrm{}^\pm `$ or $`K^0\pi ^0\mathrm{}^{}\mathrm{}^\pm `$ decay may be further subdivided into two groups.

* Diagrams involving only squarks, via which experimental upper limits bound products of $`LQ\overline{D}`$ operators.
* Diagrams involving also sleptons, which yield bounds on products of $`LQ\overline{D}`$ and $`LL\overline{E}`$ operators.

The diagrams of Fig. 4 lead to the the following effective Lagrangian for $`K^+\pi ^{}\mathrm{}^+\mathrm{}^{}`$ interactions:

$`_{K^+\pi ^{}\mathrm{}_k^{}\mathrm{}_j^+}`$ $`=`$ $`{\displaystyle \frac{M_s}{m_{\stackrel{~}{\nu }_i}^2}}\left[\lambda _{ijk}^{}\lambda _{i12}^{}\left(\overline{\mathrm{}_{j}^{}{}_{L}{}^{}}\mathrm{}_{k}^{}{}_{R}{}^{}\right)+\lambda _{ikj}\lambda _{i21}^{}\left(\overline{\mathrm{}_{j}^{}{}_{R}{}^{}}\mathrm{}_{k}^{}{}_{L}{}^{}\right)\right]K^+(p_K)\pi ^{}(p_\pi )`$ (30)
$`+`$ $`{\displaystyle \frac{f_+}{4m_{\stackrel{~}{u}_i}^2}}\lambda _{ji1}^{}\lambda _{ki2}^{}\left(p_K+p_\pi \right)^\mu \left(\overline{\mathrm{}_{j}^{}{}_{L}{}^{}}\gamma _\mu \mathrm{}_{k}^{}{}_{L}{}^{}\right)K^+(p_K)\pi ^{}(p_\pi ),`$
where $`M_s0.49`$ GeV is the constituent mass of the $`s`$ quark. In the general case, the vector matrix element of the $`K\pi `$ transition is parametrized by two momentum-dependent form factors $`f_+(t)`$ and $`f_{}(t)`$:

$`<0|\overline{s}\gamma ^\mu d|K^+(p_K)\pi ^{}(p_\pi )>=`$
$`<\pi ^+(p_\pi )|\overline{s}\gamma ^\mu d|K^+(p_K)>=(p_K+p_\pi )^\mu f_+(t)+(p_Kp_\pi )^\mu f_{}(t),`$
where $`t=(p_Kp_\pi )^2`$ is the squared momentum transfered to the lepton pair. The experimental data on semileptonic kaon decay are adequately described by a linear approximation for $`f_+`$:

$$f_+(t)=f_+(0)\left[1+\lambda _+(t/m_\pi ^2)\right],$$

where $`\lambda _+0.03`$ and $`f_{}0.31\pm 0.15`$. In the case of exact flavour $`SU(3)`$ symmetry, the following relations hold: $`f_+(0)=1`$ and $`f_{}(0)=0`$. Due to the Ademollo-Gatto theorem, the first form factor is renormalized only at second order in the $`SU(3)`$-violating interactions, and is therefore expected to be close to unity. The Ademollo-Gatto theorem is not applicable to $`f_{}`$, whose value is close to zero, as follows from experiments and the Callan-Treiman relation:

$$f_{}(m_K^2)=\frac{f_K}{f_\pi }f_+(m_K^2)0.15.$$

For our estimations, we set $`f_+(0)=1`$ and neglect $`f_{}`$. To estimate the scalar form factor, we use the relativistic quark model, which gives :

$$<0|\overline{s}d|K^+(p_K)\pi ^{}(p_\pi )>=<\pi ^+(p_\pi )|\overline{s}d|K^+(p_K)>2M_s,$$

where we keep only the leading term and drop the momentum dependence, and $`M_s0.49`$ GeV as before. In the approximation of unbroken $`SU(3)`$ symmetry, the corresponding form factors for the neutral kaons are smaller by a factor of $`\sqrt{2}`$.

We can now estimate the decay rates of the charged kaons:

$`\mathrm{\Gamma }_{K^+\pi ^+\mathrm{}_j^{}\mathrm{}_k^+}(\stackrel{~}{\nu }_i)={\displaystyle \frac{(\lambda _{ijk}\lambda _{i12}^{})^2M_s^2}{256\pi ^3m_{\stackrel{~}{\nu }_i}^4m_{K^+}^3}}{\displaystyle \underset{(m_\mathrm{}_j+m_\mathrm{}_k)^2}{\overset{(m_{K^+}m_{\pi ^+})^2}{}}}`$ $`\mathrm{\Delta }(t,m_\mathrm{}_j,m_\mathrm{}_k)\mathrm{\Delta }(t,m_{K^+},m_{\pi ^+})`$ (31)
$`\times \left(t^2m_\mathrm{}_j^2m_\mathrm{}_k^2\right){\displaystyle \frac{\text{d}t^2}{t^2}},`$
$`\mathrm{\Gamma }_{K^+\pi ^+\mathrm{}_j^{}\mathrm{}_k^+}(\stackrel{~}{u}_i)`$ $`=`$ $`{\displaystyle \frac{(\lambda _{ji1}^{}\lambda _{ki2}^{})^2}{2048\pi ^3m_{\stackrel{~}{u}_i}^4m_{K^+}^3}}{\displaystyle \underset{(m_\mathrm{}_j+m_\mathrm{}_k)^2}{\overset{(m_{K^+}m_{\pi ^+})^2}{}}}\mathrm{\Delta }(t,m_\mathrm{}_j,m_\mathrm{}_k)\mathrm{\Delta }(t,m_{K^+},m_{\pi ^+})`$ (32)
$`\times [(m_{K^+}^2+m_{\pi ^+}^2+m_\mathrm{}_j^2+m_\mathrm{}_k^2t^2)F_1(t)F_2(t)`$
$`+{\displaystyle \frac{m_\mathrm{}_j^2+m_\mathrm{}_k^2}{2}}(t^2m_\mathrm{}_j^2m_\mathrm{}_k^2)2m_{\pi ^+}^2m_{K^+}^2]`$
$`\times \left[1+\lambda _+(t^2/m_\pi ^2)\right]^2{\displaystyle \frac{\text{d}t^2}{t^2}},`$
where $`F_1(t)=q_{max}(t)+q_{min}(t)`$, $`F_2(t)=2(q_{max}^2(t)+q_{max}(t)q_{min}(t)+q_{min}^2(t))/3`$, and

$`q_{max}(t)={\displaystyle \frac{1}{4t^2}}\left[\left(m_{K^+}^2m_{\pi ^+}^2+m_\mathrm{}_j^2m_\mathrm{}_k^2\right)^2\left(\mathrm{\Delta }(t,m_\mathrm{}_j,m_\mathrm{}_k)\mathrm{\Delta }(t,m_{K^+},m_{\pi ^+})\right)^2\right],`$
$`q_{min}(t)={\displaystyle \frac{1}{4t^2}}\left[\left(m_{K^+}^2m_{\pi ^+}^2+m_\mathrm{}_j^2m_\mathrm{}_k^2\right)^2\left(\mathrm{\Delta }(t,m_\mathrm{}_j,m_\mathrm{}_k)+\mathrm{\Delta }(t,m_{K^+},m_{\pi ^+})\right)^2\right].`$
For the decay rates of the neutral kaons, the $`SU(3)`$ factor of 1/2 and the corresponding masses for the neutral mesons $`m_{K^0}`$ and $`m_{\pi ^0}`$ must be taken into account.

We obtain the following nominal numerical results for $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$ and $`K^0\pi ^0\mathrm{}^+\mathrm{}^{}`$ decays via $`\stackrel{~}{\nu }`$ exchange, and the corresponding limits on the $`\lambda \lambda ^{}`$ products (for the numerical results we have used the nominal values $`\lambda =\lambda ^{}=1`$ and $`m_{\stackrel{~}{\nu }(\stackrel{~}{u})}=100\text{ GeV}`$):

$`\mathrm{\Gamma }_{K^+\pi ^+e^+\mu ^{}}=1.38\times 10^{15}\text{ GeV},\lambda _{i21}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^24.5\times 10^6`$
$`\lambda _{i12}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^24.5\times 10^6`$
$`\mathrm{\Gamma }_{K^0\pi ^0e^+\mu ^{}}=7.71\times 10^{16}\text{ GeV},\lambda _{i21}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^27.2\times 10^6`$
$`\lambda _{i12}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^27.2\times 10^6`$
$`\mathrm{\Gamma }_{K^+\pi ^+e^{}\mu ^+}=1.38\times 10^{15}\text{ GeV},\lambda _{i12}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^21.0\times 10^6`$
$`\lambda _{i21}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^21.0\times 10^6`$
$`\mathrm{\Gamma }_{K^0\pi ^0e^{}\mu ^+}=7.71\times 10^{16}\text{ GeV},\lambda _{i12}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^27.2\times 10^6`$
$`\lambda _{i21}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^27.2\times 10^6`$
$`\mathrm{\Gamma }_{K^+\pi ^+e^+e^{}}=2.14\times 10^{15}\text{ GeV},\lambda _{i11}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^28.8\times 10^5`$
$`\lambda _{i11}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^28.8\times 10^5`$
$`\mathrm{\Gamma }_{K^0\pi ^0e^+e^{}}=1.17\times 10^{15}\text{ GeV},\lambda _{i11}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^26.8\times 10^6`$
$`\lambda _{i11}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^26.8\times 10^6`$
$`\mathrm{\Gamma }_{K^+\pi ^+\mu ^+\mu ^{}}=7.58\times 10^{16}\text{ GeV},\lambda _{i22}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^28.2\times 10^5`$
$`\lambda _{i22}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^28.2\times 10^5`$
$`\mathrm{\Gamma }_{K^0\pi ^0\mu ^+\mu ^{}}=4.38\times 10^{16}\text{ GeV},\lambda _{i22}\lambda _{i12}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^21.2\times 10^5`$
$`\lambda _{i22}\lambda _{i21}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{\nu }}}}\right)^21.2\times 10^5`$ (33)
For the decays via up-squark exchange we have the following limits:

$`\mathrm{\Gamma }_{K^+\pi ^+e^+\mu ^{}}=1.61\times 10^{16}\text{ GeV},\lambda _{2i1}^{}\lambda _{1i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^21.3\times 10^5`$
$`\mathrm{\Gamma }_{K^0\pi ^0e^+\mu ^{}}=9.06\times 10^{17}\text{ GeV},\lambda _{2i1}^{}\lambda _{1i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^22.1\times 10^5`$
$`\mathrm{\Gamma }_{K^+\pi ^+e^{}\mu ^+}=1.61\times 10^{16}\text{ GeV},\lambda _{1i1}^{}\lambda _{2i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^23.0\times 10^6`$
$`\mathrm{\Gamma }_{K^0\pi ^0e^{}\mu ^+}=9.06\times 10^{17}\text{ GeV},\lambda _{1i1}^{}\lambda _{2i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^22.1\times 10^5`$
$`\mathrm{\Gamma }_{K^+\pi ^+e^+e^{}}=2.38\times 10^{16}\text{ GeV},\lambda _{1i1}^{}\lambda _{1i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^22.7\times 10^4`$
$`\mathrm{\Gamma }_{K^0\pi ^0e^+e^{}}=1.31\times 10^{16}\text{ GeV},\lambda _{1i1}^{}\lambda _{1i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^22.0\times 10^5`$
$`\mathrm{\Gamma }_{K^+\pi ^+\mu ^+\mu ^{}}=9.33\times 10^{17}\text{ GeV},\lambda _{2i1}^{}\lambda _{2i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^22.3\times 10^4`$
$`\mathrm{\Gamma }_{K^0\pi ^0\mu ^+\mu ^{}}=5.48\times 10^{17}\text{ GeV},\lambda _{2i1}^{}\lambda _{2i2}^{}\times \left({\displaystyle \frac{100\text{ GeV}}{m_{\stackrel{~}{u}}}}\right)^23.4\times 10^5`$ (34)
We have used in our calculations the decay width $`\mathrm{\Gamma }_{exp}(K^+)=5.314\times 10^{17}`$ GeV and the present limits on $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$ and $`K_L^0\pi ^0\mathrm{}^+\mathrm{}^{}`$ decay widths (2). Note that there are different constraints from the $`K^+\pi ^+e^+\mu ^{}`$ and $`K^+\pi ^+e^{}\mu ^+`$ decays, because of the rather different experimental limits $`BR(K^+\pi ^+e^+\mu ^{})6.9\times 10^9`$ and $`BR(K^+\pi ^+e^{}\mu ^+)2.8\times 10^{11}`$. The limit obtained from $`K^0\mathrm{}^+\mathrm{}^{}`$ is typically 1-2 orders of magnitude better than that derived from $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$ decay.

As already mentioned, diagrams with non-zero sbottom-quark mixing may lead to like-sign leptons, as seen in Fig. 5. They arise from the effective lepton-number violating contact interactions

$$=\frac{\lambda _{ikp}^{}\lambda _{jqk}^{}V_{LR}}{m_{d_k}^2}\left(\overline{d_{p}^{}{}_{R}{}^{}}\nu _{i}^{}{}_{L}{}^{}\right)\left(\overline{(\mathrm{}_{j}^{}{}_{L}{}^{})^C}u_{q}^{}{}_{L}{}^{}\right)+\text{h.c.},$$
(35)
where $`V_{LR}`$ denotes left-right squark mixing matrix element.

There are two different topologies, shown in the first and second row of Fig. 5, respectively. Diagrams in the first row lead to an effective tree-level process, since the $`W`$-boson and squark masses are much bigger then the typical energy scale, which is of the order of the kaon mass. Diagrams in the second row cannot be reduced to tree-level diagrams, since one has a neutrino propagator in the loop. However, this last set of diagrams give a contribution that is typically about 2 orders of magnitude lower than the diagrams of the first row . Therefore, for the sake of simplicity, we neglect the second-row diagrams, and have derived an effective Lagrangian only for the first two diagrams of Fig. 5. Two kinds of effective interaction appear: Standard-Model-like $`K\mathrm{}\nu `$ and $`\pi \mathrm{}\nu `$ interactions and new effective interactions related to $`R`$-violating operators of the forms $`K\mathrm{}^C\nu `$ and $`\pi \mathrm{}^C\nu `$.

The effective Lagrangian for those interactions take the following forms:

$$_{K^+\mathrm{}_i^{}\overline{\nu }_i}=V_{us}\sqrt{2}G_Ff_Kp_K^\mu \left(\overline{\nu _{i}^{}{}_{L}{}^{}}\gamma _\mu \mathrm{}_{i}^{}{}_{L}{}^{}\right)K^+(p_K)$$
(36)
$$_{\pi ^+\mathrm{}_i^{}\overline{\nu }_i}=V_{ud}\sqrt{2}G_Ff_\pi p_\pi ^\mu \left(\overline{\nu _{i}^{}{}_{L}{}^{}}\gamma _\mu \mathrm{}_{i}^{}{}_{L}{}^{}\right)\pi ^+(p_\pi )$$
(37)
where $`f_K`$ and $`f_\pi =0.1307`$ GeV are the kaon and pion decay constants, respectively, $`G_F=1.1663910^5`$ GeV-2 is the Fermi constant and $`V_{us},V_{ud}`$ are CKM matrix elements, and:

$$_{K^+(\mathrm{}_j^{})^C\nu _i}(\stackrel{~}{d}_k)=\frac{\lambda _{ik2}^{}\lambda _{j1k}^{}V_{LR}}{4m_{\stackrel{~}{d}_k}^2}F_{K^+}\left(\overline{(\mathrm{}_{j}^{}{}_{L}{}^{})^C}\nu _{i}^{}{}_{L}{}^{}\right)K^+(p_K)$$
(38)
$$_{\pi ^+(\mathrm{}_j^{})^C\nu _i}(\stackrel{~}{d}_k)=\frac{\lambda _{ik1}^{}\lambda _{j1k}^{}V_{LR}}{4m_{\stackrel{~}{d}_k}^2}F_{\pi ^+}\left(\overline{(\mathrm{}_{j}^{}{}_{L}{}^{})^C}\nu _{i}^{}{}_{L}{}^{}\right)\pi ^+(p_\pi )$$
(39)
where $`F_{K^+}=m_{K^+}^2f_K/(m_s+m_u)`$, $`F_{\pi ^+}=m_{\pi ^+}^2f_\pi /(m_d+m_u)`$, $`m_s+m_u0.15`$ GeV, and $`m_d+m_u0.01`$ GeV.

The Feynman diagrams for $`K^+\pi ^{}\mathrm{}_i^+\mathrm{}_j^+`$ decay are given in terms of these effective interactions, as shown in Fig. 6. They yield the matrix elements:

$$=\frac{\lambda _{ik1}^{}\lambda _{j1k}^{}V_{LR}F_{\pi ^+}}{2m_{\stackrel{~}{d}_k}^2}V_{us}\frac{G_F}{\sqrt{2}}f_Km_\mathrm{}_i\left(\overline{(\mathrm{}_{j}^{}{}_{L}{}^{})^C}\frac{\widehat{p}_{\nu _i}}{p_{\nu _i}^2}\mathrm{}_{i}^{}{}_{R}{}^{}\right)K^+(p_K)\pi ^{}(p_\pi )$$
(40)
and

$$=\frac{\lambda _{ik2}^{}\lambda _{j1k}^{}V_{LR}F_{K^+}}{2m_{\stackrel{~}{d}_k}^2}V_{ud}\frac{G_F}{\sqrt{2}}f_\pi m_\mathrm{}_i\left(\overline{(\mathrm{}_{j}^{}{}_{L}{}^{})^C}\frac{\widehat{p}_{\nu _i}}{p_{\nu _i}^2}\mathrm{}_{i}^{}{}_{R}{}^{}\right)K^+(p_K)\pi ^{}(p_\pi ).$$
(41)
Recalling that $`V_{us}/V_{ud}0.2`$, $`F_{K^+}F_{\pi ^+}`$ and $`f_Kf_\pi `$, the main contribution comes from the last matrix element with $`i=2`$, because the matrix element is chirally suppressed in the other case. As an example, one can obtain a constraint on the product of $`\lambda _{2k2}^{}\lambda _{11k}^{}`$ and $`V_{LR}`$ from $`K^+\pi ^{}e^+\mu ^+`$ decay. In this case, we have the following numerical result 444 In  a similar constraint was found with a different choice of diagrams.:

$$\mathrm{\Gamma }(K^+\pi ^{}e^+\mu ^+)=V_{LR}^2(\lambda _{2k2}^{}\lambda _{11k}^{})^2\times \left(\frac{100\text{GeV}}{m_{\stackrel{~}{d}_k}}\right)^4\times 2.5\times 10^{28}\text{ GeV}$$
(42)
or

$$V_{LR}(\lambda _{2k2}^{}\lambda _{11k}^{})\times \left(\frac{100\text{GeV}}{m_{\stackrel{~}{d}_k}}\right)^210.$$
(43)
It is apparent that kaon decay into a pion and a like-sign lepton pair is too strongly suppressed to be useful at present: the corresponding bounds on the $`\lambda ^{}\lambda ^{}`$ product are currently of the order of 100.

## 4 Conclusions

We have discussed in this paper flavour-violating decays of kaons into charged-lepton pairs in supersymmetric theories, in both the Minimal Supersymmetric Standard Model and $`R`$-violating models. In the first case, these decays are generated by box diagrams involving charginos and neutralinos, and both the squark and the slepton mixings enter in the analysis. The process looks promising for correlating the (s)-quark and (s)-lepton mixing by a combined study of rare charged lepton and kaon decays. Despite the limits from $`\mu e\gamma `$, $`\mu `$$`e`$ conversion and $`\mathrm{\Delta }m_K`$, the kaon decay branching ratios for large $`\mathrm{tan}\beta `$ and small soft supersymmetry-breaking terms may be accessible to a future generation of experiments using new intense proton sources.

In the case of $`R`$-violating supersymmetry, such rare kaon decays may occur at tree level. In this case, $`\mu e\gamma `$ again occurs via one-loop diagrams, whilst $`\mu `$$`e`$ conversion may also occur at tree-level, but via a set of operators different from those relevant to kaon decays. In this framework, we studied the expected rates for the decays $`K\mu ^\pm e^{}`$ and $`K\pi \mu e`$, for all two- and three-body processes. Using the current experimental data, we obtained the bounds on products of $`LL\overline{E}`$ and $`LQ\overline{D}`$ operators summarized in (28) and (29). We have also noted the possibility of like-sign lepton events in the presence of non-zero $`\stackrel{~}{b}_L`$$`\stackrel{~}{b}_R`$ mixing, but for this to occur at significant rate one would need large $`R`$-violating couplings.

Our final conclusion is that lepton-flavour-violating rare kaon decays have the potential to provide important information on the issue of flavour physics. Any future observation would, in addition, help distinguish between different supersymmetric theories.

Acknowledgements: We thank Gerhard Buchalla for useful discussions.

Appendix: The rate for $`K\mu e`$ decay via box diagrams in the MSSM

The branching ratio for $`K\mu e`$ is given by:

$`BR(K\mu e)={\displaystyle \frac{2.65\lambda ^2}{2G_F\mathrm{sin}^2\theta _W}}(|K_L|^2+|K_R|^2)`$ (44)
where $`\lambda =(m_d+m_s)m_\mu /m_K^2`$, and $`K_L`$, $`K_R`$ are given by the following expressions:

$`K_L`$ $`=`$ $`K_L^c+(K_L^{n(1)}+K_L^{n(2)})`$ (45)
$`K_R`$ $`=`$ $`K_R^c+(K_R^{n(1)}+K_R^{n(2)})`$ (46)
with

$`K_L^c`$ $`=`$ $`{\displaystyle \frac{1}{4}}J_{4(A,B,X,Y)}\left({\displaystyle \frac{\lambda }{2}}C_{dAX}^{R(d)}C_{sBX}^{R(d)}C_{\mu BY}^{R(\mathrm{})}C_{eAY}^{R(\mathrm{})}+C_{dAX}^{R(d)}C_{sBX}^{L(d)}C_{\mu BY}^{L(\mathrm{})}C_{eAY}^{R(\mathrm{})}\right)`$
$`{\displaystyle \frac{1}{4}}I_{4(A,B,X,Y)}m_{\stackrel{~}{\chi }_A^{}}m_{\stackrel{~}{\chi }_B^{}}\left(C_{dAX}^{L(d)}C_{sBX}^{R(d)}C_{\mu BY}^{L(\mathrm{})}C_{eAY}^{R(\mathrm{})}+\lambda C_{dAX}^{L(d)}C_{sBX}^{L(d)}C_{\mu BY}^{R(\mathrm{})}C_{eAY}^{R(\mathrm{})}\right)`$
$`K_R^c`$ $`=`$ $`K_L^c|_{LR}`$ (48)
where $`m_{\stackrel{~}{\chi }_{A,B}}`$ and $`m_{\stackrel{~}{l}_{X,Y}}`$ denote chargino and sneutrino masses, in the chargino contribution. Moreover,

$`iJ_{4(A,B,X,Y)}`$ $`=`$ $`{\displaystyle \frac{d^4k}{(2\pi )^4}\frac{k^2}{(k^2M_{\stackrel{~}{\chi }_A}^2)(k^2M_{\stackrel{~}{\chi }_B}^2)(k^2m_{\stackrel{~}{l}_X}^2)(k^2m_{\stackrel{~}{l}_Y}^2)}}`$ (49)
$`iI_{4(A,B,X,Y)}`$ $`=`$ $`{\displaystyle \frac{d^4k}{(2\pi )^4}\frac{1}{(k^2M_{\stackrel{~}{\chi }_A}^2)(k^2M_{\stackrel{~}{\chi }_B}^2)(k^2m_{\stackrel{~}{l}_X}^2)(k^2m_{\stackrel{~}{l}_Y}^2)}}.`$ (50)
whilst the mixing coefficients $`C^{R,L}`$ appear in the fermion-sfermion-chargino interaction Lagrangian:

$`_{\mathrm{int}}`$ $`=`$ $`\overline{\mathrm{}}_i(C_{iAX}^{R(l)}P_R+C_{iAX}^{L(l)}P_L)\stackrel{~}{\chi }_A^{}\stackrel{~}{\nu }_X`$ (51)
$`+\overline{\nu }_i(C_{iAX}^{R(\nu )}P_R+C_{iAX}^{L(\nu )}P_L)\stackrel{~}{\chi }_A^+\stackrel{~}{\mathrm{}}_X`$
$`+\overline{d}_i(C_{iAX}^{R(d)}P_R+C_{iAX}^{L(d)}P_L)\stackrel{~}{\chi }_A^{}\stackrel{~}{u}_X`$
$`+\overline{u}_i(C_{iAX}^{R(u)}P_R+C_{iAX}^{L(u)}P_L)\stackrel{~}{\chi }_A^+\stackrel{~}{d}_X+h.c.`$
and their explicit expressions are given, for instance, in .

The fermion-sfermion-neutralino interaction Lagrangian is similarly written as

$$_{\mathrm{int}}=\overline{f}_i(N_{iAX}^{R(f)}P_R+N_{iAX}^{L(f)}P_L)\stackrel{~}{\chi }_A^0\stackrel{~}{f}_X$$
(52)
where $`f`$ stands for $`l,\nu ,d`$ and $`u`$. The neutralino box contributions corresponding to the permutations of the $`\mu `$ and $`e`$ in the external lines, $`K_L^n(2)`$ and $`K_R^n(2)`$, are then found to be

$`K_L^{n(1)}`$ $`=`$ $`{\displaystyle \frac{1}{4}}J_{4(A,B,X,Y)}\left({\displaystyle \frac{\lambda }{2}}N_{dAX}^{R(d)}N_{sBX}^{R(d)}N_{\mu BY}^{R(\mathrm{})}N_{eAY}^{R(\mathrm{})}+N_{dAX}^{R(d)}N_{sBX}^{L(d)}N_{\mu BY}^{L(\mathrm{})}N_{eAY}^{R(\mathrm{})}\right)`$
$`{\displaystyle \frac{1}{4}}I_{4(A,B,X,Y)}m_{\stackrel{~}{\chi }_A^0}m_{\stackrel{~}{\chi }_B^0}\left(N_{dAX}^{L(d)}N_{sBX}^{R(d)}N_{\mu BY}^{L(\mathrm{})}N_{eAY}^{R(\mathrm{})}+\lambda N_{dAX}^{L(d)}N_{sBX}^{L(d)}N_{\mu BY}^{R(\mathrm{})}N_{eAY}^{R(\mathrm{})}\right)`$
$`K_R^{n(1)}`$ $`=`$ $`K_L^{n(1)}|_{LR}`$ (54)
and

$`K_L^{n(2)}`$ $`=`$ $`{\displaystyle \frac{1}{4}}J_{4(A,B,X,Y)}\left({\displaystyle \frac{\lambda }{2}}N_{dAX}^{R(d)}N_{sBX}^{R(d)}N_{\mu BY}^{L(\mathrm{})}N_{eAY}^{L(\mathrm{})}N_{dAX}^{R(d)}N_{sBX}^{L(d)}N_{\mu BY}^{L(\mathrm{})}N_{eAY}^{R(\mathrm{})}\right)`$
$`+{\displaystyle \frac{1}{4}}I_{4(A,B,X,Y)}m_{\stackrel{~}{\chi }_A^0}m_{\stackrel{~}{\chi }_B^0}\left(N_{dAX}^{L(d)}N_{sBX}^{R(d)}N_{\mu BY}^{L(\mathrm{})}N_{eAY}^{R(\mathrm{})}+\lambda N_{dAX}^{L(d)}N_{sBX}^{L(d)}N_{\mu BY}^{L(\mathrm{})}N_{eAY}^{L(\mathrm{})}\right)`$
$`K_R^{n(2)}`$ $`=`$ $`K_L^{n(2)}|_{LR}.`$ (56)


format: text

text:
# Abundance profiles of CH3OH and H2CO toward massive young stars as tests of gas-grain chemical models

## 1 Introduction

One of the early stages of massive star formation is the “hot core” phase (Kurtz et al. 2000, Macdonald & Thompson (2000)), which is characterized by masses of $`100`$ M of molecular gas at temperatures of $`\genfrac{}{}{0pt}{}{_>}{^{}}100`$ K, leading to a rich line spectrum at submillimeter wavelengths like that of the prototypical Orion hot core. The chemical composition of hot cores is quite distinct for its high abundances of fully hydrogenated, large carbon-bearing molecules (Ohishi (1997)), as opposed to dark clouds which show predominantly unsaturated species and molecular ions (van Dishoeck & Blake (1998)). Saturated carbon chains are not expected in large quantities by steady-state gas-phase models, which led Charnley et al. (1992) and Caselli et al. (1993) to propose models where these molecules are made on dust grains, and evaporate into the gas phase when the newly-formed star heats up its surroundings. This leads to a short period ($`10^4`$ yr) when the envelope of the young star is rich in complex organic molecules. Grain surface reactions at low temperature lead to hydrogenation if atomic H is able to quantum tunnel through activation barriers as proposed by Tielens & Hagen (1982). This model makes specific predictions for the abundance ratios of CO, H2CO, CH3OH and their deuterated versions which for a reasonable choice of parameters fit the abundances observed in the Compact Ridge in Orion as shown by Charnley et al. (1997). However, it is unknown if other sources can be fit as well. Alternatively, processing by ultraviolet light or energetic particles may change the ice composition in the vicinity of young stars.

One key molecule to test models of hot core chemistry is methanol, CH3OH. At temperatures $`\genfrac{}{}{0pt}{}{_<}{^{}}100`$ K, production of CH3OH in the gas phase goes mainly via radiative association of CH$`{}_{}{}^{+}{}_{3}{}^{}`$ and H2O, an inefficient process yielding abundances of only $`10^{11}`$ relative to H2 (Lee et al. 1996). Observations of CH3OH toward dark and translucent clouds by Turner (1998) and Takakuwa et al. (1998) yield abundances of $`10^9`$, which suggests that grain chemistry operates. Toward the Orion hot core and compact ridge, CH3OH abundances of $`10^7`$ have been observed by Blake et al. (1987) and Sutton et al. (1995), and even higher abundances, $`10^6`$, have been inferred from observations of CH3OH masers (e.g., Menten et al. (1986)). Methanol has also been observed in the solid state at infrared wavelengths from the ground by, e.g., Grim et al. (1991) and Allamandola et al. (1992) and recently by Dartois et al. (1999b), at abundances up to $`10^6`$.

To study the organic chemistry in warm ($`30200`$ K) and dense ($`10^410^7`$ cm-3) circumstellar envelopes, this paper discusses observations of CH3OH and H2CO lines at submillimeter wavelengths. Although gas-grain interactions are clearly important for methanol, it is unknown how efficiently (if at all) surface reactions modify the composition of the ices. In addition, it is still open what process returns the molecules to the gas phase: thermal heating, shocks or both? Finally, the ice layer must have formed at temperatures and densities much lower than the present situation, and hence provides a unique fossil record of the conditions in the molecular cloud prior to star formation.

The sources are thirteen massive ($`L=10^310^5`$ L) stars, which are at an early stage of evolution and still embedded in $`10^210^3`$ M of dust and molecular gas. A few of the sources are hot cores by the definition of Kurtz et al. (2000), but most are in an even younger phase where most of the envelope is still at low temperatures. The physical structure of the soures has been studied by van der Tak et al. (2000), who developed detailed temperature and density profiles. The wealth of CH3OH and H2CO lines, combined with the physical models, allows the determination of abundance profiles for H2CO and CH3OH, which demonstrate that CH3OH evaporates off dust particles in the envelope on a $`10^{45}`$ yr time scale, while H2CO is predominantly formed in the gas phase. The excitation and abundance of CH3OH are therefore useful evolutionary indicators during the embedded stage of star formation. Moreover, by comparing the observed amounts of evaporated CH3OH and CO2 to a model of grain surface chemistry, we set limits on the density and the duration of the pre-protostellar phase.

## 2 Observations

### 2.1 Single-dish observations

Spectroscopy of the $`J=54`$ and $`76`$ bands of CH3OH near $`242`$ and $`338`$ GHz was performed with the 15-m James Clerk Maxwell Telescope (JCMT)111The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre, on behalf of the Particle Physics and Astronomy Research Council of the United Kingdom, the Netherlands Organization for Scientific Research and the National Research Council of Canada. on Mauna Kea, Hawaii in May and October of 1995. The beam size (FWHM) and main beam efficiency of the antenna were $`18^{\prime \prime }`$ and $`69`$% at $`242`$ GHz and $`14^{\prime \prime }`$ and $`58`$% at $`338`$ GHz. The frontends were the receivers A2 and B3i; the backend was the Digital Autocorrelation Spectrometer, covering $`500`$ MHz instantaneous bandwidth. Pointing was checked every 2 hours during the observing and was always found to be within $`5^{\prime \prime }`$. To subtract the atmospheric and instrumental background, a reference position $`180^{\prime \prime }`$ East was observed. Integration times are 30–40 minutes for each frequency setting, resulting in rms noise levels in $`T_{\mathrm{mb}}`$ per $`625`$ kHz channel of $`30`$ mK at 242 GHz to $`50`$ mK at 338 GHz. Although the absolute calibration is only correct to $`30`$%, the relative strength of lines within either frequency setting is much more accurate.

In addition to these CH3OH lines, this paper discusses observations of lines of H2CO, obtained in a similar manner, which were presented in van der Tak et al. (2000). We will also use data on both molecules for W 3 IRS5 and W 3 (H2O) from the survey by Helmich & van Dishoeck (1997), and CH3OH $`J=76`$ data from Thompson & MacDonald (1999).

### 2.2 Interferometric observations

Maps of the $`J_K=3_14_0`$ A+ (107013.852 MHz; $`E_u=28.01`$ K) and $`11_110_2`$ E (104300.461 MHz; $`E_u=157.24`$ K) lines of CH3OH were obtained for four sources in 1995-1999 with the millimeter interferometer of the Owens Valley Radio Observatory (OVRO)222The Owens Valley Millimeter Array is operated by the California Institute of Technology under funding from the U.S. National Science Foundation (AST-9981546).. This instrument consists of six 10.4 m antennas on North-South and East-West baselines, and a detailed technical description is given in Padin et al. (1991). For the sources GL 2591, NGC 7538 IRS1 and NGC 7538 IRS9, data were taken in the compact and extended configurations, while for the Southern source W 33A, a hybrid configuration with long North-South but short East-West baselines was also used to improve the beam shape. Antenna spacings range from the shadowing limit out to 90 k$`\lambda `$, corresponding to spatial frequencies of $`0.0180.44`$ arcsec-1. The observations were carried out simultaneously with continuum observations, which have been presented in van der Tak et al. (1999, 2000). We refer the reader to these papers for further observational details.

## 3 Results

### 3.1 Submillimeter spectroscopic results

Figure 1 presents the calibrated JCMT spectra, reduced using the IRAM CLASS package. All lines are due to either A- or E-type CH3OH, while some features are blends of lines. After subtraction of a linear baseline, Gaussian profiles were fitted to extract line parameters. The free parameters in these fits were the line flux $`T_{\mathrm{mb}}𝑑V`$, the FWHM line width and the central velocity. The retrieved line fluxes and widths are listed in Table 1, while the central velocities are consistent with those found for C17O and C34S in van der Tak et al. (2000).

For most sources, the line widths are also consistent with those found for C17O and C34S. However, for the two lowest-luminosity sources, GL 7009S and IRAS 20126, the methanol lines are much broader than those of C17O and C34S. The CH3OH line profiles toward these sources (Fig. 2) suggest an origin in their molecular outflows, especially in GL 7009S. The methanol emission toward young low-mass stars such as L 1157 (Bachiller et al. (1995)) and NGC 1333 IRAS 4A (Blake et al. (1995)) is also dominated by their outflows. The other sources studied here possess outflows, but those do not seem to be important for the CH3OH emission. As we will show below, the observed gas-phase CH3OH is mostly due to evaporation of grain mantles. The CH3OH line profiles thus suggest that around low-mass stars, the desorption of icy grain mantles is dominated by shocks, while radiative heating dominates for high-mass stars.

The sources exhibit clear variations in their methanol spectra, both in absolute line strength and in the relative strength between lines. Variations in absolute line strength can be due to differences in molecular abundances, the size of the emitting region or distance, while the line ratios reflect differences in the temperature and density of the emitting gas. Before modeling the observations in detail, we perform a quick-look examination using rotation diagrams, as in Blake et al. (1987) and Helmich et al. (1994). This analysis assumes that the emission is optically thin and that the population of the molecular energy levels can be described by a single temperature, the “rotational temperature”. Since the envelopes of these stars are known to have temperature and density gradients, the rotation temperature gives information where the molecule is preferentially located in the envelope. The nondetection of 13CH3OH lines limits the optical depth of the lines to $`\genfrac{}{}{0pt}{}{_<}{^{}}3`$, except in the case of NGC 6334 IRS1, where the strongest lines are probably saturated.

Table 2 presents the results of the rotation diagram analysis, which are averages over the JCMT beam. Formal errors from the least-square fit to the data and calibration errors combine into uncertainties of $`20`$% in both parameters. The results for IRAS 20126, W 3 (H2O) and NGC 6334 IRS1 are more uncertain, to $`50`$%, because the derived $`T_{\mathrm{rot}}`$ is close to the highest observed energy level. The OVRO maps presented in § 3.2 indicate that the CH3OH column densities are beam diluted by a factor of $`25`$; this information is not available for H2CO. Although the column densities of the two species are roughly correlated, with $`N`$(CH3OH) $`10\times N`$(H2CO), their excitation temperatures behave differently. While $`T_{\mathrm{rot}}`$(H2CO) is rather constant from source to source, $`6090`$ K, except toward the hot cores W 3 (H2O) and NGC 6334 IRS1, the CH3OH excitation temperatures span the full range $`10200`$ K.

The spectral line survey of Thompson & MacDonald (1999) toward W 28 A2 includes the CH3OH $`J=76`$ band, and their analysis gives $`T_{\mathrm{rot}}=49`$ K and $`N=3.8\times 10^{15}`$ cm-2, in good agreement with our results from the $`J=54`$ band. Our values of $`N`$(H2CO) toward W 3 (H2O) and NGC 6334 IRS1 are a factor of $``$10 higher than those found by Mangum & Wootten (1993), probably due to the smaller JCMT beam size.

Figure 3 compares the rotational temperatures of CH3OH and H2CO to the C2H2 excitation temperatures measured by Lahuis & van Dishoeck (2000) in absorption at $`13.7`$ $`\mu `$m for the same lines of sight. Lacking a dipole moment, C2H2 is a clean probe of temperature. The tight correlation between the CH3OH and C2H2 excitation temperatures implies that the two species trace the same gas, and that source orientation and optical depth effects do not influence the infrared data. The only exception is W 33A: this massive source becomes optically thick at $`14`$ $`\mu `$m before the warm gas is reached. In contrast, the H2CO temperature is not correlated with that of C2H2: omitting W 33A, the correlation coefficient is $`7.8`$%, versus $`98.3`$% for CH3OH. This difference indicates that in these sources, H2CO is not chemically related to CH3OH and C2H2, unlike in low-mass objects, as the high H2CO abundance in the L 1157 outflow (Bachiller & Pérez Gutiérrez (1997)) and the detection of HDCO and D2CO in IRAS 16293 (van Dishoeck et al. (1995); Ceccarelli et al. (1998)) indicate.

### 3.2 Interferometric maps of CH3OH line emission

Figure 4 presents the OVRO observations, as maps of the emission integrated over the line profile. These maps were obtained by a gridding and fast Fourier transform of the $`(u,v)`$-data with natural weighting, and deconvolution with the CLEAN algorithm. The self-calibration solutions obtained for the continuum data were used to suppress phase noise introduced by the atmosphere. Also shown are spectra at the image maxima, extracted from image cubes at full spectral resolution and coverage. In these maps and spectra, the continuum emission has not been subtracted. Only in the case of NGC 7538 IRS1, the continuum brightness is comparable to the line brightness; for the other sources, the continuum can be neglected.

The maps show single, compact emitting sources, which coincide within the errors with their counterparts in $`86230`$ GHz continuum emission (van der Tak et al. 1999, 2000). No CH3OH line emission was detected toward NGC 7538 IRS9 to the $`1\sigma `$ noise level of $`1`$ K. In the case of W 33A, a binary source at these wavelengths, the CH3OH emission is associated with the source MM1 from van der Tak et al. (2000) and with the infrared source. Although the CH3OH emission is seen to be slightly extended toward the West, no emission was detected at the position of the other continuum source, MM2.

Table 3 lists various parameters of the lines observed with OVRO. Deconvolved source sizes were obtained by fitting two-dimensional Gaussian profiles to the integrated emission maps; other parameters by fitting Gaussians to the presented spectra. The continuum emission was subtracted from the data before this fitting procedure.

The measured central velocities and line widths are consistent with the values measured with the JCMT. To check if the two telescopes indeed trace the same gas, we have calculated the line strengths expected from the JCMT rotation diagrams. For the sources where CH3OH is detected with OVRO, we predict a brightness of $`1530`$ K in the CH3OH $`107`$ GHz line and of $`1015`$ K in the $`104`$ GHz line, in good agreement with the data. For NGC 7538 IRS9, which has a low CH3OH rotation temperature, we expect $`T_B60`$ K in the $`107`$ GHz line. The lack of CH3OH emission in the OVRO data suggests that there is a gradient in $`T_{\mathrm{ex}}`$ from $`30`$ K on a $`15^{\prime \prime }`$ scale (traced by the JCMT but resolved out by OVRO) to $`100`$ K within the $`3^{\prime \prime }`$ OVRO synthesized beam. For the other sources, the bulk of the envelope may already be heated to $`100`$ K and up.

### 3.3 Maser emission

While the $`V_{\mathrm{LSR}}=53.9`$ km s-1 component in the $`107`$ GHz line profile observed toward NGC 7538 IRS1 with OVRO has a brightness and a width consistent with the JCMT data of this source, the brightness of the other velocity components is a factor of $`25`$ higher than expected from the JCMT data. This difference suggests that the observed emission is amplified by a maser process. This possibility is supported by the line profile with three narrow velocity peaks (Figure 4). Maser emission has been observed toward this source in several lower-frequency methanol lines, by Wilson et al. (1984) and Batrla et al. (1987), at velocities consistent with the two strongest components reported here. Using VLBI, Minier et al. (1998) resolved the maser emission in the $`5_16_0`$ A+ ($`6.67`$ GHz) and $`2_03_1`$ E ($`12.2`$ GHz) lines into eleven spots, with sizes of $`<10^3`$ arc seconds each. The spots at velocities $`5762`$ km s-1 lie in an outflow, but the kinematics of brightest masers, those at $`56`$ km s-1, are consistent with an origin in a disk. These two lines as well as the $`107`$ GHz line are transitions to the $`K`$ground state (so-called Class II masers), which require a strong continuum radiation field to be pumped (Cragg et al. (1992)). In the case of NGC 7538 IRS1, this radiation is provided by the H II region, consistent with the fact that among the sources studied here with OVRO, only NGC 7538 IRS1 displays maser emission, which is also the source with the strongest H II region. For the $`1_{10}1_{11}`$ line of H2CO, where Rots et al. (1981) observed maser action at similar velocities, the same pump mechanism operates, as demonstrated by Boland & de Jong (1981).

The brightness ratio of the $`6.67`$ to $`12.2`$ GHz methanol maser lines measured by Minier et al. (1998) is $`10`$. Comparing this number with the models by Sobolev et al. (1997), we find $`N`$(CH3OH)$`10^{16}`$ cm-2. This value is similar to $`N`$(H2CO) from Boland & de Jong (1981), which is a lower limit because Rots et al. (1981) did not resolve the emission spatially. The ratio $`N`$(H2CO)/$`N`$(CH3OH) for the disk of NGC 7538 IRS1 is therefore $`\genfrac{}{}{0pt}{}{_>}{^{}}1`$, which is significantly larger than in the extended envelopes.

## 4 Abundance profiles of CH3OH and H2CO

In this section, we proceed to interpret the CH3OH and H2CO observations in terms of the detailed physical models developed by van der Tak et al. (2000). These models use a power law density structure, $`n=n_0(r/r_0)^\alpha `$. Specific values for the density gradient $`\alpha `$, the size scale $`r_0`$ and the density scale $`n_0`$ for every source can be found in van der Tak et al. (2000). These values, as well as temperature profiles, have been derived from dust continuum and CS, C34S and C17O line observations at submillimeter wavelengths.

### 4.1 Constant-abundance models

Using the density and temperature profiles, we have first attempted to model the data using constant abundances of H2CO and CH3OH throughout the envelopes. This assumption is motivated by the fact that in models of pure gas-phase chemistry, the abundances of H2CO and CH3OH do not change much in the range $`20100`$ K. Molecular excitation and radiative transfer are solved simultaneously with a computer program based on the Monte Carlo method, written by Hogerheijde & van der Tak (2000). The emission in all observed lines is calculated, integrated over velocity, and convolved to the appropriate telescope beam. Comparison to observations proceeds with the $`\chi ^2`$ statistic described in van der Tak et al. (2000).

The collisional rate coefficients for CH3OH have been provided by M. Walmsley (1999, priv. comm.), and are based on the experiments by Lees & Haque (1974). A detailed description of the coefficients is given by Turner (1998). Following Johnston et al. (1992), we have set the $`\mathrm{\Delta }K=3`$ rates to $`10`$% of the $`\mathrm{\Delta }K=1`$ values.

The results of the model calculations are shown in Fig. 5, and the derived best-fit abundances are given in column 6 of Table 2. In the figure, the observed and modeled line fluxes have been converted to column densities following Helmich et al. (1994). Three types of sources can be distinguished. Most sources can be fitted with CH3OH/H2$`10^9`$, similar to the values found for dark and translucent clouds (Turner 1998, Takakuwa et al. 1998), while the “hot core”-type sources W 3(H2O), NGC 6334 IRS1 and W 28 A2 require abundances of a few times $`10^8`$. However, for three sources, no single methanol abundance gives a good fit.

### 4.2 Jump models

For the sources W 33A, GL 2591 and NGC 7538 IRS1, the lines from energy levels $`\genfrac{}{}{0pt}{}{_>}{^{}}100`$ K above ground require significantly higher abundances than the lower-excitation lines. These results suggest that the abundance of CH3OH is higher in the warm, dense gas close to the star than in the more extended, cold and tenuous gas. From the location of the break in the abundance profile, it seems likely that evaporation of grain mantles, which also occurs at $`100`$ K, plays a role. As a simple test of ice evaporation, we have run models for these three sources where the CH3OH abundance follows a step function. In these “jump” models, the CH3OH abundance is at a low level, the “base level”, far from the star, while at a threshold temperature, the abundance surges to a high value, the “top level”. The situation is sketched in Figure 6. The motivation for this model is that if methanol is present in icy grain mantles, its abundance will increase strongly when these ices evaporate. For the temperature threshold, we take $`90`$ K, which is where water ice, the most refractory and most abundant component of the grain mantles, evaporates in $`10`$ yr (Sandford & Allamandola (1993)).

Alternatively, the ice may be desorbed in a weak shock associated with the molecular outflows these sources are known to have: Jones et al. (1996) calculate that a local grain-grain velocity dispersion of $`2`$ km s-1 is sufficient to shatter ice material. Shock desorption of methanol is known to be important in low-mass protostars, for example L 1157 (Bachiller et al. (1995)) and NGC 1333 IRAS 4A (Blake et al. (1995)). The large line widths and line profiles measured here for the low-luminosity sources GL 7009S and IRAS 20126 indicate that shock desorption plays a role there; for the other sources, thermal effects probably dominate. Since CH3OH is confined to the refractory (polar) component of the ice, it is not necessary to consider the effect of a slowly rising temperature in the hot core region, as Viti & Williams (1999) did.

For the base level, we take the abundances found in the previous section, which were constrained mostly by the low-$`K`$ lines with $`E<100`$ K. We have considered jump factors of 3, 10, 30, 100 and 300 without further iteration. Between these models, jumps by factors of $`30`$ give the best match to the high-excitation lines. The results of these models are plotted as squares in Fig. 5 and listed in column 6 of Table 2.

Further constraints on the abundance profile of CH3OH in our sources may be obtained by comparing the models to the OVRO data. The points in Fig. 7 are the OVRO visibility data of W 33A and GL 2591. The observations of the $`107`$ GHz line have been integrated over velocity and binned in annuli around the source position. Superposed are model points for the constant-abundance model and for the jump model. The data do not particularly favour one model or the other, perhaps because this transition traces mostly cold gas. Collisional rate coefficients up to $`J=11`$ are eagerly awaited, so that the $`104`$ GHz line can be modeled as well.

## 5 Discussion

### 5.1 Relation to solid-state observations

The enhancement of the methanol abundance at high temperatures is very likely related to the evaporation of solid methanol, which has been observed at infrared wavelengths toward many of the sources studied here (§ 1). Ground-based data already revealed that CH3OH is not mixed with the bulk of the water ice (Skinner et al. (1992)). Recent ISO results indicate extensive ice segregation toward massive young stars; the absorption profiles of CO2 ice indicate that the methanol resides in a H2O:CH3OH:CO2=1:1:1 layer heated to $`80`$ K as shown by Ehrenfreund et al. (1998), Gerakines et al. (1999) and Dartois et al. (1999a). The column density of solid CH3OH is typically $`5`$% of that of water ice, but in two cases, W 33A and GL 7009S, it is as high as $`30`$%. It has been suggested that CH3OH ice is preferentially formed around more massive young stars, but among the sources studied here, these two do not stand out in luminosity. Instead, they are the two most embedded sources in the sample (see van der Tak et al. 2000), which may enhance the formation of solid methanol. The gas-phase abundance of CH3OH (Table 2) is much lower in GL 7009S than that in W 33A, suggesting that GL 7009S is in an earlier evolutionary state where ice evaporation affects only a small part of the envelope.

A link between solid and gaseous methanol is plausible because the sources with high gas-phase CH3OH abundances are also the ones with high fractions of annealed solid 13CO2 (Boogert et al. (2000)), and high abundances of gas-phase H2O and CO2 (Boonman et al. 2000). These molecules evaporate at $`90`$ K, like methanol. The CH3OH abundances also follow the ratios of envelope mass to stellar mass and the 45/100 $`\mu `$m colours from van der Tak et al. (2000), which confirms the picture that warmer sources have higher molecular abundances in the gas phase. These results indicate that the excitation and abundance of gaseous CH3OH can be used as evolutionary indicators during the embedded stage of massive star formation. As discussed in § 3.1, using submillimeter data to trace evolution has the advantage of being independent of source orientation or total mass, because the dust emission is optically thin at these wavelengths.

Grain mantle evaporation appears to be much less important for H2CO than for CH3OH. We have compared the H2CO data to the constant-abundance models from van der Tak et al. (2000) in similar plots as Figure 5, and found good agreement. There is no evidence for jumps in the H2CO abundance by factors $`\genfrac{}{}{0pt}{}{_>}{^{}}3`$ within the temperature range of $`20250`$ K that the observations are sensitive to. This result suggests that the formaldehyde observed in these sources is predominantly formed in the gas phase by oxidation of CH3. Gas-phase models by Lee et al. (1996) indicate an abundance of $`10^9`$, similar to the observed value, although a contribution from ice evaporation at the $`10^9`$ level cannot be excluded.

The high observed HDCO/H2CO ratios and the detections of D2CO toward the Compact Ridge in Orion (Wright et al. 1996; Turner 1990) and toward embedded low-mass objects such as IRAS 16293 (van Dishoeck et al. 1995; Ceccarelli et al. 1998) indicate H2CO formation on grains. In our sources, HDCO is not detected to H2CO/HDCO $`>10`$; in a survey of hot core-type sources, Hatchell et al. (1998) obtained DCN/HCN $`10^3`$, much lower than in embedded low-mass stars. Although the surface chemistry should qualitatively be the same, the cold ($`10`$ K) phase may last too short in our sources to build up large amounts of deuterated molecules such as those seen in the Compact Ridge and in low-mass objects.

In eleven regions of (mostly massive) star formation, Mangum & Wootten (1993) derived ratios of ortho- to para-H2CO of $`1.53`$ and took this result as evidence for grain surface formation of H2CO. However, this conclusion appears tentative, since for most of their sources, only one line of ortho-H2CO was observed. In addition, at the H2CO column densities of $`10^{14}`$ cm-2, representative of the regions studied by Mangum & Wootten and by us, the $`K=0`$ and $`K=1`$ lines of H2CO have optical depths of $`1`$, requiring careful modeling. The data and models presented in this paper are consistent with an ortho-/para-H2CO ratio of 3.

High abundances of methanol as observed in the ices, $`10^6`$ relative to hydrogen, cannot be produced in the gas phase, except maybe in shocks. This mechanism has been proposed for H2O ice by Bergin et al. (1999) and for CO2 by Charnley & Kaufman (2000); its application to CH3OH depends on the existence of a high-temperature route to form methanol, which is not yet known. Hartquist et al. (1995) proposed the reaction CH$`{}_{4}{}^{}+`$OH, but ISO-SWS observations by Boogert et al. (1998) indicate low abundances of gaseous methane in our sources, and at $`T\genfrac{}{}{0pt}{}{_>}{^{}}200`$ K, all OH should be consumed by H and H2 to form H2O. The widths of the CH3OH lines of only a few km s-1 also argue against formation in shocks.

More likely, the solid methanol is formed by reactions on or inside the ice layers. Addition of H atoms to CO molecules will lead to H2CO, and further to CH3OH. In the literature, there are three proposed sources of atoms: direct accretion from the gas phase, ultraviolet irradiation and bombardment by energetic particles. We show in Section 5.4 that the latter two mechanisms are unlikely to be important for the sources studied in this paper, and focus for now on the first.

### 5.2 Surface chemistry model

If the CH3OH ice that we see evaporating in these sources originates from H atom addition (=reduction) of CO ice on the surfaces of dust grains, it cannot have been produced under the current physical conditions. The evaporation temperature of CH3OH and H2O ice, $`90`$ K, is much higher than that of CO and O, $`20`$ K, and that of H, $`1015`$ K, depending on its surface mobility and reactivity. Any formation of CH3OH ice through surface reduction of CO must therefore have occurred before the central star heated up its envelope above $`15`$ K. The most likely phase of the cloud to form methanol ice through surface chemistry is therefore the pre-protostellar phase, when the grain temperature may have been as low as $`10`$ K and the cloud was contracting to form a dense core. However, about 1/3 of solid CO is observed inside the water ice layer, and will not evaporate until much higher temperatures are reached (Tielens et al. (1991)). If H and O atoms can be similarly trapped, solid-state reduction and oxidation may occur at temperatures well above $`20`$ K.

Tielens & Hagen (1982) and Tielens & Allamandola (1987) proposed that direct accretion of atoms and molecules from the gas phase, followed by low-temperature surface reactions, determine the composition of grain mantles. Activation barriers are offset by the long effective duration of the collision on the surface. The resulting mantle composition is determined by the relative accretion rates of H and CO onto the grains, and the relative height of the reaction barriers.

This mechanism predicts that the abundance ratios of H2CO/CH3OH and CO/CH3OH are correlated, as shown by the curve in Fig. 8, taken from Charnley et al. (1997), with the reaction probability ratio $`\varphi _\mathrm{H}`$ of CO + H to H2CO + H is $`10^3`$. Along the curve, the H/CO abundance ratio decreases down to $`10^2`$ at its tip. At the density in the central regions of our sources, $`10^6`$ cm-3, atomic H is mainly produced by cosmic-ray ionization of H2, giving a constant concentration of $`1`$ cm-3, or an abundance of $`10^6`$. The abundance of CO in these sources was measured by van der Tak et al. (2000) to be $`1\times 10^4`$, as expected if all gas-phase carbon is in CO and no significant depletion occurs. Hence, H/CO $`10^2`$, and our sources should lie right at the tip of the curve. The value $`\varphi _\mathrm{H}=10^3`$ implies that the CO$``$ H2CO reaction limits the rate of CH3OH formation, consistent with the non-detection of evaporated H2CO.

The figure also plots the observed abundance ratios, where the H2CO and CO abundances have been taken from van der Tak et al. (2000), while for CH3OH, either the constant-abundance model or the “top level” of the jump model has been plotted. The data points are seen to lie along a “sequence” which lies along the model curve, perhaps displaced toward lower CH3OH abundances by an order of magnitude. Although this behaviour is consistent with solid-state formation of H2CO and CH3OH at various H/CO ratios, the different rotation temperatures and abundance profiles of H2CO and CH3OH argue against a common chemical origin for these two species, and favour complete hydrogenation of CO into CH3OH. The figure also plots solid state observations by Keane et al. (2000), which are lower limits in the case of CO because of evaporation.

Since the abundances of CO and H2CO are almost the same in all our sources, the spread in the abundance ratios is probably due to different degrees of CH3OH evaporation and subsequent gas-phase destruction. Indeed, the abundances of gaseous CH3OH are factors of $`\genfrac{}{}{0pt}{}{_>}{^{}}10`$ lower than the values measured in the solid state. This difference is expected since methanol is destroyed in $`10^4`$ yr (Charnley et al. 1992) in reactions with ions such H$`{}_{}{}^{+}{}_{3}{}^{}`$, leading to a rich chemistry with species like CH3OCH3 and CH3OCHO. These species are indeed detected toward the “hot cores” W 3 (H2O) and NGC 6334 IRS1, weaker toward the “jump sources” GL 2591 and NGC 7538 IRS1, and absent in the low-methanol source W 3 IRS5 (Helmich & van Dishoeck 1997; van Dishoeck & van der Tak 2000).

We conclude that the observed spread in the gaseous methanol abundance is due to incomplete evaporation for the cold sources, which are the least evolved ones in our sample, and to gas-phase reactions for the hot cores, which are the most evolved. The most likely source of the methanol is grain surface chemistry in the pre-protostellar phase. The conditions in this phase cannot be derived by comparison to CO and H2CO, for which our data suggest that gas-phase processes control the chemistry. A more promising molecule is CO2.

### 5.3 The H/O competition: Density and duration of the pre-stellar phase

The ISO detection of ubiquitous solid CO2 by Gerakines et al. (1999) makes clear that oxidation of CO is a potentially important competitor to reduction. These data are shown in the bottom panel of Fig. 8, as well as gas-phase CO2 abundances from Boonman et al. (2000). Gas-grain chemistry is clearly important for CO2, as it is for CH3OH, warranting a comparative study of the two. The model by Charnley et al. (1997) did not include reactions with O, while Tielens & Hagen (1982) used O2 as the dominant form of oxygen, which is inconsistent with recent observational limits (Melnick (2000)). We have constructed a model of gas-grain chemistry based on the modified rate equation approach described in Caselli et al. (1998), and extensively tested against a Monte Carlo program. The chemical system consists of three gaseous species: H, O and CO, and eight surface species: OH, H2O, O2, CO2, HCO, H2CO, CH2OH and CH3OH. The calculations assume fast quantum tunneling of hydrogen, as proposed by Tielens & Hagen (1982). If hydrogen does not tunnel, as measurements by Katz et al. (1999) suggest, our simple system is unable to form significant mole fractions of CH3OH and CO2. In the case of reduced mobility for all species, the Eley-Rideal mechanism should be considered to model surface chemistry (Herbst (2000)), which is not included here.

The mantle composition (Table 4) depends on density through the composition of the gas phase, in particular the H/O ratio. We take the concentrations of H, O and CO from Lee et al. (1996), for the case of low metals, $`T=10`$ K, in steady state. The main result is that after reduction to HCO, CO is mostly reduced into CH3OH at low densities and mostly oxidized at higher densities. In our model, the reaction CO+O$``$CO2 has a barrier of 1000 K (see d’Hendecourt et al. (1985)), else ices of 60% CO2 would form which are not observed. Water is copiously made through the sequence O+H$``$OH and OH+H$``$H2O. The detailed composition of the grain mantles is quite uncertain because of unknown reaction rates. For instance, the fraction of H2CO depends on relative barrier height of the CO+H and H2CO+H reactions, but H2CO is never a major component of the ice layer, consistent with gas-phase and solid state observations. The results in Table 4 refer to a chemical time scale of 1000 years. Unlike inert species such as O2, H2O, CH3OH and CO2, the fractions of trace species (fractions $`\genfrac{}{}{0pt}{}{_<}{^{}}0.01`$) also depend somewhat on time, and are not used in the analysis. Species that react quickly with H maintain a constant population because they are hydrogenated on the H accretion time scale ($`1`$ day if $`n_\mathrm{H}=10^3`$ cm-3), and their molar fractions decline as the ice layer grows with time.

As an example, we will now estimate the yield of CH3OH by this scheme for the case $`n_\mathrm{H}=10^4`$ cm-3. The depletion rate of CO molecules from the gas phase is $`f_d=n_d\sigma V_{\mathrm{CO}}S`$, with $`n_d`$ and $`\sigma `$ the number density and cross section of the dust grains, $`V_{CO}`$ the thermal velocity of CO and $`S`$ the sticking coefficient, assumed to be $`10`$%. Taking silicate grains of radius $`0.1`$ $`\mu `$m, $`f_d=3\times 10^{15}`$ s-1, or $`10`$% in $`10^6`$ years. At this density, $`77`$% of this depleted CO goes into CH3OH, so that upon evaporation, the abundance ratio would be CH3OH/CO $`0.086`$ or CH3OH/H2$`10^5`$.

Assuming that the composition of ice layers is determined by surface chemistry, our model can be used to investigate the initial conditions of massive star formation by considering abundance ratios of CO, CO2 and CH3OH. The CO2/CH3OH ratio is sensitive to density, independent of time, and equals the ratio of the molar fractions of solid CO2 and CH3OH in Table 4. Time is measured by the ratio of raw material (CO) to product (CH3OH) through the fraction of CO in the solid state derived above, $`0.01(t/10^5\mathrm{yr})(n_\mathrm{H}/10^4\mathrm{cm}^3)`$. The abundance ratio in the gas phase is this fraction multiplied by the efficiency of CO$``$CH3OH conversion, equal to the molar fraction of CH3OH divided by the sum of the fractions of CO and its possible products, CO2, H2CO and CH3OH. Figure 8 plots the synthetic abundance ratios as open circles, labeled by time and density. If this model is valid, the observed CO2/CH3OH ratio constrains the density in the pre-protostellar phase to be $`10^5`$ cm-3, or higher, since some CO2 may be destroyed shortly after evaporation. At this density, hydrogenation of CO is incomplete due to the low H flux, and a significant abundance of solid H2CO is expected, consistent with the results of Keane et al. (2000). Using this limit on the density, the observed CO/CH3OH ratio constrains the time spent in the pre-protostellar phase to be $`\genfrac{}{}{0pt}{}{_<}{^{}}10^5`$ years within the grain surface chemistry scenario. The same conclusion is reached when the CO/CO2 ratio is used instead of CO/CH3OH; the uncertainty is a factor of ten due to the unknown sticking coefficient. This time scale is significantly smaller than the corresponding number for low-mass stars, where this phase lasts $`10^6`$ years as derived from the ratio of dense cores with and without stars found by Beichman et al. (1986). The same time scale of $`10^6`$ yr is obtained from our models for an abundance ratio of CO/CH3OH $`10^3`$, as observed in the compact ridge in Orion.

These models also help to understand observations of solid CO2 in other regions where no star formation is occurring. Toward the field star Elias 16 behind the Taurus dark cloud, CO2 has been detected (Whittet et al. (1998)) but CH3OH has not (Chiar et al. (1998)). The models suggest thus that the density in this region is $`\genfrac{}{}{0pt}{}{_>}{^{}}3\times 10^4`$ cm-3. The model does not explain the lack of solid CH3OH toward SgrA ($`\genfrac{}{}{0pt}{}{_<}{^{}}3`$% relative to H2O ice; Chiar et al. (2000)), where CO2 ice has been detected. For this low-density line of sight, the opposite ratio would be expected. An enhanced temperature or ultraviolet irradiation may be important here. For the disk around NGC 7538 IRS1, where $`N`$(H2CO)/$`N`$(CH3OH) $`\genfrac{}{}{0pt}{}{_>}{^{}}1`$ (§ 3.3), the models indicate that the density is $`\genfrac{}{}{0pt}{}{_>}{^{}}10^5`$ cm-3.

The model is also consistent with current observational limits on solid O2. At low densities, our models drive all O into H2O on the dust grains, but for $`10^5`$ cm-3, part of it goes into O2. However, the ratio O2/CO remains $`<1`$ on the surface, in agreement with the limit derived by Vandenbussche et al. (1999) from ISO-SWS observations of NGC 7538 IRS9. Evaporation of solid O2 plays a role, but observations of gaseous O2 with the Submillimeter Wave Astronomy Satellite give an upper limit toward massive star-forming regions of $`10^7`$ relative to H2 (Melnick (2000)), consistent with the solid-state results.

### 5.4 Alternative models

Could ultraviolet irradiation or energetic particle bombardment also produce the observed trends in the abundances? The production rate of species $`i`$ by irradiation can be written as $`dn_i/dt=\alpha _i\mathrm{\Phi }4\pi r_g^2`$, with $`\mathrm{\Phi }`$ the ultraviolet flux in photons cm-2 s-1 and $`r_g`$ the radius of the dust grains, taken to be $`10^5`$ cm. Inside dense clouds, most ultraviolet radiation is produced by cosmic-ray interaction with H2, which gives a field approximately equal to the interstellar radiation field at $`A_V=5`$ (Prasad & Tarafdar (1983)), or $`5000`$ photons cm-2 s-1. The reaction yields $`\alpha _i`$ follow from laboratory experiments such as those reported at the Leiden Observatory Laboratory database333see www.strw.leidenuniv.nl/$``$lab, in this case on a mixture of initial composition H2O:CO = 100:33. The amounts of CH3OH and CO2 produced do not depend on this ratio as long as it is $`>1`$. Using band strengths by Gerakines et al. (1996) and Kerkhof et al. (1999), we obtain $`\alpha =1.9\times 10^3`$ for CH3OH and $`\alpha =3.2\times 10^2`$ for CO2. The ratio of these numbers is in good agreement with the observed abundance ratio of $`10`$. However, the absolute values of the $`\alpha _i`$ imply production rates of $`10^8`$ s-1 per grain, compared to the accretion rate of H, O and CO of $`10^5`$ s-1. Stellar ultraviolet radiation can only affect the inner parts of the envelopes since their extinctions are $`A_V100`$ magnitudes. The absence of a strong ultraviolet field throughout the envelopes of our sources is also suggested by the observational limits from ISO-SWS on mid-infrared fine structure lines and of emission by polycyclic aromatic hydrocarbons (van Dishoeck & van der Tak  2000).

The processing of interstellar ice by cosmic rays was studied by Hudson & Moore (1999), who bombarded an H2O:CO=5:1 ice mixture with protons of energy $`0.8`$ MeV. These experiments produced abundance ratios of H2CO/CH3OH=0.6 and CO2/CH3OH=2.0, in fair agreement with the observed values. Other species formed as well, notably formic acid (HCOOH) and methane (CH4), which indeed are observed in interstellar grain mantles. However, the experiments by Hudson & Moore produced almost twice as much HCOOH and CH4 as CO2, while the observed abundances are only $`10`$% of that of CO2 (Schutte et al. (1999); Boogert et al. 1998). In addition, the particle dose in the experiments was $`2\times 10^{15}`$ cm-2, while the interstellar cosmic-ray flux is only 3 cm-2 s-1. Hence, the experiments simulate a bombardment for $`3\times 10^7`$ yr, a factor of $`1000`$ longer than the ages of the sources studied here. Stellar X-ray emission (Glassgold et al. (2000)) only acts on small scales, especially because of heavy attenuation in the envelopes of these obscured objects.

## 6 Conclusions

The chemistry of CH3OH and H2CO in thirteen regions of massive star formation is studied through single-dish (JCMT) and interferometer (OVRO) line observations at submillimeter wavelengths. Our main conclusions are:

1. The submillimeter emission lines of CH3OH toward most sources have widths of $`35`$ km s-1, consistent with those found earlier for C17O and C34S. However, in the low-luminosity sources GL 7009S and IRAS 20126, the line shapes reveal that CH3OH is present in the outflow. These results indicate that the desorption of ices in the envelopes of low-mass protostars is primarily by shocks, while thermal processes dominate in the case of massive stars.

2. Rotational temperatures of CH3OH range from $`10`$ to $`200`$ K and correlate very well with the excitation temperature of C2H2 measured in infrared absorption. This correlation suggests that both species trace the same gas, which may be outgassing of icy grain mantles. For H2CO, the range in $`T_{\mathrm{rot}}`$ is only $`6090`$ K without relation to $`T_{\mathrm{ex}}`$(C2H2), suggesting a different chemical origin. We propose that $`T_{\mathrm{rot}}`$ (CH3OH) can be used as evolutionary indicator during the embedded phase of massive star formation, independent of source optical depth or orientation.

3. Detailed non-LTE radiative transfer models of the CH3OH lines suggest a distinction of three types of sources: those with CH3OH/H2$`10^9`$, those with CH3OH/H2$`10^7`$ and those which require a jump in its abundance from $`10^9`$ to $`10^7`$. The models are consistent with the $`3^{\prime \prime }`$ size of the CH3OH $`107`$ GHz emission measured interferometrically. The location of the jump at $`T100`$ K strongly suggests that the methanol enhancement is due to evaporation of icy grain mantles. The sequence of low-methanol $``$ jump $``$ high-methanol sources corresponds to a progression in the ratio of envelope mass to stellar mass and the mean temperature of the envelope. In contrast, the observed H2CO seems primarily produced in the gas phase, since the H2CO data can be well fit with a constant abundance of a few $`\times 10^9`$ throughout the envelope. The grain surface hydrogenation of CO thus appears to be completed into CH3OH, with little H2CO left over.

4. Model calculations of gas-grain chemistry show that CO is primarily reduced (into CH3OH) at densities $`n_\mathrm{H}\genfrac{}{}{0pt}{}{_<}{^{}}10^4`$ cm-3, and primarily oxydized (into CO2) at higher densities. To keep sufficient CO on the grains, this mechanism requires temperatures of $`\genfrac{}{}{0pt}{}{_<}{^{}}15`$ K, i.e., conditions before star formation. Assuming that surface reactions proceed at the accretion rate of CO, the observed CO2 and CH3OH abundances constrain the density in the pre-protostellar phase to be $`n_\mathrm{H}\genfrac{}{}{0pt}{}{_>}{^{}}`$ a few $`10^4`$ cm-3, and the time spent in this phase to be $`\genfrac{}{}{0pt}{}{_<}{^{}}10^5`$ yr. Our surface chemistry model predicts that lines of sight through clouds with a high H/O ratio will show abundant solid methanol and less CO2. Ultraviolet photolysis and radiolysis by energetic (MeV) protons appear less efficient as ice processing mechanisms for these sources; radiolysis also overproduces HCOOH and CH4.

###### Acknowledgements.

The authors wish to thank Malcolm Walmsley for providing collisional rate coefficients for methanol, Xander Tielens, Pascale Ehrenfreund and Willem Schutte for comments on the manuscript and Ted Bergin, Eric Herbst, Friedrich Wyrowski and Wilfried Boland for useful discussions. Annemieke Boonman and Jacquie Keane kindly provided us with their results in advance of publication. We are grateful to Remo Tilanus and Fred Baas at the JCMT and Geoffrey Blake at OVRO for assistance with the observations. This research was supported by NWO grant 614-41-003 and the MURST program “Dust and Molecules in Astrophysical Environments”.