format: text

text:
# Three-Loop Radiative-Recoil Corrections to Hyperfine Splitting in Muonium

## I Introduction

Recently we initiated a program of calculating of all three-loop radiative-recoil corrections to hyperfine splitting (HFS) in muonium . Three-loop radiative-recoil corrections are enhanced by the presence of the cube of the large logarithm of the electron-muon mass ratio . All leading logarithm cubed and logarithm squared contributions of this order were calculated a long time ago  (see also reviews in ). As the first step of our program we obtained in  previously unknown single-logarithmic and nonlogarithmic radiative-recoil corrections of order $`\alpha ^2(Z\alpha )(m/M)\stackrel{~}{E}_F`$***We define the Fermi energy as
$$\stackrel{~}{E}_F=\frac{16}{3}Z^4\alpha ^2\frac{m}{M}\left(\frac{m_r}{m}\right)^3chR_{\mathrm{}},$$
(1) where $`m`$ and $`M`$ are the electron and muon masses, $`\alpha `$ is the fine structure constant, $`c`$ is the velocity of light, $`h`$ is the Planck constant, $`R_{\mathrm{}}`$ is the Rydberg constant, and $`Z`$ is the nucleus charge in terms of the electron charge ($`Z=1`$ for muonium). The Fermi energy $`\stackrel{~}{E}_F`$ does not include the muon anomalous magnetic moment $`a_\mu `$ which does not factorize in the case of recoil corrections, and should be considered on the same grounds as other corrections to hyperfine splitting. generated by graphs with two-loop polarization insertions (irreducible and reducible) in the two-photon exchange diagrams. As the next logical stage in implementing our program, we present below the calcuulation of all single-logarithmic and nonlogarithmic three-loop radiative-recoil corrections generated by diagrams with one-loop electron and muon polarization insertions in the exchanged photons. There are four gauge invariant sets of such three-loop diagrams, and we calculate all their contributions.

## II Radiative-Recoil Correction of Order $`\alpha (Z\alpha )(m/M)\stackrel{~}{E}_F`$

All four sets of diagrams considered below can be obtained from the two-photon exchange diagrams with the radiative photons in the electron or muon lines by insertions of the one-loop electron or muon polarization operators. As was discussed in , it is sufficient to calculate contributions of these diagrams in the scattering approximation. In the calculations below we use the approach developed earlier for analytic calculation of the two-loop radiative-recoil corrections of orders $`\alpha (Z\alpha )(m/M)\stackrel{~}{E}_F`$ and $`(Z^2\alpha )(Z\alpha )(m/M)\stackrel{~}{E}_F`$ in  (these corrections were also calculated numerically in ). To make this paper self-contained we first briefly remind the reader of the main steps in calculation of the corrections induced by the radiative photon insertions in the electron line.

The integral representation for the radiative corrections of order $`\alpha (Z\alpha )(m/M)\stackrel{~}{E}_F`$ generated by the graphs with radiative insertions in the electron line in Fig. 1 has the form

$$\delta E^{eline}=\alpha (Z\alpha )\stackrel{~}{E}_F\frac{1}{8\pi ^2\mu }_0^1𝑑x_0^x𝑑y\frac{d^4k}{i\pi ^2(k^2+i0)^2}$$
(2)
$`\left({\displaystyle \frac{1}{k^2+\mu ^1k_0+i0}}+{\displaystyle \frac{1}{k^2\mu ^1k_0+i0}}\right)`$

$`\{(3k_0^22𝐤^2)[{\displaystyle \frac{c_1𝐤^2+c_2(k^2)^2}{(k^2+2bk_0+a^2)^3}}+{\displaystyle \frac{c_3k^2+c_42k_0}{(k^2+2bk_0+a^2)^2}}]`$

$`3k_0[{\displaystyle \frac{c_5k^2+c_6k^22k_0}{(k^2+2bk_0+a^2)^2}}+{\displaystyle \frac{c_7k^2}{k^2+2bk_0+a^2}}]\}{\displaystyle }_1^7\delta E^{eline}_i,`$

