text:
# Precise study of the Efimov three-particle spectrum and structure functions within variational approach

## I Introduction

As known, the Efimov effect R1 ; R2 ; R3  reveals itself in a system of three particles with finite-range interaction as the appearance of the infinite number of levels near the threshold with very small energies under the condition that the two-particle subsystems have an infinitely small binding energy or an infinitely large scattering length. In the case of three identical particles, the energy ratio for neighboring highly excited near-threshold Efimov states R1 ; R4  is $`\lambda _E=\underset{n\mathrm{}}{lim}(E_n/E_{n+1})515`$. In three-nucleon systems, the conditions necessary for the effect to take place can be easily broken by taking into account the spin structure of the nuclear interaction R3 ; R6  or the Coulomb long-range repulsive potential.

Theoretical studies of the Efimov near-threshold energy spectrum were carried out both by asymptotic expansions with the use of the small parameter $`r_0/a1`$ (being the ratio of the range of forces to the two-particle scattering length) R1 ; R2 ; R3 ; R5 ; R6  and by numerical calculations R7 ; R8  based on the Faddeev equations by using separable interaction potentials. The main properties of the energy states are universal to the great extent and do not depend on a specific form of the potential. For the special choice of an interaction potential in the form of two components with essentially different radii, new additional excited energy levels appear in the three-particle system — the “trap” effect R9 . In the case of local potentials, numerical calculations of the Efimov energy spectrum were not carried out systematically because of essential difficulties due to the different scales of very small near-threshold excited states energies and the ground-state energy. This requires to develop the precise methods of calculation of the highly excited weakly bounded states of three-particle systems with local interactions.

Due to the use of Gaussian bases in the variational calculations of a few-particle systems R10 , noticeable progress was achieved in the calculations of the main properties of few-body systems R11 ; R12 ; R13 ; R14  with interactions of different nature in recent years. In the present paper, within the precise variational approach, a study of the near-threshold Efimov states of a three-particle system with Gaussian potential is carried out up to the third excited level, and all the main structure functions for various energy states are found.

## II Statement of the Problem

Consider a system of three identical particles

$$H=\frac{\mathrm{}^2}{2m}\underset{i=1}{\overset{3}{}}\mathrm{\Delta }_iV_0\underset{i>j=1}{\overset{3}{}}\mathrm{exp}\left(\frac{r_{ij}^2}{r_0^2}\right)$$
(1)
with a two-particle interaction taken in the Gaussian form for simplicity, where $`V_0`$ is the intensity and $`r_0`$ is the radius of forces. Let us study the properties of the energy spectrum and the main structure functions of the system in symmetric (with respect to permutations of particles) states with zero total angular momentum in the range of parameters where the Efimov effect takes place. It is convenient to use the dimensionless variables and to measure all the distances in $`r_0`$ and energies in $`\mathrm{}^2/mr_0^2`$. In this case, the dimensionless Hamiltonian

$$H=\frac{1}{2}\underset{i=1}{\overset{3}{}}\mathrm{\Delta }_ig\underset{i>j=1}{\overset{3}{}}\mathrm{exp}\left(r_{ij}^2\right)$$
(2)
contains only one combination of physical parameters, the coupling constant $`g=mr_0^2V_0/\mathrm{}^2`$, which determines the properties of the spectrum and wave functions of the system.

The solution of the problem on eigenvalues and eigenfunctions for bound states is carried out in the framework of the variational Galerkin method in the Gaussian representation with the use of special optimization schemes (see R12 ; R13 ; R14 ) which enhance the convergence and ensure a high accuracy of calculations with the least dimensions of variational bases. Within this approach, the variational wave function of symmetric states of the three-particle system with zero total angular momentum can be presented as the expansion in Gaussian functions depending on relative distances $`r_{ij}=\left|𝐫_i𝐫_j\right|`$,

$$\mathrm{\Psi }(r_{12},r_{13},r_{23})=\underset{k=1}{\overset{𝒦}{}}D_k𝐒|\psi _k\underset{k=1}{\overset{𝒦}{}}D_k𝐒\mathrm{exp}\left(a_kr_{12}^2b_kr_{13}^2c_kr_{23}^2\right),$$
(3)
with a set of variational parameters $`a_k`$, $`b_k`$, $`c_k`$. For the bound states, the linear coefficients $`D_k`$ are found in the variational approach from the system of linear algebraic equations of a problem on eigenvalues (which is equivalent to the Schrödinger equation in the chosen Gaussian basis):

$$\underset{n=1}{\overset{𝒦}{}}D_n\left(𝐒\psi _k\left|H\right|𝐒\psi _nE𝐒\psi _k|𝐒\psi _n\right)=0,k=1,2,\mathrm{},𝒦.$$
(4)
Note that all the matrix elements in Eq. (4) for Hamiltonian (2) can be calculated in explicit form. The chosen method of calculation has proved its high efficiency and accuracy in a number of problems of a few particles with interactions of different nature. But, in the case of the near-threshold Efimov spectrum, we are faced with an additional difficulty due to very small values of the energies of these levels. Since the ratio of the Efimov neighboring energy levels $`\lambda _E=\underset{n\mathrm{}}{lim}(E_n/E_{n+1})515`$ is rather large, there exists a real difficulty for all the methods (including the method based on the Faddeev equations R9  and variational methods as well) in calculating the energies with essentially different scales and the corresponding wave functions with oscillations at essentially different distances.

In the present work, the main structure properties of the ground and three excited states of a system of three particles are calculated with the coupling constants near $`gg_{cr}`$, where a two-particle bound state appears and the Efimov effect reveals itself. The numerical analysis of highly excited near-threshold Efimov levels appears to be a nontrivial problem for all the methods including the variational method because of the necessity to prepare the “initial” configurations for the wave functions of highly excited near-threshold bound states to start the minimization of energy in the nonlinear variational parameters $`a_k`$, $`b_k`$, $`c_k`$.

## III The Problem of Initial Configurations

An advantage of the variational method lies in the possibility to calculate (in the approximation of a given number of the basis functions) the whole spectrum of low-energy bound states for a given Hamiltonian or to calculate the energy and wave function of the chosen ground or excited state separately. The further minimization in the nonlinear parameters $`a_k`$, $`b_k`$, $`c_k`$ enables one to improve the variational estimation for the energy and to achieve the best result at a fixed dimension of the basis in (3). It is clear that if the initial configuration gives, in fact, the bound (ground or excited) state, we have to approach the exact solution for this state by minimizing the nonlinear parameters $`a_k`$, $`b_k`$, $`c_k`$ and expanding the basis dimension. However, at a finite and even rather large number $`𝒦`$ of the basis functions in (3) and with a random set of the nonlinear parameters $`a_k`$, $`b_k`$, $`c_k`$ (we call it a random “initial configuration”), nobody can be sure not only in the fact that the studied energy level is close to the exact value, but, even more, that it is below the threshold of a decomposition into subsystems at all. In the latter case, as a rule, the further minimization in nonlinear parameters should be terminated on approaching the rated level to the threshold (from the side of the continuum spectrum), because the wave function can have nothing in common with a bound state. Such problems reveal themselves starting from the second excited state which exists in a narrow interval near $`gg_{cr}=2.684005`$ (the critical constant for the chosen Gaussian potential). There is no such problem for the ground state of three particles, whose initial configuration can be simply prepared even with one Gaussian function at sufficienly large $`g`$. The mentioned problem is also easily solved for the first excited state if one takes several Gaussian functions, among which at least two or three ones depend on the parameters $`a_k`$, $`b_k`$, $`c_k`$ taken from the ground-state function, while the rest have to form a cluster function R13 . This enables one to rather simply construct the necessary configuration (even in the case where the initial configuration corresponds to an energy slightly exceeding the two-particle threshold). For the ground and first excited states, the problem is simplified also due to the fact that these levels exist at all the interaction constants $`g`$ greater than certain critical values (see Table 1). Therefore, it is sufficient to enhance the attractive potential, to form the initial configurations at a sufficiently large $`g`$, and then, by gradually decreasing $`g`$ (the evolution in coupling constant), to achieve the necessary value of $`g`$ by minimizing the energy in the parameters $`a_k`$, $`b_k`$, $`c_k`$. We used both variants to form the initial configurations for the ground and first excited states.

For the following excited near-threshold Efimov states with extremely small binding energies existing within a narrow interval of $`g`$ near the critical two-particle constant $`g_{cr}`$, the problem of initial configurations requires special methods for its solution. One way is to use a wider class of Hamiltonians, in particular, those with an additional interaction giving a possibility to form relatively easily the configurations of highly excited levels. Then the above interaction is eliminated step by step (in order to return to the initial problem) with a simultaneous change in the parameters $`a_k`$, $`b_k`$, $`c_k`$ within the variational procedure in order to keep the energies of the formed configurations below the threshold. The Coulomb attractive interaction seems to be suitable for this purpose giving rise to the infinite number of levels for the three-particle system. An additional advantage of this interaction is the fact that all matrix elements in the Gaussian basis for the Coulomb potential are known to be calculated in explicit form.

Another way to prepare the initial configurations is to deal with a Hamiltonian with different masses of particles. It is known that the Efimov effect takes place in three-particle systems with different masses R1 ; R2 ; R3 . Moreover, the ratio of the energies of neighboring levels depends significantly on the ratio of masses (see also R15 ). From the analysis of the power asymptotics of a three-particle wave function in the momentum representation within the approach based on the Faddeev equations, it can be shown that the following secular equation is the condition for the solvability of three-particle equations:

$$L(1;2,3)L(2;1,3)+L(1;3,2)L(3;1,2)+L(2;3,1)L(3;2,1)+$$

$$L(1;3,2)L(2;1,3)L(3;2,1)+L(1;2,3)L(2;3,1)L(3;1,2)=1;$$
(5)
where

$`L(i;j,k)N({\displaystyle \frac{\mu _{jk}}{\mu _i}};{\displaystyle \frac{\mu _{jk}}{m_k}});N(a;b){\displaystyle \frac{\mathrm{sinh}(s_0\theta /2)}{b\sqrt{a}s_0\mathrm{cosh}(s_0\pi /2)}};\mathrm{cos}\theta {\displaystyle \frac{ab^2}{a+b^2}};`$ (6)
$`\mu _{jk}{\displaystyle \frac{m_jm_k}{m_j+m_k}};\mu _i{\displaystyle \frac{m_i(m_j+m_k)}{m_i+m_j+m_k}}.`$ (7)
In the case of equal masses, the secular equation (5) becomes essentially simpler R1 :

$$N(\frac{3}{4};\frac{1}{2})=\frac{1}{2},$$
(8)
or

$$\frac{8}{\sqrt{3}}\frac{\mathrm{sinh}(\pi s_0/6)}{s_0\mathrm{cosh}(\pi s_0/2)}=1,$$
(9)
and one has for the power $`s_0`$ of the asymptotics of the Faddeev equations solution (the Danilov-Minlos-Faddeev-Efimov constant), $`s_0=1.0062378`$. This constant determines, in particular, the ratio of Efimov neighboring energy levels (the scale invariance multiplier),

$$\lambda _E=\underset{n\mathrm{}}{lim}(E_n/E_{n+1})=\mathrm{exp}\left(\frac{2\pi }{s_0}\right)515.035.$$
(10)
Fig. 1 shows the dependence of the ratio $`\lambda _E=\mathrm{exp}(2\pi /s_0)`$ of the energies of Efimov neighboring levels on the ratio of the masses of particles. It is seen clearly that the case of equal mases is the most unsuitable for the Efimov effect to reveal itself (the scale multiplier $`\lambda _E`$ is the greatest) and, at the same time, it is the most difficult case for numerical calculations. But, in the case of essentially different masses, one may have the ratio of energies not too far from 1 instead of (10) (for any ratio of masses, of course, we have $`\lambda _E>1`$ ). In particular, in the limit $`m_1/m_3=m_2/m_30`$ (the one-center problem), relation (5) yields the limiting secular equation for $`s_0`$:

