format: text

text:
# Kaon-baryon coupling constants in the QCD sum rule approach

## I Introduction

To understand kaon-nuclear physics, it is important to know the hadronic coupling strengths involving the kaons. Among them, $`g_{KN\mathrm{\Lambda }}`$ and $`g_{KN\mathrm{\Sigma }}`$ are the most relevant coupling constants. In contrast to $`g_{\pi NN}`$, however, the determination of these kaon couplings has some difficulties both in the experimental side and in the theoretical side, e.g. see .

Among other theoretical approaches, QCD sum rule method  has been used to extract these kaon couplings. However, compared to the large number of works devoted to $`g_{\pi NN}`$, there have been only few QCD sum rule estimates on $`g_{KN\mathrm{\Lambda }}`$ and $`g_{KN\mathrm{\Sigma }}`$, for which there are still ambiguities in among the calculations. Thus the results are quite different from each other. More detailed analyses are needed both experimentally and theoretically to understand this discrepancy, and to understand kaon-nuclear physics.

In Ref. , the OPE was calculated only up to the leading term coming from the quark condensate and to leading order in $`m_s`$ in the sum rule structure proportional to $`\text{/}qi\gamma _5`$. However, the next leading term, dimension 5 $`\overline{q}g_s\sigma Gq`$ may contribute to the OPE side with considerable amount as in nucleon mass sum rule . In addition, operators of dimension 7 may also be important in the OPE side as a further power correction. Thus, in this paper we re-analyze our QCD sum rule calculation including higher dimensional condensates, such as $`\overline{q}g_s\sigma Gq`$ and $`\overline{q}q\frac{\alpha _s}{\pi }G^2`$, and study the contribution of these condensates on the previous results.

In Sec. II we present our sum rules for $`g_{KN\mathrm{\Lambda }}`$ and $`g_{KN\mathrm{\Sigma }}`$, and Sec. III we discuss some uncertainties in our sum rules and summarize our results.

## II QCD sum rules for $`g_{KN\mathrm{\Lambda }}`$ and $`g_{KN\mathrm{\Sigma }}`$

We will closely follow the procedures given in Refs.. Consider the three point function constructed of the two baryon currents $`\eta _B`$, $`\eta _B^{}`$ and the pseudoscalar meson current $`j_5`$.

$`A(p,p^{},q)={\displaystyle 𝑑x𝑑y0|T(\eta _B^{}(x)j_5(y)\overline{\eta }_B(0))|0e^{i(p^{}xqy)}}.`$ (1)
In order to obtain $`g_{KN\mathrm{\Lambda }}`$, we will use the following currents for the nucleon and the $`\mathrm{\Lambda }`$ .