where the dimensionless integration momentum $`k`$ is measured in units of the electron mass, the small parameter $`\mu `$ is defined as half the ratio of the electron and muon masses ($`\mu =m/(2M))`$, the auxiliary functions of the Feynman parameters $`a(x,y)`$ and $`b(x,y)`$ are defined by the relationships

$$a^2=\frac{x^2}{y(1y)},b=\frac{1x}{1y},$$
(3)
and explicit expressions for the coefficient functions $`c_i`$ are collected in Table I. We used the Yennie gauge for the radiative photons in the derivation of the explicit expressions for these coefficient functions $`c_i`$.

After the Wick rotation and transition to four-dimensional spherical coordinates, the expression in eq.(2) acquires the form

$$\delta E^{eline}=\frac{\alpha (Z\alpha )}{\pi ^2}\frac{m}{M}\stackrel{~}{E}_F\frac{1}{8\pi \mu ^2}_0^\pi 𝑑\theta \mathrm{sin}^2\theta _0^{\mathrm{}}𝑑k^2𝒟(k,\theta )$$
(4)
$`{\displaystyle \frac{k^2}{k^2+\mu ^2\mathrm{cos}^2\theta }}{\displaystyle \frac{1}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }},`$

where

$$𝒟(k,\theta )=_0^1dx_0^xdy\{(2+\mathrm{cos}^2\theta )[(c_1\mathrm{sin}^2\theta +c_2k^2)\left(\frac{}{a^2}\right)^2+2c_3\frac{}{a^2}](k^2+a^2)$$
(5)
$`8bc_4(2+\mathrm{cos}^2\theta )\mathrm{cos}^2\theta {\displaystyle \frac{}{a^2}}+12b\mathrm{cos}^2\theta (c_5{\displaystyle \frac{}{a^2}}c_7)12c_6\mathrm{cos}^2\theta {\displaystyle \frac{}{a^2}}(k^2+a^2)\}.`$

We introduced derivatives with respect to $`a^2`$ in order to reduce the powers in the denominators before integration over angles.

Next we separate the $`\mu `$-dependent and $`\mu `$-independent terms in the integrand with the help of the identity

$$\frac{1}{(k^2+\mu ^2\mathrm{cos}^2\theta )[(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta ]}$$
(6)
$`={\displaystyle \frac{\mu ^2}{(k^2+a^2)^24\mu ^2b^2k^4}}\left[{\displaystyle \frac{1}{\mu ^2k^2+\mathrm{cos}^2\theta }}{\displaystyle \frac{4b^2k^2}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }}\right]`$

$`{\displaystyle \frac{\mu ^2}{(k^2+a^2)^2}}\left[{\displaystyle \frac{1}{\mu ^2k^2+\mathrm{cos}^2\theta }}{\displaystyle \frac{4b^2k^2}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }}\right],`$

where, at the last step, we omitted the term proportional to $`\mu ^2`$ in the denominator before the square bracket since, as can be shown, this term generates recoil corrections which are at least quadratic in the recoil factor $`\mu `$. Then the integral for the radiative corrections in eq.(4) becomes a sum of $`\mu `$-dependent and $`\mu `$-independent integrals

$$\delta E^{eline}\frac{\alpha (Z\alpha )}{\pi ^2}\frac{m}{M}\stackrel{~}{E}_F\frac{1}{8\pi }_0^\pi 𝑑\theta \mathrm{sin}^2\theta _0^{\mathrm{}}𝑑k^2𝒟(k,\theta )$$
(7)
$`{\displaystyle \frac{k^2}{(k^2+a^2)^2}}\left[{\displaystyle \frac{1}{\mu ^2k^2+\mathrm{cos}^2\theta }}{\displaystyle \frac{4b^2k^2}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }}\right]=\delta E_{\mu i}^{eline}+\delta E_{ci}^{eline},`$

which we will further call $`\mu `$-integrals and $`c`$-integrals, respectively.

It is convenient to write the angular integrals in eq.(7) in terms of the standard functions $`\mathrm{\Phi }_n^i`$ defined as

$$\mathrm{\Phi }_n(k)\frac{1}{\pi \mu ^2}_0^\pi 𝑑\theta \mathrm{sin}^2\theta \mathrm{cos}^{2n}\theta \frac{(k^2+a^2)^24\mu ^2b^2k^4}{(k^2+\mu ^2\mathrm{cos}^2\theta )[(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta ]}$$
(8)
$`={\displaystyle \frac{1}{\pi }}{\displaystyle _0^\pi }𝑑\theta \mathrm{sin}^2\theta \mathrm{cos}^{2n}\theta \left[{\displaystyle \frac{1}{\mu ^2k^2+\mathrm{cos}^2\theta }}{\displaystyle \frac{4b^2k^2}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }}\right]`$

$`=\mathrm{\Phi }_n^S(k)+\mathrm{\Phi }_n^\mu (k)+\mathrm{\Phi }_n^C(k).`$

Explicit expressions for these functions and their properties are collected in Appendix A (see also ).

The $`c`$-integrals with the functions $`\mathrm{\Phi }_n^C(k)`$ do not contain any free parameters and generate only contributions linear in $`m/M`$ to HFS. These integrals are pure numbers, which admit analytic calculation. The $`\mu `$-integrals with the functions $`\mathrm{\Phi }_n^S(k)`$, $`\mathrm{\Phi }_n^\mu (k)`$ parametrically depend on $`\mu `$ and generate both nonrecoil and recoil contributions to HFS. Contributions of a fixed order in the small mass ratio (which are often enhanced by the large logarithms of the mass ratio) can be extracted from the $`\mu `$-integrals with the help of an auxiliary parameter $`\sigma `$ chosen such that it satisfies the inequality $`1\sigma \mu ^1`$. The parameter $`\sigma `$ is used to separate the momentum integration into two regions, a region of small momenta $`0k\sigma `$, and a region of large momenta $`\sigma k<\mathrm{}`$. In the region of small momenta one uses the condition $`\mu k1`$ to simplify the integrand, and in the region of large momenta the same goal is achieved with the help of the condition $`k1`$. Note that for $`k\sigma `$ both conditions on the integration momenta are valid simultaneously, so in the sum of the low-momenta and high-momenta integrals all $`\sigma `$-dependent terms cancel and one obtains a $`\sigma `$-independent result for the total momentum integral (for a more detailed discussion of this method, see, e.g., ). All contributions, up to and including the corrections quadratic in the small mass ratio $`m/M`$, were analytically calculated earlier in this framework  (see also ).

## III Diagrams with Radiative Photons in the Electron Line and Electron Vacuum Polarization

The general expression for the contribution to HFS arising from the diagrams in Fig. 2The graphs with the crossed exchanged photons are not shown explicitly in this figure and similar figures below. is obtained from the integral in eq.(4) by insertion in the integrand of the doubled one-loop electron polarization $`(\alpha /\pi )k^2I_e(k)`$

$$2\left(\frac{\alpha }{\pi }\right)k^2I_e(k)=2\left(\frac{\alpha }{\pi }\right)k^2_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)},$$
(9)
where the additional multiplicity factor 2 corresponds to the fact that we can insert the vacuum polarization in each of the exchanged photons.

Then all contributions to HFS generated by the diagrams in Fig. 2 are given by the integral

$$\delta E^{ee}=\frac{\alpha ^2(Z\alpha )}{\pi ^3}\frac{m}{M}\stackrel{~}{E}_F\frac{1}{4\pi \mu ^2}_0^\pi 𝑑\theta \mathrm{sin}^2\theta _0^{\mathrm{}}𝑑k^2𝒟(k,\theta )$$
(10)
$`{\displaystyle \frac{k^4}{k^2+\mu ^2\mathrm{cos}^2\theta }}{\displaystyle \frac{1}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }}{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{4+k^2(1v^2)}}.`$

To obtain the radiative-recoil corrections of order $`\alpha ^2(Z\alpha )(m/M)\stackrel{~}{E}_F`$ contained in this integral we follow the route described in the previous section. First we write the integral in eq.(10) as a sum of seven $`\mu `$\- and seven $`c`$-integrals. These integrals are collected in Table II and Table III, where the dimensionless contributions $`\delta ϵ`$ to the energy shifts are defined as $`\delta E_i=\delta ϵ_i\alpha ^2(Z\alpha )(m/M)\stackrel{~}{E}_F/\pi ^3`$.

The $`c`$-integrals in Table III automatically contain only contributions linear in the recoil factor $`m/M`$, and can be immediately calculated numerically. The respective results are again presented in Table III, and the total contribution of all $`c`$-integrals is

$$\delta E_c^{ee}=\delta ϵ_c^{ee}\frac{\alpha ^2(Z\alpha )}{\pi ^3}\frac{m}{M}\stackrel{~}{E}_F=6.96182(3)\frac{\alpha ^2(Z\alpha )}{\pi ^3}\frac{m}{M}\stackrel{~}{E}_F.$$
(11)
The situation with the $`\mu `$-integrals is more complicated. Besides the recoil contributions linear in the small mass ratio $`m/M`$, they contain both nonrecoil contributions and the recoil contributions of higher order in the mass ratio. We would like to remove the already known nonrecoil contributions, to extract in analytic form all coefficients before the logarithmically enhanced terms linear in the recoil factor, and to calculate numerically the nonlogarthmic term which is linear in the recoil factor. To this end, we need to throw away the recoil contributions of higher orders which are also contained in the $`\mu `$-integrals and which are too small from the phenomenological point of view. If preserved, these higher order recoil contributions result only in the loss of accuracy in the numerical integrations. It is easy to see that the integrals with the function $`\mathrm{\Phi }_0^s(k)`$ generate only already known nonrecoil contributions and can be safely omitted for our present goals. We extract the contributions linear (and logarithmically enhanced) in the recoil factor from the $`\mu `$-integrals by separating the integration region with the help of an auxiliary parameter $`\sigma `$ as described in the previous section. All logarithms of the mass ratio originate from the high momentum parts of the $`\mu `$-integrals, where we can use the high energy asymptotic expansion of the polarization operator. Thus we obtain in the analytic form all coefficients before the logarithms of the mass ratio. The radiative-recoil corrections of order $`\alpha (Z\alpha )(m/M)\stackrel{~}{E}_F`$ in eq.(10), generated by the insertions of the radiative photons in the electron line, are linear in the large logarithm of the mass ratio. Hence the corrections of order $`\alpha ^2(Z\alpha )(m/M)\stackrel{~}{E}_F`$ are quadratic in this large logarithm. The coefficient before the logarithm squared, which we obtain in this way, coincides with the one obtained earlier , and the analytically obtained coefficient before the single-logarithmic term, as well as the nonlogarithmic contribution, are new. The results of the calculation of the $`\mu `$-integrals are collected in Tables IV and V.

As an illustration of our methods let us obtain the first three entries in Table IV. We start with the integral in the first line in Table II. It contains a nonrecoil contribution

$$\delta ϵ_{\mu 1}^{ee}(nonrecoil)=_0^1𝑑x_0^x𝑑yc_1_0^{\mathrm{}}𝑑k^2\frac{k^4}{(k^2+a^2)^3}\mathrm{\Phi }_0^s(k)_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}$$
(12)
$`={\displaystyle \frac{1}{\mu }}{\displaystyle _0^1}𝑑x{\displaystyle _0^x}𝑑yc_1{\displaystyle _0^{\mathrm{}}}𝑑k^2{\displaystyle \frac{k^3}{(k^2+a^2)^3}}{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{4+k^2(1v^2)}}.`$

This integral, as well as all other nonrecoil contributions of order $`\alpha ^2(Z\alpha )E_F`$, was calculated analytically in , and numerically in , and we will not consider it here. We write the recoil part of the integral in the first line in Table II as a sum of low-momentum and high-momentum integralsHere, as well as in similar cases below, we allow ourselves a slightly confusing notation, using the same symbol $`\delta ϵ_{\mu 1}^{ee}`$ both for the total contribution of the $`\mu `$-integral and for its recoil part. We hope this will not lead to any misunderstanding since in this paper we are interested only in the recoil corrections.

$$\delta ϵ_{\mu 1}^{ee}=\frac{1}{2}_0^1𝑑x_0^x𝑑yc_1_0^{\mathrm{}}𝑑k^2\frac{k^4}{(k^2+a^2)^3}\left[2\mathrm{\Phi }_0^\mu (k)\mathrm{\Phi }_1^\mu (k)\mathrm{\Phi }_2^\mu (k)\right]_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}$$
(13)
$`=\delta ϵ_{\mu 1}^{ee>}+\delta ϵ_{\mu 1}^{ee<}.`$

First we estimate the high momentum part $`\delta ϵ_{\mu 1}^{ee>}`$ considering the leading term in the integrand in the high momentum limit (more technical details on similar estimates can be found in )

$$\delta ϵ_{\mu 1}^{ee>}=\frac{1}{2}_0^1𝑑x_0^x𝑑yc_1_{\sigma ^2}^{\mathrm{}}𝑑k^2\frac{k^4}{(k^2+a^2)^3}\left[2\mathrm{\Phi }_0^\mu (k)\mathrm{\Phi }_1^\mu (k)\mathrm{\Phi }_2^\mu (k)\right]_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}$$
(14)
$`{\displaystyle _\sigma ^{\mathrm{}}}{\displaystyle \frac{dk^2}{k^4}}\mathrm{ln}^2k^2\left[2\mathrm{\Phi }_0^\mu (k)\mathrm{\Phi }_1^\mu (k)\mathrm{\Phi }_2^\mu (k)\right]{\displaystyle \frac{1}{\sigma ^2}}\mathrm{ln}^2\sigma .`$

We see that the high momentum contribution is suppressed like $`1/\sigma ^2`$. We have already extracted the recoil factor from $`\delta ϵ_{\mu 1}^{ee}`$ explicitly, so now we are looking only for such contributions which are not additionally suppressed. Hence, in the leading order in the recoil parameter we can omit the high-momentum contrbution to $`\delta ϵ_{\mu 1}^{ee}`$, and in our approximation the total contribution to the energy shift is given by the low-momentum integral

$$\delta ϵ_{\mu 1}^{ee<}=\frac{1}{2}_0^1𝑑x_0^x𝑑yc_1_0^{\sigma ^2}𝑑k^2\frac{k^4}{(k^2+a^2)^3}\left[2\mathrm{\Phi }_0^\mu (k)\mathrm{\Phi }_1^\mu (k)\mathrm{\Phi }_2^\mu (k)\right]_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}.$$
(15)
To extract the leading nonvanishing contribution to this integral at $`\mu 0`$ we substitute in the integrand the leading terms in the small $`\mu `$ expansion of the functions $`\mathrm{\Phi }_i^\mu (k)`$ (see Appendix A) and obtain

$$\delta ϵ_{\mu 1}^{ee<}\frac{21}{16}_0^1𝑑x_0^x𝑑yc_1_0^{\sigma ^2}𝑑k^2\frac{k^4}{(k^2+a^2)^3}_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}.$$
(16)
The momentum integral is ultravioletly convergent, so the magnitude of the integral with an infinite upper limit differs from the integral in eq.(16) only by inverse powers of $`\sigma `$, which we omit anyway. Hence, the total $`\mu 1`$ contribution to the energy shift is given by the integral

$$\delta ϵ_{\mu 1}^{ee}=\frac{21}{16}_0^1𝑑x_0^x𝑑yc_1_0^{\mathrm{}}𝑑k^2\frac{k^4}{(k^2+a^2)^3}_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)},$$
(17)
which is just the first entry in Table IV.

Let us turn now to the next integrals in Table II. Each of the integrals $`\delta ϵ_{\mu 2}^{ee}`$ and $`\delta ϵ_{\mu 3}^{ee}`$ generates spurious logarithm cubed contributions which cancel in the sum of these terms. It is therefore convenient to rearrange these terms in such way that spurious logarithm cubed terms do not arise at all. To this end we write the coefficient functions $`c_2`$ in the form