$$\left\{1+8\left[1+\mathrm{cosh}^1\left(\frac{s_0\pi }{2}\right)\right]\left[\mathrm{sinh}\left(\frac{s_0\pi }{4}\right)/s_0\right]^2\right\}\mathrm{cosh}^2\left(\frac{s_0\pi }{2}\right)=1.$$
(11)
This results in $`s_0=1.139759`$, and the corresponding ratio of energies $`\lambda _E=247.826`$. If, vice versa, $`m_1/m_3=m_2/m_3\mathrm{}`$ (the two-center problem), the value of $`s_0`$ increases infinitely as $`s_0k\sqrt{m_1/m_3}`$, where the coefficient $`k=0.40103`$ is found from the secular equation

$$\frac{\mathrm{exp}(k\sqrt{2})}{k\sqrt{2}}=1.$$
(12)
The corresponding ratio of energy levels goes to 1 as $`\lambda _E1+\frac{2\pi }{k}\sqrt{\frac{m_3}{m_1}}`$ at $`m_3/m_10`$. It should be noted that, due to the essential dependence of the parameter $`\lambda _E`$ on the ratio of masses, it is more probable to observe the Efimov effect in a real physical system just with different masses of particles whose pairwise interaction is of the resonance type. Close to such systems might be three-cluster nuclei consisting of a rather heavy magic core and two neutrons. In particular, if the nucleus 40Ca serves to be the core, it follows from (5) that $`\lambda _E262`$ (i.e., it is twice smaller than that for equal masses). If, to say, we consider the system of mesoatoms of deuterium and tritium interacting with an electron, then, due to the small mass of an electron in comparison with the masses of mesoatoms, one has a large value of $`s_0`$ from (5) enough to result in $`\lambda _E1.26`$. This value is not essentially greater than 1 and much less than that in the case of equal masses. This means that the energies of neighboring Efimov levels would be of the same order of magnitude if they were observed in such a system.

Returning to the problem of formation of initial configurations, we note that the configuration for a given (the $`n`$-th excited) state can be prepared rather easily at different masses. Then one has to restore step by step the equality of masses by minimization of the energy of this level in nonlinear parameters. In this way, we can get the necessary variational wave function corresponding to the configuration of a bound state.

To obtain the initial configurations for variational wave functions, we paid the main attention to the method of evolution in the coupling constant $`g`$ at a fixed (zero) energy, which is possible because the interaction potential in (1) has a definite sign. The method, in essence, consists in a variational procedure with respect to the critical coupling constant $`g`$ to be determined for the given ground or excited state at a fixed (zero) energy. It is the general statement that the variational method can be formulated with respect to the coupling constant $`g`$ at a fixed energy provided that the interaction potential is negative definite. Taking an arbitrary initial configuration with a sufficient number of the basis functions with some parameters $`a_k`$, $`b_k`$, $`c_k`$, we find the corresponding coupling constant $`g`$ of the chosen excited state by variational procedure and advance step by step towards small $`g`$ by varying the parameters $`a_k`$, $`b_k`$, $`c_k`$ (or sometimes, if necessary, by spreading the basis). Thus, we approach the two-particle critical constant from the upper side, where a two-particle ground state appears. When we are close enough to the region $`gg_{cr}`$, where the Efimov excited state under consideration must exist, we return to the commonly used variational procedure in energy using the obtained set of the parameters $`a_k`$, $`b_k`$, $`c_k`$ and start to increase the parameter $`g`$. If we were close enough to the critical constant $`g_{cr}`$, the three-particle energy level crosses the two-particle threshold and goes down below the threshold. Thus, we obtain the required initial configuration. Such a possibility is connected with the fact that the two-particle energy threshold behaves itself by the square law near the critical coupling constant,

$$E_{(2)}\mathrm{\hspace{0.17em}0.107}(gg_{cr})^2,$$
(13)
while three-particle levels reveal a linear behavior R1 ,

$$E_nk_n(gg_{left}^{(n)}),$$
(14)
where $`g_{left}^{(n)}`$ is the critical constant for the $`n`$-th level appearance on the left from the critical constant, and $`k_n=E_n/g`$ (at the critical point $`g=g_{left}^{(n)}`$) is the corresponding angular coefficient (see Table 1). Using the above method, we succeeded in preparing the initial configurations for the second and third excited states which exist near the critical constant $`g_{cr}`$ and belong to the Efimov series of levels. Note that the greater the number of the excited near-threshold Efimov level, the greater are the difficulties in preparing the initial configuration. It is not only due to the necessity to increase the dimension of the variational basis, but also because the nonlinear variational parameters differ from one another by several orders of magnitude. Moreover, there appears a considerable quantity of nonlinear parameters such that they are roughly by $`\lambda _E`$ times smaller than those typical of the previous energy level. This makes it necessary to significantly decrease the steps in the process of variation in nonlinear parameters and decelerates the whole procedure of preparing the initial configuration. Because the energy of already the third excited level is less by nine orders of magnitude in comparison with the ground state energy and the dimension of the basis which is necessary to ensure a reliable accuracy achieves several hundreds of functions, the calculation of levels higher than the third one becomes a serious problem.

## IV Asymmetry of the Efimov Energy Levels

It is convenient to depict the calculated energies of the ground and excited states of a three-particle system as the square roots of the moduli of three-particle energies minus the two-particle energy threshold versus the reciprocal scattering length (instead of the coupling constant $`g`$), because such dependences are model-independent R1 ; R2 ; R3 ; R5  to the great extent. Since the orders of magnitude of the values (of both energies and reciprocal scattering lengths) are essentially different for different levels (the ratio for neighboring levels being about $`\sqrt{\lambda _E}=\mathrm{exp}(\pi /s_0)=22.692`$ ), Fig. 2 shows each level in its own scale in order to emphasize the universal properties of highly excited states. Namely, the horizontal axis depicts $`a^1\lambda _E^{(n2)/2}`$, and the vertical axis shows $`(E_{(2)}E_n)^{1/2}\lambda _E^{(n2)/2}`$. The ground and first excited states exist at any $`g`$ greater than the corresponding critical values (see Table 1). The reciprocal scattering length $`a^1`$ is approximately proportional to $`gg_{cr}`$ near the two-particle critical constant $`gg_{cr}`$. Note that the separation of three-particle levels from the zero energy occurs by law (14), though this linear behavior takes place only in a very narrow interval near the corresponding three-particle critical coupling constant, and the angular coefficient $`k_n`$ decreases essentially with increase in the number of the excited state. The second and third excited states appear with increase in $`g`$ (in the region $`g<g_{cr}`$) and then disappear (in the region $`g>g_{cr}`$) at the two-particle threshold (the corresponding values of the coupling constants $`g`$ and the reciprocal scattering lengths are given in the second and third columns in Table 1). A rather good coincidence of the curves for the second ($`n=2`$) and third ($`n=3`$) excited levels depicted in Fig. 2 testifies to that the asymptotic formulae R1 ; R2  are valid already starting from the second excited state and indicates the fact that highly excited Efimov states possess the scaling property. Although we restricted ourselves by the calculations of only three excited levels, it should be assumed that the rest of levels is to be determined by the asymptotic relations, and thus the properties of the whole spectrum are completely described. A significant asymmetry of the curves obtained in our calculations for the second and third excited states relative to the critical concentration point $`g_{cr}`$ of the Efimov spectrum testifies to the asymmetry of all the rest levels. Each of them extends to the right side from the critical point much farther than to the left. This essential asymmetry of the levels relative to the point $`g=g_{cr}`$ (or $`a^1=0`$) is connected with the following fact. Though the potential of the effective long-range interaction $`R^2`$ R1 ; R2 ; R3  ranges up to distances of the order of $`\left|a\right|`$, it remains attractive at essentially larger distances at $`a>0`$. Whereas, at $`a<0`$, it becomes practically zero out of the region of about $`\left|a\right|`$ (see also R15 ). The above-mentioned asymmetry is clearly seen in Fig. 3, where the second and third excited states are depicted in natural scale versus the coupling constant $`g`$. This figure contains a fragment of the near-threshold area in larger scale in order to make the third level visible, since it looks almost like a dot and cannot be distinguished from the two-particle threshold within the main figure. Note that the largest binding energy of each of the Efimov levels (minus the energy of the two-particle threshold) happens at the constants rather far from the critical concentration point of the spectrum (see Fig. 2 and Fig. 3), though the infinite number of levels appears only at the limiting point $`g_{cr}`$. The three-particle levels cross the two-particle threshold (with increase in $`g`$) at very small angles at the points in Fig. 3, where the solid lines of these levels become the dashed line of the two-particle threshold. After crossing the two-particle threshold, the three-particle levels may exist over the threshold as virtual states (see R16 ). It is worth to note that, with the essential enlargement of the coupling constant $`g`$ up to $`gg_{cr}^{(2)}=16.3`$, the second excited level appears below the two-particle threshold, the third one does at $`gg_{cr}^{(3)}=30.8`$, the fourth does at $`gg_{cr}^{(4)}=32.55`$, the fifth does at $`gg_{cr}^{(5)}=50.89`$, and so on. Their binding energies increase with the coupling constant. The angular coefficients determining the angle between the curve of the $`n`$-th level and the two-particle energy threshold (near the critical point of the corresponding level appearance),

$$E_nE_{(2)}\stackrel{~}{k}_n\left(gg_{cr}^{(n)}\right),$$
(15)
are as follows: $`\stackrel{~}{k}_20.0045`$, $`\stackrel{~}{k}_30.1`$, $`\stackrel{~}{k}_40.01`$, and $`\stackrel{~}{k}_50.03`$. The fact that, for the third and fourth levels, $`g_{cr}^{(3)}`$ and $`g_{cr}^{(4)}`$ are close one to another, as well as that the angular coefficients $`\stackrel{~}{k}_3`$ and $`\stackrel{~}{k}_4`$ are of different order of magnitude, testifies to a different nature of these levels.

## V Structure Functions of the Efimov States

The wave functions of the three-particle ground and excited states calculated within the precise variational approach enable us to calculate directly such structure functions of the system as density distributions, formfactors, pair correlation functions, momentum distributions, etc. It is a rather simple problem to calculate any such an average since the wave functions of each state have the simple form (3) of a superposition of Gaussian functions with the already known linear ($`D_k`$) and nonlinear ($`a_k`$, $`b_k`$, $`c_k`$) parameters.

Fig. 4 presents the one-particle density distributions,

$$n(r)=\frac{1}{3}\underset{i=1}{\overset{3}{}}\delta \left(𝐫\left(𝐫_i𝐑_{c.m.}\right)\right),$$
(16)
for the ground and excited states at the critical two-particle coupling constant $`g_{cr}`$. The density distributions (16) are shown in logarithmic scale on both axes. This is made, on the one hand, in order that the typical “halo” structure appearing for different states at distances of different scales can be shown in the same figure. On the other hand, on the chosen logarithmic scale, the density distributions reveal clearly the almost periodic dependence of the wave functions of Efimov states on the global radius logarithm R1  $`\mathrm{sin}\left(\left|s_0\right|\mathrm{ln}\left(R\sqrt{m\left|E\right|/\mathrm{}^2}\right)+\mathrm{\Delta }\right)`$. Moreover, this asymptotic behavior is valid from the distances of about the radius of forces (in our case, of about 1) to the distances of about $`\mathrm{}\left(m\left|E\right|\right)^{1/2}`$. The less the binding energy of an excited state, the longer is the extension of this asymptotic behavior. The density distribution of the first excited state is seen to change sharply its behavior between the short and long distances. The long-range “halo” is situated around the more condensed central core similar to the density distribution of the ground state. The density distribution of the second excited state changes its behavior two times having a “halo” in the form of two concentric spherical layers put one into another. They are of essentially different radii and densities and have the central part of the same type as in the case of the first excited state. The density distribution of the third excited state changes its behavior three times (three concentric layers around the central core), and so on. The ratio of the radii of neighboring layers is about $`\sqrt{\lambda _E}22.7`$. Since the density distributions are normalized by $`n(r)𝑑𝐫=1`$, they decrease at short distances when the number of the excited state increases, because the size of the system is growing.