$`\eta _N`$ $`=`$ $`ϵ_{abc}(u_a^TC\gamma _\mu u_b)\gamma _5\gamma ^\mu d_c,`$ (2)
$`\eta _\mathrm{\Lambda }`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{3}}}ϵ_{abc}\left[(u_a^TC\gamma _\mu s_b)\gamma _5\gamma ^\mu d_c(d_a^TC\gamma _\mu s_b)\gamma _5\gamma ^\mu u_c\right],`$ (3)
where u and d are the up and down quark fields ($`a,b`$ and $`c`$ are color indices), $`T`$ denotes the transpose in Dirac space, and $`C`$ is the charge conjugation matrix. For the $`K^{}`$ we choose the current

$`j_K^{}=\overline{s}i\gamma _5u.`$ (4)
The general expression for $`A(p,p^{},q)`$ has the following form

$`A(p,p^{},q)`$ $`=`$ $`F_1(p^2,p^2,q^2)i\gamma _5+F_2(p^2,p^2,q^2)\text{/}qi\gamma _5`$ (5)
$`+`$ $`F_3(p^2,p^2,q^2)\text{/}Pi\gamma _5+F_4(p^2,p^2,q^2)\sigma ^{\mu \nu }\gamma _5q_\mu p_\nu ^{},`$ (6)
where $`q=p^{}p`$ and $`P=\frac{p+p^{}}{2}`$. Recently, in Ref.  it was reported that in the case of $`g_{\pi NN}`$ the $`\sigma ^{\mu \nu }\gamma _5`$ structure is independent of the effective models employed in the phenomenological side and further provides the $`\pi NN`$ coupling with less uncertainties from QCD parameters. Motivated by this result $`g_{KN\mathrm{\Lambda }}`$ and $`g_{KN\mathrm{\Sigma }}`$ were calculated from this structure in Refs. . In this paper, however, we construct the sum rule for only the $`\text{/}qi\gamma _5`$ structure as before, and compare this with our previous one.

On the phenomenological side, keeping the first two terms we have

$`\lambda _N\lambda _\mathrm{\Lambda }{\displaystyle \frac{M_B}{(p^2M_N^2)(p^2M_\mathrm{\Lambda }^2)}}(\text{/}qi\gamma _5)g_{KN\mathrm{\Lambda }}{\displaystyle \frac{1}{q^2m_K^2}}{\displaystyle \frac{f_Km_K^2}{2m_q}}`$ (7)
$`+`$ $`\lambda _N\lambda _\mathrm{\Lambda }^{}{\displaystyle \frac{M_B^{}}{(p^2M_N^2)(p^2M_\mathrm{\Lambda }^{}^2)}}(\text{/}qi\gamma _5)g_{KN\mathrm{\Lambda }^{}}{\displaystyle \frac{1}{q^2m_K^2}}{\displaystyle \frac{f_Km_K^2}{2m_q}}`$ (8)
$`+`$ $`\mathrm{higher}\mathrm{resonances},`$ (9)
where $`M_B=\frac{1}{2}(M_N+M_\mathrm{\Lambda })`$, and $`M_B^{}=\frac{1}{2}(M_NM_\mathrm{\Lambda }^{})`$. Here $`\mathrm{\Lambda }^{}`$ means the $`\mathrm{\Lambda }`$(1405), and we introduce (–) sign for the $`\mathrm{\Lambda }`$ (1405) mass because it is a negative parity state. However, this is not relevant in the following calculation. $`\lambda _N`$, $`\lambda _\mathrm{\Lambda }`$ and $`\lambda _\mathrm{\Lambda }^{}`$ are the coupling strengths of the baryons to their currents. $`m_q`$ is the average of the quark masses, $`f_K`$ the kaon decay constant and $`m_K`$ the kaon mass. We take $`f_K`$ = 0.160 GeV and $`m_s`$ = 0.150 GeV.

As for the OPE side, the new contribution from the quark-gluon condensates is given by

$`\sqrt{{\displaystyle \frac{2}{3}}}{\displaystyle \frac{7}{2^43\pi ^2}}\mathrm{ln}(p^2)(\overline{q}g_s\sigma Gq+\overline{s}g_s\sigma Gs),`$ (10)
and from dimension 7 ops.

$`+\sqrt{{\displaystyle \frac{2}{3}}}{\displaystyle \frac{5}{2^33^2}}{\displaystyle \frac{1}{p^2}}(\overline{q}q+\overline{s}s){\displaystyle \frac{\alpha _s}{\pi }}G^2,`$ (11)
where we take the limit $`p^2p^2`$ and let $`\overline{u}u=\overline{d}d\overline{q}q`$, $`\overline{u}g_s\sigma Gu=\overline{d}g_s\sigma Gd\overline{q}g_s\sigma Gq`$. Here we collect only the terms which contribute to the $`\text{/}q/q^2`$ structure such as Figs. 1 and 2. Using the standard values for $`\overline{s}g_s\sigma Gs=0.8\overline{q}g_s\sigma Gq`$ and $`\overline{q}g_s\sigma Gq=m_0^2\overline{q}q=0.8\overline{q}q`$  the sum rule after Borel transformation to $`p^2=p^2`$ becomes

$`\lambda _N\lambda _\mathrm{\Lambda }{\displaystyle \frac{M_B}{M_\mathrm{\Lambda }^2M_N^2}}\left(e^{M_N^2/M^2}e^{M_\mathrm{\Lambda }^2/M^2}\right)g_{KN\mathrm{\Lambda }}{\displaystyle \frac{f_Km_K^2}{2m_q}}+A\left(e^{M_N^2/M^2}e^{M_\mathrm{\Lambda }^{}^2/M^2}\right)=`$ (12)
$`\sqrt{{\displaystyle \frac{2}{3}}}\left({\displaystyle \frac{33}{40\pi ^2}}E_1M^4+({\displaystyle \frac{11m_s^2}{60\pi ^2}}{\displaystyle \frac{21}{100\pi ^2}})E_0M^2+({\displaystyle \frac{m_s}{3}}\overline{s}s+{\displaystyle \frac{1}{8}}{\displaystyle \frac{\alpha _s}{\pi }}G^2)\right)\overline{q}q.`$ (13)
Here, A is the unknown constant coming from $`\lambda _\mathrm{\Lambda }^{}g_{KN\mathrm{\Lambda }^{}}`$, and

$`E_i=1{\displaystyle \underset{k=0}{\overset{i}{}}}{\displaystyle \frac{s_0^k}{k!(M^2)^k}}e^{\frac{s_0}{M^2}},`$ (14)
where $`s_0`$ is a continuum threshold. One should be cautious, however, that there may be non-accounted terms, which can not be neglected by using this simple Borel transformation.

For $`\lambda _N`$ and $`\lambda _\mathrm{\Lambda }`$, we use the values obtained from the following baryon sum rules for the $`N`$ and $`\mathrm{\Lambda }`$ :

$`E_2^NM^6+bE_0^NM^2+{\displaystyle \frac{4}{3}}a^2=2(2\pi )^4\lambda _N^2e^{M_N^2/M^2},`$ (15)
$`E_2^\mathrm{\Lambda }M^6+{\displaystyle \frac{2}{3}}am_s(13\gamma )E_0^\mathrm{\Lambda }M^2+bE_0^\mathrm{\Lambda }M^2+{\displaystyle \frac{4}{9}}a^2(3+4\gamma )=2(2\pi )^4\lambda _\mathrm{\Lambda }^2e^{M_\mathrm{\Lambda }^2/M^2},`$ (16)
where $`a(2\pi )^2\overline{q}q`$, $`b\pi ^2\frac{\alpha _s}{\pi }G^2`$, and $`\gamma \overline{s}s/\overline{q}q10.2`$. We use different thresholds for $`\lambda _N`$ and $`\lambda _\mathrm{\Lambda }`$ in Eqs. (15) and (16). We take $`s_N`$ = (1.440 GeV)2 for the nucleon sum rule and $`s_\mathrm{\Lambda }`$ = (1.405 GeV)2 for the $`\mathrm{\Lambda }`$ sum rule considering the next excited nucleon and $`\mathrm{\Lambda }`$ state, respectively.

$`g_{KN\mathrm{\Lambda }}`$ , however, does not display a plateau as a function of the Borel mass. This is because there is no usual power correction term like ($`\frac{a}{M^2},\frac{b}{M^4}`$, and so on) in the r.h.s. of Eq. (13) even including up to dimension 7 operators. We need more higher dimensional operators to get those terms. Thus, in this case we prefer to use a best-fit method. Eq. (13) has the following form :

$`g_{KN\mathrm{\Lambda }}f_1(M^2)+Af_2(M^2)=f_3(M^2).`$ (17)
Then, we get $`g_{KN\mathrm{\Lambda }}`$ and the unknown constant A by minimizing $`(g_{KN\mathrm{\Lambda }}f_1+Af_2f_3)^2`$ with a fixed $`s_0`$ and an appropriate Borel interval:

$`{\displaystyle _{M_{min}^2}^{M_{max}^2}}(g_{KN\mathrm{\Lambda }}f_1+Af_2f_3)^2𝑑M^2=\mathrm{minimum}.`$ (18)
We fix the continuum threshold $`s_0`$ = 2.074 GeV2 taking into account the next term from the N(1440), i.e., $`N(1440)\mathrm{\Lambda }`$, in the phenomenological side.

The Borel interval $`M^2`$ is restricted by the following conditions : OPE convergence and pole dominance. The lower limit of $`M^2`$, $`M_{min}^2`$ is determined as the value at which the contribution of the highest dimensional operators is less than 10$`\%`$ of total OPE. The upper limit $`M_{max}^2`$ is determined by restricting the continuum contribution to be less than 50$`\%`$. Then, we get

$`|g_{KN\mathrm{\Lambda }}|`$ $`=`$ $`2.49,`$ (19)
$`|A|`$ $`=`$ $`0.00174\mathrm{GeV}^7,`$ (20)
and the Borel interval (0.478, 1.068) GeV2 for basic inputs (i.e., $`\overline{q}q`$ = – (0.230 GeV)3, $`\frac{\alpha _s}{\pi }G^2`$ = 0.012 GeV4, $`m_s`$ = 0.150 GeV, and $`m_0^2`$ = 0.8 GeV2). Here we denote the absolute value because we can not determine signs of the coupling strengths ($`\lambda _N`$, $`\lambda _\mathrm{\Lambda }`$ and $`\lambda _\mathrm{\Lambda }^{}`$) in the baryon sum rules. We also calculate the average deviation $`\overline{\delta }`$ $``$ $`_i^N\left|1RHS(M_i^2)/LHS(M_i^2)\right|/N`$ = 8.8 $`\times `$ 10-2 to test the reliability of our fitting, and it shows that the deviation is less than 10 %.

Table I shows variations of $`g_{KN\mathrm{\Lambda }}`$ for other input parameters, which are coming from the uncertainty of the basic inputs. For example, the first line in Table I shows that $`|g_{KN\mathrm{\Lambda }}|`$ = 3.11 (or 1.93) if we change the quark condensate to $`\overline{q}q`$ = –(0.210 GeV)3 (or –(0.250 GeV)3) while other basic inputs are fixed. In the last line we take $`m_0^2`$ = 0.6 GeV2 from the lowest value of the standard QCD sum rule estimate , and 1.4 GeV2 which was evaluated in the instanton vacuum in Ref. . Total variation is about $`\pm `$ 1.25 on the above $`g_{KN\mathrm{\Lambda }}`$ value. On the other hand, the unknown constant $`|A|`$ varies from 0.00120 to 0.00203 GeV7.

Next, consider $`g_{KN\mathrm{\Sigma }}`$. The current of $`\mathrm{\Sigma }^{}`$ is obtained by making an SU(3) rotation from the nucleon current

$`\eta _\mathrm{\Sigma }=\sqrt{2}ϵ_{abc}\left[(u_a^TC\gamma _\mu s_b)\gamma _5\gamma ^\mu d_c+(d_a^TC\gamma _\mu s_b)\gamma _5\gamma ^\mu u_c\right].`$ (21)
In this case the contribution of the quark-gluon condensate is given by

$`\sqrt{2}{\displaystyle \frac{1}{2^43\pi ^2}}\mathrm{ln}(p^2)(\overline{q}g_s\sigma Gq+\overline{s}g_s\sigma Gs),`$ (22)
and from dimension 7 ops.

$`+\sqrt{2}{\displaystyle \frac{1}{2^33^2}}{\displaystyle \frac{1}{p^2}}(\overline{q}q+\overline{s}s){\displaystyle \frac{\alpha _s}{\pi }}G^2.`$ (23)
Then, within the same approximation as before we get the following sum rule.

$`\lambda _N\lambda _\mathrm{\Sigma }{\displaystyle \frac{M_B}{M_\mathrm{\Sigma }^2M_N^2}}\left(e^{M_N^2/M^2}e^{M_\mathrm{\Sigma }^2/M^2}\right)g_{KN\mathrm{\Sigma }}{\displaystyle \frac{f_Km_K^2}{2m_q}}+B\left(e^{M_N^{}^2/M^2}e^{M_\mathrm{\Sigma }^2/M^2}\right)=`$ (24)
$`+\sqrt{2}\left({\displaystyle \frac{3}{40\pi ^2}}E_1M^4+({\displaystyle \frac{m_s^2}{60\pi ^2}}+{\displaystyle \frac{3}{100\pi ^2}})E_0M^2{\displaystyle \frac{1}{40}}{\displaystyle \frac{\alpha _s}{\pi }}G^2\right)\overline{q}q,`$ (25)
where $`M_B=\frac{1}{2}(M_N+M_\mathrm{\Sigma })`$ and $`N^{}`$ is N(1440). B is the unknown constant coming from $`\lambda _N^{}g_{KN^{}\mathrm{\Sigma }}`$. Again for $`\lambda _\mathrm{\Sigma }`$, we take the value from the following sum rule for the $`\mathrm{\Sigma }`$:

$`E_2^\mathrm{\Sigma }M^62am_s(1+\gamma )E_0^\mathrm{\Sigma }M^2+bE_0^\mathrm{\Sigma }M^2+{\displaystyle \frac{4}{3}}a^2=2(2\pi )^4\lambda _\mathrm{\Sigma }^2e^{M_\mathrm{\Sigma }^2/M^2}.`$ (26)
We fix the continuum threshold $`s_\mathrm{\Sigma }`$ = (1.660 GeV)2 considering the next $`\mathrm{\Sigma }`$ state, $`\mathrm{\Sigma }`$ (1660).

Using the continuum threshold $`s_0`$ = 2.356 GeV2 taking into account the next term from the N(1535), i.e., $`N(1535)\mathrm{\Sigma }`$, in the phenomenological side we get

$`|g_{KN\mathrm{\Sigma }}|`$ $`=`$ $`0.395,`$ (27)
$`|B|`$ $`=`$ $`0.00148\mathrm{GeV}^7`$ (28)
for the same basic inputs. The Borel interval is (0.488, 1.584) GeV2 and the average deviation of the fit $`\overline{\delta }`$ is 9.7 % in this case. We present the variation of $`g_{KN\mathrm{\Sigma }}`$ on other parameters in Table II. The total variation is about $`\pm `$ 0.377. On the other hand, $`|B|`$ varies from 0.00117 to 0.00184 GeV7.

## III Discussion

SU(3) symmetry, using de Swart’s convention, predicts

$`g_{KN\mathrm{\Lambda }}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}(32\alpha _D)g_{\pi NN},`$ (29)
$`g_{KN\mathrm{\Sigma }}`$ $`=`$ $`+(2\alpha _D1)g_{\pi NN},`$ (30)
where $`\alpha _D`$ is the fraction of the D type coupling, $`\alpha _D=\frac{D}{D+F}`$. In Table III we compare our results with previous QCD sum rule estimates  and an SU(3) symmetry prediction, where we denote the error-bar allowing for SU(3) symmetry breaking at the 20 % level. Here we take $`\alpha _D`$ from a recent analysis of hyperon semi-leptonic decay data by Ratcliffe, $`\alpha _D`$=0.64, and $`g_{\pi NN}`$ from an analysis of the $`np`$ data by Ericson et al., $`g_{\pi NN}`$=13.43. A comparison to fitting analyses of experimental data is also provided. SU(3) symmetry predicts $`|g_{KN\mathrm{\Lambda }}/g_{KN\mathrm{\Sigma }}|`$ = 3.55 taking $`\alpha _D`$ = 0.64, while our results show that this ratio is 6.30 using the basic inputs, and the order of SU(3) symmetry breaking is rather huge.

Let us remark on $`g_{\pi NN}`$ which was calculated in Ref.  using the three-point function method. After including dimension 5 and 7 condensates as in the previous section the sum rule becomes

$`\lambda _N^2{\displaystyle \frac{e^{M_N^2/M^2}}{M^2}}M_Ng_{\pi NN}{\displaystyle \frac{f_\pi m_\pi ^2}{\sqrt{2}m_q}}+C\left(e^{M_N^{}^2/M^2}e^{M_N^2/M^2}\right)=`$ (31)
$`\left({\displaystyle \frac{1}{\pi ^2}}E_1M^4{\displaystyle \frac{1}{5\pi ^2}}E_0M^2+{\displaystyle \frac{1}{9}}{\displaystyle \frac{\alpha _s}{\pi }}G^2\right)\overline{q}q,`$ (32)
where C is the unknown constant from $`\lambda _N^{}g_{\pi NN^{}}`$ and $`f_\pi `$ = 0.133 GeV. The contribution of the quark-gluon condensates in the OPE side is important as in the $`g_{KN\mathrm{\Lambda }}`$ and $`g_{KN\mathrm{\Sigma }}`$ sum rules. In this case we use the PCAC relation $`f_\pi ^2m_\pi ^2=4m_q\overline{q}q`$ first, then the quark condensate becomes an overall factor on both sides. However, the coupling strength $`\lambda _N`$ is still related to the quark condensate as shown in Eq. (15).

Following the same method in the previous section, and using $`\overline{q}q`$ = –(0.230 GeV)3, $`\frac{\alpha _s}{\pi }G^2`$ = 0.012 GeV4, and $`s_0`$ = 2.074 GeV2 as a pure continuum threshold we get

$`|g_{\pi NN}|`$ $`=`$ $`3.65\pm 2.31,`$ (33)
$`|C|`$ $`=`$ $`0.00261\pm 0.00091\mathrm{GeV}^7,`$ (34)
the Borel interval (0.460, 1.110) GeV2, and the average deviation of the fit $`\overline{\delta }`$ 9.3 % at the central value. Here the uncertainty comes from using different input parameters, i.e. $`\overline{q}q`$ = –(0.210 GeV)3 (or –(0.250 GeV)3) , $`\frac{\alpha _s}{\pi }G^2`$ = 0.015 GeV4, and $`m_0^2`$ = 0.6 (or 1.4 GeV2) as before. In this case the most error bar comes from uncertainties of the quark condensate, i.e. from the coupling strength $`\lambda _N`$.

Now, let us discuss some uncertainties in our sum rules. In Eqs. (13), (25) and (32) the contribution of the quark-gluon condensate is about 25 %, 40 %, and 20 %, respectively, of the leading term at $`M^2`$ = 1 GeV2. Thus the accurate value of this condensate is one of important factors in our sum rules, and a more precise estimate may be needed (e.g., see Ref.).

As we mentioned before, we need more higher dimensional operators to get some power correction terms in our sum rules. Their contribution will be much smaller than that of dimension 7 operators at the relevant Borel region around $`M^2`$ $``$ 1 GeV2. However, those operators may contribute because the lower limit of the Borel interval for each coupling constant is much less than 1 GeV2 in our sum rules.

We find that the coupling constants become 2 or 3 times larger than the previous ones if we take the coupling strengths ($`\lambda _N`$, $`\lambda _\mathrm{\Lambda }`$ and $`\lambda _\mathrm{\Sigma }`$) from the chiral-odd baryon sum rules . For example, we get 7.04, 0.890, and 14.49 for $`|g_{KN\mathrm{\Lambda }}|`$, $`|g_{KN\mathrm{\Sigma }}|`$, and $`|g_{\pi NN}|`$, respectively, for the basic inputs. Because the coupling strengths from each baryon sum rule (the chiral-even and chiral-odd) are not the same in the whole Borel region and the discrepancy between the coupling strengths is larger in the low Borel region, we get quite different coupling constants. Of course, it should be judged by the stability of the sum rule whether one chooses the coupling strengths from the chiral-even sum rules or those from the chiral-odd sum rules.

As a final remark, in the case of $`g_{\pi NN}`$ it was shown that there is a higher pseudoscalar resonance contamination from the $`\pi `$ (1300) and $`\pi `$ (1800) in the three-point function method . Maybe there is a similar contamination from the K(1460) and K(1830) on the kaon-baryon couplings. Although the masses of the K(1460) and K(1830) are quite uncertain and these states need further experimental confirmation, we can briefly estimate the contribution of the K(1460) as done in Ref. . Using the parameters from recent works , we get

$`\left[{\displaystyle \frac{f_Mm_M^2}{Q^2+m_M^2}}\right]_{Q^2=1\mathrm{GeV}^2}=21.3\mathrm{and}2.2\mathrm{MeV}`$ (35)
for the kaon and K(1460), respectively. Here $`f_M`$ is the decay constant and $`m_M`$ is the meson mass. We take $`f_K`$ = 108 MeV, $`f_{K(1460)}`$ = 3.3 MeV and $`m_K`$ = 496 MeV, $`m_{K(1460)}`$ = 1.45 GeV in Ref. . Comparing the values in Eq. (35) to those for the pion and $`\pi `$ (1300), i.e. 1.7 and 0.4 MeV, the contamination from the excited kaon state on the kaon-baryon couplings seems smaller than that from the excited pion state on $`g_{\pi NN}`$.

In summary, including higher dimensional condensates we re-analyze our previous QCD sum rule estimate on $`g_{KN\mathrm{\Lambda }}`$ and $`g_{KN\mathrm{\Sigma }}`$ in the $`\text{/}qi\gamma _5`$ structure. The contribution of dimension 5 quark-gluon condensates is comparable to that of the leading term, and the present result is much different from the previous one.

###### Acknowledgements.

The author thanks Prof. Su H. Lee for valuable discussions and comments. This work is supported by Research Fellowship of the Japan Society for the Promotion of Science (JSPS).


text:
# Variational–Wavelet Approach to RMS Envelope Equations

## 1 Introduction

In this paper we consider the applications of a new numerical-analytical technique which is based on the methods of local nonlinear Fourier analysis or Wavelet analysis to the nonlinear beam/accelerator physics problems related to root-mean-square (rms) envelope dynamics . Such approach may be useful in all models in which it is possible and reasonable to reduce all complicated problems related with statistical distributions to the problems described by systems of nonlinear ordinary/partial differential equations. In this paper we consider approach based on the second moments of the distribution functions for the calculation of evolution of rms envelope of a beam.

The rms envelope equations are the most useful for analysis of the beam self–forces (space–charge) effects and also allow to consider both transverse and longitudinal dynamics of space-charge-dominated relativistic high–bright ness axisymmetric/asymmetric beams, which under short laser pulse–driven radio-frequency photoinjectors have fast transition from nonrelativistic to relativistic regime -.

From the formal point of view we may consider rms envelope equations after straightforward transformations to standard Cauchy form as a system of nonlinear differential equations which are not more than rational (in dynamical variables). Such rational type of nonlinearities allow us to consider some extension of results from -, which are based on application of wavelet analysis technique to variational formulation of initial nonlinear problem.

Wavelet analysis is a relatively novel set of mathematical methods, which gives us a possibility to work with well-localized bases in functional spaces and give for the general type of operators (differential, integral, pseudodifferential) in such bases the maximum sparse forms.

An example of such type of basis is demonstrated on Fig. 1.

Our approach in this paper is based on the generalization  of variational-wavelet approach from -, which allows us to consider not only polynomial but rational type of nonlinearities.

So, our variational-multiresolution approach gives us possibility to construct explicit numerical-analytical solution for the following systems of nonlinear differential equations

$$\dot{z}=R(z,t)\text{or}Q(z,t)\dot{z}=P(z,t),$$
(1)
where $`z(t)=(z_1(t),\mathrm{},z_n(t))`$ is the vector of dynamical variables $`z_i(t)`$,

$`R(z,t)`$ is not more than rational function of z,

$`P(z,t),Q(z,t)`$ are not more than polynomial functions of z and P,Q,R have arbitrary dependence of time.

The solution has the following form

$$z(t)=z_N^{slow}(t)+\underset{jN}{}z_j(\omega _jt),\omega _j2^j$$
(2)
which corresponds to the full multiresolution expansion in all time scales. Formula (2) gives us expansion into a slow part $`z_N^{slow}`$ and fast oscillating parts for arbitrary N. So, we may move from coarse scales of resolution to the finest one for obtaining more detailed information about our dynamical process. The first term in the RHS of equation (2) corresponds on the global level of function space decomposition to resolution space and the second one to detail space. In this way we give contribution to our full solution from each scale of resolution or each time scale (detailed description we give in part 3.2 and numerical illustration in part 7 below). The same is correct for the contribution to power spectral density (energy spectrum): we can take into account contributions from each level/scale of resolution.

In part 2 we describe the different forms of rms equations. Starting in part 3.1 from variational formulation of initial dynamical problem we construct via multiresolution analysis (3.2) explicit representation for all dynamical variables in the base of compactly supported (Daubechies) wavelets. Our solutions (3.3) are parametrized by solutions of a number of reduced algebraical problems one from which is nonlinear with the same degree of nonlinearity and the rest are the linear problems which correspond to particular method of calculation of scalar products of functions from wavelet bases and their derivatives. Then we consider further extension of our previous results. In part 4 we consider modification of our construction to the periodic case, in part 5 we consider generalization of our approach to variational formulation in the biorthogonal bases of compactly supported wavelets and in part 6 to the case of variable coefficients. In part 7 we consider results of numerical calculations.

## 2 RMS Equations

Below we consider a number of different forms of RMS envelope equations, which are from the formal point of view not more than nonlinear differential equations with rational nonlinearities and variable coefficients. Let $`f(x_1,x_2)`$ be the distribution function which gives full information about noninteracting ensemble of beam particles regarding to trace space or transverse phase coordinates $`(x_1,x_2)`$. Then (n,m) moments are:

$$x_1^nx_2^mf(x_1,x_2)dx_1dx_2$$
(3)
The (0,0) moment gives normalization condition on the distribution. The (1,0) and (0,1) moments vanish when a beam is aligned to its axis. Then we may extract the first nontrivial bit of ‘dynamical information’ from the second moments

$`\sigma _{x_1}^2`$ $`=`$ $`<x_1^2>={\displaystyle x_1^2f(x_1,x_2)dx_1dx_2}`$
$`\sigma _{x_2}^2`$ $`=`$ $`<x_2^2>={\displaystyle x_2^2f(x_1,x_2)dx_1dx_2}`$ (4)
$`\sigma _{x_1x_2}^2`$ $`=`$ $`<x_1x_2>={\displaystyle x_1x_2f(x_1,x_2)dx_1dx_2}`$
RMS emittance ellipse is given by

$$\epsilon _{x,rms}^2=<x_1^2><x_2^2><x_1x_2>^2$$
(5)
Expressions for twiss parameters are also based on the second moments.

We will consider the following particular cases of rms envelope equations, which described evolution of the moments (4) (- for full designation):

for asymmetric beams we have the system of two envelope equations of the second order for $`\sigma _{x_1}`$ and $`\sigma _{x_2}`$:

$`\sigma _{x_1}^{^{\prime \prime }}+\sigma _{x_1}^{^{}}{\displaystyle \frac{\gamma ^{}}{\gamma }}+\mathrm{\Omega }_{x_1}^2\left({\displaystyle \frac{\gamma ^{}}{\gamma }}\right)^2\sigma _{x_1}`$ $`=`$ $`{\displaystyle \frac{I}{I_0(\sigma _{x_1}+\sigma _{x_2})\gamma ^3}}+{\displaystyle \frac{\epsilon _{nx_1}^2}{\sigma _{x_1}^3\gamma ^2}},`$
$`\sigma _{x_2}^{^{\prime \prime }}+\sigma _{x_2}^{^{}}{\displaystyle \frac{\gamma ^{}}{\gamma }}+\mathrm{\Omega }_{x_2}^2\left({\displaystyle \frac{\gamma ^{}}{\gamma }}\right)^2\sigma _{x_2}`$ $`=`$ $`{\displaystyle \frac{I}{I_0(\sigma _{x_1}+\sigma _{x_2})\gamma ^3}}+{\displaystyle \frac{\epsilon _{nx_2}^2}{\sigma _{x_2}^3\gamma ^2}}`$ (6)
the envelope equation for an axisymmetric beam is

$$\sigma ^{^{\prime \prime }}+\sigma ^{^{}}\frac{\gamma ^{}}{\gamma }+\mathrm{\Omega }^2\left(\frac{\gamma ^{}}{\gamma }\right)^2\sigma =\frac{k_s}{\sigma \gamma ^3}+\frac{\epsilon _{n,th}^2}{\sigma ^3\gamma ^2}$$
(7)
Also we have related Lawson’s equation for evolution of the rms envelope in the paraxial limit, which governs evolution of cylindrical symmetric envelope under external linear focusing channel of strenghts $`K_r`$:

$$\sigma ^{^{\prime \prime }}+\sigma ^{^{}}\left(\frac{\gamma ^{}}{\beta ^2\gamma }\right)+K_r\sigma =\frac{k_s}{\sigma \beta ^3\gamma ^3}+\frac{\epsilon _n^2}{\sigma ^3\beta ^2\gamma ^2},$$

where

$$K_rF_r/r\beta ^2\gamma mc^2,\beta \nu _b/c=\sqrt{1\gamma ^2}$$

After transformations to Cauchy form we can see that all this equations from the formal point of view are not more than ordinary differential equations with rational nonlinearities and variable coefficients and correspond to the form (1) (also,we may consider regimes in which $`\gamma `$, $`\gamma ^{}`$ are not fixed functions/constants but satisfy some additional differential constraint/equation, but this case does not change our general approach).

## 3 Rational Dynamics

The first main part of our consideration is some variational approach to this problem, which reduces initial problem to the problem of solution of functional equations at the first stage and some algebraical problems at the second stage. We have the solution in a compactly supported wavelet basis. Multiresolution expansion is the second main part of our construction. The solution is parameterized by solutions of two reduced algebraical problems, one is nonlinear and the second are some linear problems, which are obtained from one of the next wavelet constructions: the method of Connection Coefficients (CC), Stationary Subdivision Schemes (SSS).

### 3.1  Variational Method

Our problems may be formulated as the systems of ordinary differential equations

$`Q_i(x){\displaystyle \frac{\mathrm{d}x_i}{\mathrm{d}t}}=P_i(x,t),x=(x_1,\mathrm{},x_n),`$ (8)
$`i=1,\mathrm{},n,\underset{i}{\mathrm{max}}degP_i=p,\underset{i}{\mathrm{max}}degQ_i=q`$
with fixed initial conditions $`x_i(0)`$, where $`P_i,Q_i`$ are not more than polynomial functions of dynamical variables $`x_j`$ and have arbitrary dependence of time. Because of time dilation we can consider only next time interval: $`0t1`$. Let us consider a set of functions

$`\mathrm{\Phi }_i(t)=x_i{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}(Q_iy_i)+P_iy_i`$ (9)
and a set of functionals

$`F_i(x)={\displaystyle _0^1}\mathrm{\Phi }_i(t)𝑑tQ_ix_iy_i_0^1,`$ (10)
where $`y_i(t)(y_i(0)=0)`$ are dual (variational) variables. It is obvious that the initial system and the system

$$F_i(x)=0$$
(11)
are equivalent. Of course, we consider such $`Q_i(x)`$ which do not lead to the singular problem with $`Q_i(x)`$, when $`t=0`$ or $`t=1`$, i.e. $`Q_i(x(0)),Q_i(x(1))\mathrm{}`$.

In part 5 we consider more general approach, which is based on possibility taking into account underlying symplectic structure and on more useful and flexible analytical approach, related to bilinear structure of initial functional. Now we consider formal expansions for $`x_i,y_i`$:

$`x_i(t)=x_i(0)+{\displaystyle \underset{k}{}}\lambda _i^k\phi _k(t)y_j(t)={\displaystyle \underset{r}{}}\eta _j^r\phi _r(t),`$ (12)
where $`\phi _k(t)`$ are useful basis functions of some functional space ($`L^2,L^p`$, Sobolev, etc) corresponding to concrete problem and because of initial conditions we need only $`\phi _k(0)=0`$.

$$\lambda =\{\lambda _i\}=\{\lambda _i^r\}=(\lambda _i^1,\lambda _i^2,\mathrm{},\lambda _i^N),r=1,\mathrm{},N,i=1,\mathrm{},n,$$
(13)
where the lower index i corresponds to expansion of dynamical variable with index i, i.e. $`x_i`$ and the upper index $`r`$ corresponds to the numbers of terms in the expansion of dynamical variables in the formal series. Then we put (12) into the functional equations (11) and as result we have the following reduced algebraical system of equations on the set of unknown coefficients $`\lambda _i^k`$ of expansions (12):

$`L(Q_{ij},\lambda ,\alpha _I)=M(P_{ij},\lambda ,\beta _J),`$ (14)
where operators L and M are algebraization of RHS and LHS of initial problem (8), where $`\lambda `$ (13) are unknowns of reduced system of algebraical equations (RSAE)(14).

$`Q_{ij}`$ are coefficients (with possible time dependence) of LHS of initial system of differential equations (8) and as consequence are coefficients of RSAE.

$`P_{ij}`$ are coefficients (with possible time dependence) of RHS of initial system of differential equations (8) and as consequence are coefficients of RSAE.

$`I=(i_1,\mathrm{},i_{q+2}),J=(j_1,\mathrm{},j_{p+1})`$ are multiindexes, by which are labelled $`\alpha _I`$ and $`\beta _I`$ — other coefficients of RSAE (14):

$$\beta _J=\{\beta _{j_1\mathrm{}j_{p+1}}\}=\underset{1j_kp+1}{}\phi _{j_k},$$
(15)
where p is the degree of polinomial operator P (8)

$$\alpha _I=\{\alpha _{i_1}\mathrm{}\alpha _{i_{q+2}}\}=\underset{i_1,\mathrm{},i_{q+2}}{}\phi _{i_1}\mathrm{}\dot{\phi _{i_s}}\mathrm{}\phi _{i_{q+2}},$$
(16)
where q is the degree of polynomial operator Q (8), $`i_{\mathrm{}}=(1,\mathrm{},q+2)`$, $`\dot{\phi _{i_s}}=\mathrm{d}\phi _{i_s}/\mathrm{d}t`$.

Now, when we solve RSAE (14) and determine unknown coefficients from formal expansion (12) we therefore obtain the solution of our initial problem. It should be noted if we consider only truncated expansion (12) with N terms then we have from (14) the system of $`N\times n`$ algebraical equations with degree $`\mathrm{}=max\{p,q\}`$ and the degree of this algebraical system coincides with degree of initial differential system. So, we have the solution of the initial nonlinear (rational) problem in the form

$`x_i(t)=x_i(0)+{\displaystyle \underset{k=1}{\overset{N}{}}}\lambda _i^kX_k(t),`$ (17)
where coefficients $`\lambda _i^k`$ are roots of the corresponding reduced algebraical (polynomial) problem RSAE (14). Consequently, we have a parametrization of solution of initial problem by solution of reduced algebraical problem (14). The first main problem is a problem of computations of coefficients $`\alpha _I`$ (16), $`\beta _J`$ (15) of reduced algebraical system. As we will see, these problems may be explicitly solved in wavelet approach.

Next we consider the construction of explicit time solution for our problem. The obtained solutions are given in the form (17), where $`X_k(t)`$ are basis functions and $`\lambda _k^i`$ are roots of reduced system of equations. In our first wavelet case $`X_k(t)`$ are obtained via multiresolution expansions and represented by compactly supported wavelets and $`\lambda _k^i`$ are the roots of corresponding general polynomial system (14) with coefficients, which are given by CC or SSS constructions. According to the variational method to give the reduction from differential to algebraical system of equations we need compute the objects $`\alpha _I`$ and $`\beta _J`$.

### 3.2 Wavelet Framework

Our constructions are based on multiresolution approach. Because affine group of translation and dilations is inside the approach, this method resembles the action of a microscope. We have contribution to final result from each scale of resolution from the whole infinite scale of spaces. More exactly, the closed subspace $`V_j(j𝐙)`$ corresponds to level j of resolution, or to scale j. We consider a r-regular multiresolution analysis of $`L^2(𝐑^n)`$ (of course, we may consider any different functional space) which is a sequence of increasing closed subspaces $`V_j`$:

$$\mathrm{}V_2V_1V_0V_1V_2\mathrm{}$$
(18)
satisfying the following properties:

$`{\displaystyle \underset{j𝐙}{}}V_j=0,\overline{{\displaystyle \underset{j𝐙}{}}}V_j=L^2(𝐑^n),`$
$`f(x)V_j<=>f(2x)V_{j+1},`$
$`f(x)V_0<=>f(xk)V_0,,k𝐙^n.`$ (19)
There exists a function $`\phi V_0`$ such that {$`\phi _{0,k}(x)=\phi (xk),k𝐙^n`$} forms a Riesz basis for $`V_0`$.

The function $`\phi `$ is regular and localized: $`\phi `$ is $`C^{r1},\phi ^{(r1)}`$ is almost everywhere differentiable and for almost every $`x𝐑^n`$, for every integer $`\alpha r`$ and for all integer p there exists constant $`C_p`$ such that

$$^\alpha \phi (x)C_p(1+|x|)^p$$
(20)
Let $`\phi (x)`$ be a scaling function, $`\psi (x)`$ is a wavelet function and $`\phi _i(x)=\phi (xi)`$. Scaling relations that define $`\phi ,\psi `$ are

$`\phi (x)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{N1}{}}}a_k\phi (2xk)={\displaystyle \underset{k=0}{\overset{N1}{}}}a_k\phi _k(2x),`$ (21)
$`\psi (x)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{N2}{}}}(1)^ka_{k+1}\phi (2x+k).`$ (22)
Let indices $`\mathrm{},j`$ represent translation and scaling, respectively and

$$\phi _{jl}(x)=2^{j/2}\phi (2^jx\mathrm{})$$
(23)
then the set $`\{\phi _{j,k}\},k𝐙^n`$ forms a Riesz basis for $`V_j`$. The wavelet function $`\psi `$ is used to encode the details between two successive levels of approximation. Let $`W_j`$ be the orthonormal complement of $`V_j`$ with respect to $`V_{j+1}`$:

$$V_{j+1}=V_jW_j.$$
(24)
Then just as $`V_j`$ is spanned by dilation and translations of the scaling function, so are $`W_j`$ spanned by translations and dilation of the mother wavelet $`\psi _{jk}(x)`$, where

$$\psi _{jk}(x)=2^{j/2}\psi (2^jxk).$$
(25)
All expansions which we used are based on the following properties:

$`\{\psi _{jk}\},j,k𝐙\text{is a Hilbertian basis of }L^2(𝐑)`$
$`\{\phi _{jk}\}_{j0,k𝐙}\text{is an orthonormal basis for}L^2(𝐑),`$
$`L^2(𝐑)=\overline{V_0{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}W_j},`$ (26)
or
$`\{\phi _{0,k},\psi _{j,k}\}_{j0,k𝐙}\text{is an orthonormal basis for}L^2(𝐑).`$
### 3.3  Wavelet Computations

Now we give construction for computations of objects (17),(18) in the wavelet case. We use compactly supported wavelet basis: orthonormal basis for functions in $`L^2(𝐑)`$.

Let be $`f:𝐑𝐂`$ and the wavelet expansion is

$`f(x)={\displaystyle \underset{\mathrm{}𝐙}{}}c_{\mathrm{}}\phi _{\mathrm{}}(x)+{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k𝐙}{}}c_{jk}\psi _{jk}(x)`$ (27)
If in formulae (29) $`c_{jk}=0`$ for $`jJ`$, then $`f(x)`$ has an alternative expansion in terms of dilated scaling functions only $`f(x)=\underset{\mathrm{}𝐙}{}c_J\mathrm{}\phi _J\mathrm{}(x)`$. This is a finite wavelet expansion, it can be written solely in terms of translated scaling functions. Also we have the shortest possible support: scaling function $`DN`$ (where $`N`$ is even integer) will have support $`[0,N1]`$ and $`N/2`$ vanishing moments. There exists $`\lambda >0`$ such that $`DN`$ has $`\lambda N`$ continuous derivatives; for small $`N,\lambda 0.55`$. To solve our second associated linear problem we need to evaluate derivatives of $`f(x)`$ in terms of $`\phi (x)`$. Let be $`\phi _{\mathrm{}}^n=\mathrm{d}^n\phi _{\mathrm{}}(x)/\mathrm{d}x^n`$. We consider computation of the wavelet - Galerkin integrals. Let $`f^d(x)`$ be d-derivative of function $`f(x)`$, then we have $`f^d(x)=_{\mathrm{}}c_l\phi _{\mathrm{}}^d(x)`$, and values $`\phi _{\mathrm{}}^d(x)`$ can be expanded in terms of $`\phi (x)`$

$`\varphi _{\mathrm{}}^d(x)`$ $`=`$ $`{\displaystyle \underset{m}{}}\lambda _m\phi _m(x),`$ (28)
$`\lambda _m`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\phi _{\mathrm{}}^d(x)\phi _m(x)dx,`$
where $`\lambda _m`$ are wavelet-Galerkin integrals. The coefficients $`\lambda _m`$ are 2-term connection coefficients. In general we need to find $`(d_i0)`$