$$c_2=\frac{4(12y)}{y(1y)}x(1+x+2\mathrm{ln}x)+\frac{4}{(1y)^3}[(1x)(12y/x)$$
(18)
$`+2(4+4y/x)\mathrm{ln}x+x(1+x+2\mathrm{ln}x)(45y+2y^2)]c_{2a}+c_{2b}.`$

We have chosen $`c_{2a}`$ essentially as the singular part of the function $`c_2`$ at $`y0`$ and multiplied this singular part by the factor $`(12y)/(1y)`$ in order to simplify integration over $`y`$. Since $`(12y)dy/(y(1y))=da^2/a^2`$ this makes integration over $`y`$ trivial. The spurious logarithm cubed term originates only from the integral of the function $`c_{2a}`$. On the other hand, the integrals over $`x`$ and $`y`$ with the function $`c_{2b}`$ remain finite even if we omit $`a^2`$ in the denominator $`(k^2+a^2)^3`$, which leads to significant simplifications in the high-momentum part of the integral $`\delta ϵ_{\mu 2}^{ee}`$.

Following the same logic, we represent the coefficient functions $`c_3`$ also as a sum of two functions $`c_{3a}`$ and $`c_{3b}`$

$$c_3=\frac{12y}{y(1y)}(16x2x^2)+\frac{1}{(1y)^2}[2(1y)(16x2x^2)+(1y)(8+2x)$$
(19)
$`(1/x)(26(1x)+2x(1x)6(1x^2)y/x6x(1y)+16\mathrm{ln}x)]c_{3a}+c_{3b}.`$

Now the sum of the recoil contributions to HFS generated by the $`\mu 2`$ and $`\mu 3`$ integrals in Table II may be written in the form

$$\delta ϵ_{\mu 2}^{ee}+\delta ϵ_{\mu 3}^{ee}=\delta ϵ_{\mu 2a+\mu 3a}^{ee}+\delta ϵ_{\mu 2b}^{ee}+\delta ϵ_{\mu 3b}^{ee}.$$
(20)
Let us consider calculation of the $`\delta ϵ_{\mu 2a+\mu 3a}^{ee}`$ contribution

$$\delta ϵ_{\mu 2a+\mu 3a}^{ee}=\frac{1}{2}_0^1𝑑x_0^x𝑑y_0^{\mathrm{}}𝑑k^2\left[\frac{c_{2a}k^6}{(k^2+a^2)^3}\frac{c_{3a}k^4}{(k^2+a^2)^2}\right]$$
(21)
$`{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{4+k^2(1v^2)}}\left[2\mathrm{\Phi }_0^\mu (k)+\mathrm{\Phi }_1^\mu (k)\right]=ϵ_{\mu 2a+\mu 3a}^{ee>}+ϵ_{\mu 2a+\mu 3a}^{ee<}.`$

As mentioned above, calculation of the $`y`$ integral becomes trivial after the change of variables $`ya^2`$. In the high-momentum integral we next expand over the inverse powers of $`k^2`$, and use the asymptotic expansions of the arising functions $`V_{lnm}`$ (these functions were defined and their properties were discussed in )

$$ϵ_{\mu 2a+\mu 3a}^{ee>}=\frac{1}{2}_{\sigma ^2}^{\mathrm{}}\frac{dk^2}{k^2}\left[\left(\frac{16}{3}\mathrm{ln}k\frac{2\pi ^2}{3}+10\right)\left(\frac{16}{3}\mathrm{ln}k+\frac{20}{3}\right)\right]\left(\frac{2}{3}\mathrm{ln}k\frac{5}{9}\right)$$
(22)
$`\left[2\mathrm{\Phi }_0^\mu (k)+\mathrm{\Phi }_1^\mu (k)\right]\left({\displaystyle \frac{\pi ^2}{6}}+{\displaystyle \frac{5}{6}}\right)\left\{{\displaystyle \frac{4}{3}}\left[2V_{110}+V_{111}\right]{\displaystyle \frac{10}{9}}\left[2V_{100}+V_{101}\right]\right\}`$

$`\left({\displaystyle \frac{\pi ^2}{3}}{\displaystyle \frac{5}{3}}\right)\mathrm{ln}^2{\displaystyle \frac{M}{m}}+\left({\displaystyle \frac{4\pi ^2}{9}}{\displaystyle \frac{20}{9}}\right)\mathrm{ln}{\displaystyle \frac{M}{m}}+{\displaystyle \frac{\pi ^4}{18}}{\displaystyle \frac{5\pi ^2}{18}}`$

$`+\left({\displaystyle \frac{\pi ^2}{3}}+{\displaystyle \frac{5}{3}}\right)\mathrm{ln}^2\sigma +\left({\displaystyle \frac{5\pi ^2}{9}}{\displaystyle \frac{25}{9}}\right)\mathrm{ln}\sigma .`$

In the low-momentum part of the integral we again preserve only the leading terms in the small $`\mu `$ expansion of the functions $`\mathrm{\Phi }_i^\mu (k)`$

$$ϵ_{\mu 2a+\mu 3a}^{ee<}=\frac{3}{4}_0^1𝑑x_0^x𝑑y_0^{\sigma ^2}𝑑k^2\left[\frac{c_{2a}k^6}{(k^2+a^2)^3}+\frac{c_{3a}k^4}{(k^2+a^2)^2}\right]_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}$$
(23)
and integrate explicitly over $`y`$

$$ϵ_{\mu 2a+\mu 3a}^{ee<}=\left(\frac{\pi ^2}{3}\frac{5}{3}\right)\mathrm{ln}^2\sigma +\left(\frac{5\pi ^2}{9}+\frac{25}{9}\right)\mathrm{ln}\sigma +\frac{14(\pi ^25)}{27}$$
(24)
$`+{\displaystyle \frac{3}{4}}{\displaystyle _0^1}dx{\displaystyle _0^{\mathrm{}}}dk^2\{(1+6x+2x^2)[\mathrm{ln}k^2\mathrm{ln}(k^2+a_1^2){\displaystyle \frac{a_1^2}{k^2+a_1^2}}]`$

$`4x(1+x+2\mathrm{ln}x)[\mathrm{ln}k^2+\mathrm{ln}(k^2+a_1^2)+{\displaystyle \frac{2a_1^2}{k^2+a_1^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{a_1^4}{(k^2+a_1^2)^2}}]\}`$

$`{\displaystyle _0^1}dv{\displaystyle \frac{v^2(1v^2/3)}{4+k^2(1v^2)}}\},`$

where $`a_1^2=x/(1x)`$. Note that we again extended the momentum integration to infinity in the finite integrals for the low-momentum contribution. Finally, summing the high- and low-momentum contributions we obtain the $`\sigma `$-independent result

$$ϵ_{\mu 2a+\mu 3a}^{ee}=\left(\frac{\pi ^2}{3}\frac{5}{3}\right)\mathrm{ln}^2\frac{M}{m}+\left(\frac{4\pi ^2}{9}\frac{20}{9}\right)\mathrm{ln}\frac{M}{m}+\frac{\pi ^4}{18}\frac{5\pi ^2}{18}+\frac{14(\pi ^25)}{27}$$
(25)
$`+{\displaystyle \frac{3}{4}}{\displaystyle _0^1}dx{\displaystyle _0^{\mathrm{}}}dk^2\{(1+6x+2x^2)[\mathrm{ln}k^2\mathrm{ln}(k^2+a_1^2){\displaystyle \frac{a_1^2}{k^2+a_1^2}}]`$

$`4x(1+x+2\mathrm{ln}x)[\mathrm{ln}k^2+\mathrm{ln}(k^2+a_1^2)+{\displaystyle \frac{2a_1^2}{k^2+a_1^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{a_1^4}{(k^2+a_1^2)^2}}]\}`$

$`{\displaystyle _0^1}dv{\displaystyle \frac{v^2(1v^2/3)}{4+k^2(1v^2)}}\},`$

which is convenient for numerical calculations.

As our final example let us derive an expression for $`\delta ϵ_{\mu 2b}^{ee}`$ in the third line in Table IV. We start with the general expression for this contribution from Table II

$$\delta ϵ_{\mu 2b}^{ee}=\frac{1}{2}_0^1𝑑x_0^x𝑑y_0^{\mathrm{}}𝑑k^2\frac{c_{2b}k^6}{(k^2+a^2)^3}_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}\left[2\mathrm{\Phi }_0^\mu (k)+\mathrm{\Phi }_1^\mu (k)\right].$$
(26)
As we already mentioned above in this case we are free to omit $`a^2`$ in comparison with $`k^2`$ in the denominator in the high-momentum part of this integral

$$ϵ_{\mu 2b}^{ee>}=\frac{1}{2}_0^1𝑑x_0^x𝑑y_{\sigma ^2}^{\mathrm{}}𝑑k^2c_{2b}_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}\left[2\mathrm{\Phi }_0^\mu (k)+\mathrm{\Phi }_1^\mu (k)\right]$$
(27)
$`{\displaystyle \frac{1}{2}}{\displaystyle _\sigma ^{\mathrm{}}}{\displaystyle \frac{dk^2}{k^2}}\left({\displaystyle \frac{10\pi ^2}{3}}{\displaystyle \frac{271}{9}}\right)\left({\displaystyle \frac{2}{3}}\mathrm{ln}k{\displaystyle \frac{5}{9}}\right)\left[2\mathrm{\Phi }_0^\mu (k)+\mathrm{\Phi }_1^\mu (k)\right]`$

$`\left({\displaystyle \frac{5\pi ^2}{6}}{\displaystyle \frac{271}{36}}\right)\left\{{\displaystyle \frac{4}{3}}\left[2V_{110}+V_{111}\right]{\displaystyle \frac{10}{9}}\left[2V_{100}+V_{101}\right]\right\}`$

$`\left({\displaystyle \frac{5\pi ^2}{3}}+{\displaystyle \frac{271}{18}}\right)\mathrm{ln}^2{\displaystyle \frac{M}{m}}+\left({\displaystyle \frac{20\pi ^2}{9}}+{\displaystyle \frac{542}{27}}\right)\mathrm{ln}{\displaystyle \frac{M}{m}}{\displaystyle \frac{5\pi ^4}{18}}+{\displaystyle \frac{271\pi ^2}{108}}`$

$`+\left({\displaystyle \frac{5\pi ^2}{3}}{\displaystyle \frac{271}{18}}\right)\mathrm{ln}^2\sigma +\left({\displaystyle \frac{25\pi ^2}{9}}+{\displaystyle \frac{1355}{54}}\right)\mathrm{ln}\sigma .`$