The fact that the ratio of the sizes of neighboring states of the system is about $`\sqrt{\lambda _E}`$ is a consequence of the universal model-independent properties of Efimov states and is confirmed once more by calculations of the r.m.s. radii (see Table 1 with the calculated $`R_{rms}`$ at the critical value $`g_{cr}`$). Note that the smallest value of $`R_{rms}`$ for a given Efimov level is achieved at the constant $`g`$ corresponding to the largest binding energy, rather far to the right from the critical point $`g_{cr}`$ (see Fig. 2). When the coupling constant approaches the value, where the Efimov levels appear at the threshold ($`E0`$), the calculated mean square radius

$$R_{rms}r^2^{1/2}=\left(r^2n(r)𝑑𝐫\right)^{1/2}$$
(17)
reveals the behavior

$$r^2^{1/2}=\frac{c_1}{\sqrt{E}}+c_2+\mathrm{}.$$
(18)
That is, the radius goes to infinity, in accordance with the general physical principles, as $`r^2^{1/2}(E)^{1/2}`$ at $`E0`$ (we recall that the energy level near the threshold depends linearly on the coupling constant according to (14)). It is important to notice that the growth of the r.m.s. radii at $`E0`$ is connected mainly with the extension of the region where the asymptotic behavior of the wave function $`\mathrm{sin}\left(\left|s_0\right|\mathrm{ln}\left(R\sqrt{m\left|E\right|/\mathrm{}^2}\right)+\mathrm{\Delta }\right)`$ is valid, i.e., up to the distances of order $`\mathrm{}\left(m\left|E\right|\right)^{1/2}`$. At the same time, the wave function at short and intermediate distances does not change (except for the normalization factor depending on the total size of the system). As a result, in all the intervals of $`g`$ where a given Efimov level exists, the behavior of the wave function and, in particular, the density distribution is not changed at short and intermediate distances (but only the normalization does), which is confirmed by Fig. 5. This figure shows the density distributions of the second excited state at several different values of the coupling constant, from the point $`g=2.682195`$ near the level appearance to the point $`g_{right}^{(2)}=2.7244`$ of its disappearance. Other typical constants chosen by us are as follows: the two-particle critical constant $`g_{cr}`$ (the point of the Efimov spectrum concentration) and the constant $`g=2.702`$, where the energy of the second excited state (minus the two-particle threshold energy) is the greatest.

Specific sharp changes in the density profiles near the distances, where the corresponding wave functions pass through zero, manifest themselves in formfactors,

$$F(q)=\mathrm{exp}(i𝐪𝐫)n(r)𝑑𝐫,$$
(19)
in the form of specific dips of finite depth. Fig. 6 shows the formfactors of the ground state (no dips), the first (one dip), second (two dips), and third (three dips) excited states. A logarithmic scale commonly used for formfactors on the vertical axis is accomplished with the logarithmic scale in the squared momentum transfer in order to show clearly the almost periodic repetition of the dips connected with regularities in the “halo” structure of the density distributions of excited states. Since the greater distances in a density distribution correspond to a less momentum transfer in the corresponding formfactor and vice versa, a new dip in the formfactor of the next excited state appears at a less momentum transfer, whereas the previous dips remain practically at their positions.

Fig. 7 shows the pair correlation functions

$$g_2(r)=\frac{1}{3}\underset{i>j=1}{\overset{3}{}}\delta \left(𝐫\left(𝐫_i𝐫_j\right)\right)$$
(20)
of the ground state and three excited states. We multiply the pair correlation functions by $`r^2`$ in order to show clearly the following interesting effect. With increase in the number of the excited state, we observe the spreading of the area, where the pair correlation function is proportional to the reciprocal squared distance between a pair of particles. That is, the two-particle subsystem in the system of three particles has the density distribution very similar to the squared two-particle wave function. This gives rise to the effective long-range interaction R1 ; R2 ; R3  $`R^2`$ in the system of three particles in the region from distances of order of the radius of forces $`Rr_0`$ to the distances of about $`R\mathrm{}\left(m\left|E\right|\right)^{1/2}`$, where the wave function starts its exponential decrease due to a finite binding energy.

The momentum distributions

$$D(p)=\frac{1}{3}\underset{i=1}{\overset{3}{}}\delta \left(𝐩𝐩_i\right)$$
(21)
for the ground and three excited states are presented in Fig. 8. It is found that the momentum distributions vanish at approximately the same limiting momentum $`p_0`$ (in our case, $`10^0<p_0<10^1`$). This is due to the fact that the largest momenta in the considered three-particle system are typical of the shortest distances in the coordinate representation (of order of the radius of forces). There, the wave function has the sharpest behavior, but it is practically unchangeable under growing the number of an excited state. With increase in the number of a state, the wave function extends to larger and larger distances almost without change in the central part (see, for illustration, Fig. 4 depicting the density distributions). Note that, with increase in the number $`n`$ of the state, the momentum distribution $`D(p)`$ increases essentially at small momenta proportionally to $`\left(\lambda _E^{3/2}\right)^n`$. Indeed, since the size of the system increases by $`\sqrt{\lambda _E}`$ times for each of the next excited states and its volume increases by $`\left(\lambda _E\right)^{3/2}`$ times, the region of momenta where the momentum distribution is mainly concentrated decreases at the same extent. Taking into account $`D(p)𝑑𝐩=1`$, we get that the momentum distribution density grows by $`\left(\lambda _E\right)^{3/2}1.17\times 10^4`$ times at small $`p`$ with every increment of the state number. Fig. 8 also demonstrates the extension (towards zero momentum) of the area of the power dependence of the momentum distribution $`D(p)p^4`$, which looks like a linear decreasing function on the logarithmic scale. Such a behavior is connected with the fact that the greater the number $`n`$ of the state, the greater is the role of the kinetic energy entering the Green’s function in the formation of the power dependence of the wave function in the momentum representation. Really, in the formal operator representation of the Schrödinger equation solution

$$|\mathrm{\Psi }_n=(H_0E_n)^1V|\mathrm{\Psi }_n,$$
(22)
the energy $`E_n0`$ with increase in $`n`$. In this case, instead of the Green’s function, we have only the kinetic energy in the denominator. The expression $`V|\mathrm{\Psi }_n`$ varies slightly, because it is nonzero only within the region of the short-range interaction, where the wave function is changed insignificantly at growing $`n`$. After the integration with the $`\delta `$ \- function, the squared modulus of the wave function in $`D(p)`$ leads to the value in the denominator which is proportional to the squared kinetic energy of one particle, by giving rise to $`D(p)p^4`$.

## VI Conclusions

We have developed a variational approach with the use of the Gaussian basis and precise optimization schemes and have studied the spectrum and wave functions of three particles with Gaussian pairwise interaction near the critical coupling constant, where the Efimov effect takes place. A special technique is developed for preparing the initial configurations of highly excited near-threshold weakly bounded energy states using both the deformation of the initial Hamiltonian (by changing the interaction or the mass ratio) and the method of evolution in the coupling constant at a fixed (in particular, zero) energy. For several calculated levels from the Efimov series, it is shown that they have an essential asymmetry with respect to the critical point and the largest binding energy (minus the threshold energy) significantly far to the right from the critical coupling constant which is the point of the spectrum concentration. The scale invariance is confirmed for the Efimov levels starting from the second excited state of three particles.

For the first time, we have carried out the study of the specific behavior of the density distributions, formfactors, pair correlation functions, and momentum distributions for the ground and three excited states of the Efimov series. The details of the density distributions are found, and the halo-type structure is revealed for the excited states. It is shown that the corresponding formfactors have some dips of finite depth which are arranged periodically in the momentum transfer on the logarithmic scale. The behavior of pair correlation functions and momentum distribution profiles is found, and the main conclusions following from the asymptotic formulae for the Efimov spectrum are confirmed. The main specific properties of the structure functions of three particles within the area of parameters where the Efimov effect reveals itself would take place with other attractive interaction potentials as well.

The precise calculations of the spectrum and the wave functions of highly excited near-threshold levels confirm the high accuracy of the variational method in the Gaussian representation combined with the efficient technique for the maximum optimization of the basis.


format: text

text:
# Proton Zemach radius from measurements of the hyperfine splitting of hydrogen and muonic hydrogen

## I Introduction

The hyperfine splitting of the ground state of the hydrogen atom is among the most accurately measured quantities H-hfs1 ; H-hfs2 :

$$\mathrm{\Delta }E_{exp}^{\mathrm{HFS}}=\mathrm{1\hspace{0.17em}420\hspace{0.17em}405\hspace{0.17em}751.\hspace{0.17em}7667}\pm 0.0009\mathrm{Hz},$$
(1)
The relative experimental uncertainty in Eq. 1 does not exceed $`10^{12}`$. The theoretical predictions for $`\mathrm{\Delta }E^{\mathrm{HFS}}`$, based on QED, are less accurate. This is partly due to the computational difficulties which increase very fast for the higher order terms in the perturbative expansion in powers of $`\alpha `$ and $`(Z\alpha )`$ sapir ; eides , and to the limited precision of the fundamental constants involved (the Rydberg known to $`10^{11}`$, electron to proton mass ratio and $`\alpha `$ known to parts of $`10^8`$ eides ). The main uncertainties come, however, from the insufficient knowledge of the structure of protons. Because of this, the comparison of theoretical results with the experimental value of Eq. 1 can not test QED beyond the contribution of proton polarizability effects of the order of a few ppm. To perform a more precise test of QED, therefore, either additional information on the electromagnetic structure of protons should be used, or the comparison should be done between theoretical and experimental results on the hyperfine splitting in hydrogen-like bound states of charged point-like leptons. Muonium is most appropriate for the latter; indeed, the recent measurements of the hyperfine splitting of the ground state of muonium muonium-hfs-exp  have been shown to agree with theory up to $`\mathrm{0.5\hspace{0.17em}10}^7`$ (see eides  and references therein), so that the correctness of QED results about the hyperfine splitting of hydrogen-like atoms has been experimentally confirmed at least with the same precision. As for the use of information on the electromagnetic structure of protons from other sources in attempts to reduce the theoretical uncertainty of $`\mathrm{\Delta }E^{\mathrm{HFS}}`$, until recently there was no theoretical progress in this direction, and the only realistic idea for years seemed to be to estimate the proton structure contributions from complementary measurements of the hyperfine splitting in muonic hydrogen pla ; wein . We are now going to critically analyze this idea in the context of recent theoretical results on the polarizability of protons faustov-HHFS3 ; cherednikova , of the development of new precision spectroscopy instrumentation kottmann-2000 , and of our better understanding of the dominating proton structure contributions to the hyperfine splitting in hydrogen-like atoms. This will lead us to the alternative point of view to look at the hyperfine splitting measurements in hydrogen and muonic hydrogen as measurements of the Zemach radius of the proton by assuming that all QED predictions are credible. This way we shall obtain a “first experimental” value of the proton Zemach radius from the hyperfine splitting of hydrogen which may be compared to theoretical values based on different proton form factor fits.