$`\mathrm{\Lambda }_{\mathrm{}_1\mathrm{}_2\mathrm{}\mathrm{}_n}^{d_1d_2\mathrm{}d_n}={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \phi _\mathrm{}_i^{d_i}(x)dx}`$ (29)
For Riccati case we need to evaluate two and three connection coefficients

$`\mathrm{\Lambda }_{\mathrm{}}^{d_1d_2}={\displaystyle _{\mathrm{}}^{\mathrm{}}}\phi ^{d_1}(x)\phi _{\mathrm{}}^{d_2}(x)𝑑x,\mathrm{\Lambda }^{d_1d_2d_3}={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\phi ^{d_1}(x)\phi _{\mathrm{}}^{d_2}(x)\phi _m^{d_3}(x)𝑑x`$ (30)
According to CC method  we use the next construction. When $`N`$ in scaling equation is a finite even positive integer the function $`\phi (x)`$ has compact support contained in $`[0,N1]`$. For a fixed triple $`(d_1,d_2,d_3)`$ only some $`\mathrm{\Lambda }_\mathrm{}m^{d_1d_2d_3}`$ are nonzero: $`2N\mathrm{}N2,2NmN2,|\mathrm{}m|N2`$. There are $`M=3N^29N+7`$ such pairs $`(\mathrm{},m)`$. Let $`\mathrm{\Lambda }^{d_1d_2d_3}`$ be an M-vector, whose components are numbers $`\mathrm{\Lambda }_\mathrm{}m^{d_1d_2d_3}`$. Then we have the first reduced algebraical system : $`\mathrm{\Lambda }`$ satisfy the system of equations $`(d=d_1+d_2+d_3)`$