Dealing with the low-momentum part of the integral in eq.(26) we would like to extract analytically leading logarithms of $`\sigma `$. This would allow us to get rid of any trace of the parameter $`\sigma `$ in the final expression for this integral, making it much more suitable for further numerical calculations. The double integral over $`x`$ and $`y`$, with the integrand including the function $`c_{2b}`$, cannot be calculated analytically as easily as the integral with the function $`ϵ_{\mu 2a+\mu 3a}^{ee}`$ above. In order to overcome this difficulty we use the identity

$$\frac{k^4}{(k^2+a^2)^3}=\frac{1}{k^2+4}+\left[\frac{k^4}{(k^2+a^2)^3}\frac{1}{k^2+4}\right]$$
(28)
in the low-momentum part of the integral in eq.(26). Analytic calculation of the integral with the first term on the right hand side (RHS) is simple. On the other hand, all logarithms of $`\sigma `$ are supplied by this integral since the second term on the RHS in eq.(28) decreases at large $`k`$ faster than the first term. Then we obtain the low-momentum contribution to the integral in eq.(26) in the form

$$\delta ϵ_{\mu 2b}^{ee<}=\left(\frac{5\pi ^2}{3}+\frac{271}{18}\right)\mathrm{ln}^2\sigma +\left(\frac{25\pi ^2}{9}\frac{1355}{54}\right)\mathrm{ln}\sigma \frac{1}{54}\left(30\pi ^2271\right)\left(8\mathrm{ln}23\right)$$
(29)
$`{\displaystyle \frac{3}{4}}{\displaystyle _0^1}𝑑x{\displaystyle _0^x}𝑑yc_{2b}{\displaystyle _0^{\mathrm{}}}𝑑k^2\left[{\displaystyle \frac{k^6}{(k^2+a^2)^3}}{\displaystyle \frac{k^2}{k^2+4}}\right]{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{4+k^2(1v^2)}}.`$

Finally, the total expression for the $`\mu 2b`$ integral has the form

$$\delta ϵ_{\mu 2b}^{ee}=\left(\frac{5\pi ^2}{3}+\frac{271}{18}\right)\mathrm{ln}^2\frac{M}{m}+\left(\frac{20\pi ^2}{9}+\frac{542}{27}\right)\mathrm{ln}\frac{M}{m}$$
(30)
$`{\displaystyle \frac{5\pi ^4}{18}}+{\displaystyle \frac{271\pi ^2}{108}}{\displaystyle \frac{1}{54}}\left(30\pi ^2271\right)\left(8\mathrm{ln}23\right)`$

$`{\displaystyle \frac{3}{4}}{\displaystyle _0^1}𝑑x{\displaystyle _0^x}𝑑yc_{2b}{\displaystyle _0^{\mathrm{}}}𝑑k^2\left[{\displaystyle \frac{k^6}{(k^2+a^2)^3}}{\displaystyle \frac{k^2}{k^2+4}}\right]{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{4+k^2(1v^2)}}.`$

Calculation of the remaining $`\mu `$-integrals goes along the similar lines, and does not require any additional new tricks. Calculating the nonlogarithmic contributions numerically and collecting all $`\mu `$-integrals from Table IV, we obtain the total contribution of all $`\mu `$-integrals

$$\delta ϵ_\mu ^{ee}=\frac{5}{2}\mathrm{ln}^2\frac{M}{m}+\frac{22}{3}\mathrm{ln}\frac{M}{m}+4.45606(1).$$
(31)
The sum of $`\mu `$\- and $`c`$-integrals in eq.(31) and eq.(11) gives all radiative-recoil corrections of order $`\alpha ^2(Z\alpha )(m/M)\stackrel{~}{E}_F`$ generated by the diagrams in Fig. 2

$$\delta ϵ^{ee}=\delta ϵ_\mu ^{ee}+\delta ϵ_c^{ee}=\frac{5}{2}\mathrm{ln}^2\frac{M}{m}+\frac{22}{3}\mathrm{ln}\frac{M}{m}+11.41788(3).$$
(32)
## IV Diagrams with Radiative Photons in the Electron Line and Muon Vacuum Polarization

Let us consider now the diagrams in Fig. 3. The only difference between these diagrams and the diagrams in Fig. 2 from the previous section is that they contain the muon vacuum polarization insertion§§§We ascribe an extra factor $`Z`$ to each photon emission by the heavy line, but being somewhat inconsequential, do not write any $`Z`$ factor in the muon vacuum polarization.

$$2\left(\frac{\alpha }{\pi }\right)k^2I_\mu (k)=2\left(\frac{\alpha }{\pi }\right)k^2_0^1𝑑v\frac{v^2(1v^2/3)}{\mu ^2+k^2(1v^2)},$$
(33)
instead of the electron vacuum polarization, and the respective contribution to HFS has the form (compare eq.(10))

$$\delta E^{e\mu }=\frac{\alpha ^2(Z\alpha )}{\pi ^3}\frac{m}{M}\stackrel{~}{E}_F\frac{1}{4\pi \mu ^2}_0^\pi 𝑑\theta \mathrm{sin}^2\theta _0^{\mathrm{}}𝑑k^2𝒟(k,\theta )$$
(34)
$`{\displaystyle \frac{k^4}{k^2+\mu ^2\mathrm{cos}^2\theta }}{\displaystyle \frac{1}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }}{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{\mu ^2+k^2(1v^2)}}.`$

We transform the integrand in the sum of $`\mu `$\- and $`c`$-integrals exactly as in the previous section, the only difference being that each integrand now contains as a factor the muon vacuum polarization insertion from eq.(33) instead of the electron vacuum polarization insertion. This immediately leads to significant simplification of further calculations in comparison to the case of the electron polarization. The muon polarization insertion in eq.(34) is suppressed as $`k^2\mu ^2`$ at integration momenta much smaller than the muon mass. But the characteristic integration momenta in the $`c`$-integrals are of order one in our dimensionless units, and hence now all $`c`$-integralsdo not generate corrections linear in the recoil factor $`m/M`$. Low-momenta parts of the $`\mu `$-integrals, where integration goes over momenta $`k\sigma 1/\mu `$, are also suppressed as $`\sigma ^2\mu ^2`$. Only the high-momenta parts of the $`\mu `$-integrals, where integration effectively goes over the momenta comparable to $`1/\mu `$, generate contributions linear in $`m/M=2\mu `$. Such contributions are present only in the $`\mu `$-integrals connected with the coefficient functions $`c_2`$, $`c_3`$, $`c_6`$, and $`c_7`$. Note that the integrals with the coefficient functions $`c_2`$ and $`c_3`$ contain the function $`\mathrm{\Phi }_0^s(k)=1/(\mu k)`$ in the integrand. In the case of the electron vacuum polarization (and absence of any polarization at all) the integrals with this function generated nonrecoil corrections only. This essentially happened because the characteristic integration momenta were of order one, and in this region the contribution of the function $`1/(\mu k)`$ is obviously enhanced as $`1/\mu `$. In the present case, characteristic integration momenta are about $`1/(\mu k)`$, and the integrals with the function $`\mathrm{\Phi }_0^s(k)`$ generate recoil corrections on a par with the integrals with the other functions $`\mathrm{\Phi }(k)`$.

Integration over the Feynman parameters $`x`$ and $`y`$ in the high-momenta integrals simplifies and may be performed analytically, and we obtain all one-dimensional momenta integrals which generate radiative-recoil corrections in the case of muon polarization insertions

$$\delta ϵ_{\mu 2}^{e\mu >}2_{\sigma ^2}^{\mathrm{}}𝑑k^2\left[\left(\frac{4}{3}\mathrm{ln}k\frac{\pi ^2}{6}+\frac{5}{2}\right)+\left(\frac{5\pi ^2}{6}\frac{271}{36}\right)\right]\left[2\mathrm{\Phi }_0^s(k)+2\mathrm{\Phi }_0^\mu (k)+\mathrm{\Phi }_1^\mu (k)\right]$$
(35)
$`{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{\mu ^2+k^2(1v^2)}},`$

$$\delta ϵ_{\mu 3}^{e\mu >}2_{\sigma ^2}^{\mathrm{}}𝑑k^2\left[\left(\frac{4}{3}\mathrm{ln}k\frac{5}{3}\right)+\left(\frac{2\pi ^2}{3}+\frac{89}{18}\right)\right]\left[2\mathrm{\Phi }_0^s(k)+2\mathrm{\Phi }_0^\mu (k)+\mathrm{\Phi }_1^\mu (k)\right]$$
(36)
$`{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{\mu ^2+k^2(1v^2)}},`$

$$\delta ϵ_{\mu 6}^{e\mu >}\left(\pi ^2\frac{21}{2}\right)_\sigma ^{\mathrm{}}𝑑k^2\mathrm{\Phi }_1^\mu (k)_0^1𝑑v\frac{v^2(1v^2/3)}{\mu ^2+k^2(1v^2)},$$
(37)
$$\delta ϵ_{\mu 7}^{e\mu >}\left(\pi ^2+\frac{15}{2}\right)_\sigma ^{\mathrm{}}𝑑k^2\mathrm{\Phi }_1^\mu (k)_0^1𝑑v\frac{v^2(1v^2/3)}{\mu ^2+k^2(1v^2)}.$$
(38)
The sum of all these contributions has the form

$$\delta ϵ_\mu ^{e\mu >}_{\sigma ^2}^{\mathrm{}}𝑑k^2\left\{\frac{7}{2}\left[2\mathrm{\Phi }_0^s(k)+2\mathrm{\Phi }_0^\mu (k)+\mathrm{\Phi }_1^\mu (k)\right]+3\mathrm{\Phi }_1^\mu (k)\right\}_0^1𝑑v\frac{v^2(1v^2/3)}{\mu ^2+k^2(1v^2)},$$
(39)
or explicitly