## II Hyperfine splitting of hydrogen atom ground state

To analyze the various sources of uncertainty in the theoretical value of $`\mathrm{\Delta }E^{\mathrm{HFS}}`$ we put it in the traditional form bodwin :

$$\mathrm{\Delta }E_{th}^{\mathrm{HFS}}=E^\mathrm{F}\left(1+\delta ^{\mathrm{QED}}+\delta ^{\mathrm{str}}\right)$$
(2)
where $`E^\mathrm{F}`$ is the Fermi splitting fermi  expressed in terms of the electron and proton masses $`m_e,m_p`$ and the dipole magnetic moment of the proton $`\mu _p`$:

$$E^\mathrm{F}=\frac{8}{3}\alpha ^4c^2\frac{m_e^2m_p^2}{(m_e+m_p)^3}\mu _p,$$
(3)
while $`\delta ^{\mathrm{QED}}`$ and $`\delta ^{\mathrm{str}}`$ are correction terms related to higher order QED effects and to proton electromagnetic structure due to strong interactions. (Eq. 2 is only correct in the leading order, since higher orders QED and structure effects mix up karshen-mix .) Up to terms of order $`O(\alpha ^3)`$ (without distinguishing $`\alpha `$ from $`(Z\alpha )`$), $`\delta ^{\mathrm{QED}}`$ is given by sapir ; eides :

$$\delta ^{\mathrm{QED}}=a_e+\frac{3}{2}\alpha ^2+\alpha ^2(\mathrm{log}2\frac{5}{2})\frac{8\alpha ^3}{3\pi }\mathrm{log}\alpha (\mathrm{log}\alpha \mathrm{log}4+\frac{281}{480})+18.984\times \frac{\alpha ^3}{\pi }+\mathrm{}$$
(4)
where $`a_e`$ is the anomalous magnetic moment of the electron. Note that the expression for $`\delta ^{\mathrm{QED}}`$ does not involve the mass ratio $`m_e/m_p`$; all terms which depend on proton mass or come from strong interactions are included in $`\delta ^{\mathrm{str}}`$. In turn, $`\delta ^{\mathrm{str}}`$ splits into a “static” part $`\delta ^{\mathrm{rigid}}`$ that accounts for the elastic electromagnetic form factors of the proton and can be calculated using data from elastic scattering experiments, a part $`\delta ^{\mathrm{pol}}`$ that comes from the internal dynamics of the proton and could only be evaluated using data on inelastic processes with protons, and a part $`\delta ^{\mathrm{hvp}}`$ describing the strong interaction effects outside the proton, such as hadron vacuum polarization: $`\delta ^{\mathrm{str}}=\delta ^{\mathrm{rigid}}+\delta ^{\mathrm{pol}}+\delta ^{\mathrm{hvp}}`$. Two types of “static” proton structure corrections are incorporated in $`\delta ^{\mathrm{rigid}}`$, associated with the spatial distribution of the charge and magnetic moment within the proton and with recoil effects, respectively: $`\delta ^{\mathrm{rigid}}=\delta ^{\mathrm{Zemach}}+\delta ^{\mathrm{recoil}}`$. The former has been calculated by Zemach zemach  and may be put in the form cherednikova :

$$\delta ^{\mathrm{Zemach}}=\frac{2\alpha m_{ep}}{\pi ^2}\frac{d^3p}{(\text{p}^2+\alpha ^2m_{ep}^2)^2}\left(\frac{1}{\mu _p}G_E(\text{p}^2)G_M(\text{p}^2)1\right)=2\alpha m_{ep}R_p+O(\alpha ^3)$$
(5)
where $`m_{ep}=m_em_p/(m_e+m_p)`$ and $`R_p`$ is the first moment of the convolution of the proton charge and magnetic moment distributions, also known as Zemach radius of the proton. The recoil correction $`\delta ^{\mathrm{recoil}}`$ denotes the contribution of all terms which depend on the ratio $`m_e/m_p`$; for sake of simplicity we skip here the rather lengthy explicit expression of $`\delta ^{\mathrm{recoil}}`$ which may be found in bodwin . The proton polarizability correction $`\delta ^{\mathrm{pol}}`$ and the hadron vacuum polarization correction have been evaluated recently faustov-HHFS3 ; faustov-HHFS2 , using the available data on the proton polarized structure functions and on electron-positron annihilation into hadrons.

Eq. 2 may now be put in a more detailed form:

$$\mathrm{\Delta }E_{th}^{\mathrm{HFS}}=E^\mathrm{F}\left(1+\delta ^{\mathrm{QED}}+\delta ^{\mathrm{Zemach}}+\delta ^{\mathrm{recoil}}+\delta ^{\mathrm{pol}}+\delta ^{\mathrm{hvp}}\right)$$
(6)
It is important for what follows to review the order of magnitude of the various terms in the right-hand side of Eq. 6 and briefly discuss the uncertainty of each of them. The uncertainty of $`E^\mathrm{F}`$ is due to the uncertainty of $`m_e`$, $`m_p`$, $`\mu _p`$ and $`\alpha `$ and does not exceed 0.01 ppm partprop . $`\delta ^{\mathrm{QED}}`$ is dominated by the anomalous magnetic moment of the electron; the uncertainty of the term comes from the uncalculated terms of relative order $`O(\alpha ^4)`$ and higher and is estimated not to exceed 0.001 ppm eides . $`\delta ^{\mathrm{rigid}}`$ is close to 40 ppm (bodwin  and references therein); its uncertainty – 2-3 ppm – is due to the uncertainty of the actual values of the proton radii and is of the order of 5%. The recoil term $`\delta ^{\mathrm{recoil}}5.68`$ ppm sapir  adds little uncertainty to $`\delta ^{\mathrm{rigid}}`$. As for the proton polarizability correction $`\delta ^{\mathrm{pol}}`$, until recently there existed only the upper limit $`\delta ^{\mathrm{pol}}<4`$ ppm sapir , obtained by ascribing the whole discrepancy between $`\mathrm{\Delta }E_{exp}^{\mathrm{HFS}}`$ and $`\mathrm{\Delta }E_{th}^{\mathrm{HFS}}`$ to the contribution from $`\delta ^{\mathrm{pol}}`$. The present value of the proton polarizability correction, calculated on the ground of experimental data on the polarized structure function faustov-HHFS3 , is 1.6 $`\pm `$ 0.6 ppm. The hadron vacuum polarization correction $`\delta ^{\mathrm{hvp}}10^8`$ faustov-HVPmu  is much too small. The overall uncertainty of $`\mathrm{\Delta }E_{th}^{\mathrm{HFS}}`$ is therefore of the order of 2-3 ppm and is entirely due to proton structure effects. All discussed quantities are summarized in the leftmost two columns of Table 1; the righmost two columns of the Table contain the numerical values of the corresponding quantities for muonic hydrogen and will be discussed in next section.

As already pointed out, the challenging comparison of $`\mathrm{\Delta }E_{th}^{\mathrm{HFS}}`$ and $`\mathrm{\Delta }E_{exp}^{\mathrm{HFS}}`$ for hydrogen can not be regarded as a test of the QED contributions $`\delta ^{\mathrm{QED}}`$ of order $`O(\alpha ^3)`$ and higher because of the significant overall theoretical uncertainty. We adopt instead an alternative point of view to assume that the theoretical values of $`\delta ^{\mathrm{QED}}`$, $`\delta ^{\mathrm{recoil}}`$, $`\delta ^{\mathrm{hvp}}`$ and $`\delta ^{\mathrm{pol}}`$ are accurate and use the experimental data to determine the Zemach radius of the proton $`R_p`$ as:

$$R_p=\left(\mathrm{\Delta }E_{exp}^{\mathrm{HFS}}/E^\mathrm{F}1\delta ^{\mathrm{QED}}\delta ^{\mathrm{recoil}}\delta ^{\mathrm{pol}}\delta ^{\mathrm{hvp}}\right)/(2m_{ep}\alpha )$$
(7)
The above assumption is justifiable since all four correction terms are objects of QED, the only difference of $`\delta ^{\mathrm{hvp}}`$ and $`\delta ^{\mathrm{pol}}`$ from the former two being that their evaluation requires the use of additional phenomenological information beyond first principles. From Eq. 7 we get the experimental value $`R_p=1.040(16)`$ fm, where the uncertainty $`\pm 0.016`$ fm comes from the theoretical uncertainty of $`\delta ^{\mathrm{pol}}`$ faustov-HHFS3 .