$`A\mathrm{\Lambda }^{d_1d_2d_3}=2^{1d}\mathrm{\Lambda }^{d_1d_2d_3},A_{\mathrm{},m;q,r}={\displaystyle \underset{p}{}}a_pa_{q2\mathrm{}+p}a_{r2m+p}`$ (31)
By moment equations we have created a system of $`M+d+1`$ equations in $`M`$ unknowns. It has rank $`M`$ and we can obtain unique solution by combination of LU decomposition and QR algorithm. The second reduced algebraical system gives us the 2-term connection coefficients:

$`A\mathrm{\Lambda }^{d_1d_2}=2^{1d}\mathrm{\Lambda }^{d_1d_2},d=d_1+d_2,A_{\mathrm{},q}={\displaystyle \underset{p}{}}a_pa_{q2\mathrm{}+p}`$ (32)
For nonquadratic case we have analogously additional linear problems for objects (31). Solving these linear problems we obtain the coefficients of nonlinear algebraical system (16) and after that we obtain the coefficients of wavelet expansion (19). As a result we obtained the explicit time solution of our problem in the base of compactly supported wavelets. We use for modelling D6, D8, D10 functions and programs RADAU and DOPRI for testing.

In the following we consider extension of this approach to the case of periodic boundary conditions, the case of presence of arbitrary variable coefficients and more flexible biorthogonal wavelet approach.

## 4 Variational Wavelet Approach for Periodic Trajectories

We start with extension of our approach to the case of periodic trajectories. The equations of motion corresponding to our problems may be formulated as a particular case of the general system of ordinary differential equations $`dx_i/dt=f_i(x_j,t)`$, $`(i,j=1,\mathrm{},n)`$, $`0t1`$, where $`f_i`$ are not more than rational functions of dynamical variables $`x_j`$ and have arbitrary dependence of time but with periodic boundary conditions. According to our variational approach we have the solution in the following form

$`x_i(t)=x_i(0)+{\displaystyle \underset{k}{}}\lambda _i^k\phi _k(t),x_i(0)=x_i(1),`$ (33)
where $`\lambda _i^k`$ are again the roots of reduced algebraical systems of equations with the same degree of nonlinearity and $`\phi _k(t)`$ corresponds to useful type of wavelet bases (frames). It should be noted that coefficients of reduced algebraical system are the solutions of additional linear problem and also depend on particular type of wavelet construction and type of bases.

This linear problem is our second reduced algebraical problem. We need to find in general situation objects

$`\mathrm{\Lambda }_{\mathrm{}_1\mathrm{}_2\mathrm{}\mathrm{}_n}^{d_1d_2\mathrm{}d_n}={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \phi _\mathrm{}_i^{d_i}(x)\mathrm{d}x},`$ (34)
but now in the case of periodic boundary conditions. Now we consider the procedure of their calculations in case of periodic boundary conditions in the base of periodic wavelet functions on the interval  and corresponding expansion (35) inside our variational approach. Periodization procedure gives us

$`\widehat{\phi }_{j,k}(x)`$ $``$ $`{\displaystyle \underset{\mathrm{}Z}{}}\phi _{j,k}(x\mathrm{})`$ (35)
$`\widehat{\psi }_{j,k}(x)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}Z}{}}\psi _{j,k}(x\mathrm{})`$
So, $`\widehat{\phi },\widehat{\psi }`$ are periodic functions on the interval . Because $`\phi _{j,k}=\phi _{j,k^{}}`$ if $`k=k^{}\mathrm{mod}(2^j)`$, we may consider only $`0k2^j`$ and as consequence our multiresolution has the form $`{\displaystyle \underset{j0}{}}\widehat{V}_j=L^2[0,1]`$ with $`\widehat{V}_j=\mathrm{span}\{\widehat{\phi }_{j,k}\}_{k=0}^{2j1}`$ . Integration by parts and periodicity gives useful relations between objects (36) in particular quadratic case $`(d=d_1+d_2)`$:

$`\mathrm{\Lambda }_{k_1,k_2}^{d_1,d_2}=(1)^{d_1}\mathrm{\Lambda }_{k_1,k_2}^{0,d_2+d_1},\mathrm{\Lambda }_{k_1,k_2}^{0,d}=\mathrm{\Lambda }_{0,k_2k_1}^{0,d}\mathrm{\Lambda }_{k_2k_1}^d`$ (36)
So, any 2-tuple can be represent by $`\mathrm{\Lambda }_k^d`$. Then our second additional linear problem is reduced to the eigenvalue problem for $`\{\mathrm{\Lambda }_k^d\}_{0k2^j}`$ by creating a system of $`2^j`$ homogeneous relations in $`\mathrm{\Lambda }_k^d`$ and inhomogeneous equations. So, if we have dilation equation in the form $`\phi (x)=\sqrt{2}_{kZ}h_k\phi (2xk)`$, then we have the following homogeneous relations

$$\mathrm{\Lambda }_k^d=2^d\underset{m=0}{\overset{N1}{}}\underset{\mathrm{}=0}{\overset{N1}{}}h_mh_{\mathrm{}}\mathrm{\Lambda }_{\mathrm{}+2km}^d,$$
(37)
or in such form $`A\lambda ^d=2^d\lambda ^d`$, where $`\lambda ^d=\{\mathrm{\Lambda }_k^d\}_{0k2^j}`$. Inhomogeneous equations are:

$$\underset{\mathrm{}}{}M_{\mathrm{}}^d\mathrm{\Lambda }_{\mathrm{}}^d=d!2^{j/2},$$
(38)
where objects $`M_{\mathrm{}}^d(|\mathrm{}|N2)`$ can be computed by recursive procedure

$$M_{\mathrm{}}^d=2^{j(2d+1)/2}\stackrel{~}{M_{\mathrm{}}^d},\stackrel{~}{M_{\mathrm{}}^k}=<x^k,\phi _{0,\mathrm{}}>=\underset{j=0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{k}{j}\right)n^{kj}M_0^j,\stackrel{~}{M_0^{\mathrm{}}}=1.$$
(39)
So, we reduced our last problem to standard linear algebraical problem. Then we use the same methods as in part 3.3. As a result we obtained for closed trajectories of orbital dynamics the explicit time solution (35) in the base of periodized wavelets (37).

## 5 Variational Approach in Biorthogonal <br>Wavelet Bases

Now we consider further generalization of our variational wavelet approach. Because integrand of variational functionals is represented by bilinear form (scalar product) it seems more reasonable to consider wavelet constructions  which take into account all advantages of this structure.

The action functional for loops in the phase space is

$$F(\gamma )=_\gamma p𝑑q_0^1H(t,\gamma (t))𝑑t$$
(40)
The critical points of $`F`$ are those loops $`\gamma `$, which solve the Hamiltonian equations associated with the Hamiltonian $`H`$ and hence are periodic orbits.

Let us consider the loop space $`\mathrm{\Omega }=C^{\mathrm{}}(S^1,R^{2n})`$, where $`S^1=R/𝐙`$, of smooth loops in $`R^{2n}`$. Let us define a function $`\mathrm{\Phi }:\mathrm{\Omega }R`$ by setting

$$\mathrm{\Phi }(x)=_0^1\frac{1}{2}<J\dot{x},x>dt_0^1H(x(t))𝑑t,x\mathrm{\Omega }$$
(41)
The critical points of $`\mathrm{\Phi }`$ are the periodic solutions of $`\dot{x}=X_H(x)`$. Computing the derivative at $`x\mathrm{\Omega }`$ in the direction of $`y\mathrm{\Omega }`$, we find

$`\mathrm{\Phi }^{}(x)(y)={\displaystyle \frac{d}{dϵ}}\mathrm{\Phi }(x+ϵy)|_{ϵ=0}={\displaystyle _0^1}<J\dot{x}H(x),y>dt`$ (42)
Consequently, $`\mathrm{\Phi }^{}(x)(y)=0`$ for all $`y\mathrm{\Omega }`$ iff the loop $`x`$ satisfies the equation

$$J\dot{x}(t)H(x(t))=0,$$
(43)
i.e. $`x(t)`$ is a solution of the Hamiltonian equations, which also satisfies $`x(0)=x(1)`$, i.e. periodic of period 1.

But now we need to take into account underlying bilinear structure via wavelets.

We started with two hierarchical sequences of approximations spaces :

$`\mathrm{}V_2V_1V_0V_1V_2\mathrm{},\mathrm{}\stackrel{~}{V}_2\stackrel{~}{V}_1\stackrel{~}{V}_0\stackrel{~}{V}_1\stackrel{~}{V}_2\mathrm{},`$
and as usually, $`W_0`$ is complement to $`V_0`$ in $`V_1`$, but now not necessarily orthogonal complement. New orthogonality conditions have now the following form:

$$\stackrel{~}{W}_0V_0,W_0\stackrel{~}{V}_0,V_j\stackrel{~}{W}_j,\stackrel{~}{V}_jW_j$$
(44)
translates of $`\psi `$ $`\mathrm{span}`$ $`W_0`$, translates of $`\stackrel{~}{\psi }\mathrm{span}\stackrel{~}{W}_0`$. Biorthogonality conditions are

$$<\psi _{jk},\stackrel{~}{\psi }_{j^{}k^{}}>=_{\mathrm{}}^{\mathrm{}}\psi _{jk}(x)\stackrel{~}{\psi }_{j^{}k^{}}(x)dx=\delta _{kk^{}}\delta _{jj^{}},$$
(45)
where $`\psi _{jk}(x)=2^{j/2}\psi (2^jxk)`$. Functions $`\phi (x),\stackrel{~}{\phi }(xk)`$ form dual pair:

$$<\phi (xk),\stackrel{~}{\phi }(x\mathrm{})>=\delta _{kl},<\phi (xk),\stackrel{~}{\psi }(x\mathrm{})>=0\text{for}k,\mathrm{}.$$
(46)
Functions $`\phi ,\stackrel{~}{\phi }`$ generate a multiresolution analysis. $`\phi (xk)`$, $`\psi (xk)`$ are synthesis functions, $`\stackrel{~}{\phi }(x\mathrm{})`$, $`\stackrel{~}{\psi }(x\mathrm{})`$ are analysis functions. Synthesis functions are biorthogonal to analysis functions. Scaling spaces are orthogonal to dual wavelet spaces. Two multiresolutions are intertwining $`V_j+W_j=V_{j+1},\stackrel{~}{V}_j+\stackrel{~}{W}_j=\stackrel{~}{V}_{j+1}`$. These are direct sums but not orthogonal sums.

So, our representation for solution has now the form

$$f(t)=\underset{j,k}{}\stackrel{~}{b}_{jk}\psi _{jk}(t),$$
(47)
where synthesis wavelets are used to synthesize the function. But $`\stackrel{~}{b}_{jk}`$ come from inner products with analysis wavelets. Biorthogonality yields

$$\stackrel{~}{b}_\mathrm{}m=f(t)\stackrel{~}{\psi }_\mathrm{}m(t)dt.$$
(48)
So, now we can introduce this more complicated construction into our variational approach. We have modification only on the level of computing coefficients of reduced nonlinear algebraical system. This new construction is more flexible. Biorthogonal point of view is more stable under the action of large class of operators while orthogonal (one scale for multiresolution) is fragile, all computations are much more simpler and we accelerate the rate of convergence. In all types of (Hamiltonian) calculation, which are based on some bilinear structures (symplectic or Poissonian structures, bilinear form of integrand in variational integral) this framework leads to greater success.

## 6 Variable Coefficients

In the case when we have situation when our problem is described by a system of nonlinear (rational) differential equations, we need to consider extension of our previous approach which can take into account any type of variable coefficients (periodic, regular or singular). We can produce such approach if we add in our construction additional refinement equation, which encoded all information about variable coefficients . According to our variational approach we need to compute integrals of the form

$$_Db_{ij}(t)(\phi _1)^{d_1}(2^mtk_1)(\phi _2)^{d_2}(2^mtk_2)dx,$$
(49)
where now $`b_{ij}(t)`$ are arbitrary functions of time, where trial functions $`\phi _1,\phi _2`$ satisfy a refinement equations:

$$\phi _i(t)=\underset{k𝐙}{}a_{ik}\phi _i(2tk)$$
(50)
If we consider all computations in the class of compactly supported wavelets then only a finite number of coefficients do not vanish. To approximate the non-constant coefficients, we need choose a different refinable function $`\phi _3`$ along with some local approximation scheme

$$(B_{\mathrm{}}f)(x):=\underset{\alpha 𝐙}{}F_{\mathrm{},k}(f)\phi _3(2^{\mathrm{}}tk),$$
(51)
where $`F_{\mathrm{},k}`$ are suitable functionals supported in a small neighborhood of $`2^{\mathrm{}}k`$ and then replace $`b_{ij}`$ in (49) by $`B_{\mathrm{}}b_{ij}(t)`$. In particular case one can take a characteristic function and can thus approximate non-smooth coefficients locally. To guarantee sufficient accuracy of the resulting approximation to (49) it is important to have the flexibility of choosing $`\phi _3`$ different from $`\phi _1,\phi _2`$. In the case when D is some domain, we can write

$$b_{ij}(t)_D=\underset{0k2^{\mathrm{}}}{}b_{ij}(t)\chi _D(2^{\mathrm{}}tk),$$
(52)
where $`\chi _D`$ is characteristic function of D. So, if we take $`\phi _4=\chi _D`$, which is again a refinable function, then the problem of computation of (49) is reduced to the problem of calculation of integral

$`H(k_1,k_2,k_3,k_4)=H(k)=`$ (53)
$`{\displaystyle _{𝐑^s}}\phi _4(2^jtk_1)\phi _3(2^{\mathrm{}}tk_2)\phi _1^{d_1}(2^rtk_3)\phi _2^{d_2}(2^stk_4)dx`$
The key point is that these integrals also satisfy some sort of refinement equation :

$$2^{|\mu |}H(k)=\underset{\mathrm{}𝐙}{}b_{2k\mathrm{}}H(\mathrm{}),\mu =d_1+d_2.$$
(54)
This equation can be interpreted as the problem of computing an eigenvector. Thus, we reduced the problem of extension of our method to the case of variable coefficients to the same standard algebraical problem as in the preceding sections. So, the general scheme is the same one and we have only one more additional linear algebraic problem by which we in the same way can parameterize the solutions of corresponding problem.

## 7 Numerical Calculations

In this part we consider numerical illustrations of previous analytical approach. Our numerical calculations are based on compactly supported Daubechies wavelets and related wavelet families.

On Fig. 2 we present according to formulae (2) contributions to approximation of our dynamical evolution (top row on the Fig. 3) starting from the coarse approximation, corresponding to scale $`2^0`$ (bottom row) to the finest one corresponding to the scales from $`2^1`$ to $`2^5`$ or from slow to fast components (5 frequencies) as details for approximation. Then on Fig. 3, from bottom to top, we demonstrate the summation of contributions from corresponding levels of resolution given on Fig. 2 and as result we restore via 5 scales (frequencies) approximation our dynamical process(top row on Fig. 3 ). In this particular model case we considered for approximation simple two frequencies harmonic process. But the same situation we have on the Fig. 5 and Fig. 6 in case when we added to previous 2-frequencies harmonic process the noise as perturbation. Again, our dynamical process under investigation (top row of Fig. 6) is recovered via 5 scales contributions (Fig. 5) to approximations (Fig. 6). The same decomposition/approximation we produce also on the level of power spectral density in the process without noise (Fig. 4) and with noise (Fig. 7). On Fig. 8 we demonstrate the family of localized contributions to beam motion, which we also may consider for such type of approximation.

It should be noted that complexity of such algorithms are minimal regarding other possible. Of course, we may use different multiresolution analysis schemes, which are based on different families of generating wavelets and apply such schemes of numerical–analytical calculations to any dynamical process which may be described by systems of ordinary/partial differential equations with rational nonlinearities .

## Acknowledgments

We would like to thank Professor James B. Rosenzweig and Mrs. Melinda Laraneta for nice hospitality, help, support and discussions before and during Workshop and all participants for interesting discussions.


text:
# On mixing times for stratified walks on the 𝑑-cube

## 1  Introduction.

The stratified random walk (SRW) on the $`d`$-cube $`Q_d`$ is the Markov chain whose state space is the set of vertices of the $`d`$-cube and whose transition probabilities are defined thus:

Given a set of non-zero probabilities $`p=(p_0,p_1,\mathrm{},p_{d1})`$, from any vertex with $`k`$ 1’s, the process moves either to any neighboring vertex with $`k+1`$ 1’s with probability $`\frac{p_k}{d}`$; or to any neighboring vertex with $`k1`$ 1’s with probability $`\frac{p_{k1}}{d}`$; or to itself with the remaining probability. The simple random walk on the $`d`$-cube corresponds to the choice $`p_k=1`$ for all $`k`$.

Vaguely speaking, the mixing time of a Markov chain is the time it takes the chain to have its distribution close to the stationary distribution under some measure of closeness. Chung and Graham studied the SRW on the $`d`$-cube in , mainly with algebraic methods, and found bounds for the mixing times under total variation and relative pointwise distances. Here we use non-algebraic methods, the electric and coupling approaches, in order to study the same SRW and get exact results for maximal commute times and bounds for cover times and mixing times under total variation distance. We take advantage of the fact that there seems to be some inequality or another linking hitting times, commute times, cover times and any definition of mixing time with any other under any measure of closeness (see Aldous and Fill  and Lovász and Winkler ).

## 2 The electric approach

On a connected undirected graph $`G=(V,E)`$ such that the edge between vertices $`i`$ and $`j`$ is given a resistance $`r_{ij}`$ (or equivalently, a conductance $`C_{ij}=1/r_{ij}`$), we can define the random walk on $`G`$ as the Markov chain $`𝐗=\{𝐗(n)\}_{n0}`$ that from its current vertex $`v`$ jumps to the neighboring vertex $`w`$ with probability $`p_{vw}=C_{vw}/C(v)`$, where $`C(v)=_{w:wv}C_{vw}`$, and $`wv`$ means that $`w`$ is a neighbor of $`v`$. There may be a conductance $`C_{zz}`$ from a vertex $`z`$ to itself, giving rise to a transition probability from $`z`$ to itself. Some notation: $`𝔼_aT_b`$ and $`𝔼_aC`$ denote the expected value, starting from the vertex $`a`$, of respectively, the hitting time $`T_b`$ of the vertex $`b`$ and the cover time $`C`$, i. e., the number of jumps needed to visit all the states in $`V`$; $`R_{ab}`$ is the effective resistance, as computed by means of Ohm’s law, between vertices $`a`$ and $`b`$.

A Markov chain is reversible if $`\pi _i(i,j)=\pi _j(j,i)`$ for all $`i,j`$, where $`\{\pi _.\}`$ is the stationary distribution and $`(,)`$ are the transition probabilities. Such a reversible Markov chain can be described as a random walk on a graph if we define conductances thus:

$$C_{ij}=\pi _i(i,j).$$
(2.1)
We will be interested in finding a closed form expression for the commute time $`𝔼_0T_d+𝔼_dT_0`$ between the origin, denoted by 0, and its opposite vertex, denoted by $`d`$.

Notice first that the transition matrix for $`𝐗=\{𝐗(n),n0\}`$, the SRW on the $`d`$-cube, is doubly stochastic and therefore its stationary distribution is uniform. If we now collapse all vertices in the cube with the same number of 1’s into a single vertex, and we look at the SRW on this collapsed graph, we obtain a new reversible Markov chain $`𝐒=\{𝐒(n),n0\}`$, a birth-and-death chain in fact, on the state space $`\{0,1,\mathrm{}d\}`$, with transition probabilities

$`(k,k+1)`$ $`=`$ $`{\displaystyle \frac{dk}{d}}p_k,`$ (2.2)
$`(k,k1)`$ $`=`$ $`{\displaystyle \frac{k}{d}}p_{k1},`$ (2.3)
$`(k,k)`$ $`=`$ $`1(k,k+1)(k,k1).`$ (2.4)
It is plain to see that the stationary distribution of this new chain is the Binomial with parameters $`d`$ and $`\frac{1}{2}`$. It is also clear that the commute time between vertices 0 and $`d`$ is the same for both $`𝐗`$ and $`𝐒`$. For the latter we use the electric machinery described above, namely, we think of a linear electric circuit from $`0`$ to $`d`$ with conductances given by (2.1) for $`0id`$, $`j=i1,i,i+1`$, and where $`\pi _i=\left({\displaystyle \genfrac{}{}{0pt}{}{d}{i}}\right){\displaystyle \frac{1}{2^d}}`$.

It is well known (at least since Chandra et al. proved it in ) that

$$𝔼_aT_b+𝔼_aT_b=R_{ab}\underset{z}{}C(z),$$
(2.5)
where $`R_{ab}`$ is the effective resistance between vertices $`a`$ and $`b`$.

If this formula is applied to a reversible chain whose conductances are given as in (2.1), then it is clear that

$$C(z)=\pi _z$$

and therefore the summation in (2.5) equals 1. We get then this compact formula for the commute time:

$$𝔼_aT_b+𝔼_aT_b=R_{ab},$$
(2.6)
where the effective resistance is computed with the individual resistors having resistances

$$r_{ij}=\frac{1}{C_{ij}}=\frac{1}{\pi _i(i,j)}.$$

In our particular case of the collapsed chain, because it is a linear circuit, the effective resistance $`R_{0d}`$ equals the sum of all the individual resistances $`r_{i,i+1}`$, so that (2.6) yields

$$𝔼_0T_d+𝔼_dT_0=R_{0d}=2^d\underset{k=0}{\overset{d1}{}}\frac{1}{p_k\left(\genfrac{}{}{0pt}{}{d1}{k}\right)}.$$
(2.7)
Because of the particular nature of the chain under consideration, it is clear that $`𝔼_0T_d+𝔼_dT_0`$ equals the maximal commute time ($`\tau ^{}`$ in the terminology of Aldous ) between any two vertices.

(i) For simple random walk, formula (2.7) is simplified by taking all $`p_k=1`$. This particular formula was obtained in  with a more direct argument, and it was argued there that

$$\underset{k=0}{\overset{d1}{}}\frac{1}{\left(\genfrac{}{}{0pt}{}{d1}{k}\right)}=2+o(1).$$

An application of Matthews’ result (see ), linking maximum and minimum expected hitting times with expected cover times, yields immediately that the expected cover time is $`𝔼_vC=\mathrm{\Theta }(|V|\mathrm{log}|V|),`$ which is the asymptotic value of the lower bound for cover times of walks on a graph $`G=(V,E)`$ (see ). Thus we could say this SRW is a “rapidly covered” walk.

(ii) The so-called Aldous cube (see ) corresponds to the choice $`p_k=\frac{k}{d1}`$. This walk takes place in the “punctured cube” that excludes the origin. Formula (2.7) thus, must exclude $`k=0`$ in this case, for which we still get a closed-form expression for the commute time between vertex $`d`$, all of whose coordinates are 1, and vertex $`s`$ which consists of the collapse of all vertices with a single 1:

$$𝔼_sT_d+𝔼_dT_s=2^d\underset{k=1}{\overset{d1}{}}\frac{1}{\left(\genfrac{}{}{0pt}{}{d2}{k1}\right)}.$$
(2.8)
The same argument used in (i) tells us that the summation in (2.8) equals $`2+o(1)`$ and, once again, Matthews’ result tells us that the walk on the Aldous cube has a cover time of order $`|V|\mathrm{log}|V|`$.

(iii) The choice $`p_k=\frac{1}{\left(\genfrac{}{}{0pt}{}{d1}{k}\right)}`$ would be in the terminology of  a “slow walk”: the commute time is seen to be exactly equal to $`|V|\mathrm{log}_2|V|`$ and thus the expected cover time is $`O(|V|\mathrm{log}^2|V|)`$. In general, the SRW will be rapidly covered if and only if