$$\delta ϵ_\mu ^{e\mu >}_{\sigma ^2}^{\mathrm{}}dk^2\{\frac{7}{2}[\frac{2}{\mu k}(\sqrt{1+\mu ^2k^2}\mu k)+(\mu k\sqrt{1+\mu ^2k^2}+\mu ^2k^2+\frac{1}{2})]$$
(40)
$`+3(\mu k\sqrt{1+\mu ^2k^2}+\mu ^2k^2+{\displaystyle \frac{1}{2}})\}{\displaystyle }_0^1dv{\displaystyle \frac{v^2(1v^2/3)}{\mu ^2+k^2(1v^2)}}.`$

We can extend the momentum integration region in this integral to zero, since contributions from small momenta are additionally suppressed by powers of $`\sigma \mu `$. It is also natural to rescale the integration momentum $`k\mu k`$. Then the expression for the radiative-recoil corrections has the form

$$\delta ϵ_\mu ^{e\mu }=_0^{\mathrm{}}dk\{7[2(\sqrt{1+k^2}k)+k(k\sqrt{1+k^2}+k^2+\frac{1}{2})]$$
(41)
$`+6k(k\sqrt{1+k^2}+k^2+{\displaystyle \frac{1}{2}})\}{\displaystyle }_0^1dv{\displaystyle \frac{v^2(1v^2/3)}{1+k^2(1v^2)}}.`$

The integral in the first line of this equation is proportional to the integral for the muon polarization contribution of order $`\alpha (Z\alpha )(m/M)\stackrel{~}{E}_F`$ in , and the contribution of the integral in the second line may be calculated in the same way. Finally we obtain

$$\delta ϵ_\mu ^{e\mu }=\frac{5\pi ^2}{12}+\frac{1}{18}=4.0567796\mathrm{}.$$
(42)
## V Diagrams with Radiative Photons in the Muon Line and Muon Vacuum Polarization

The expression in eq.(2), for the radiative corrections generated by the diagrams with the radiative photon insertions in the electron line Fig. 1, was obtained without expansion in the mass ratio. Hence, after the substitutions $`mM`$ and $`\alpha Z^2\alpha `$ it goes into the expression for the corrections generated by the diagrams in Fig. 4We use the convention that the heavy line charge is $`Ze`$.

$$\delta E^{\mu line}=\frac{(Z^2\alpha )(Z\alpha )}{\pi ^2}\frac{m}{M}E_F\frac{1}{4}_0^1𝑑x_0^x𝑑y\frac{d^4k}{i\pi ^2(k^2+i0)^2}$$
(43)
$`\left({\displaystyle \frac{1}{k^2+4\mu k_0+i0}}+{\displaystyle \frac{1}{k^24\mu k_0+i0}}\right)`$

$`\{(3k_0^22𝐤^2)[{\displaystyle \frac{c_1𝐤^2+c_2(k^2)^2}{(k^2+2bk_0+a^2)^3}}+{\displaystyle \frac{c_3k^2+c_42k_0}{(k^2+2bk_0+a^2)^2}}]`$

$`3k_0[{\displaystyle \frac{c_5k^2+c_6k^22k_0}{(k^2+2bk_0+a^2)^2}}+{\displaystyle \frac{c_7k^2}{k^2+2bk_0+a^2}}]\}.`$

The dimensionless integration momentum here is measured in muon mass units, and apparent incomplete symmetry with the expression in eq.(2) is due to asymmetry of the parameter $`\mu =m/(2M)`$.

Unlike the case of the radiative photon insertions in the electron line in Section II, the integral in eq.(43) does not generate nonrecoil contributions . This happens because the radiative photon insertions in the muon line suppress the contribution from the integration momenta of order of the electron mass which were responsible for the nonrecoil correction by an extra factor $`\mu ^2`$. Hence, for calculation of the leading order recoil corrections we can safely allow $`\mu 0`$ in eq.(43)

$$\delta E^{\mu line}\frac{(Z^2\alpha )(Z\alpha )}{\pi ^2}\frac{m}{M}E_F\frac{1}{2}_0^1𝑑x_0^x𝑑y\frac{d^4k}{i\pi ^2(k^2+i0)^3}$$
(44)
$`\{(3k_0^22𝐤^2)[{\displaystyle \frac{c_1𝐤^2+c_2(k^2)^2}{(k^2+2bk_0+a^2)^3}}+{\displaystyle \frac{c_3k^2+c_42k_0}{(k^2+2bk_0+a^2)^2}}]`$

$`3k_0[{\displaystyle \frac{c_5k^2+c_6k^22k_0}{(k^2+2bk_0+a^2)^2}}+{\displaystyle \frac{c_7k^2}{k^2+2bk_0+a^2}}]\}.`$

After the Wick rotation and transition to the four-dimensional spherical coordinates we have

$$\delta E^{\mu line}=\frac{(Z^2\alpha )(Z\alpha )}{\pi ^2}\frac{m}{M}E_F\frac{1}{2\pi }_0^\pi 𝑑\theta \mathrm{sin}^2\theta _0^{\mathrm{}}𝑑k^2𝒟(k,\theta )\frac{1}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta },$$
(45)
where the differential operator $`𝒟(k,\theta )`$ was defined in eq.(5). There are no $`\mu `$-dependent terms in eq.(45), and this expression for the energy shift is similar to the $`c`$-integrals in eq.(7). It can be formally obtained from the expression in eq.(7) by the substitution

$`{\displaystyle \frac{b^2k^4}{(k^2+a^2)^2}}Z^2.`$

All corrections to HFS described by the integral in eq.(45) were analytically calculated in . Contributions to the energy shift generated by the diagrams in Fig. 5 are obtained from the expression in eq.(45) by insertion in the integrand of the doubled muon vacuum polarization

$$2\left(\frac{\alpha }{\pi }\right)k^2I_\mu (k)=2\left(\frac{\alpha }{\pi }\right)k^2_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}.$$
(46)
This expression is apparently different from the respective expression in eq.(33) because now our dimensionless momenta are measured in terms of the muon mass. Finally, the integral for the contribution to HFS generated by the diagrams in Fig. 5 has the form

$$\delta E^{\mu \mu }=\frac{\alpha (Z^2\alpha )(Z\alpha )}{\pi ^3}\frac{m}{M}E_F\frac{1}{\pi }_0^\pi 𝑑\theta \mathrm{sin}^2\theta _0^{\mathrm{}}𝑑k^2𝒟(k,\theta )$$
(47)
$`{\displaystyle \frac{k^2}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }}{\displaystyle _0^1}𝑑v{\displaystyle \frac{v^2(1v^2/3)}{4+k^2(1v^2)}}=\left({\displaystyle \underset{1}{\overset{7}{}}}\delta ϵ_i^{\mu \mu }\right){\displaystyle \frac{\alpha (Z^2\alpha )(Z\alpha )}{\pi ^3}}{\displaystyle \frac{m}{M}}E_F.`$

All seven contributions are collected in Table VI, and the total contribution to HFS generated by the diagrams in Fig. 5 is

$$\delta ϵ^{\mu \mu }=1.80176(2).$$
(48)
## VI Diagrams with Radiative Photons in the Muon Line and Electron Vacuum Polarization

Let us turn now to the diagrams in Fig. 6. The only difference between these diagrams and the diagrams in Fig. 5 is that now we have electron polarization insertions instead of the muon polarization insertions. The respective analytic expression is obtained from the one in eq.(47) by substitution of the electron polarization instead of the muon polarization

$$2k^2_0^1𝑑v\frac{v^2(1v^2/3)}{4+k^2(1v^2)}2k^2_0^1𝑑v\frac{v^2(1v^2/3)}{4(m/M)^2+k^2(1v^2)}.$$
(49)
The typical integration momenta in eq.(45) are of order of the muon mass and, since we are calculating to linear accuracy in $`m/M`$, it is sufficient to substitute instead of the electron polarization operator the leading terms in its expansion over $`m/M`$

$$2k^2_0^1𝑑v\frac{v^2(1v^2/3)}{4(m/M)^2+k^2(1v^2)}\frac{2}{3}\mathrm{ln}k^2+\frac{4}{3}\mathrm{ln}\frac{M}{m}\frac{10}{9}.$$
(50)
Then the total contribution to HFS generated by the diagrams in Fig. 6 is given by the integral

$$\delta E^{\mu e}=\frac{\alpha (Z^2\alpha )(Z\alpha )}{\pi ^3}\frac{m}{M}E_F\frac{1}{2\pi }_0^\pi 𝑑\theta \mathrm{sin}^2\theta _0^{\mathrm{}}𝑑k^2𝒟(k,\theta )$$
(51)
$`{\displaystyle \frac{1}{(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta }}\left({\displaystyle \frac{2}{3}}\mathrm{ln}k^2+{\displaystyle \frac{4}{3}}\mathrm{ln}{\displaystyle \frac{M}{m}}{\displaystyle \frac{10}{9}}\right)`$

$`=\left(\delta ϵ_c^{\mu e}+\delta ϵ_s^{\mu e}\right){\displaystyle \frac{\alpha (Z^2\alpha )(Z\alpha )}{\pi ^3}}{\displaystyle \frac{m}{M}}E_F,`$

where the contribution $`\delta ϵ_c^{\mu e}=_iϵ_{ci}^{\mu e}`$ corresponds to the first term in the last brackets, and the contribution $`\delta ϵ_s^{\mu e}`$ corresponds to the second and third terms in the same brackets. We calculate $`\delta ϵ_c^{\mu e}`$ numerically and collect all seven contributions to this integral in Table VII. The sum of all these contributions is

$$\delta ϵ_c^{\mu e}=12.94447(4).$$
(52)
The integral for $`\delta ϵ_s^{\mu e}`$ is proportional to the integral for the radiative-recoil corrections generated by the radiative photon insertions in the muon line, which is known analytically . Hence, we can immediately put down an analytic result for $`\delta ϵ_s^{\mu e}`$

$$\delta ϵ_s^{\mu e}=\left(\frac{4}{3}\mathrm{ln}\frac{M}{m}\frac{10}{9}\right)\left(\frac{9}{2}\zeta (3)3\pi ^2\mathrm{ln}2+\frac{39}{8}\right).$$
(53)
Then the total contribution to HFS generated by the diagrams in Fig. 6 is equal to

$$\delta ϵ^{\mu e}=\left(6\zeta (3)4\pi ^2\mathrm{ln}2+\frac{13}{2}\right)\mathrm{ln}\frac{M}{m}+24.32115(4).$$
(54)
## VII Discussion of Results