The Zemach radius of the proton is defined in terms of an integral of the charge and magnetic form factors of the proton $`G_E(k)`$ and $`G_E(k)`$ over space-like transfer momenta $`k,k^2=𝐤^2`$ (see Eq. 5, or equivalently, by the first moment of the convolution of the charge and magnetic moment distributions $`\rho _E(r)`$ and $`\rho _M(r)`$ in coordinate space zemach :

$$R_p=d^3rrd^3r^{}\rho _E(𝐫𝐫^{})\rho _M(𝐫^{})$$
(8)
The directly observable quantity which is most sensitive to the Zemach radius of the proton is the hyperfine splitting of bound systems involving protons (compare to the Lamb shift which is related to the proton r.m.s. charge radius pachucki-FSZmH .) The experimental value of $`R_p`$ sets important restrictions on the theoretical models of proton electromagnetic structure and, in particular, on the parametrization of proton form factors, in terms of which the theoretical values are calculated. Consider as an example the values of $`R_p`$ calculated from Eq. 5 using a few popular approximations of the proton form factors. Numerical calculations give $`R_p`$=1.02 fm for the dipole fit, and $`R_p`$=1.067 fm for the fit of simon ; unfortunately, no information on the uncertainty of the parameters of the fit is available, and no conclusions could be made on the compatibility of these values with the ”experimental” one. Both fits are consistent with $`\mu _pG_E(𝐤^2)/G_M(𝐤^2)1`$ for $`𝐤^2<5`$ GeV2/c2. To account for the recent experimental results on the form factor ratio gayou , we also evaluated $`R_p`$ by using the Simon’s fit for either the charge or magnetic form factor and expressing the other one using the relation $`\mu _pG_E(𝐤^2)/G_M(𝐤^2)=10.13(𝐤^20.04)`$ gayou , and got $`R_p`$=1.060 fm and $`R_p`$=1.073 fm respectively. Though preliminary, these estimates show that the current theoretical uncertainty of $`R_p`$ significantly exceeds the experimental one, and that the experimental results on the proton Zemach radius may be used as a test for the quality of models of the proton in the limit of low transfer momenta.

## III Hyperfine splitting of muonic hydrogen atom ground state

Muonic hydrogen is the only other hydrogen-like atom in which the hyperfine splitting of the ground state could be measured with high precision. Due to the large muon mass $`m_\mu /m_e\mathrm{2\hspace{0.17em}10}^2`$, the binding energy of the ground state of muonic hydrogen is of the order of 200 Ry, and the radius of the muon orbit is $`a_0/200`$ so that the energy levels of muonic hydrogen are more “sensitive” to the details of the proton structure than the levels of normal hydrogen. The expressions in Eqs. 2,3 and 5 for the hyperfine splitting of a hydrogen-like atom and the various contributions to it apply for muonic hydrogen as well; however, the energy scale and the relative size of the various terms differ significantly from the hydrogen case (see the rightmost two columns of Table 1.) The Fermi splitting now is 183 meV that corresponds to a hyperfine transition wavelength of 6.1 $`\mu `$m. The explicit form of the higher order terms in $`\delta ^{\mathrm{QED}}`$ may differ from Eq. 4 since different momenta are expected to give the main contribution in loop integrations, while proton structure modifies the proton voerteces. We do not know of any published explicit expression of $`\delta ^{\mathrm{QED}}`$ for the ground state of muonic hydrogen. To our opinion the lack of interest in the topic is due to the unclear perspectives of an experimental verification of the theoretical results. Indeed, the analogous contributions to the Lamb shift of muonic hydrogen were evaluated with a very high accuracy pachucki-LSmH ; kinoshita1  as soon as the proposal for the experiment kottmann  was about to be put forward. We are therefore convinced that the uncertainty of $`\delta ^{\mathrm{QED}}`$ may be brought down to 0.1 ppm if necessary. The evaluation of the recoil terms $`\delta ^{\mathrm{recoil}}`$ may not be that easy since the mass ratio $`m_\mu /m_p0.11`$ is much larger compared to hydrogen, and terms of order $`O(\alpha (m_\mu /m_p)^n),n=2,3,4`$, are all expected to contribute by more than 1 ppm. Most appropriate in this case might be an essentially two-body approach based on the quasipotential equation todorov ; faustov-QP ; we assume that the theoretical uncertainty may be brought down below $`10^6`$ this way, and put 1 ppm in the rightmost column of Table 1. Since we do not know of any published result on $`\delta ^{\mathrm{recoil}}`$, we take as estimate of the magnitude of the recoil effects the formula for the leading recoil term in muonium of Ref. sapir  and get $`\delta ^{\mathrm{recoil}}(3\alpha /\pi )(m_\mu /m_p)\mathrm{ln}(m_\mu /m_p)\mathrm{1.7\; 10}^3`$. The proton polarizability correction has been evaluated using the same methods as for hydrogen: $`\delta ^{\mathrm{pol}}=(4.6\pm 0.8)\times 10^4`$ cherednikova . The term $`\delta ^{\mathrm{hvp}}`$ describing the hadron vacuum polarization was shown to contribute by approximately 20 ppm faustov-vp , unlike hydrogen where it does not give any considerable contribution in hydrogen.

With the content of Table 1 in mind, we are now ready to discuss the information that would be provided by measurements of the hyperfine splitting of the ground state of muonic hydrogen atoms.

The proton structure correction $`\delta ^{\mathrm{str}}`$ in muonic hydrogen is enhanced (compared to hydrogen) by a factor of $`\mathrm{2\; 10}^2`$. Therefore, a measurement of $`\mathrm{\Delta }E^{\mathrm{HFS}}`$ in $`(\mu ^{}p)_{1s}`$ can not be a good test of QED since QED effects are overshadowed by the proton structure corrections. Further on, in both hydrogen and muonic hydrogen, the proton structure corrections $`\delta ^{\mathrm{str}}`$ is dominated by two independent terms: the Zemach term $`\delta ^{\mathrm{rigid}}`$ and the polarizability term $`\delta ^{\mathrm{pol}}`$. While the Zemach term is directly related to a well defined physical parameter - the Zemach radius of the proton $`R_p`$ (see Eq. 5), $`\delta ^{\mathrm{pol}}`$ is expressed in terms of the form factors and polarized structure functions of the proton in an indirect and case-dependent way and is not associated with a single parameter. Compared to hydrogen, both these terms scale approximately as $`(m_\mu /m_e)`$. This all brings us to the conclusion that opposite to what was believed by some authors pla ; wein , the measurements of $`\mathrm{\Delta }E^{\mathrm{HFS}}`$ in hydrogen and muonic hydrogen atoms are not complementary in a sense which would let us extract the values of two universal parameters of the proton, characterizing its charge and magnetic distribution and polarizability. However, if assuming that all terms in the right-hand side of Eq. 7 are evaluated correctly by theory, these measurements may be regarded as repeated experimental determination of the Zemach radius of the proton. While the discutible point in this assumption is the credibility ot the theoretical evaluation of $`\delta ^{\mathrm{pol}}`$ (if neglecting $`\delta ^{\mathrm{hvp}}`$), the repeated measurements of $`R_p`$ in hydrogen and muonic hydrogen are the best way to verify it: compatible values of $`R_p`$ extracted from the hyperfine splitting in hydrogen and muonic hydrogen will confirm the reliability of the theoretical values of $`\delta ^{\mathrm{pol}}`$ and vice versa.

The accuracy of $`R_p`$ depends on the uncertainty of $`\delta ^{\mathrm{pol}}`$; a measurement of the hyperfine splitting of the ground state of muonic hydrogen based on the results of cherednikova  would give the value of $`R_p`$ accurate to 1%. This would be more precise than the value obtained in the previous section because of the smaller relative uncertainty of the theoretical uncertainty of $`\delta ^{\mathrm{pol}}`$ (see Table 1), however, things may change with the more refined theoretical results to come in the future. As already mentioned, such an accuracy would fairly allow to filter the numerous theoretical estimates of $`R_p`$ and detect a deviation of $`G_E/G_M`$ from 1 by distinguishing the values of $`R_p`$ obtained with and without account of the JLab experimental results gayou . It would be preferable for this purpose to have the value of $`R_p`$ accurate to 0.5% or better, that requires in turn that the theoretical uncertainty of $`\delta ^{\mathrm{pol}}`$ be brought below $`\mathrm{3\hspace{0.17em}10}^5`$ and that the experimental error of $`\mathrm{\Delta }E_{exp}^{\mathrm{HFS}}`$ not exceed 30 ppm.

The muonic hydrogen Lamb shift experiment, currently in progress at PSI kottmann , may provide at a later stage as a by-product the hyperfine splitting of the $`2S`$-state of the muonic hydrogen atom with a relative accuracy of the order of $`\mathrm{0.5\hspace{0.17em}10}^3`$. According to pachucki-LSmH , the Zemach correction $`\delta ^{\mathrm{rigid}}`$ for the $`2S`$-state is again $`\mathrm{8\hspace{0.17em}10}^3`$. Unless the (yet uncalculated) polarization correction $`\delta ^{\mathrm{pol}}`$ happens to be anomalously large in this case, this measurement would therefore provide the value of $`R_p`$ with an accuracy of 5-10% kottmann  which is below the expected accuracy of the measurements in the 1S state. The significant improvement of the accuracy of the proton rms charge radius expected from PSI experiment will not help imcreasing the accuracy of $`R_p`$ either because $`R_p`$ is independent of the proton rms charge radius. In order to determine the Zemach radius of the proton from the hyperfine splitting of muonic hydrogen atoms in the $`1s`$ or $`2s`$ states with accuracy 1% or better, the experimental error of $`\mathrm{\Delta }E_{\mathrm{exp}}^{\mathrm{HFS}}`$ should not exceed 50 ppm – a requirement that is not met by the PSI Lamb shift experiment kottmann . An alternative experimental method, based on the response of the muon transfer rate from hydrogen to oxygen to the population of the $`(\mu ^{}p)_{1S}`$ para-state werth , was proposed recently adamczak . The method takes advantage of the recent progress in the development of tunable lasers in the far infrared range around 6.1 $`\mu `$m kottmann-2000 . The efficiency of the method has been demonstrated by means of Monte Carlo simulations. The experimental error limits have not been discussed; the main source of experimental uncertainty is expected to be the Doppler broadening of the transition lines.

## IV Hyperfine splitting of muonic hydrogen molecular ion

In conclusion, we would like to briefly outline an alternative possibility for determining the proton Zemach radius from a measurement of the hyperfine splitting of the muonic molecular ion $`p\mu ^{}p`$. The hydrogen muonic molecular ions $`(p\mu ^{}p)`$ are formed in the excited ortho-state with orbital momentum $`J=1`$ in collisions of muonic hydrogen atoms with hydrogen molecules ponomarev :

$$(\mu ^{}p)+H_2[(p\mu p)_Jpee]^{}$$
(9)
The formation rate $`\lambda _{p\mu p}`$ is proportional to the hydrogen target number density $`\phi `$: $`\lambda _{p\mu p}=\phi /\phi _0\times \mathrm{2.2\hspace{0.17em}10}^6\mathrm{s}^1`$ npa384 . At high densities $`\lambda _{p\mu p}`$ exceeds the muon weak decay rate $`\lambda _0=\mathrm{0.45\hspace{0.17em}10}^6\mathrm{s}^1`$ so that the muons spend most of their life time bound in muonic molecular ions. The spin-orbit and spin-spin interactions of the protons and the muon split the ortho-state with $`J=1`$ into 5 hyperfine state labeled with the quantun numbers $`(sF)`$ of the total spin $`\text{s}=\text{s}_{p1}+\text{s}_{p2}+\text{s}_\mu `$ and the total angular momentum $`\text{F}=\text{s}+\text{J}`$ npa384 ; aissing ; the para-state $`J=0`$ has no hyperfine structure. The energy separation between the ortho-levels with $`s`$=1/2 and $`s`$=3/2 is $``$138 meV, while the separation within each of these groups is two orders of magnitude smaller.

A measurement of the splitting between the hyperfine states with $`s`$=1/2 and $`s`$=3/2 might be based on the strong dependence of the rate of spontaneous ortho-para transitions $`\lambda _{op}=\lambda _{op}(s)`$ on $`s`$: $`\lambda _{op}(1/2)=\mathrm{7.2\hspace{0.17em}10}^4\mathrm{s}^1`$, $`\lambda _{op}(3/2)=\mathrm{0.8\hspace{0.17em}10}^4\mathrm{s}^1`$. Since only the ortho-states with $`s`$=1/2 are initially populated at high densities, the observable ortho-para transition rate is close to $`\lambda _{op}(1/2)`$. The idea of the experimental method would be to use a tunable laser to stimulate transitions from the $`s`$=1/2 to the $`s`$=3/2 hyperfine states, for which $`\lambda _{op}(3/2)`$ is an order of magnitude smaller, which would result in a resonance drop of the observable rate of ortho-para transitions. Though such an experiment would require the development of tunable far IR narrow band lasers in the range of 8.8 $`\mu `$m, it has some significant advantages compared with measurement in gaseous hydrogen at room tempertures – suppressed Doppler broadening, small muon stopping volume, high energy density of the laser beam, high rate of laser- stimulated ortho-para transitions – which deserve a careful consideration in future. As for the theory, the accuracy of the currently available results on the hyperfine structure of $`(p\mu ^{}p)_{J=1}`$ is limited to the leading order Breit and Zemach corrections and to the muon anomalous magnetic moment yaf ; aissing ; therefore, the theoretical uncertainty is currently about $`10^4`$. In order to determine $`R_p`$ to 1% from $`p\mu ^{}p`$ hyperfine splitting measurements, next order QED effects have to be taken into account together with proton polarizability and molecular ion finite size effects melezhik ; shima  so that the theoretical uncertainty be reduced by an order of magnitude to meet an experimental error below 50 ppm. This is a challenging task since $`(p\mu ^{}p)`$ is a bound system of particles with comparable masses in which three-body relativistic dynamics may show up yet in the next-to-leading order.

## V Conclusions

By assuming that the numerical results for the various terms in the theoretical expression for the hyperfine splitting of the ground state of the hydrogen atom are correct (within the limits of the claimed accuracy), we have determined the value of the proton Zemach radius $`R_p`$. We have also demonstrated that comparison of the experimental value of $`R_p`$ with theoretical calculation is a sensitive test of the quality of the approximation of the proton form factors at low momentum transfer and of the possible deviation of the ratio of the charge and magnetic form factors from 1. To verify the credibility of the theoretical evaluation of proton polarizability effects, we consider the possibility of measuring the hyperfine splitting of the ground state of muonic hydrogen atoms or, in a remote perspective, of muonic hydrogen molecular ions.

## VI Acknowledgment

The work on the paper was supported by NATO collaborative linkage grant PST.CLG.978454. D.B. also recognizes the partial support of grant Phy-1001 of the Bulgarian National Fund for Scientific Research. The authors are grateful to A.P. Martynenko and A. Le Padellec for helpful discussions.


text:
# HIGH-ENERGY EMISSION FROM PULSARS: THE POLAR CAP SCENARIO

## INTRODUCTION

Pulsar astronomy is at a fascinating juncture with many exciting new results from the Chandra and XMM X-ray telescopes, and the prospects of next generation gamma-ray experiments, led by space missions INTEGRAL (just launched), GLAST and AGILE, and a number of atmospheric Cerenkov telescopes coming on line in the next few years. These will enable a continuation of ground-breaking discoveries. Pulsars are expected to be detected by GLAST (Gamma-Ray Large Area Space Telescope: http://www-glast.stanford.edu) in profusion, some that are radio-selected, like most of the present EGRET/Comptel pulsars, and perhaps even more that are detected via independent pulsation searches. This theoretical review summarizes relevant characteristics of the polar cap model, emphasizing distinctions from the competing outer gap model. These features include acceleration properties, the X-ray to gamma-ray spectral shape, high energy cutoffs, pulse profiles and flux observabilities in different wavebands, and how these characteristics generally depend on pulsar period and period derivative. The polar cap scenario exhibits definitive signatures that will be readily tested by the detections of GLAST and other experiments, thereby establishing cogent observational diagnostics.

This polar cap review complements the outer gap review by Cheng (2003). The key distinguishing characteristic of the two models that is usually cited is that the acceleration zone for the polar cap model is confined to within a few stellar radii of the pulsar surface, while the region of acceleration in the outer gap model is proximate to the light cylinder. However, in terms of the physical manifestations that are most directly probed by astronomical observations, the presence or absence of an extremely strong magnetic field is the conspicuous feature. In polar cap models, the strong field permits single photon pair creation that attenuates super-GeV photons in Crab and Vela-like pulsars, whereas pair creation in outer gap models is mediated through the more familiar two-photon process involving surface thermal X-rays as targets. It is a principal contention of this paper that the two models will not generate coincident predictions for a large population of pulsars (contrasting isolated sources) that sample a significant range of periods and period derivatives. Since such model discrimination will probably be rendered by the gamma-ray observations of GLAST, here the focus is on magnetospheric cascade spectral properties. The reader is referred to the review by Michel (2003), and Mestel (2000) for discussions of pulsar electrodynamics. Recent polar cap model reviews include those of Harding (2001) and Rudak, et al. (2002). The issues of radio pulsar death lines and radio quiescence at high fields are beyond the scope of this paper; discussions can be found in Baring (2001a) and Zhang (2003).

## POLAR CAP MODEL BASICS

The critical physics ingredient for the polar cap model is the presence of a strong magnetic field in the acceleration and cascade emission region. Such high fields follow from the contention that the induced electric fields that seed particle acceleration in the oblique rotators exists at low altitudes near the neutron star surface at the magnetic poles. Since the earliest polar cap pulsar models of Sturrock (1971) and Ruderman and Sutherland (1975; hereafter RS75), there has been a potpourri of variations and updates, with the primary division being whether or not there is free emission of particles from the neutron star surface. Whether the surface temperature $`T`$ of the neutron star exceeds the ion, $`T_\mathrm{i}`$ and electron, $`T_\mathrm{e}`$, thermal emission temperatures controls the nature of the acceleration zone.

If $`T<T_\mathrm{i}`$, ions will be trapped in the neutron star crust (RS75, Usov and Melrose 1995) and a vacuum gap will develop at the surface forming the locale of the region of particle acceleration and radiation. If $`T>T_\mathrm{e}`$, free emission of particles of either sign of charge will occur. The flow of particles is then limited only by space charge, and since such particle flow all along each open field line is unable to supply the corotation charge that is required to short out the electric field component $`E_{}`$ along the magnetic field lines, an accelerating potential will develop (Arons and Scharlemann, 1979; Muslimov and Tsygan, 1992). In space charge-limited flow (SCLF) models, the accelerating $`E_{}`$ is screened at a height where the particles radiate $`\gamma `$-rays that produce pairs. This so-called pair formation front (e.g. Arons, 1983; Harding and Muslimov, 1998; Muslimov and Harding, 2003) can occur at high altitudes above the polar cap, depending on the colatitude of the field line (discussed below in the section on Slot Gaps), a property that may prove necessary to explain the spectral cutoffs in the some or most of the EGRET pulsars. Many surface temperatures have now been measured for canonical X-ray pulsars in the range $`T10^510^6`$ K (Becker and Trümper, 1997), though higher values are obtained ($`T4\times 10^67\times 10^6`$ K: see Perna et al., 2001) in observations of anomalous X-ray pulsars and soft gamma repeaters, so both vacuum gaps and SCLF models need to be considered depending on the source.

In the strong electric fields, the acceleration of primary electrons is rapid and ceases when one of two types of radiative cooling becomes significant, thereby establishing the maximum Lorentz factor $`\gamma _e`$ of these particles, and generating a quasi-monoenergetic primary distribution prior to cascading. The cooling mechanisms are curvature radiation (present in models from the earliest days of pulsar theory; e.g. Sturrock, 1971) induced by the non-uniform magnetic field, and resonant (magnetic) inverse Compton scattering of thermal X-rays from the stellar surface (e.g. Sturner and Dermer, 1994), a more recent incorporation. Both are strong functions of the magnetic field strength and either the electron’s Lorentz factor or the field geometry. Curvature radiation-initiated cascades generally have $`\gamma _e10^7`$ (e.g. Daugherty and Harding, 1989; see also Harding and Muslimov, 1998), while inverse-Compton seeded pair cascades yield $`\gamma _e3\times 10^5`$$`\mathrm{\hspace{0.17em}10}^6`$ (e.g. Sturner, 1995; see also Harding and Muslimov, 1998).

The primary photons propagate through the magnetosphere until they achieve sufficient angles $`\theta _{\mathrm{kB}}`$ with respect to the magnetic field to permit the creation of pairs via $`\gamma e^+e^{}`$ above the threshold energy of $`\mathrm{\hspace{0.17em}2}m_ec^2/\mathrm{sin}\theta _{\mathrm{kB}}`$. This propagation is impacted by general relativistic distortions of photon trajectories and field structure (e.g. Gonthier and Harding, 1994; Harding, et al., 1997), as is the magnitude of the field in the local inertial frame, so that curved spacetime properties significantly modify the rates of pair creation. For small polar cap sizes, corresponding to longer pulsar periods, the primary photons fail to acquire sufficient angles $`\theta _{\mathrm{kB}}`$ at low to moderate altitudes prior to the decline of the dipole field, thereby permitting the photons to escape unattenuated; pair creation is quenched at high altitudes since the rate is a strongly increasing function of $`B\mathrm{sin}\theta _{\mathrm{kB}}`$ (e.g. Tsai and Erber, 1974). It is this effect that is primarily responsible for the existence of a theoretical death line for radio pulsars (Sturrock, et al. 1976) at longer periods. The inability of simpler curvature radiation-seeded cascades to account for emission from pulsars of the longest periods has lead to recent refinements (Zhang, et al., 2000; Hibschman and Arons, 2001; Harding, et al., 2002) that incorporate the influence of non-resonant inverse Compton scattering of surface X-rays in seeding cascades and subsequent pair creation for significantly smaller polar cap sizes and small period derivatives.

The first generation of pair creation initiates the pair cascade, with pairs generally being created in excited transverse (to the field) momentum states, the so-called Landau levels. Rapid de-excitation via cyclotron and synchrotron radiation then follows (at least in pulsars with Crab-like or lower fields), and these secondary photons travel to higher altitudes and create further pairs and successive generations of photons in a cascade. The typical number of generations is around 3–4, and the total number of pairs per primary electron being of the order of $`\mathrm{\hspace{0.17em}10}^3`$$`\mathrm{\hspace{0.17em}10}^4`$ (e.g. Daugherty and Harding, 1982). The cumulative product is an emission spectrum that comprises a curvature/inverse Compton continuum that is cut off at hard gamma-ray energies by pair creation, with the addition of several synchrotron components at successively lower energies, terminating when the magnetosphere becomes transparent to $`\gamma e^+e^{}`$ at significant altitudes.

A notable exception to this cascade scenario arises in highly-magnetized pulsars such as PSR 1509-58, when the surface polar field $`B_0`$ exceeds around $`\mathrm{\hspace{0.17em}6}\times 10^{12}`$Gauss. Pairs are then produced in the zeroth (ground state) Landau level (e.g. Baring and Harding, 2001; hereafter BH01), so that cyclotron/synchrotron emission is prohibited. Cascading is then effectively quenched at low altitudes and the pair yield diminished. While Zhang and Harding (2000) suggested a possible alleviation of this cascade suppression via Landau level excitation of higher generation pairs by Compton scatterings with surface X-rays from the surface, Baring and Harding (2001) determined that the population of excited Landau states relative to that in the ground state is small, since such excitation can only achieved for photon energies exceeding the cyclotron resonance (e.g. Gonthier et al., 2000).

Another profound alteration to the nature of cascades arises in high field pulsars, due to the action of magnetic photon splitting, $`\gamma \gamma \gamma `$, a third-order quantum electrodynamical process in which a single photon splits into two lower-energy photons (e.g. Adler, 1971; Baring and Harding, 1997). The rate of splitting, like that of magnetic pair creation $`\gamma e^+e^{}`$, is generally a rapidly increasing function of field strength (for fields $`B\begin{array}{c}<\hfill \\ \hfill \end{array}10^{14}`$Gauss), photon energy and photon propagation angle with respect to the field. Since splitting possesses no energy threshold, it can dominate the first order process of pair creation if B is sufficiently high, typically above $`10^{13}`$Gauss (Harding, et al., 1997; hereafter HBG97). This leads to an alternative channel for cascade quenching, with gamma-rays being reprocessed without yielding pairs so that synchrotron generations are suppressed. The operation of splitting produces distinctive spectral bumps and polarization signals in the EGRET/Comptel band (HBG97); splitting-influenced pulsar cascades are also addressed in depth in Baring and Harding (2001).

## SPECTRAL AND TEMPORAL SIGNATURES

Spectral and temporal properties can define the most potent means for distinguishing between polar cap and outer gap scenarios. Yet such diagnostics are limited in power because of the significant number of parameters in each of the models. This degree of flexibility renders each model capable of approximately describing the major characteristics presented in a given dataset on an individual pulsar: a classic example is the comparison of polar cap models (e.g. Daugherty and Harding, 1996; Dyks and Rudak, 2000) with outer gap predictions (Romani, 1996) for the Vela pulsar. Such a predictive redundancy of two disparate models for individual sources is a limitation that is lifted when global characteristics of the pulsar population are considered. This provides motive for emphasizing global trends and properties for the entire pulsar population in this paper. The reasonable premise underpinning this tack is that the two scenarios cannot generate identical trends for luminosities, spectral indices, maximum gamma-ray energies, non-thermal X-ray indices, pulse profiles, polarization signatures, etc. as functions of period, period derivative, and viewing perspective for dozens or hundreds of well-measured pulsars. Since such large gamma-ray pulsar databases will be afforded by the GLAST mission, such an approach is both pertinent and timely.

### Trends and Distinguishing Characteristics

In the polar cap model, there are two distinguishing spectral features that serve as indicators of the field strength, namely structure sampling the cyclotron resonance that generally arises in the soft or hard X-ray bands for Crab-like and Vela-like pulsars, and the super-GeV cutoff due to attenuation by $`\gamma e^\pm `$ pair creation. In principal, these features are smeared out by distribution of the emission of magnetospheric altitudes and colatitudes; in practice this spectral structure degradation impacts the cyclotron feature much more than it does the high energy cutoff. Moreover, the cyclotronic structure is blueshifted by the minimum Lorentz factor of pairs generated in a synchrotron cascade: in most cases this further broadens the structure rendering field diagnostics difficult if not impossible. In addition, structure in the primary curvature spectrum due to cooling breaks can confuse the situation in the hard X-ray band.

Monte Carlo cascade simulations (e.g. Daugherty and Harding, 1982) generate minimum Lorentz factors $`\gamma _{\text{MIN}}`$ of around 50–100 for secondary pairs in models of Vela-like pulsars. The synchrotron photon energy for this Lorentz factor is $`\epsilon _{\text{MIN}}\gamma _{\text{MIN}}B/B_{\mathrm{cr}}`$, noting that a factor of $`\mathrm{\hspace{0.17em}1}/\gamma _{\text{MIN}}`$ is introduced to account for the cascade beaming angle. This energy, the energy of cascade cessation, is typically in hard X-rays for Crab-like pulsars, and generally $`\gamma _{\text{MIN}}`$ depends on $`B_0`$, $`P`$, etc. in more or less the same manner (Baring and Harding, 2000) that $`\epsilon _{\text{MAX}}`$ does in Eq. (1) below, since the same pair creation physics applies to both. The result is an $`\epsilon _{\text{MIN}}`$ with a fairly weak field dependence. Note, however, that the reprocessing that leads to the establishment of $`\epsilon _{\text{MIN}}`$ is initiated at slightly lower altitudes (corresponding to more cascade generations) than the ultimate attenuation that defines $`\epsilon _{\text{MAX}}`$, thereby complicating parametric specifications of $`\epsilon _{\text{MIN}}`$. The spectrum possesses a break at $`\epsilon _{\text{MIN}}`$, below which it assumes the flat $`\epsilon ^{2/3}`$ form that signifies curvature or synchrotron emission from quasi-monoenergetic pairs. For the Vela pulsar (see Pavlov et al., 2001 for a spectrum), when such a flat slope is extrapolated from the gamma-ray band, the model flux in the optical band grossly underpredicts that observed, implying that another component is present below X-ray wavelengths. Note that when the surface field becomes sufficiently high, ground state pair creation suppresses the synchrotron component and the $`\epsilon _{\text{MIN}}`$ feature disappears.

The second major spectral feature is the maximum energy of emission, which is controlled by attenuation due to magnetic pair creation $`\gamma e^\pm `$ during photon propagation through the pulsar magnetosphere. Such attenuation provides a characteristic super-exponential turnover (e.g. Daugherty and Harding, 1996) that contrasts that expected in outer gap models (e.g. see Thompson, 2001, for a comparison). Pair creation occurs at the threshold $`\epsilon \mathrm{sin}\theta _{\mathrm{kB}}=2`$ for high fields, i.e. $`B\begin{array}{c}>\hfill \\ \hfill \end{array}0.1B_{\mathrm{cr}}`$, and above threshold at $`\epsilon \mathrm{sin}\theta _{\mathrm{kB}}0.2B_{\mathrm{cr}}/B`$ for lower fields (e.g. see Daugherty and Harding, 1983). Here, $`\theta _{\mathrm{kB}}`$ is the angle of photon propagation relative to B, and hereafter photon energies $`\epsilon `$ are expressed in units of $`m_ec^2`$. Hence, the mean free path for photon attenuation in curved fields is $`\lambda _{\mathrm{pp}}\rho _c/\epsilon \mathrm{max}\{2,\mathrm{\hspace{0.33em}0.2}/B\}`$, i.e. when $`\epsilon \mathrm{sin}\theta _{\mathrm{kB}}`$ crosses above threshold during propagation. The radius of field curvature is $`\rho _c=[Prc/2\pi ]^{1/2}`$ for a pulsar period $`P`$. The approximate dependence of pair creation cutoff energies $`\epsilon _{\text{MAX}}`$ on $`B_0`$, $`R_0`$ and pulsar period $`P`$ (in seconds) can be summarized in the relation (Harding, 2001; Baring, 2001a)

$$\epsilon _{\text{MAX}}0.4\sqrt{P}\left(\frac{r}{R_0}\right)^{1/2}\mathrm{max}\{1,\frac{0.1B_{\mathrm{cr}}}{B_0}\left(\frac{r}{R_0}\right)^3\}\text{GeV}.$$
(1)
Accurate numerical determinations derived from the codes developed in HBG97 and Baring and Harding (2001), are plotted in Figure 1; these include the effects of general relativity on spacetime curvature, field enhancement and photon energy in non-rotating systems. At fields $`B_0\begin{array}{c}>\hfill \\ \hfill \end{array}0.7B_{\mathrm{cr}}`$ photon splitting acts to further reduce $`\epsilon _{\text{MAX}}`$, as discussed in Baring and Harding (2001); the operation of splitting proved necessary to account for the turnover inferred from Comptel data and EGRET upper limits to PSR 1509-58 (HBG97). For magnetars, pulsars with fields above $`\mathrm{\hspace{0.17em}4}\times 10^{13}`$Gauss, photon splitting and pair creation should prohibit any emission above $`100`$MeV, though prominent signals below 100 MeV are possible (Baring, 2001b) due to the efficiency of resonant Compton scattering.

There is clearly a strong anti-correlation between the maximum energy and the surface magnetic field, which seems to be augmented by an apparent decline of emission altitude with $`B_0`$. Such a trend is a distinctive characteristic that can be probed by GLAST; there appears to be no prediction of such a trend in outer gap models. The maximum energy is generally in the 1–10 GeV band for normal young pulsars, but can be much lower (e.g. HBG97; BH01) for highly magnetized ones, and also much higher for millisecond pulsars so that sub TeV-band (i.e. $`50`$–100 GeV) signals are possible (Bulik, et al., 2000) for polar cap models via synchrotron/curvature cascades if the field is low enough. It also should be remarked that the cutoff energy depends on pulse phase, with slightly greater values achieved between the pulse peaks in the case of Vela modelling (Daugherty and Harding, 1996); such a property matches the EGRET observations (Kanbach et al., 1994). Furthermore, there is an asymmetry in $`\epsilon _{\text{MAX}}`$ values between the peaks due to geometrical effects that is discussed in the next Section below. In passing it is noted that it is a generic property of polar cap models that small rotator obliquities ($`<45`$ degrees between rotation and dipole axes) are usually needed to simultaneously produce the observed $`\epsilon _{\text{MAX}}`$ values and double peaked pulse profiles of phase separation appropriate for gamma-ray pulsars (Harding and Daugherty, 1998).

The gamma-ray spectral index $`\alpha `$ extending below $`\epsilon _{\text{MAX}}`$ depends on the details of the cascading, and its behavior can be broadly summarized as follows. The primary photon index from cooling curvature radiation is $`\mathrm{\hspace{0.17em}5}/3`$, and should be realized in millisecond pulsars (Bulik, et al., 2000), where the fields are sufficiently low that synchrotron components contribute only below 1 MeV. For Vela-like fields, successive generations of cascading sequentially steepen the index (Wei, et al., 1997; Harding and Daugherty, 1998), saturating at $`\alpha =2`$. While most gamma-ray pulsars have $`\alpha <2`$ and can be simply described via such cascading, the Crab possesses a spectrum steeper than $`\alpha =2`$ that must depend on a spatial convolution in some subtle and obscure manner. When the local field in the emission region exceeds around $`\mathrm{\hspace{0.17em}6}\times 10^{12}`$Gauss, such as should be the case for PSR 1509-58, pair creation occurs only in the ground state and quenches cascading, yielding just a bare, flat curvature spectrum ($`\alpha =5/3`$). The contribution from resonant Compton scattering can also be significant in high field pulsars, depending on the proximity of the acceleration region to the surface and the surface X-ray temperature; its extremely flat spectrum is discussed below in the context of magnetars.

The presence of strong fields virtually guarantees a strong polarization signal in polar cap models, and when these couple with spectral structure and temporal information, particularly powerful observational diagnostics are possible. This may be fruitful at the lower end of the cascade continuum in Vela-like objects, but it is an especially valuable tool for highly-magnetized pulsars since the attenuation cutoffs fall in the Comptel band, and should exhibit strong and distinctive polarization signatures. Hard gamma-ray experiments like GLAST are generally not afforded the opportunity to act as polarimeters, being limited by multiple scattering in trackers above 300 MeV. Medium energy gamma-ray experiments, on the other hand, are ideally suited to polarization studies (via their sampling of Compton scattering kinematics). Gamma-ray polarimetry is no longer a distant dream, given the prospects (Lei, et al., 1997) that the recently-launched INTEGRAL mission will detect polarization at the 10% level from the Crab pulsar (at 200–600 keV), and also in a handful of other sources. Polarimetric capability in the hard X-ray and soft gamma-ray bands is a high priority for next-generation advanced Compton telescopes (e.g. see Kanbach et al., 2000, and the Web pages for the MEGA \[http://www.gamma.mpe-garching.mpg.de/MEGA/mega.html\] and ACT \[http://gamma.nrl.navy.mil/ngram/\] consortia).

### Observabilities

Spectral signatures provide the second layer of detail in observational diagnostics of gamma-ray pulsars; the primary layer is clearly the flux of a source at earth, i.e. its observability. This essentially represents the normalization of the model spectrum convolved with the pulsar’s distance. In recent years there has been an evolution in the understanding of what controls the detectability, making way for a further dimension of discrimination between polar cap and outer gap models in the GLAST era, when luminosity trends will not be crimped by small population statistics. The traditional approach of the EGRET community was to use the dipole spin-down power $`\dot{E}B_0^2/P^4`$ as an indicator of a pulsar’s observability, which is in fact mirrored by the population of field X-ray pulsars (Becker and Trümper, 1997); recent Chandra detections of pulsars in globular clusters indicate (Grindlay et al., 2002) a weaker luminosity dependence on $`B_0/P^2`$. While theoretically motivated (see just below), the spin-down luminosity choice did not match the trend subsequently established by EGRET: that there is clearly an almost linear correlation between the inferred luminosity $`L_\gamma `$ of gamma-ray pulsars and $`B_0/P^2`$, i.e. the voltage drop across the open field lines for standard polar caps. There is a very modest scatter in this correlation, largely due to the uncertainty in determining source distance by folding radio dispersion measures into the Taylor and Cordes (1993) galactic electron model (just recently updated in Cordes and Lazio, 2002).

It is interesting to note that such a linear correlation with $`B_0/P^2`$ of the inferred luminosity of gamma-ray pulsars was predicted (Harding, 1981) in the context of the polar cap model. At the time, only 2 gamma-ray pulsars were known, strengthening the impact of the polar cap model on the understanding of magnetospheric emission in such pulsars. The contention of a $`LB_0/P^2`$ dependence is based on the premises that the radiative luminosity is proportional to the Goldreich-Julian current, and that the cascade emission is initiated by pairs of Lorentz factor that is almost independent of pulsar $`B_0`$ and period. Subsequent predictions by competing analyses/models (e.g. Sturner and Dermer, 1994; Romani and Yadigaroglu, 1995; Cheng et al., 1998; Rudak and Dyks, 1999) and revisions (Zhang and Harding, 2000) all post-dated the EGRET database. The current status is that the polar cap expectations (Sturner and Dermer, 1994; Zhang and Harding, 2000) match the data slightly more accurately than their outer gap counterparts (Romani and Yadigaroglu, 1995; Cheng et al., 1998), with each group of researchers offering different $`B_0`$ and $`P`$ dependences for the luminosity (see Harding, 2001 for a review). This situation is presently limited by small number statistics, however in the GLAST era such correlations will be established on a firm basis.

Motivations for considering their luminosity dependence are not confined to model discrimination and refinement; assumed luminosity “laws” can dictate period selection in pulsation searches. This is a salient issue for GLAST, since it will be capable of blind period searches on gamma-ray sources with no radio counterparts. The period dependence is the most critical element to the $`L_\gamma (P,\dot{P})`$ relationship. While EGRET observed most pulsars high up on a $`B_0^2/P^4/d_{\text{PSR}}^2`$ rank-ordered list (where $`d_{\text{PSR}}`$ is the pulsar distance), certain gamma-ray pulsars (notably the longer period pulsars PSR 0656+14 and PSR 1055-52) are surprisingly low in spin-down luminosity, and millisecond pulsars proved extraordinarily difficult to detect, till the detection of PSR 0218+4232 (see Kuiper et al., 2000). Clearly, a gamma-ray luminosity dependence $`L_\gamma (P,\dot{P})`$ that differs from the spin-down one will dramatically modify the observability criterion, particularly if the period dependence is substantially different. Accurate determination of $`L_\gamma (P,\dot{P})`$ will influence the weight on different period ranges that will be applied to GLAST source data in pulsation searches. This will be particularly germane to cases where the source is near the galactic plane and pulsations are not evident at the frequencies of known radio pulsars within the GLAST source localization.

The spectral shape also affects the observability (Baring and Harding, 2000), a more subtle influence. This is a consequence of how the luminosity is distributed in the gamma-ray band, specifically that portion that emerges above the threshold sensitivity for a specific gamma-ray detector. The driving parameter for such an apportionment is the maximum emergent energy $`\epsilon _{\text{MAX}}`$, and to a lesser extent $`\epsilon _{\text{MIN}}`$ and the gamma-ray spectral index $`\alpha `$, since the spectra are generally flat enough for the bulk of the luminosity to emerge at the highest energies. These parameters control the normalization of the pulsar gamma-ray power-law for a given luminosity. An appropriate definition of a detector’s observability $`𝒪(\epsilon _{\text{TH}})`$ is the integral flux above an effective instrumental energy threshold $`\epsilon _{\text{TH}}`$. For pulsars with $`\epsilon _{\text{MIN}}\epsilon _{\text{TH}}`$, the usual case for GLAST considerations, this scales as the luminosity divided by the spectral normalization, yielding $`𝒪(\epsilon _{\text{TH}})L_\gamma (P,\dot{P})\epsilon _{\text{MAX}}^{(\alpha 2)}/d_{\text{PSR}}^2`$ (Baring and Harding, 2000; Baring, 2001a). Modest modifications to this dependence are possible, in particular if $`\epsilon _{\text{MIN}}\begin{array}{c}>\hfill \\ \hfill \end{array}\epsilon _{\text{TH}}`$ (see Baring 2001a for a discussion). The essential feature is that the observability is substantially different from that inferred from the spin-down formula, assuming that $`\epsilon _{\text{MAX}}`$ scales with $`B_0`$ and $`P`$ approximately as the low field alternative offered in Eq. (1). Accordingly, observabilities predicted for GLAST pulsation searches should follow a dependence (Baring, 2001a) somewhere in between $`B_0/P^2`$ and $`B_0^2/P^{5/2}`$. Using the latter possibility, Baring and Harding (2000) generated a revised rank-ordered listing that indicated a dramatic rearrangement from the traditional EGRET ordering. A depiction of this re-ordering in $`𝒪(\epsilon _{\text{TH}})`$ versus $`\dot{E}/d_{\text{PSR}}^2`$ space is given in Figure 2, where each axis can be used to define a rank ordering.

The limited scatter around a linear dependence in the figure indicates that there are generally only modest changes to the rankings for most pulsars when adopting the observability as an updated criterion for detectability. The most notable changes (see Baring, 2001a) included the much higher ranking of the “outlier” longer period gamma-ray pulsars PSR 0656+14 and PSR 1055-52, and the marked lowering of millisecond pulsars (PSR 1939+2134, PSR 0437-4715, PSR 1744-1134, etc.) in the ranks, specifically out of the top 40. Both of these reflect the weaker dependence of the revised observability on $`P`$. These refinements of rank-orderings mute questions of why PSR 1055-52 was seen by EGRET. Further revisions are in progress, and largely focus on the details of spectral shape. Only modest influences are expected from new pulsar distance determinations resulting from the updated model for the galactic electron distribution, NE2001 (Cordes and Lazio, 2002).

## RECENT DEVELOPMENTS

## Slot Gaps

During the Compton Gamma-Ray Observatory era, polar cap modelling focussed mostly on magnetospheric emission processes in order to take advantage of the rapid increase in gamma-ray pulsar data afforded by EGRET. In the last few years, considerable effort has been invested by Harding and Muslimov in the exploration of the acceleration region in order to determine the altitudes and particle energies associated with acceleration in the polar cap potentials (e.g. see Harding and Muslimov 1998, Muslimov and Harding 2003). Such details are essential to forming a self-consistent polar cap model for outward cascade emission as well as for polar cap reheating and subsequent X-ray reprocessing of the cascade power. In particular, a major question was raised by the modeling of Vela data: Daugherty and Harding (1996) ascertained that high altitudes ($`23`$ stellar radii) for the emission region were required to simultaneously explain the pulse phase separation and the $`5`$GeV maximum energy of emission. The altitude of the pair formation front (PFF) was a free parameter in their model, and a concerted effort in studying pulsar electrodynamics was needed to ascertain whether high altitudes could naturally be expected.

The recent work of Muslimov and Harding (2003) has affirmed the altitude choice of Daugherty and Harding (1996). In an extension of their sequence of papers on PFF development and location under the influence of general relativistic frame-dragging effects, they have explored colatitudinal variations of the PFF height. Their principal result is that there is a large variation in height with colatitude $`\theta `$ along the surface of the polar cap. This arises because the accelerating potential drops smoothly to zero at the rim of the polar cap, thereby extending the acceleration region to high altitudes in the vicinity of the rim. Pair quenching of the smaller electric field is precipitated more gradually due to the prolongation of acceleration in this confined region, referred to as the slot gap (after Arons, 1983). The result is that the dominant emission region spans a range of altitudes from relatively near the surface at the magnetic axis up to several stellar radii above the surface in the vicinity of the rim. Accordingly, hard X-ray and gamma-ray pulsar emission contains both core and conal components. In addition, the broad distribution of emission altitudes could aid smoothing of the super-exponential turnovers imposed by pair creation attenuation near $`\epsilon _{\text{MAX}}`$. The angular width of the slot gap scales approximately as $`(P/B)^{1/2}`$ for $`B_0\begin{array}{c}<\hfill \\ \hfill \end{array}4\times 10^{12}`$Gauss. When this is incorporated in solid angle modifications to the emergent gamma-ray luminosity, the slot gap model of Muslimov and Harding (2003) impressively accounts for the inferred luminosities in most cases; Geminga and the millisecond pulsar PSR 0218+4232 provide exceptions that require refinements to the model.

## Pulse Asymmetry

The GLAST experiment will afford phase-resolved spectroscopy at unprecedented statistical significance. It will provide well-defined pulse profiles in much smaller energy bins than were possible with EGRET or Comptel. Such developments enable new diagnostic capabilities. One effect that can be probed by this advance is that of pulse asymmetry, which has been studied in detail by Dyks and Rudak (2002). As the magnetosphere rotates, the photon trajectories slip across the curved field lines slightly. This relative motion differs between the leading and trailing rims of the polar cap, translating to different shapes for the two peaks in the pulse profile. Since pair creation rates are critically sensitive to the angle $`\theta _{\mathrm{kB}}`$ between the photon path and the local field line, one expects significant correlation of pulse asymmetries with photon energy. Dyks and Rudak find that near the pair creation turnover at $`\epsilon _{\text{MAX}}`$, the trailing peak dominates since $`\theta _{\mathrm{kB}}`$ is larger on average for the leading rim. Equivalently, $`\epsilon _{\text{MAX}}`$ is lower for the leading pulse. Concomitantly, the pair reprocessing is enhanced for the leading peak so that it becomes the more prominent of the two peaks at sub-GeV energies. While $`\epsilon _{\text{MAX}}`$ is an increasing function of $`P`$ for the leading rim, non-monotonic dependence of $`\epsilon _{\text{MAX}}`$ on period is exhibited by the trailing rim, providing constraints that will aid isolation of the rotator obliquity for an assumed polar cap size. Note that although Dyks and Rudak explored these asymmetry effects in flat spacetime, they are also present in the curved spacetime magnetospheric models of Muslimov and Harding (2003). Observe also that this rotational effect is a strong function of pulse period, being much more pronounced for fast rotators, as expected. Fortunately, these are the brighter portion of pulsars so such distinctive properties will provide powerful probes of polar cap models in the GLAST era.

## POWERFUL OBSERVATIONAL DIAGNOSTICS

In conclusion, this paper has identified key properties that the polar cap (PC) model has that are dependent on high $`B`$ physics, and that are palpably distinguishable from outer gap model characteristics; these can be summarized as follows. In the PC model there is no pulsed TeV emission in garden-variety pulsars, though millisecond pulsars should exhibit sub-TeV emission. The attenuation of gamma-rays due to one photon pair creation generates super-exponential cutoffs in the 10 MeV – 10 GeV band when the surface polar field is in the range $`\mathrm{\hspace{0.17em}10}^{11}`$Gauss $`\begin{array}{c}<\hfill \\ \hfill \end{array}B_0\begin{array}{c}<\hfill \\ \hfill \end{array}4\times 10^{13}`$Gauss; these have been posited as a key pulsar diagnostic for GLAST. Moreover, $`\epsilon _{\text{MAX}}`$ should decline with increasing $`B_0`$ if $`B_0\begin{array}{c}<\hfill \\ \hfill \end{array}10^{14}`$Gauss. At the same time, an increase of $`\epsilon _{\text{MAX}}`$ with emission altitude implies a correlation between $`\epsilon _{\text{MAX}}`$ and pulse separation. Recent work has also suggested that there should exist a significant pulse asymmetry near $`\epsilon _{\text{MAX}}`$ and for lower cascade photon energies. In the special case of magnetars, pulsars with fields above $`\mathrm{\hspace{0.17em}4}\times 10^{13}`$Gauss, attenuation dues to photon splitting and pair creation should prohibit any emission above $`100`$MeV, though prominent sub-100 MeV signals are possible due to the efficiency of resonant Compton scattering, and can be probed by GLAST. The hard X-ray/soft $`\gamma `$-ray spectral slope is expected to steepen slightly as $`B_0`$ rises when $`\mathrm{\hspace{0.17em}3}\times 10^{11}`$Gauss $`\begin{array}{c}<\hfill \\ \hfill \end{array}B_0\begin{array}{c}<\hfill \\ \hfill \end{array}10^{13}`$Gauss, though flat curvature radiation spectra with $`\alpha =5/3`$ are anticipated for both highly-magnetized and millisecond pulsars, i.e. at both ends of the pulsar field range. Most of these features will be probed by new gamma-ray missions INTEGRAL, AGILE, GLAST and ground-based air Čerenkov telescopes (HESS, Veritas, MAGIC, CANGAROO-III), as well as current X-ray missions CHANDRA, XMM, and RXTE. In addition, strong and distinctive polarization signatures are expected in pulsar spectra, which may be explored for the first time by INTEGRAL; positive detections anticipated for the Crab pulsar would break new ground in for the astrophysics of gamma-ray pulsars.

E-mail address of Matthew G. Baring: baring@rice.edu