$$\underset{k=0}{\overset{d1}{}}\frac{1}{p_k\left(\genfrac{}{}{0pt}{}{d1}{k}\right)}=c+o(1),$$

for some constant $`c`$.

Remark. A formula as compact as (2.7) could be easily obtained through the commute time formula (2.6). It does not seem that it could be obtained that easily, by just adding the individual hitting times $`𝔼_iT_{i+1}`$. (A procedure that is done, for instance, in , , , and in the next section).

## 3 The coupling approach

In order to asess the rate of convergence of the SRW on the cube $`Q_d`$ to the uniform stationary distribution $`\pi `$, we will bound the mixing time $`\tau `$ defined as

$$\tau =\mathrm{min}\{t:d(t^{})\frac{1}{2e},\text{ for all }t^{}>t\},$$

where

$$d(t)=\underset{𝐱Q_d}{\mathrm{max}}P_𝐱(𝐗(t)=)\pi (),$$

and $`\theta _1\theta _2`$ is the variation distance between probability distributions $`\theta _1`$ and $`\theta _2`$, one of whose alternative definitions is (see Aldous and Fill ), chapter 2):

$$\theta _1\theta _2=\mathrm{min}(V_1V_2),$$

where the minimum is taken over random pairs $`(V_1,V_2)`$ such that $`V_m`$ has distribution $`\theta _m,m=1,2.`$

The bound for the mixing time is achieved using a coupling argument that goes as follows: let $`\{𝐗(t),t0\}`$ and $`\{𝐘(t),t0\}`$ be two versions of the SRW on $`Q_d`$ such that $`𝐗(0)=𝐱`$ and $`𝐘(0)\pi `$. Then

$$_𝐱(𝐗(t)=)\pi ()(𝐗(t)𝐘(t)).$$
(3.1)
A coupling between the processes $`𝐗`$ and $`𝐘`$ is a bivariate process such that its marginals have the distributions of the original processes and such that once the bivariate process enters the diagonal, it stays there forever. If we denote by

$$T_𝐱=inf\{t;𝐗(t)=𝐘(t)\}$$

the coupling time, i. e., the hitting time of the diagonal, then (3.1) translates as

$$_𝐱(𝐗(t)=)\pi ()(T_𝐱>t),$$
(3.2)
and therefore,

$$d(t)\underset{𝐱Q_d}{\mathrm{max}}(T_𝐱>t).$$
(3.3)
If we can find a coupling such that $`𝔼T_𝐱=O(f(d))`$, for all $`xQ_d`$ and for a certain function $`f`$ of the dimension $`d`$, then we will also have $`\tau =O(f(d))`$. Indeed, if we take $`t=2ef(d)`$, then (3.3) and Markov’s inequality imply that $`d(t)1/2e`$ and the definition of $`\tau `$ implies $`\tau =O(f(d))`$.

We will split $`T_𝐱`$ as $`T_𝐱=T_𝐱^1+T_𝐱^2`$, where $`T_𝐱^1`$ is a coupling time for the birth-and-death process $`𝐒`$, and $`T_𝐱^2`$ is another coupling time for the whole process X, once the bivariate $`𝐒`$ process enters the diagonal, and we will bound the values of $`𝔼T_𝐱^1`$ and $`𝔼T_𝐱^2`$.

More formally, for any $`𝐱,𝐲Q_d`$ define

$`s(𝐱)={\displaystyle \underset{i=1}{\overset{d}{}}}x_i`$ (3.4)
$`d(𝐱,𝐲)={\displaystyle \underset{i=1}{\overset{d}{}}}|x_iy_i|.`$ (3.5)
Define also for the birth-and-death process $`𝐒(t)=s(𝐗(t))`$ its own mixing time:

$$\tau ^{(S)}=inf\{t;d_S(t)\frac{1}{2e}\},$$

where $`d_S(t)=\mathrm{max}_i_i(𝐒(t)=)\pi _S()`$, and $`\pi _S`$ is the stationary distribution of $`𝐒`$. Notice that $`s(𝐘(0))\pi _S`$ since $`𝐘(0)\pi `$.

Now we will prove that $`\tau ^{(S)}=O(f_S(d)),`$ for a certain function $`f_S`$ of the expected hitting times of the “central states”, and that this bound implies an analogous bound for $`𝔼T_𝐱^1`$. Indeed, as shown by Aldous , we can bound $`\tau ^{(S)}`$ by a more convenient stopping time

$$\tau ^{(S)}K_2\tau _1^{(3)}$$
(3.6)
where $`\tau _1^{(3)}=\mathrm{min}_\mu \mathrm{max}_i\mathrm{min}_{U_i}𝔼_iU_i`$ and the innermost minimum is taken over stopping times $`U_i`$ such that $`_i(S(U_i))=\mu ()`$. In particular,

$`\tau _1^{(3)}`$ $``$ $`\underset{b}{\mathrm{min}}\underset{i}{\mathrm{max}}\underset{U_i^b}{\mathrm{min}}𝔼_iU_i^b`$ (3.7)
$``$ $`\underset{b}{\mathrm{min}}\underset{i}{\mathrm{max}}𝔼_iT_b`$ (3.8)
$`=`$ $`\mathrm{max}(𝔼_0T_{d/2},𝔼_dT_{d/2})`$ (3.9)
where the innermost minimum in (3.7) is taken over stopping times $`U_i^b`$ such that $`_i(S(U_i)=b)=1`$. Expression (3.9) follows from (3.8) since we are dealing with birth and death chains. Therefore, combining (3.6) and (3.9) we have

$$\tau ^{(S)}K_2\mathrm{max}(𝔼_0T_{d/2},𝔼_dT_{d/2}):=f_S(d).$$
(3.10)
In general, a coupling implies an inequality like (3.2). However, the inequality becomes an equality for a certain maximal (non-Markovian) coupling, described by Griffeath . Let $`T_𝐱^1`$ be the coupling time for the maximal coupling between $`s(𝐗(t))`$ and $`s(𝐘(t))`$ such that

$$P_𝐱(S(𝐗(t))=)\pi _S()=(T_𝐱^1>t).$$

Then

$$d_S(t)=\underset{𝐱Q_d}{\mathrm{max}}P_𝐱(S(𝐗(t))=)\pi _S()=\underset{𝐱Q_d}{\mathrm{max}}(T_𝐱^1>t).$$

By the definition of $`\tau ^{(S)}`$ it is clear that $`P(T_𝐱^1>\tau ^{(S)})\frac{1}{2e}`$. Moreover, by the “submultiplicativity” property (see Aldous and Fill , chapter 2)

$$d(s+t)2d(s)d(t),s,t0,$$
(3.11)
we have that

$$P(T_𝐱^1>k\tau ^{(S)})\frac{1}{2e^k},k1.$$
(3.12)
Thus