Collecting the results in eq.(32), eq.(42), eq.(48), and eq.(54) we obtain the three-loop single–logarithmic and nonlogarithmic corrections generated by the one-loop electron and muon polarization insertions in the exchanged photons

$$\delta E=\delta ϵ^e\frac{\alpha ^2(Z\alpha )}{\pi ^3}\frac{m}{M}\stackrel{~}{E}_F+\delta ϵ^\mu \frac{\alpha (Z^2\alpha )(Z\alpha )}{\pi ^3}\frac{m}{M}\stackrel{~}{E}_F,$$
(55)
where

$$\delta ϵ^e=\frac{5}{2}\mathrm{ln}^2\frac{M}{m}+\frac{22}{3}\mathrm{ln}\frac{M}{m}+7.36110(3),$$
(56)
$$\delta ϵ^\mu =\left(6\zeta (3)4\pi ^2\mathrm{ln}2+\frac{13}{2}\right)\mathrm{ln}\frac{M}{m}+22.51939(5).$$
(57)
If we recall that in the real muonium $`Z=1`$ and that the logarithm squared term was already calculated earlier  then the new three-loop single-logarithmic and nonlogarithmic corrections obtained above may be written as

$$\delta E_f=\left[\left(6\zeta (3)4\pi ^2\mathrm{ln}2+\frac{83}{6}\right)\mathrm{ln}\frac{M}{m}+29.88049(6)\right]\frac{\alpha ^3}{\pi ^3}\frac{m}{M}\stackrel{~}{E}_F.$$
(58)
Combining this result with the result of our earlier paper  we obtain all three-loop single–logarithmic and nonlogarithmic corrections generated by the electron and muon polarization insertions in the exchanged photons

$$\delta E_{tot}=\left[\left(4\pi ^2\mathrm{ln}2+\frac{67}{12}\right)\mathrm{ln}\frac{M}{m}+9.59318(6)\right]\frac{\alpha ^3}{\pi ^3}\frac{m}{M}\stackrel{~}{E}_F,$$
(59)
or, numerically

$$\delta E=0.0288\text{kHz}.$$
(60)
This correction has the same scale as some other corrections to hyperfine splitting in muonium calculated recently . All these results improve the accuracy of the theory of hyperfine splitting and affect the value of the electron-muon mass ratio derived from the experimental data  on hyperfine splitting (see, e.g., reviews in ). We postpone discussion of the phenomenological implications of the result above until the completion of the calculations of the remaining three-loop radiative recoil corrections.

###### Acknowledgements.

This work was supported in part by the NSF grant PHY-0138210. Work of V.A. Shelyuto was also supported by the Russian Foundation for Basic Research.

## A Standard auxiliary functions

Calculations described in this paper are greatly facilitated by the use of the auxiliary functions $`\mathrm{\Phi }_n(k)`$ ($`n=0,1,2,3`$) defined by the relationship

$$\mathrm{\Phi }_n(k)\frac{1}{\pi \mu ^2}_0^\pi 𝑑\theta \mathrm{sin}^2\theta \mathrm{cos}^{2n}\theta \frac{(k^2+a^2)^24\mu ^2b^2k^4}{(k^2+\mu ^2\mathrm{cos}^2\theta )[(k^2+a^2)^2+4b^2k^2\mathrm{cos}^2\theta ]}.$$
(A1)
The integral over angles may be explicitly calculated, and the result of the integration is conveniently written as a sum

$$\mathrm{\Phi }_n(k)\mathrm{\Phi }_n^S(k)+\mathrm{\Phi }_n^\mu (k)+\mathrm{\Phi }_n^C(k),$$
(A2)
where

$$\mathrm{\Phi }_n^S(k)=\frac{\delta _{n0}}{\mu k},$$
(A3)
$$\mathrm{\Phi }_0^\mu (k)=W(\xi _\mu )\frac{1}{\sqrt{\xi _\mu }},$$
(A4)
$$\mathrm{\Phi }_1^\mu (k)=\xi _\mu W(\xi _\mu )+\frac{1}{2},$$
(A5)
$$\mathrm{\Phi }_2^\mu (k)=\xi _\mu \left(\xi _\mu W(\xi _\mu )\frac{1}{2}\right)+\frac{1}{8},$$
(A6)
$$\mathrm{\Phi }_3^\mu (k)=\xi _\mu \left[\xi _\mu \left(\xi _\mu W(\xi _\mu )\frac{1}{2}\right)+\frac{1}{8}\right]+\frac{1}{16},$$
(A7)
$$\mathrm{\Phi }_0^C(k)=W(\xi _C),$$
(A8)
$$\mathrm{\Phi }_1^C(k)=\xi _CW(\xi _C)\frac{1}{2},$$
(A9)
$$\mathrm{\Phi }_2^C(k)=\xi _C\left(\xi _CW(\xi _C)\frac{1}{2}\right)\frac{1}{8},$$
(A10)
$$\mathrm{\Phi }_3^C(k)=\xi _C\left[\xi _C\left(\xi _CW(\xi _C)\frac{1}{2}\right)+\frac{1}{8}\right]\frac{1}{16}.$$
(A11)
The standard function $`W(\xi )`$ has the form

$$W(\xi )=\sqrt{1+\frac{1}{\xi }}1.$$
(A12)
and

$$\xi _\mu =\mu ^2k^2,\xi _C=\frac{(k^2+a^2)^2}{4b^2k^2}.$$
(A13)
One may easily obtain asymptotic expressions for the function $`W(\xi )`$

$$\underset{\xi 0}{lim}W(\xi )\frac{1}{\sqrt{\xi }},$$
(A14)
$`\underset{\xi \mathrm{}}{lim}W(\xi ){\displaystyle \frac{1}{2\xi }}.`$

High- and low-momentum asymptotic expressions for the functions $`\mathrm{\Phi }_i(k)`$ may also be easily calculated. Let us cite low-momentum expansions, which were used in the main text for calculation of the contributions to the hyperfine splitting of relative order $`\mu ^2`$

$$\mathrm{\Phi }_0^\mu (k)1+\frac{\mu k}{2},$$
(A15)
$$\mathrm{\Phi }_1^\mu (k)\frac{1}{2}\mu k+(\mu k)^2,$$
(A16)
$$\mathrm{\Phi }_2^\mu (k)\frac{1}{8}\frac{(\mu k)^2}{2}.$$
(A17)


text:
# Cosmic ray spectrum and anisotropies from the knee to the second knee

## 1 Introduction

A puzzling feature of the cosmic ray (CR) spectrum is the so-called knee, i.e. a steepening taking place at an energy $`E_{knee}3\times 10^{15}`$ eV. Several models have been proposed so far in order to explain this feature, although none of them has managed to become broadly accepted. Some proposals focus on a possible crossover between different acceleration mechanisms below and above the knee , or exploit the possibility of a change in the particle acceleration efficiency . Other scenarios include the effects of nuclear photodisintegration processes at the sources , the recent explosion of a single source , and increased leakage from the Galaxy due to a change in the confinement efficiency of CRs by galactic magnetic fields . However, considering the information on the spectrum alone it is difficult to discriminate among the predictions worked out from the proposed models, and so other kind of observations, namely the CR mass composition and the anisotropy measurements, are decisive in testing the different proposals.

For instance, the scenario that considers nuclear photodisintegration at the sources predicts a CR composition that becomes lighter above the knee, while other (rigidity dependent) explanations predict the opposite trend. Unfortunately, the observational status is far from being clear in this respect , with some observations  suggesting that the CR composition could turn lighter while others finding that the heavier components become dominant .

On the other hand, anisotropy observations around and above the knee appear in this context as a valuable means of imposing further restrictions on the possible scenarios. In particular, scenarios in which the knee is due to a change in the acceleration efficiency at the sources , or to the effects of nuclear photodisintegration in them , are not expected to lead to any particular change in the behavior of the observed anisotropies in correlation with the observed change in the spectral index at the knee energy, while scenarios in which the knee is linked to a change in the efficiency of CR confinement in the Galaxy  can predict such correlation. Moreover, to account for the observed energy dependence of the anisotropy amplitude and phase of first harmonic in the whole range of energies between the knee and the ankle is a highly non-trivial requirement for any of these scenarios.

The aim of this work is to perform a more detailed comparison between the anisotropy observations and the predictions within the scenario in which the knee is due to a crossover from a regime in which the CR transport is dominated by transverse diffusion to one dominated by drift effects , recently reconsidered in detail in refs. , and hereafter referred as the diffusion/drift model. This analysis will be done for the whole range of energies in which drift effects are relevant ($`E>E_{knee}`$) and where CRs are dominated by the galactic component ($`E<E_{ankle}5\times 10^{18}`$ eV). The CR transport in the Galaxy is determined by solving numerically a diffusion equation that takes account of the regular galactic magnetic field as well as a random turbulent field with a Kolmogorov spectrum of inhomogeneities, while for the extragalactic component an isotropic flux is adopted. Then, under different assumptions concerning the galactic field model considered, the particular set of parameters adopted within the given model (i.e. the field amplitudes and scale heights assumed for both the regular and random field components), the source distribution and the spectrum of the extragalactic (isotropic) flux, the results are contrasted to the observed CR energy spectrum, the anisotropy amplitude and phase of first harmonic measured by different experiments. It is found that the results show the right tendency to naturally account for the observed features, and hence these observations provide further support to this model.

## 2 The diffusion/drift model

To understand the main physical ingredients of the diffusion/drift model considered in refs. , let us assume the CR transport in the Galaxy to be described by a steady-state diffusion equation $`𝐉=Q`$, where $`Q`$ is the source and the CR current is related to the CR density $`N`$ through

$$𝐉=D_{}_{}ND_{}_{}N+D_A𝐛\times N,$$
(1)
with $`𝐛`$ being the unit vector in the direction of the regular magnetic field $`𝐁_0`$, i.e. $`𝐛𝐁_0/|𝐁_0|`$, and $`_{}=𝐛(𝐛)`$, while $`_{}=_{}`$. The components of the diffusion tensor are $`D_{}`$ (along the direction of $`𝐁_0`$), $`D_{}`$ in the perpendicular direction, while $`D_A`$ is associated to the antisymmetric part and determines the drift effects. Assuming for simplicity that the regular magnetic field is directed in the azimuthal direction and that in a first approach both the Galaxy and the CR sources can be considered to have cylindrical symmetry, one finds that $`D_{}`$ plays no role in the diffusion equation, which involves then only $`D_{}`$ and $`D_A`$. These diffusion coefficients are obtained from the approximate expressions given in ref. . Indeed, defining the turbulence level $`\eta B_{rand}^2/(B_{rand}^2+B_0^2)`$, where $`B_{rand}`$ is the rms amplitude of the random component of the magnetic field, the relevant diffusion coefficients are well described in the range of interest of this work by

$$D_{}8.7\times 10^{27}e^{3.24\eta }\left(1\eta \right)^{1/6}\left(\frac{r_L}{\mathrm{pc}}\right)^{1/3}\frac{\mathrm{cm}^2}{\mathrm{s}}$$
(2)
and

$$D_A3\times 10^{28}\left(1+0.049\left(\frac{\eta }{1\eta }\right)^{3/2}\right)^1\left(\frac{r_L}{\mathrm{pc}}\right)\frac{\mathrm{cm}^2}{\mathrm{s}},$$
(3)
where $`r_L`$ is the CR Larmor radius. The spatial distribution of sources will be specified below in Section 3.

The change in the spectrum from the diffusion to the drift dominated regimes can be simply understood from the fact that, as turns out from Eqs.(2) and (3), the energy dependence of the diffusion coefficients is given by $`D_{}(𝐱)D_{}^0(𝐱)(E/E_0)^{1/3}`$ and $`D_A(𝐱)D_A^0(𝐱)E/E_0`$ 111In general, for a spectrum of random magnetic field inhomogeneities satisfying $`\mathrm{d}B^2/\mathrm{d}kk^a`$ (where $`a=5/3`$ corresponds to the Kolmogorov case, or $`a=3/2`$ for a Kraichnan  hydromagnetic spectrum), one has $`D_{}E^{2a}`$.. Hence, one finds that at low energies ($`E<ZE_{knee}`$, where $`Z`$ is the CR charge) the transverse diffusion is the dominant process affecting the transport of CRs and leading to $`\mathrm{d}N/\mathrm{d}E(D_{}^0/D_{})\mathrm{d}Q/\mathrm{d}EE^{\beta 1/3}`$, where $`\beta `$ is the spectral index of the source, while at high energies ($`EZE_{knee}`$) one has instead $`\mathrm{d}N/\mathrm{d}E(D_A^0/D_A)\mathrm{d}Q/\mathrm{d}EE^{\beta 1}`$. Thus, we see that for a constant source spectral index $`\beta 2.4`$, one has the correct spectral slope below the knee ($`\mathrm{d}N/\mathrm{d}EE^\alpha `$, with $`\alpha 2.7`$). Each CR component of charge $`Z`$ starts to be affected by drifts at an energy $`EZE_{knee}`$, and its spectrum progressively steepens, with the spectral index finally changing by $`\mathrm{\Delta }\alpha 2/3`$ in a decade of energy. The envelope of the total spectrum obtained by adding together the different nuclear components nicely fits the change from a spectrum $`E^{2.7}`$ below the knee, to one $`E^3`$ above it. Moreover, since all the lighter components are strongly suppressed above $`10^{17}`$ eV, the dominant iron component will progressively steepen its spectrum until the overall spectrum becomes $`E^{2.72/3}`$ above a few $`\times 10^{17}`$ eV, hence also reproducing the behavior observed at the so-called second knee, namely a second steepening taking place at $`E_{sk}4\times 10^{17}`$ eV, in which the spectral slope is close to $`\alpha 3.3`$. Furthermore, we have shown  that this scenario predicts a CR mass composition which is compatible with those experiments that report a heavier composition above the knee. Concerning the anisotropy observations, we have also shown  that a few $`\%`$ anisotropy with an excess from a direction near the galactic center and a deficit near the galactic anticenter direction around $`10^{18}`$ eV  could be accounted for by diffusion and drift currents affecting the galactic CR component. In this way, this simple scenario that just takes into account well established properties of the propagation of charged particles in regular and turbulent magnetic fields can provide an explanation of quite diverse galactic CR observations without requiring any additional assumptions.

## 3 The cosmic ray spectrum and anisotropies

The contribution to the anisotropy from galactic CRs of a given charge $`Z_i`$ will be

$$\delta _i=\frac{3𝐉_𝐢}{cN_i},$$
(4)
where $`𝐉_𝐢`$ is the CR current corresponding to the given component, and is hence given by an expression analogous to Eq.(1) but involving the densities $`N_i`$ corresponding to the CRs with charge $`Z_i`$. The total anisotropy will then be given by $`\delta =_if_i\delta _i`$, where $`f_iN_i/N`$ are the fractional abundances of all CR species (with $`N`$ the total, galactic plus extragalactic, CR density). If we assume that the extragalactic flux is isotropic, it can be shown that this flux, which presumably gives the dominant contribution above the ankle, will not be enhanced by the diffusion process and hence will remain isotropic . This implies that the presence of the assumed isotropic extragalactic component just reduces the fractions $`f_i`$ of the galactic components, having the effect of suppressing the growth of the total anisotropy amplitude as the energy of the ankle is approached.

Notice that under the assumptions of cylindrical symmetry and that the regular magnetic field is in the azimuthal direction, one has that $`_{}N=N`$ and also that the CR current will be perpendicular to the regular magnetic field (i.e. lying on the $`rz`$ plane). Moreover, the contribution to the CR current arising from the transverse diffusion will be in the direction of $`N`$, while the drift part will be orthogonal to $`N`$.

As we mentioned before, at low energies the diffusion is completely determined by $`D_{}`$, and hence the contribution to the anisotropy from CRs of charge $`Z_i`$ is $`\delta _iD_{}N_i/N_i`$, increasing with energy roughly as $`\delta _iE^{1/3}`$. However, at higher energies there is a crossover to the drift dominated regime, leading to the behavior $`\delta _iE`$. At even higher energies, above $`10^{18}`$ eV, the CR density gradually starts to be dominated by the extragalactic component, and hence the anisotropy will be determined by the intrinsic anisotropy of the extragalactic component (which for simplicity was assumed to be negligible in this work).

The field models considered in this work are similar to those studied in ref.  (except for a minor modification of the vertical profile function). We will consider for simplicity that the regular magnetic field is directed in the azimuthal direction and that the system possesses azimuthal symmetry. The propagation region of the galactic CRs will be taken as a cylinder of radius $`R=20`$ kpc and height $`2H`$, with a halo size $`H=10`$ kpc. Adopting cylindrical coordinates $`(r,\varphi ,z)`$ throughout, the regular galactic magnetic field can be expressed as $`𝐁_\mathrm{𝟎}=(B_0^{disk}+B_0^{halo})\widehat{\varphi }`$, where $`B_0^{disk}`$ and $`B_0^{halo}`$ are $`\varphi `$-independent functions that correspond to the disk and halo regular field components, respectively. The disk component is assumed to be given by

$$B_0^{disk}(r,z)=B_d\sqrt{\frac{1+(r_{obs}/r_c)^2}{1+(r/r_c)^2}}\mathrm{sin}\left(\frac{\pi r}{4\mathrm{kpc}}\right)\mathrm{th}^2\left(\frac{r}{1\mathrm{kpc}}\right)\frac{1}{\mathrm{cosh}(z/z_d)},$$
(5)
where $`r_c=4`$ kpc is a core radius which smoothes out near the galactic center the overall $`1/r`$ behavior usually considered, $`z_d`$ is the disk scale height, and the value of $`B_d`$ is chosen such that $`B_0^{disk}(r_{obs},z_{obs})=1`$ $`\mu `$G at our observation point ($`r_{obs}=8.5`$ kpc, $`z_{obs}=0`$), the minus sign arising from the fact that the local galactic magnetic field is nearly directed in the $`\widehat{\varphi }`$ direction. The th$`{}_{}{}^{2}(r)`$ factor has been introduced just to insure that there are no singularities in the $`r=0`$ axis, since otherwise a singular drift current along that axis would be artificially produced. Notice also that this regular field component has reversals in its direction for $`r=4,`$ 8, 12 and 16 kpc.

The halo component can be described either by a field structure with radial reversals and also symmetric with respect to the galactic plane (denoted by $`R`$-$`S`$), or by a model without radial reversals (denoted $`NR`$) which may correspond to an independently generated halo field. In the latter case four different configurations arise from choosing the relative orientation of the halo field relative to the disk one, and from the symmetry with respect to the galactic plane ($`S`$ models having $`B(r,z)=B(r,z)`$, while $`A`$ models having $`B(r,z)=B(r,z)`$). In this work, we will consider only the $`NR`$-$`A+`$ and $`NR`$-$`S+`$ cases, the $`+`$ sign denoting that sign$`(B_0^{halo}B_0^{disk})>0`$ at our galactocentric radius for $`z>0`$, since a detailed discussion comparing the different field configurations has already been worked out previously .

The halo component is then taken as

$$B_0^{halo}(r,z)=B_h\sqrt{\frac{1+(r_{obs}/r_c)^2}{1+(r/r_c)^2}}\mathrm{th}^2\left(\frac{r}{1\mathrm{kpc}}\right)\frac{1}{\mathrm{cosh}(z/z_h)},$$
(6)
where the normalization factor is assumed to be $`B_h=1\mu `$G, and a halo scale height $`z_h`$ is introduced in the vertical profile function. The function $``$ is just unity for the symmetric halo case, and it is taken as $`(z)=\mathrm{tanh}(z/z_d)`$ in the case of a halo model antisymmetric with respect to the galactic plane. The Faraday rotation measures of pulsars and extragalactic radio sources might indicate that the global field structure is best described by a bi-symmetric spiral field in the disk with reversed direction from arm to arm, and an azimuthal field in the halo with reversed directions below and above the galactic plane, which could correspond to an A0 dynamo field configuration . In that case, the $`NR`$-$`A+`$ field model may reflect appropriately the main features of the large-scale galactic field structure. However, according to present observational data, the origin of the galactic field in the halo could also be related to the origin of the disk field, for instance if a galactic wind is present, somehow extending the field disk properties into the halo. Hence, the $`NR`$-$`S+`$ field model, as well as other variants discussed in ref. , are also plausible.