$`𝔼T_𝐱^1`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}𝔼\left(T_𝐱^1\mathrm{𝟏}((k1)\tau ^{(S)}<T_𝐱^1k\tau ^{(S)})\right)`$
$``$ $`\tau ^{(S)}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}kP((k1)\tau ^{(S)}<T_𝐱^1k\tau ^{(S)})`$
$``$ $`\tau ^{(S)}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}kP((k1)\tau ^{(S)}<T_𝐱^1)`$
$``$ $`\tau ^{(S)}\left(1+{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}k{\displaystyle \frac{1}{2e^{k1}}}\right).`$
Since the series in the right hand side converges we have

$$𝔼T_𝐱^1=O(f_S(d)).$$

Once the bivariate $`𝐒`$ process hits the diagonal

$$𝐃=\{(𝐱,𝐲)Q_d\times Q_d;\underset{i=1}{\overset{d}{}}x_i=\underset{i=1}{\overset{d}{}}y_i\},$$
(3.13)
we devise one obvious coupling that forces the bivariate $`𝐗`$ process to stay in $`𝐃`$ and such that the distance defined in (3.5) between the marginal processes does not decrease. In words: we select one coordinate at random; if the marginal processes coincide in that coordinate, we allow them to evolve together; otherwise we select another coordinate in order to force two new coincidences. Formally, for each ($`𝐗(t),𝐘(t))𝐃`$, let $`I_1,I_2`$ and $`I_3`$ be the partition of $`\{0,1,\mathrm{},d\}`$ such that

$`I_1`$ $`=`$ $`\{i;X_i(t)=Y_i(t)\}`$
$`I_2`$ $`=`$ $`\{i;X_i(t)=0,Y_i(t)=1\}`$
$`I_3`$ $`=`$ $`\{i;X_i(t)=1,Y_i(t)=0\}`$
Given $`(𝐗(t),𝐘(t))𝐃`$, choose $`i`$ u.a.r. from $`\{0,1,\mathrm{},d\}`$.

* If $`iI_1`$;

1. If $`X_i(t)=1`$ then make $`X_i(t+1)=Y_i(t+1)=0`$ with probability $`p_{s(𝐗(t))1}`$; otherwise $`X_i(t+1)=Y_i(t+1)=1`$.
    2. If $`X_i(t)=0`$ then make $`X_i(t+1)=Y_i(t+1)=1`$ with probability $`p_{s(𝐗(t))}`$; otherwise $`X_i(t+1)=Y_i(t+1)=0`$.
* If $`iI_2`$;

1. Select $`jI_3`$ u.a.r.;
    2. Make $`X_i(t+1)=Y_j(t+1)=1`$ with probability $`p_{s(𝐗(t))}`$; otherwise $`X_i(t+1)=Y_j(t+1)=0`$.
* If $`iI_3`$;

1. Select $`jI_2`$ u.a.r.;
    2. Make $`X_i(t+1)=Y_j(t+1)=0`$ with probability $`p_{s(𝐗(t))1}`$; otherwise $`X_i(t+1)=Y_j(t+1)=1`$.

Then, it is easy to check that $`(𝐗(t+1),𝐘(t+1))𝐃`$ and $`d(𝐗(t+1),𝐘(t+1))d(𝐗(t),𝐘(t))`$. Moreover, noticing that $`|I_2|=|I_3|=d(𝐗(t),𝐘(t))/2`$, we have

$$\left(d(𝐗(t+1),𝐘(t+1))=d(𝐗(t),𝐘(t))2|𝐗(t),𝐘(t)\right)=\frac{d(𝐗(t),𝐘(t))}{2d}(p_{s(𝐗(t))}+p_{s(𝐗(t))1})$$
(3.14)
$$\left(d(𝐗(t+1),𝐘(t+1))=d(𝐗(t),𝐘(t))|𝐗(t),𝐘(t)\right)=1\frac{d(𝐗(t),𝐘(t))}{2d}(p_{s(𝐗(t))}+p_{s(𝐗(t))1}).$$
(3.15)
In this case, it is straightforward to compute

$`m(i,s(𝐗(t)))`$ $`=`$ $`i𝔼\left[d(𝐗(t+1),𝐘(t+1))|d(𝐗(t),𝐘(t))=i,𝐗(t),𝐘(t)\right]`$ (3.16)
$`=`$ $`{\displaystyle \frac{i}{d}}(p_{s(𝐗(t))}+p_{s(𝐗(t))1}).`$
Let $`T_𝐱^2`$ be the coupling time for the second coupling just described. That is, let $`T_𝐱^2=inf\{t>T_𝐱^1:d(𝐗(t),𝐘(t))=0\}`$. Then, as a consequence of the optional sampling theorem for martingales we have the following comparison lemma (cf. Aldous and Fill , Chapter 2).

###### Lemma 3.17

$$𝔼\left[T_𝐱^2|d(𝐗(T_𝐱^1),𝐘(T_𝐱^1))=L,(𝐗(T_𝐱^1),𝐘(T_𝐱^1))𝐃,s(𝐗(T_𝐱^1))=s\right]\underset{i=1}{\overset{L}{}}\frac{d}{i(p_s+p_{s1})}$$
(3.18)
for all $`s=0,1,\mathrm{},d`$.

Proof. Define $`(𝐗^{}(t),𝐘^{}(t))=(𝐗(t+T_𝐱^1),𝐘(t+T_𝐱^1)`$ for all $`t0`$. Define $`Z(t)=d(𝐗^{}(t),𝐘^{}(t))`$ and $`_t=\sigma (𝐗^{}(t),𝐘^{}(t))`$. Then, it follows from (3.16) that

$$m(i,s)=iE\left[Z(1)|Z(0)=i,s(𝐗^{}(0))=s\right].$$
(3.19)
Also, for all $`s\{1,\mathrm{},d1\}`$, $`0<m(1,s)m(2,s)\mathrm{}m(d,s)`$. Fix $`s\{1,\mathrm{},d1\}`$ and write

$$h(i)=\underset{j=1}{\overset{i}{}}\frac{1}{m(i,s)}$$
(3.20)
and extend $`h`$ by linear interpolation for all real $`0xd`$. Then $`h`$ is concave and for all $`i1`$

$`𝔼\left[h(Z(1))|Z(0)=i,s(𝐗^{}(0))=s\right]`$ $``$ $`h(im(i,s))`$
$``$ $`h(i)m(i,s)h^{}(i)`$
$`=`$ $`h(i)1,`$
where the first inequality follows from the concavity of $`h`$ and $`h^{}`$ is the first derivative of $`h`$. Now, defining $`\stackrel{~}{h}`$ such that

$$h(i)=1+\underset{j}{}\left[h(Z(1))|Z(0)=i,s(𝐗^{}(0))=s\right]h(j)+\stackrel{~}{h}(i)$$
(3.21)
and

$$M(t)=t+h(Z(t))+\underset{u=0}{\overset{t1}{}}\stackrel{~}{h}(Z(u))$$
(3.22)
we have that $`M`$ is an $`_t`$-martingale and applying the optional sampling theorem to the stopping time $`T_0=inf\{t;Z(t)=0\}`$ we have

$$𝔼\left[M(T_0)|Z(0)=i,s(𝐗^{}(0))=s\right]=𝔼\left[M(0)|Z(0)=i,s(𝐗^{}(0))=s\right]=h(i).$$
(3.23)
Noticing that $`M(T_0)T_0`$ and $`T_0=T_𝐱^2`$, we obtain the desired result $``$

Since $`s(𝐗(t))`$ is distributed as $`\pi _S`$, we can write:

$$𝔼\left[T_𝐱^2|d(𝐗(T_𝐱^1),𝐘(T_𝐱^1))=L,(𝐗(T_𝐱^1),𝐘(T_𝐱^1))𝐃\right]\underset{s=0}{\overset{d}{}}\pi _S(s)\underset{i=1}{\overset{L}{}}\frac{d}{i(p_s+p_{s1})}:=g(d).$$
(3.24)
Putting the pieces together, we have found a coupling time $`T_𝐱`$ for the whole process such that

$$𝔼T_𝐱f_S(d)+g(d).$$

The task now is to find explicit bounds for $`f_S(d)`$ and $`g(d)`$ for particular workable cases.

To avoid unnecessary complications, we will assume $`d=2m`$, and compute only the hitting times for the $`𝐒`$ process of the type $`E_0T_m`$. Hitting times in birth-and-death processes assume the following closed-form (see  for an electrical derivation):

$$𝔼_kT_{k+1}=\frac{1}{\pi _kP(k,k+1)}\underset{i=0}{\overset{k}{}}\pi _i,0kd1,$$

and in our case this expression turns into

$$𝔼_kT_{k+1}=\frac{\underset{i=0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{d}{i}\right)}{\left(\genfrac{}{}{0pt}{}{d1}{k}\right)p_k}.$$

Therefore

$$𝔼_0T_m=\underset{k=0}{\overset{m1}{}}\frac{\underset{i=0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{2m}{i}\right)}{\left(\genfrac{}{}{0pt}{}{2m1}{k}\right)p_k}.$$
(3.25)
(i) In case all $`p_k=1`$, we have the simple random walk on the cube, and it turns out there is an even more compact expression of (3.25), namely:

$$𝔼_0T_m=\underset{k=0}{\overset{m1}{}}\frac{\underset{i=0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{2m}{i}\right)}{\left(\genfrac{}{}{0pt}{}{2m1}{k}\right)}=m\underset{k=0}{\overset{m1}{}}\frac{1}{2k+1},$$
(3.26)
as was proved in , and the right hand side of (3.26) equals $`m[H(2m){\displaystyle \frac{1}{2}}H(m)]`$, where $`H(n)=1+\frac{1}{2}+\mathrm{}+\frac{1}{n}`$, allowing us to conclude immediately that in this case $`𝔼_0T_m=𝔼_0T_{d/2}\frac{d}{4}\mathrm{log}d+\frac{d}{4}\mathrm{log}2`$. Also, we have

$$𝔼[T_𝐱^2|d(𝐗(T_𝐱^1),𝐘(T_𝐱^1))=L,(𝐗(T_𝐱^1),𝐘(T_𝐱^1))𝐃]\frac{d}{2p}\underset{i=1}{\overset{L}{}}\frac{1}{i}\frac{d}{2p}O(\mathrm{log}L).$$
(3.27)
Thus in this case both $`f_S(d)`$ and $`g(d)`$, and a fortiori $`𝔼T_𝐱`$ and $`\tau `$, are $`O(d\mathrm{log}d)`$.

(ii) For the Aldous cube, $`p_k=\frac{k}{d1}`$, and (3.25) becomes (recall this cube excludes the origin):

$$𝔼_1T_m=\underset{k=1}{\overset{m1}{}}\frac{\underset{i=0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{2m}{i}\right)}{\left(\genfrac{}{}{0pt}{}{2m2}{k1}\right)}=\underset{k=0}{\overset{m2}{}}\frac{\underset{i=0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{2m}{i}\right)}{\left(\genfrac{}{}{0pt}{}{2m2}{k}\right)}+\underset{k=0}{\overset{m2}{}}\frac{\left(\genfrac{}{}{0pt}{}{2m}{k+1}\right)}{\left(\genfrac{}{}{0pt}{}{2m2}{k}\right)}.$$
(3.28)
After some algebra, it can be shown that the second summand in (3.28) equals $`(2m1)[H(2m1){\displaystyle \frac{1}{m}}]`$, and the first summand can be bounded by twice the expression in (3.26), on account of the fact that $`\left({\displaystyle \genfrac{}{}{0pt}{}{2m1}{k}}\right)2\left({\displaystyle \genfrac{}{}{0pt}{}{2m2}{k}}\right),`$ for $`0km1`$. Therefore, we can write

$$𝔼_1T_{d/2}\frac{3}{2}d\mathrm{log}d+\text{ smaller  terms},$$

thus improving by a factor of $`\frac{1}{2}`$ the computation of the same hitting time in .

Also, we have

$$𝔼[T_𝐱^2|d(𝐗(T_𝐱^1),𝐘(T_𝐱^1))=L,(𝐗(T_𝐱^1),𝐘(T_𝐱^1))𝐃,s(𝐗(T_𝐱^1))=s]\underset{i=1}{\overset{L}{}}\frac{d(d1)}{i(2s1)}$$
(3.29)
Thus, in this case

$`𝔼[T_𝐱^2|d(𝐗(T_𝐱^1),𝐘(T_𝐱^1))=d,(𝐗(T_𝐱^1),𝐘(T_𝐱^1))𝐃]`$ (3.30)
$``$ $`{\displaystyle \underset{s=1}{\overset{d}{}}}\pi _S(s){\displaystyle \underset{i=1}{\overset{d}{}}}{\displaystyle \frac{d(d1)}{i(2s1)}}`$
$``$ $`\mathrm{\Phi }(\sqrt{d}/3){\displaystyle \underset{i=1}{\overset{d}{}}}{\displaystyle \frac{d(d1)}{i}}+(1\mathrm{\Phi }(\sqrt{d}/3)){\displaystyle \underset{i=1}{\overset{d}{}}}{\displaystyle \frac{d(d1)}{i(2d/31)}}`$
$``$ $`e^{d/9}d(d1)\mathrm{log}d+(1e^{d/9})d\mathrm{log}d.`$
And so $`\tau =O(d\mathrm{log}d)`$ also in this case.

(iii) Slower walks. Consider the case when the probability $`p_k`$ grows exponentially in $`k`$, more specifically

$$p_k=\left(\frac{k+1}{n+1}\right)^\alpha $$
(3.31)
with $`\alpha >1`$. In this case, it seems that (3.25) is useless to get a closed expression for $`𝔼_0T_{d/2}`$. However, Graham and Chung  provide the following bound

$$𝔼_iT_{d/2}c_0(\alpha )d^\alpha ,\text{ for all }dd_0(\alpha ),0id$$
(3.32)
where $`c_0(\alpha )`$ and $`d_0(\alpha )`$ are constants depending only on $`\alpha `$. Moreover, (3.24) becomes

$`g(d)`$ $`=`$ $`{\displaystyle \underset{s=0}{\overset{d}{}}}\pi _S(s){\displaystyle \underset{i=1}{\overset{d}{}}}{\displaystyle \frac{d(d+1)^\alpha }{i((s+1)^\alpha +s^\alpha )}}`$ (3.33)
$`=`$ $`d(d+1)^\alpha {\displaystyle \underset{s=0}{\overset{d}{}}}{\displaystyle \frac{\pi _S(s)}{(s+1)^\alpha +s^\alpha }}{\displaystyle \underset{i=1}{\overset{d}{}}}{\displaystyle \frac{1}{i}}`$
$``$ $`d(d+1)^\alpha \mathrm{log}d{\displaystyle \underset{s=0}{\overset{d}{}}}{\displaystyle \frac{\pi _S(s)}{(s+1)^\alpha +s^\alpha }}`$
$``$ $`d(d+1)^\alpha \mathrm{log}d{\displaystyle \underset{s=0}{\overset{d}{}}}{\displaystyle \frac{\pi _S(s)}{(s+1)^\alpha }}`$
$`=`$ $`d(d+1)^\alpha \mathrm{log}d𝔼\left[{\displaystyle \frac{1}{(1+X)^\alpha }}\right],`$
where $`X`$ is a Binomial ($`d,\frac{1}{2}`$) random variable. Jensen’s inequality and the same argument that lead to (3.30) show that $`𝔼[(1+X)^\alpha ]O(d^\alpha )`$ and (3.33) can be bound by $`O(d\mathrm{log}d)`$. This fact together with (3.32) allows us to conclude that $`\tau =O(d^\alpha )`$ in this case, thus improving on the rate of the mixing time provided in  by a factor of $`\mathrm{log}d`$.

Acknowledgments.  This paper was initiated when both authors were visiting scholars at the Statistics Department of the University of California at Berkeley. The authors wish to thank their hosts, especially Dr. David Aldous with whom they shared many fruitful discussions. This work was partially supported by FAPESP 99/00260-3.


text:
# 1 Introduction

## 1 Introduction

Various observational facts, notably the dynamics of spiral galaxies, clusters and superclusters, and their theoretical explanations in modern cosmology and astrophysics can be reconciled only by postulating the existence of the non-luminous matter at amounts exceeding the total matter density provided by baryons. Thus, if gravity does not undergo some drastic changes at distances larger than few kpc (a highly improbable option from various points of view), the explanation of the rotational curves of spiral galaxies requires the existence of halos, composed of dark matter objects.

One of the most intriguing options from the point of view of particle physics is the possibility that dark matter is composed of massive particles of yet unknown identity. On the one hand these (presumably) weakly interacting massive particles (WIMPs) are the subject of intense experimental searches at numerous underground detectors . On the other hand, some of the extensions of the Standard Model naturally incorporate such stable particles. For example, in supersymmetric extensions of SM the lightest supersymmetric particle can be stable and neutral, thus providing a viable dark matter candidate. Extensive theoretical studies of stable SUSY particles, in particular neutralinos, have been performed over the last two decades  for wide ranging values of the supersymmetric masses and couplings. We believe, however, that unless supersymmetry is experimentally verified and the stability of the lightest superpartner is established in collider experiments, other possibilities for WIMPs should not be discarded.

In this letter we adopt a completely different, yet fully justified, phenomenological approach to the physics of WIMPs. We do not make any particular assumptions about the origin and possible identity of these particles, apart from the assumption that their masses are larger than a few GeV. We then assess what can be inferred about the properties of WIMPs directly from experimental limits implied by various searches. An additional incentive to address WIMP physics in a model independent way stems from recently discovered problems with subgalactic structure formation within the non-interacting cold dark matter scenario . A “generalized” WIMP with sufficient self-interaction may avoid these problems .

Among important low-energy properties of WIMPs are their electromagnetic form factors. It has been known for a long time that the possibility of charged WIMPs is strongly disfavored  and stringent limits exist in the case of a fractional charge . Here we assume that WIMPs are neutral and discuss other electro-magnetic form factors such as magnetic and electric dipole moments, anapole moment and polarizability. We put strong constraints on the possible electromagnetic form factors of WIMPs using the result of direct experimental searches aimed at the detection of the WIMP-nucleus recoil signal. These constraints depend on the number density and average velocities of WIMPs in the halo. The limits on fluxes of energetic neutrinos, originating from the annihilation of WIMPs in the center of the Sun and the Earth, can also be used to obtain indirect constraints on the form factors. These constraints are less certain, as they require additional assumptions about annihilation rates and neutrino branching ratios.

This paper is organized as follows. The calculations of the WIMP-nuclei cross sections are performed in the next section. In section 3 we present direct limits on the form factors, using the results of various underground experiments. In section 4 we obtain indirect limits on the form factors, using limits on neutrino fluxes. Our conclusions are reserved for section 5.

## 2 Electromagnetic form factors of WIMPs

The interaction of any compact object and a slowly varying electro-magnetic field can be parametrized via a set of electro-magnetic form factors. In our case, the object is a WIMP with spin $`S`$, and its interactions, linear in external fields $`𝐁`$ and $`𝐄`$, can be written as

$`H=\mu 𝐁{\displaystyle \frac{𝐒}{S}}d𝐄{\displaystyle \frac{𝐒}{S}}a𝐣{\displaystyle \frac{𝐒}{S}}`$ (1)
$`{\displaystyle \frac{1}{4S(2S1)}}[S_iS_j+S_jS_i{\displaystyle \frac{2}{3}}\delta _{ij}S(S+1)]\left(Q{\displaystyle \frac{}{x_i}}E_j+M{\displaystyle \frac{}{x_i}}B_j\right)+\mathrm{}.`$
Here $`\mu `$ and $`d`$ are magnetic and electric dipole moments, and $`M`$ and $`Q`$ are magnetic and electric quadrupole moments. The anapole moment $`a`$ is the form factor which describes the contact interaction with the external current density $`𝐣`$ .

If the spin of the WIMP is zero, these moments do not exist and the interaction with electromagnetic field is given by the charge radius $`r_D`$ of the WIMP and the polarizabilities

$$H=\frac{1}{6}er_D^2\frac{}{x_i}E_i\frac{1}{2}\chi _EE^2\frac{1}{2}\chi _BB^2\chi _{EB}𝐄𝐁+\mathrm{}$$
(2)
If time reversal is a good symmetry, then $`d`$, $`M`$ and $`\chi _{EB}`$ necessarily vanish. Conserved parity also forbids the existence of $`a`$. Here we do not assume that these symmetries are preserved a priori, and we allow all of the form factors to exist.

In what follows we calculate the cross section of the elastic scattering of a WIMP off a nucleus due to the existence of the electromagnetic form factors. We start with the scattering due to the magnetic moment which is arguably the most interesting case.

The calculation can be done in the non-relativistic limit, as $`v/c10^3`$ on average for WIMPs in the halo. We assume that the possible size of the form factors is of the order of, or smaller than, the respective characteristic values for nuclei, $`\mu <\mu _n`$, $`Q<e\mathrm{fm}^3`$, and so on. It can then easily be checked that this range of velocities also ensures the applicability of the Born approximation. In our units, $`\mathrm{}=1`$, $`c=1`$, $`e^2=\alpha =1/137`$ and $`\mu _n=e/2m_p`$.

The cross-section for the scattering of a Wimp through the interaction of its magnetic moment with the magnetic moment of the nucleus is

$$\sigma =\frac{32\pi }{3}M_R^2(N,D)\mu _D^2\mu _N^2\frac{\left(S+1\right)\left(I+1\right)}{3SI},$$
(3)
where $`M_R(N,D)`$ is the reduced mass of the WIMP-nucleus system, $`\mu _D`$ and $`\mu _N`$ are the magnetic moments of the WIMP and the nucleus, respectively, and $`I`$ is the spin of the nucleus. The cross-section for the interaction between the electric dipole moment of the WIMP and the charge of the nucleus is

$$\sigma =8\pi Z^2\left(\frac{d}{e}\right)^2\left(\frac{\alpha }{v}\right)^2\frac{S+1}{3S}\mathrm{ln}\frac{q_{max}}{q_{min}}.$$
(4)
It is logarithmically divergent at large distances. In practice, $`q_{max}`$ can be chosen to be of the order of $`M_R(N,D)v`$ and $`q_{min}`$ is determined either by the inverse of the impact parameters at which the nuclear charge is screened by the electrons or by the lowest momentum transfer that can be detected in the experiment, whichever is larger. The cross-section for the interaction between the nuclear electric current and the anapole of the WIMP is given by

$$\sigma =\frac{1}{2\pi }Z^2\alpha ^2\frac{M_R^2(N,D)}{4m_p^2}\left(\frac{a}{\mu _n}\right)^2v^2\frac{S+1}{3S},$$
(5)
where $`\mu _n`$ is the nuclear magneton, and $`m_p`$ is the proton mass. The contribution of the nuclear spin current is negligibly small. The cross-section for scattering through the interaction between the electric quadrupole moment of the WIMP and the electric charge of the nucleus takes the form

$$\sigma =\frac{8\pi }{9}M_R^2(N,D)Z^2\alpha ^2\left(\frac{Q}{e}\right)^2\frac{\left(S+1\right)\left(2S+3\right)}{10S\left(2S1\right)}.$$
(6)
In Eqs. (4)-(6) the normalization of the WIMP spin dependence is such that it is equal to 1 for the lowest spin for which these moments may exist, $`(S+1)/3S=1`$ for $`S=1/2`$ and $`(S+1)(2S+3)/10S(2S1)=1`$ for $`S=1`$. The cross-section for scattering of a WIMP because of its charge radius in the electric field of the nucleus is given by

$$\sigma =\frac{4\pi }{9}M_R^2(N,D)Z^2\alpha ^2r_D^4.$$
(7)
Finally, the cross-section for the scattering of a WIMP through its polarizability in the electric field of the nucleus, is

$$\sigma \frac{144\pi }{25}M_R^2(N,D)Z^4\alpha ^2\frac{\chi _E^2}{r_0^2}.$$
(8)
This cross-section diverges for a point-like nucleus. We therefore took the charge distribution of the nucleus to be homogeneous in a sphere with radius $`r_01.2\mathrm{fm}\sqrt[3]{A}`$.

Scattering due to the magnetic quadrupole moments and magnetic and mixed polarizabilities, $`\chi _{EB}`$, was also considered. These effects turned out to be additionally suppressed and do not lead to interesting limits. We therefore do not quote them here.

## 3 Direct constraints

We illustrate the bounds on the electro-magnetic form factors of WIMPs which are implied by the limit on the recoil signal with results of the NaI detector of the DAMA experiment , and Ge detector of the CDMS collaboration . Other experiments aimed at the direct detection of dark matter can be used as well  to get limits which would typically be a few times weaker (see also Ref.  for more complete listing of experiments). Traditionally, the results of the direct searches are presented in the form of a normalized bound on the “WIMP–proton”cross-section in the spin dependent case, and a normalized “WIMP–nucleon”cross-section in the spin independent case. Such normalized bounds are useful to compare results for different types of nuclei, to compare results from different experiments, and to compare theoretical predictions with experimental data, but for our purpose of determining bounds on the WIMP electro-magnetic form factors, it is necessary to reverse this normalization. The bounds, quoted in  are, in some sense, specifically formulated for neutralino-nucleon scattering, which necessarily implies a nuclear theory input when switching from nucleus-neutralino to nucleon-neutralino description. It is clear that in the case of the scattering due to the WIMP electro-magnetic form factors nuclear theory is not involved, and much of the nuclear physics uncertainty is not there. Thus, in order to obtain bounds on the form factors, we recover the nucleus-WIMP cross section from the nucleon-WIMP by the use of the equation

$$\sigma =\frac{M_R^2(N,D)}{M_R^2(p,D)}\frac{I(N)}{I(p,n)}\sigma _{p,n}.$$
(9)
Here $`p`$ and $`n`$ stand for proton and neutron, respectively, and the ratio of the reduced masses takes into account a kinematic factor. In the spin dependent case, $`I(N)/I(p)`$ is a spin factor, which takes the values $`I(N)/I(p)=0.009`$ for Iodine, and $`I(N)/I(p)=0.055`$ for Sodium . In the spin independent case, $`I(N)/I(n)=(AZ)^2`$ for Ref.  and $`I(N)/I(n)=A^2`$ for Ref. . For simplicity, we set all nuclear form factors equal to one, which is a better approximation for Sodium than for Iodine.111For a more rigorous treatment electric and magnetic form factors of nuclei should be included, for many of which the experimental results are available. A more accurate analysis could be done by the experimental collaborations using their original data. The ensuing bounds on the WIMP electro-magnetic form factors are listed in Table(1).

The bounds presented in this table are obtained with the standard assumptions that the mass density of WIMPs in our halo is 0.3 Gev/cm3 and that their velocities are of the order of $`200`$ km/sec. The data obtained with germanium detectors are sensitive to all electro-magnetic form factors discussed here except the magnetic dipole moment, as the most abundant germanium isotope has no spin. Let us assume that the recently reported annual modulation of the DAMA signal, consistent with $`m_D50`$ GeV , and the absence of the signal, resulting in strong limits on spin-independent cross section, quoted by CDMS , are both correct. Then the presence of a WIMP magnetic moment at the level of few$`\times \mathrm{\hspace{0.17em}10}^5`$$`1\times 10^4`$ compared to the size of the nuclear magneton can reconcile the two experimental results.

It is instructive to consider how these bounds are relaxed in the scaling regime, when the WIMP mass is larger than nuclear masses and their abundance is kept as a free parameter. To this end we introduce the parameter $`\delta _D`$, characterizing the fraction of WIMPs of species $`D`$ in the halo dark matter density ($`\delta _D=1`$ in Table 1). The counting rate in a detector is proportional to the number density of the WIMP particles and therefore should scale with $`M_D`$ and $`\delta _D`$ according to the following formula:

$$R_{N,D}(\delta _D,M_D)=R_N^0\delta _D\left(\frac{M_D^0}{M_D}\right),$$
(10)
where $`M_D^0`$ is the lowest mass of WIMPs at which the scaling sets in and $`R_N^0`$ is the counting rate, corresponding to $`M_D=M_D^0`$ and $`\delta _D=1`$. It is easy to see from the figures in Refs. that $`M_D=100`$ GeV is close to the mass where the scaling behaviour begins , and therefore we can simply rescale the bounds given in Table 1 for an arbitrary $`\delta _D`$ and arbitrary and heavy $`M_D`$, by multiplying them by $`\sqrt{\delta _DM_D/100\mathrm{G}\mathrm{e}\mathrm{V}}`$.

## 4 Indirect constraints

The elastic scattering of WIMPs in the solar (and earth) media leads to gravitational trapping. Due to the scattering, WIMPs dissipate energy and eventually they accumulate in the center of the sun (and earth). Their density grows until the accumulation rate is counter-balanced by the pair annihilation rate, and an equilibrium is reached. In the process of annihilation, among other SM products, energetic neutrinos can be produced. The interaction of these neutrinos in the earth in turn leads to the production of energetic muons, which can be detected with neutrino telescopes.

Thus limits on the flux of energetic neutrinos can be used to set bounds on the electro-magnetic form factors of WIMPs. For WIMPS interacting through their electro-magnetic form factors, the trapping rate in the center of the sun and earth is determined from the elastic WIMP-nucleus cross sections that we calculated in Eqs.(3)-(8).

We use the results of Ref. , which follows the original treatment in Refs. , and estimate the neutrino flux at the surface of the earth due to WIMP annihilation in the sun as

$$\varphi _\nu (560cm^2s^1)N_{eff}\delta _D\sigma _{p,36}\frac{\text{GeV}^2}{M_D^2}.$$
(11)
In this formula $`N_{eff}`$ is the average number of neutrinos per annihilation event and $`\sigma _{p,36}`$ is the WIMP-proton elastic cross section in units of $`10^{36}`$ cm2. It is assumed that the equilibrium number of WIMPs in the sun has been reached, so that the annihilation and capture rates are equal.

The experimental limit on the flux of energetic upward muons (for example obtained by the Kamioka, MACRO and Baksan experiments, ), is $`\varphi _\mu <1.4\times 10^{14}`$ cm-2s-1, and the probability that a neutrino directed towards the detector produces a muon at the detector is $`P(100\mathrm{G}\mathrm{e}\mathrm{V})10^7`$ . Using Eq.(11), we deduce the following indirect limits on the electro-magnetic form factors of WIMPs with $`m_D=100`$ GeV,

$`\left|\mu _D/\mu _n\right|\sqrt{\left(S+1\right)/3S}`$ $`<`$ $`510^6`$
$`\left|d/e\right|\sqrt{\left(S+1\right)/3S}`$ $`<`$ $`210^{22}\mathrm{cm}`$
$`\left|Q/e\right|\sqrt{\left(S+1\right)\left(2S+3\right)/10S\left(2S1\right)}`$ $`<`$ $`10^8\mathrm{fm}^2`$ (12)
$`\left|a/\mu _n\right|\sqrt{\left(S+1\right)/3S}`$ $`<`$ $`10^2\mathrm{fm}`$
$`r_D^2`$ $`<`$ $`1.410^8\mathrm{fm}^2`$
$`\left|\chi \right|`$ $`<`$ $`5.610^7\mathrm{fm}^3,`$
where $`N_{eff}`$ is taken to be 0.1 and $`\delta _D=1`$. Note that for higher WIMP masses the absorption of neutrinos in the core of the sun becomes important .

If the WIMP-nucleus elastic scattering cross-sections in the earth are comparable to the WIMP-proton cross-section, than the neutrino signal from the center of the earth is approximately four orders of magnitude weaker than the signal from the sun . However, the limit on $`\chi _E`$ from the neutrino flux from the center of the earth is almost one order of magnitude better than the limit derived from the neutrino flux from the sun, as the corresponding WIMP-nucleus elastic scattering cross-sections are enhanced by a large factor, $`Z^4510^5`$, for heavy (Fe) nuclei (see Eq.(8)) comprising the core of the earth.

Although indirect bounds are definitely stronger, they do require a number of additional assumptions. In particular, these bounds require that annihilation actually occurs, which in practice may not be true. For example, if a WIMP caries a conserved quantum number, and there are no “anti-WIMPs” - such as heavy Dirac neutrinos in the absence of antineutrinos - the annihilation is forbidden.

## 5 Discussion

We have shown that limits on recoil signals from underground detectors can be translated into bounds on the electro-magnetic form factors of WIMPs. Similar bounds can be inferred from the limits on energetic neutrino fluxes from the center of the earth and the sun, generated by annihilating WIMPs. Our analysis shows that all electromagnetic form factors have to be very small when expressed in characteristic nuclear units. In particular, the magnetic moment of a WIMP has to be five orders of magnitude smaller than the magnetic moment of a proton.

Unfortunately, all quoted bounds are trivially satisfied for the case of neutralino dark matter. Indeed, because it is a spin-1/2 particle, a neutralino cannot possess a quadrupole moment. In addition, because it is a Majorana fermion, a neutralino cannot have a magnetic or electric dipole moment either ($`\overline{\chi }\sigma _{\mu \nu }\chi 0`$). Finally, an anapole moment and a polarizability can be generated only radiatively. This results only in a small correction to the dominant neutralino-nucleus elastic cross section induced by squark exchange .

A magnetic moment can arise naturally if WIMPs are Dirac fermions. If we assume for a moment that the WIMP is a heavy Dirac neutrino, a large magnetic moment can be induced due to virtual $`W`$-charged lepton exchange. The value of the induced magnetic moment grows quadratically with the size of the Yukawa coupling until the point where the Yukawa interaction becomes essentially strong. In latter case a natural estimate for the size of the magnetic moment is $`e/m_\nu `$, with $`m_\nu `$ presumably of the order of 1 TeV or heavier. In this case, the interaction through the magnetic moment can be very important and, for light nuclei, can even dominate $`Z`$-boson exchange. It is interesting to note that the electric dipole moment of a WIMP is also very constrained. The limit of $`\mathrm{5\hspace{0.17em}10}^{21}`$ e cm is stronger than any other EDM constraint, apart from EDMs of electrons, neutrons and heavy atoms. The method we used to obtain this limit is in the spirit of the original paper by Purcell and Ramsey , where the first bound on the neutron EDM was inferred from neutron-nucleus scattering.

A model independent approach to the physics of WIMPs is partly motivated by recently discovered problems within the cold dark matter model. While large-scale structure can be described reasonably well, the numerical simulations of galactic substructures appear to be in conflict with observational data . A recently proposed cure, invoking self-interacting dark matter , is incompatible with the very restrictive neutralino model. This fact gives further impetus to our approach to study WIMPs in a manner as model independent as possible. The limits on the electro-magnetic form factors obtained in this work show that the amount of self-interaction, which could be induced by these form factors, is not sufficient to generate a cross section of WIMP-WIMP scattering at the level of $`10^{25}`$ cm2 as required by the self-interacting WIMP scenario.

The authors would like to thank Keith Olive, Arkady Vainshtein and Viktor Zacek for usefull discussions. This work was supported in part by the Department of Energy under Grant No. DE-FG-02-94-ER-40823.