For the random component, its intensity is taken as

$$B_{rand}(r,z)=B_r\sqrt{\frac{1+(r_{obs}/r_c)^2}{1+(r/r_c)^2}}\frac{1}{\mathrm{cosh}(z/z_r)},$$
(7)
where $`z_r`$ is the corresponding scale height and $`B_r`$ the local rms amplitude.

The source term in the diffusion equation is assumed to take the form $`Q(r,z)=𝒬(r)\theta (h_s|z|)`$, with $`h_s=200`$ pc, and where the radial profile of the source distribution is considered to be either localized at a given radius $`r_s`$ (i.e. $`𝒬(r)=\delta (rr_s)`$), or as being constant within an inner and an outer radius.

Since we deal here with a rigidity-dependent scenario, the results obtained depend on $`E/Z`$ and a convolution of the contribution of each component of charge $`Z`$ has to be made summing over all components in order to get the all-particle prediction. The source spectrum of a given CR component of charge $`Z`$ is assumed to have an energy dependence d$`Q_Z/\mathrm{d}EE^{\beta _Z}`$, where the normalization factors and spectral indices were obtained from the experimental data at low energies below the knee , which provide the local observed spectral indices $`\alpha _Z=\beta _Z+1/3`$, and extrapolated to higher energies without considering upper cutoffs. Also the relative abundances of the different elements were obtained from the same fits to the low energy observed CR densities . In addition, an isotropic extragalactic component given by

$$\left(\frac{\mathrm{dN}}{\mathrm{dE}}\right)_{XG}=1.7\times 10^{33}\left(\frac{E}{10^{19}\mathrm{eV}}\right)^{2.4}\mathrm{m}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{eV}^1$$
(8)
is added to the galactic flux as well. The energy dependence assumed for the extragalactic CR contribution is similar to that inferred for the production at the sources (i.e. $`\beta 2.4`$), since diffusion effects cannot enhance this flux. The amplitude of the extragalactic component is normalized to fit the observed fluxes for $`E10^{19}`$ eV.

In figure 1, the CR spectra and anisotropy amplitudes obtained for both the $`NR`$-$`A+`$ and $`NR`$-$`S+`$ field models are compared to the observations from several experiments 222The anisotropy experimental data points are obtained from the amplitude of the first harmonic in right ascension $`A`$, which is related to the anisotropy $`\delta `$ through $`A=\delta \mathrm{cos}d\mathrm{cos}\lambda `$, where $`d`$ is declination at which the observations are performed and $`\lambda `$ is the latitude of the direction where the flux is maximum. The plotted points correspond to $`A/\mathrm{cos}d`$ for the different experiments, and thus are actually a lower limit for $`\delta `$, as $`\mathrm{cos}\lambda `$ cannot be recovered from right ascension harmonic analysis.. The set of field parameters ($`B_r,z_d,z_h,z_r`$) used for each model is indicated in the figure in $`\mu `$G and kpc respectively. The effects of varying either the field model parameters or the source distribution have already been discussed in detail in ref. . Very briefly, a turbulence increase, which can be performed either by increasing the random field rms amplitude $`B_r`$ or by increasing the random field scale height $`z_r`$, or also by decreasing the vertical scale heights $`z_d`$ and $`z_h`$ of the regular components, tends to suppress the drift effects, thus yielding larger fluxes and shifting the bent in the spectrum to larger energies.

On the other hand, the drift effects are stronger the farther the source is located, since the drift currents can remove the propagating CRs from the galactic plane all the way from the source to the observer. Moreover, extended source distributions tend to produce flatter spectra due to the smaller associated density gradients. The plots of figure 1 also show that, when considering an extended ring of sources with constant strength, the dominant observed flux is produced near the border which is closest to the observation point.

From figure 1 we observe that, assuming quite plausible field model parameters, the diffusion/drift scenario is capable of reproducing remarkably well not only the observed CR spectrum, but also the anisotropy measurements. Notice from Eqs.(1) and (4) that $`|N|`$ plays a crucial role in the enhancement of the anisotropy amplitude, and thus it turns out that close by sources produce larger anisotropies as compared to farther ones. The effect of considering extended source distributions may be thought as given by the sum of the results produced by localized sources, weighted by the relative contribution of each one to the total CR density observed. Hence, it can happen that just a small region within the extended distribution gives the overwhelming contributions to the total CR flux and anisotropy observed. However, it can also happen that different regions in the source distribution contribute with CR drift currents directed along very different directions and comparable in magnitude, giving some regions of suppressed anisotropy, while producing anisotropy enhancements in others. Anyhow, the results shown for localized sources should be indicative of the actual trend that would be obtained for different source distributions, as follows from figure 1.

Although some results for the $`NR`$-$`S+`$ case exhibit somewhat low anisotropies, it should be noticed that, even in the case of localized sources, we are actually dealing with rings around the galactic center, since non-central point-like sources would break the azimuthal symmetry assumed. Removing this simplifying assumption and solving the CR diffusion equation for a point-like nearby source, a further increase in the anisotropy amplitude might be achieved. Indeed, it has been recognized that the effect of discrete nearby sources, as for instance pulsars and supernova remnants such as Vela, Loop III and Geminga , could actually give the dominant contribution to the observed anisotropies at low energies below the knee. The precise estimate of this effect requires however the knowledge of the spatial distribution, power, age, and evolution of the sources. Hence, an accurate fit of our results to the experimental data below the knee is not really necessary, but it is reassuring to observe that the predictions within this scenario tend to exhibit a quite acceptable agreement with observations above the knee.

In performing the sum over all components in order to obtain the all-particle CR flux, we have taken into account the contribution of all nuclear species from hydrogen to nickel. We have actually also explored the possible contribution of heavier elements up to uranium, following ref. . While the spectral indices of the elements with $`Z28`$ are measured at low energies (below $`10^{14}`$ eV), and then the inferred source index $`\beta _Z`$ is simply extrapolated to higher energies assuming that it remains constant, there is no information available on the spectral indices corresponding to ultra-heavy elements (i.e., with $`Z>28`$). It was proposed in ref.  to obtain them just extrapolating the possible $`Z`$-dependence of the observed spectral indices of the lighter elements, using a parametrization given by $`\alpha _Z=A+BZ^C`$ to fit the data. While this three-parameter expression favored a non-linear extrapolation (with $`C=1.51\pm 0.13`$), also the case of a linear extrapolation (i.e, $`C1`$) was considered and found to fit well the measured spectral indices . Assuming the linear extrapolation of spectral indices, we found that the contribution of ultra-heavy elements is completely negligible within the diffusion/drift scenario. These ultra-heavy elements would only have a noticeable effect with the extreme non-linear extrapolations (for which for instance the spectral index for uranium results $`\alpha _U1.9`$), but in our model this could not lead to a sufficient suppression of the spectrum above $`10^{17}`$ eV, unless one adopts an upper rigidity cutoff in the source near $`EZ\times 10^{17}`$ eV. In ref.  instead, an ad hoc steepening at the knee of $`\mathrm{\Delta }\alpha 2`$ was adopted, and the resulting suppression was actually too large above the second knee. One could also mention that the most natural expectation would be to have an universal rigidity source spectrum, and hence the differences in the spectral indices observed at low energies most likely would reflect the effects of spallation processes affecting heavy nuclei. In this case no particular enhancements of ultra-heavy nuclei above the knee should be expected.

The anisotropy amplitudes of different galactic CR components are plotted in figure 2 as functions of the energy, for the $`r_s=36`$ kpc extended source ring, the $`NR`$-$`A+`$ case and the same field parameters as before. As expected within the diffusion/drift scenario, the onset of the drift-dominated regime is a rigidity-dependent phenomenon that increases the anisotropy of all galactic CR components. It should be pointed out that, although light elements (particularly protons and He nuclei) are strongly suppressed at high energies due to leakage effects, they still contribute significantly to the total anisotropy. For instance, He nuclei are found to contribute around $`3040\%`$ of the total galactic CR anisotropy in the whole energy range between $`10^{14}`$ and $`10^{18}`$ eV for the case shown in figure 2, even though its abundance falls below $`10\%`$ at $`10^{18}`$ eV. It would certainly be of great interest to resolve observationally the anisotropies of separate CR components around and above the knee, since that could help discriminate among some of the models proposed so far. Indeed, KASCADE may actually be capable of performing this kind of observations in the near future333We thank Ralph Engel for bringing this issue to our attention..

In figure 3 we display the anisotropy phase of first harmonic in right ascension vs. energy, again compared to the experimental observations compiled in refs. , for the same cases considered in figure 1. At low energies, the transverse diffusion dominates and the direction of maximum intensity corresponds to the galactic center, as expected . At higher energies the phase of the anisotropy depends actually on the detailed geometry of the regular field adopted, but in the models considered the differences are not large. However, as shown in ref. , for other magnetic field models the anisotropy maximum may point towards the galactic center near $`10^{18}`$eV.

As a summary, in this work we considered the diffusion/drift model to study the propagation of galactic CRs. In agreement with previous results , we find that this scenario explains remarkably well the observed spectrum from below the knee up to the ankle, and in particular the changes in spectral slope at the knee and at the second knee. Furthermore, from the very nature of this scenario, the anisotropy amplitude naturally experiences a crossover from the transverse diffusion dominated regime to the drift dominated one, which shows up around the knee region. Indeed, this is a particular feature that very clearly differentiates this diffusion/drift scenario from other proposals in which such a change in the energy dependence of the anisotropy amplitude is not to be expected. Recalling some well known statements claiming a kind of correlation between the spectrum steepening at the knee and the occurrence of an anisotropy enhancement , one could hence expect anisotropy considerations to favor this kind of scenarios against other proposals. Hence, more detailed and precise measurements concerning anisotropy observations will certainly be decisive in order to test this diffusion/drift scenario, as well as the other models proposed so far, and could also contribute significantly to our knowledge on the galactic field structure and the galactic CR source distribution.

## Acknowledgments

Work supported by CONICET and Fundación Antorchas, Argentina.