format: text

text:
# High Resolution X-ray Spectroscopy and the Nature of Warm Absorbers in AGN

## 1 The Unanswered Questions about Warm Absorbers: How Chandra and XMM Can Help

Although absorption features due to highly-ionized species were first clearly seen in AGN spectra taken with ROSAT almost ten years ago (Turner et al. 1993), most of the fundamental facts about them remain obscure today. This is true at even the most basic, empirical level.

### 1.1 The ionization distribution

For example, we do not even know which ionization states are present in the absorbing gas. Because increasing levels of ionization tend to move atomic features higher in energy and then remove them altogether, the character of ionization states possible to see in X-ray spectra is strongly influenced by the specifics of instrumentation. Consequently, our efforts to determine how much matter there is in which ionization state suffer from extremely strong instrumental selection effects.

In ROSAT-data, the spectral resolution and signal/noise were so poor that the only feature strong enough to be seen was the O K-edge, and the ionization state of O could be determined only roughly. Labelling ionization level by the ionization parameter $`\xi L_{ion}/nr^2`$ (where $`L_{ion}`$ is conventionally defined as the luminosity between 1 and 1000 Ryd, $`n`$ is the H-nucleus density, and $`r`$ is the distance to the source of radiation), this sort of data was capable only of identifying the presence of matter with $`\xi <30`$, and perhaps distinguishing ionization parameters at the highest end of this range from those an order of magnitude or more smaller.

Unfortunately, ROSAT was was totally incapable of detecting any matter more strongly ionized because higer ionization destroys OVIII. Other elements’ K-edges could not be seen in ROSAT data either because their energy was outside its limited range (e.g., Fe) or because their abundance was too small to permit detection with such crude data. To put this (semi)-quantitatively, in gas with solar abundances most of whose O has not been stripped, O K-edges were detectable by ROSATwhen the column density $`N>10^{21}`$ cm-2. Even if the ROSAT energy range extended to the K-edge of ionized Fe at 7–9 keV, a column density almost three orders of magnitude larger would have been necessary in order for a significant edge to have been formed. In this sense, ROSAT-level data was roughly 1000 times less sensitive to highly-ionized than to weakly-ionized gas.

Data of the sort produced by ASCA was better, but not greatly. Although its resolution and signal/noise were good enough to permit estimates of the depths of the OVII and OVIII edges separately (e.g., Reynolds et al. 1995), it could see no other edge features reliably. The next-deepest K-edges (Si, S, Ar, etc.) are so much smaller in opacity that they would be detectable in ASCA data only for much greater absorbing columns. In principle, ASCA’s nominal energy range included the Fe K-edge, but the effective area there was so small as to make searching for this edge impractical. Thus, like ROSAT, ASCA could detect matter with $`\xi <30`$, but could improve upon ROSAT by permitting $`\xi 30`$ to be distinguished from, say, $`\xi =3`$. Again like ROSAT, it was effectively blind to more highly-ionized matter.

The advent of the grating spectrometers on Chandra and XMM has already dramatically changed this frustrating situation, and will undoubtedly improve matters further in the near future. Their tremendously finer spectral resolution permits study of line transitions, not just edges, while their much greater signal/noise allows one to search for much weaker features. With data of this quality, it is possible to search for species expected to be common up to $`\xi 1000`$ even when the column density of that gas is as small as $`<10^{21}`$ cm-2—and we are starting to find them. Several different recent warm absorber studies employing grating data report that more than one ionization component is detected, with $`\xi `$ ranging anywhere from a few tens to $`1000`$ (Kaastra et al. 2000, Collinge et al. 2001, Kaspi et al. 2001).

Thus, in the next few years we may reasonably hope to make much progress on the most basic empirical question of warm absorber studies: measuring the distribution $`dN/d\xi `$.

### 1.2 Location and origin

After identifying the absorber’s ionization state, the next obvious question is where to place it in the AGN system. Proposed distances have ranged all the way from inside or near the broad-line region at 0.01 – 0.1 pc from the nucleus (Reynolds et al. 1995, George et al. 1998, Netzer et al. 2001) to many tens of pc away in the narrow-line region (Bottorff et al. 2000). Some have argued that there is absorbing material across this entire range of radii (Otani et al. 1996, Morales et al. 2000). Others (e.g., Krolik & Kriss 1995) have argued that the absorbing gas should be identified with the scattering region posited by Seyfert galaxy unification schemes, perhaps $`1`$ pc from the nucleus when $`L_{ion}10^{44}`$ erg s-1.

The state of affairs regarding the absorbing matter’s origin is even worse. A few of the diverse suggestions that have been floated include: evaporation off “bloated stars” in the broad-line region (Netzer 1996); gas evaporated off the torus obscuring the nucleus (Krolik & Kriss 1995); and a wind driven off the accretion disk (Elvis 2000; Bottorff et al. 2000). The absorbing gas’s ultimate fate has been left wholly unconsidered.

This, too, is a topic on which we may hope for enlightenment from new spectroscopy data. As has been remarked in numerous papers (Reynolds et al. 1995, Krolik & Kriss 2001, Netzer et al. 2001), variability in absorption features can strongly constrain the location of the gas responsible. The virtue of the new data sets is that, by making so many features observable, our odds of seeing variability in one or more of them are greatly improved. Moreover, because the new data are less subject to ionization selection effects (as discussed in the previous subsection), the location bounds placed by variability can be applied more generally.

### 1.3 Kinematics

Independent of where the absorbing gas is located, we’d also like to know which way it’s moving. For this issue, high-resolution spectroscopy is essential. If the absorbing gas possesses ions with ionization potentials $`1`$ keV (as is necessary to have numerous soft X-ray features), its temperature cannot be much greater than $`10^7`$ K. Transonic motions at this temperature would then be $`300`$ km s-1 or slower. It immediately follows that spectral resolution $`\mathrm{\Delta }E/E>1000`$ is a prerequisite for beginning to detect motions at this speed.

Because resolution of this caliber is exactly what Chandra and XMM have now delivered, we are beginning to learn about this question. Kaspi et al. (2001), for example, have shown that the absorber in NGC 3783 is moving toward us (relative to the nucleus) at $`500`$ km s-1.

### 1.4 Geometry

Although the depth of an edge may tell us the total column density of ions along the line of sight, it does not tell us whether they are all clumped together, or spread out smoothly, or clustered at a variety of locations. And, of course, absorption tells us nothing about the existence of matter along other directions.

Other constraints are necessary to make progress answering questions of this sort. For example, resonance line emission components offset toward the red from absorption components are the classic signature of outflows. The relative equivalent width of the emission and absorption components reveals the ratio between the mean optical depth on lines of sight other than ours to the ratio on the line toward us. Similarly, the relative velocity widths further constrain the gas’s geometric disposition. Making use of the emission/absorption equivalent width indicator, Kaspi et al. (2001) have already suggested that the absorbing gas in NGC 3783, while occupying a large part of solid angle around the nucleus, does not completely cover it.

Clumpiness in the gas can be analyzed in any of several ways, given quality spectroscopy: Emission lines that are either collisionally-excited or created by recombination yield measurements of the gas’s emission measure; combined with independent density estimates, one can then infer the total emitting volume and compare it to the available volume in the region. Alternatively, spectroscopic data that requires gas at several different ionization levels co-existing in more or less the same location would strongly suggest a clumpy, inhomogeneous environment.

### 1.5 Dynamics

Given our ignorance about the warm absorber’s kinematics, location, and origin, it should come as no surprise that almost nothing is known with confidence about its dynamics. We have as yet advanced hardly at all past listing the “usual suspects”: gravity, thermal presure gradients (Balsara & Krolik 1993), radiation pressure (Chelouche & Netzer 2001, Morales & Fabian 2002), and magnetic fields (Bottorff et al. 2000).

## 2 Warm Photoionized Equilibria

To make further progress, it is helpful to consider a few constraints imposed by the basic facts of photoionization physics. We begin with the most elementary of these, a relationship that links the ionization parameter, the distance, the gas column density, and its volume filling factor. From the definition of $`\xi `$ and the fact that the column density $`N=n\mathrm{\Delta }r`$, where $`\mathrm{\Delta }r`$ is the total distance along the line of sight occupied by the absorbing gas, we immediately see (as shown by Turner et al. 1993) that

$$r=\frac{L_{ion}}{N\xi }\frac{\mathrm{\Delta }r}{r}.$$
(1)
Because $`\mathrm{\Delta }r/r`$ cannot possibly be greater than unity, the inferred ionization parameter and column density of an absorber, when combined with the luminosity of the ionizing source, place an upper bound on its distance from that source:

$$rr_{max}\frac{L_{ion}}{N\xi }.$$
(2)
For typical numbers, $`r_{max}30L_{ion,44}N_{22}^1(\xi /100)^1`$ pc. Here we have scaled the ionizing luminosity to $`10^{44}`$ erg s-1 and the column density to $`10^{22}`$ cm-2. For fixed $`r_{max}`$, the more highly clumped the gas is, the smaller the radius at which it might be found.

A second constraint is provided by the nature of photoionization equilibrium. If the gas’s cooling time is shorter than the time required for it flow through the absorbing region, we can expect that it is in energy balance with the radiation. If that is the case, its temperature is a unique function of its ionization parameter $`\xi `$.

If we wish to link photoionization physics with fluid dynamics, it is more convenient to write the ionization parameter in a related, but significantly different form: $`\mathrm{\Xi }\xi /(4\pi ckT)`$. The factor of $`4\pi ck`$ makes $`\mathrm{\Xi }`$ renders it dimensionless; more significantly, by dividing by the temperature, a given $`\mathrm{\Xi }`$ at any particular radius corresponds to a certain pressure. Viewed from the opposite perspective, if the boundary condition is a fixed pressure and the gas is in energy balance, $`\mathrm{\Xi }`$ is determined, but $`T`$ is uniquely determined only when $`T(\mathrm{\Xi })`$ is single-valued. Although $`T(\xi )`$ is, in general, unique, when the ionizing continuum has a form typical of AGN, $`T(\mathrm{\Xi })`$, in general, is multi-valued in the range of temperatures between $`5\times 10^4`$ K and $`>10^6`$ K (Krolik, McKee & Tarter 1981; fig. 1). As a result, the boundary conditions may not suffice to determine the temperature; its history (and possibly the influence of thermal instabilities) is instead the final deciding factor. The fact that the thermal balance may be multi-valued is central to the study of warm absorbers because the multiple-valued solutions occur in exactly the range of conditions capable of making soft X-ray features.

When a detailed calculation is made of the thermal equilibrium curve, it is necessary to make several guesses about matters that we don’t truly know and yet have noticeable effect on the result. The distribution of elemental abundances is the first of these. Generally-speaking, solar abundances are assumed; other choices can alter the detailed shape of the curve, but within quite wide bounds have at most modest effect on the shape of $`T(\mathrm{\Xi })`$. More important to the thermal balance curve is the shape of the ionizing continuum. In low-redshift objects, we can observe the continuum at energies just below 1 Rydberg and again above $`500`$ eV, but the most important segment of the ionizing continuum is precisely in the band from 13 eV to $`500`$ eV that we can never directly see. The best we can do is smoothly interpolate from either side. Fortunately, as shown in Figure 1, different plausible guesses do not result in qualitatively different results. Both equilibrium calculations shown in the figure were based on continuum shapes meant to apply to the Seyfert galaxy NGC 3783 (one in Kaspi et al. 2001, the other in Krolik & Kriss 2001). Although rather different, both predict the same general sort of behavior—a nearly vertical rise in temperature from $`5\times 10^4`$ K to almost $`10^6`$ K. The value of $`\mathrm{\Xi }`$ where this takes place is called $`\mathrm{\Xi }_c`$.

In fact, it is worthwhile to focus greater attention on the “vertical branch” of the photoionized thermal equilibrium curve because it is significant for several reasons. The first, as just mentioned, is that the range of $`\xi `$ inferred from the ionization stages seen in warm absorber spectra is $`10`$ – 1000; this is almost exactly the range traversed on the vertical branch. In other words, if warm absorbers are in a state of photoionization equilibrium and thermal balance, they “live” on the vertical branch.

Second, the vertical branch has a rather special thermodynamic property: it is marginally stable to isobaric thermal perturbations. Were $`dT/d\mathrm{\Xi }>0`$, constant pressure perturbations would be stable because increasing the temperature from an initial point on the equilibrium curve must result in net cooling; conversely, were $`dT/d\mathrm{\Xi }<0`$, the same sort of perturbation grows exponentially because the excursion causes net heating. However, on the vertical branch $`d\mathrm{\Xi }/dT=0`$: in this case, if the temperature is perturbed while the pressure is fixed, the gas remains in a state of zero net radiative energy gain or loss. Consequently, there is no temperature “feedback”, either positive or negative. Perturbations are neither strongly damped nor strongly amplified—they simply do whatever the outside agency creating them demands. Although it is true that the curve of thermal balance can never be precisely vertical, to the extent that it is nearly so, the net cooling or heating that is engendered by the perturbations is small, and the magnitude of growth permitted in a region of weak instability is capped at a low level as a new segment of the equilibrium curve is quickly reached.

Third, far from being a special case, the value of $`\mathrm{\Xi }`$ at which the vertical branch occurs is actually a preferred value of the pressure-based ionization parameter. Whenever large amounts of cool material are exposed to a strong ionizing flux (as might readibly be imagined near an AGN), the low-density portion of this cool matter finds itself rapidly ionized and heated. If the heating time is short compared to the flow time, the ambient pressure rises equally rapidly. At fixed distance from the source of ionizing radiation, this means the value of $`\mathrm{\Xi }`$ surrounding the cool matter rapidly falls. Nothing stops this fall in $`\mathrm{\Xi }`$ until it reaches $`\mathrm{\Xi }_c`$, where the net energy exchange via radiation switches from heating to cooling. Similarly, if the pressure overshoots to the net-cooling side of the vertical branch, matter cools rapidly, losing pressure, until balance is achieved—at exactly $`\mathrm{\Xi }_c`$.

This portion of the argument may be summed up very simply: Given that, by definition, warm absorbers are exposed to AGN fluxes, it makes sense to suppose that they are photoionized. Let us suppose for the moment that the gas is in photoionization equilibrium and radiative balance. Then the particular ions observed (OVII,VIII; NeX; MgXII; SiXIV; etc.) indicate that the gas lies on the vertical branch of the $`T(\mathrm{\Xi })`$ curve. Because that branch allows a wide range of temperatures to coexist at a single pressure, it follows that the material is likely to contain multiple sub-regions whose temperature could be anywhere within the range found on the vertical branch of the equilibrium curve. Moreover, pressure balance tends to force photoionized gases in the presence of a copious mass source toward $`\mathrm{\Xi }_c`$. That is, it is a natural consequence of photoionization physics that we should see gas with $`\mathrm{\Xi }=\mathrm{\Xi }_c`$ and that it should be inhomogeneous in temperature.

## 3 Global Picture

These special thermodynamic properties of the vertical branch find a natural application in the context of AGN warm absorbers. These are frequently seen in type 1 Seyfert galaxies, a variety of AGN in which we have excellent reason to believe that extremely optically thick obscuring matter wraps around the nucleus with roughly toroidal geometry (Antonucci 1993, Krolik 1999). Where the ionizing radiation of the nucleus strikes the inner edge of the obscuring stuff, exactly the process described in the previous paragraph can be expected to take place: copious amounts of matter should be ionized and heated, maintaining $`\mathrm{\Xi }\mathrm{\Xi }_c`$. Until the temperature rises quite high, the heating time is far shorter than the flow time: as shown in detail in Krolik & Kriss (2001), if the flow speed is of the magnitude observed ($`500`$ km s-1), this criterion is met for $`T<4\times 10^6r_{pc}^{1/2}`$ K, for $`r_{pc}`$ the distance from the central source in parsecs.

As the newly-liberated material is heated, it will expand in order to maintain pressure balance. One might then ask, “Will the expanding parcels eventually be confined? Or will they spread out and fill the entire volume of the region?” Suppose that the matter we see by its “warm absorber” features is clumpy. If the volume outside the clumps has lower pressure, the clumps will expand at roughly their internal sound speeds. The ratio in column density between the low-pressure background and the high-pressure absorbing lumps may then be estimated as

$$N_{bkgd}/N_{abs}(r/\mathrm{\Delta }r_{abs})(c_s/v_{flow}),$$
(3)
where the absorber has radial thickness $`\mathrm{\Delta }r_{abs}`$ and sound speed $`c_s`$. Because the sound speed at $`10^5`$ K (the typical temperature inferred for warm absorber regions) is $`50`$ km s-1, it is clear that the column density of the background is automatically close to the column density of the absorber if there is strong clumping. That is, one must expect the background to be a significant part of the absorbing structure. Put yet another way, it is probably best to consider the entire system as volume-filling, albeit one that is likely to be strongly inhomogeneous.

Although the ionization parameter $`\mathrm{\Xi }`$ is tightly constrained by photoionized thermodynamics, the temperature, at least within the range permitted at $`\mathrm{\Xi }\mathrm{\Xi }_c`$, is hardly constrained at all. We should expect, then, that the evaporating matter in this region might exist anywhere in the temperature range $`5\times 10^4`$$`10^6`$ K, i.e., $`30<\xi <1000`$. The issue is not “What is the favored value of $`T`$ or $`\xi `$?”, but “What is the distribution of $`T`$ and $`\xi `$ within this broad range?” Small accidents having to do with irregularities in the source of matter (the obscuring torus?) or “bumps” in the outflow of heated gas could affect this distribution, so different objects could easily differ in detail.

At fixed pressure, the hottest regions have the least density and so, per unit mass, occupy the largest volume. For this reason, one might expect most of the volume to be occupied by relatively hot gas, but it is likely studded by smaller regions of higher density and lower temperature. Each of these sub-regions can easily change its temperature over time.

## 4 Implications

This conceptual picture has numerous implications for the character of warm absorbers. The most basic, of course, is that the features we see should originate in ions found in photoionized equilibrium over a span of a factor of 20 or so in temperature. There are also other, more quantitative implications.

### 4.1 Location

Returning to equations 1 and 2, we can rewrite them as

$$rr_{max}(\xi _{max})=\frac{L_{ion}}{\xi _{max}N(\xi _{max})},$$
(4)
where $`\xi _{max}`$ is the highest $`\xi `$ (corresponding to the highest $`T`$) attainable in a flow time. Because the gas at $`\xi _{max}`$ is the hottest that we can expect to find in the absorber region, it should occupy the largest part of the volume. That is, for the hottest gas $`\mathrm{\Delta }r/r1`$, while cooler gas is strongly clumped. Scaled in terms of typical numbers, the predicted distance

$$r3\frac{L_{ion,44}}{N_{22}}\text{ pc},$$
(5)
where $`\xi _{max}`$ has been set at 1000.

A few parsecs is, of course, the scale of the inner torus (as seen, for example, in near-infrared imaging—Thatte et al. 1997, Marco & Alloin 2000—or H2O maser spots—Greenhill et al. 1996). It is likewise the expected scale of the polarizing reflection region that permits us to see type 1 Seyfert nuclei in type 2 Seyfert galaxies.

### 4.2 Dynamics

With an estimate of the distance to the nucleus, the strength of gravity may be more quantitatively estimated. Making use of equation 5, we find that the free-fall speed in the absorber region is $`30(L/L_E)^{1/2}N_{22}^{1/2}`$ km s-1. Here $`L_E`$ is the Eddington luminosity and we have supposed that most of the bolometric luminosity is in the ionizing band. Because the sound speed in the hottest phase of the absorber is $`100`$ km s-1, we can expect that thermal pressures alone may be able to drive an outflow unless $`L/L_E1`$. Whether radiation or magnetic forces supplement the thermal pressure gradient remains an open question, but the fact that we see blue-shifted absorption lines indicating an outflow should now be no surprise.

### 4.3 Geometry

If the warm absorber is formed by this sort of photoionized evaporation off the inner edge of the obscuring torus, the torus itself blocks the warm gas from expanding to the side. A biconical expansion and outflow results. For this reason, we might expect the warm absorber to occupy fully only a part of solid angle around the nucleus. Moreover, depending on our viewing angle, parts of the obscuring matter might block our view of selected portions of the far side of the outflow. Thus, emission line equivalent widths should indicate substantial, but incomplete, solid angle coverage—as is tentatively seen in NGC 3783 (Kaspi et al. 2001).

### 4.4 Variability and location

This picture also makes specific predictions about variability in absorption features. The column densities we see are, of course, controlled by ionization balance in the line-of-sight gas. The governing equation can be written (in abridged form) as

$`{\displaystyle \frac{d\mathrm{ln}n_i}{dt}}={\displaystyle \frac{n_{i1}}{n_i}}F_x^{}\sigma _{ion}^{}+`$ $`n_e\alpha _{rec}^{\prime \prime }{\displaystyle \frac{n_{i+1}}{n_i}}+\mathrm{}`$ (6)
$`F_x\sigma _{ion}n_e\alpha _{rec}\mathrm{}`$ (7)
where the density of ionization stage $`i`$ is $`n_i`$, $`i\pm 1`$ correspond to the stages once more or less ionized, $`F_x`$ is the photon flux at the relevant ionization edge, $`\sigma _{ion}`$ is the appropriately averaged ionization cross section, and the primes and double primes denote the quantities relevant to the $`i+1`$ and $`i1`$ ionization stages. When $`n_{i+1}n_i`$ (as is often true for the particular ionization stages seen), two terms dominate the rest, so that

$$\frac{d\mathrm{ln}n_i}{dt}n_e\alpha _{rec}^{\prime \prime }n_{i+1}/n_iF_x\sigma _{ion}.$$
(8)
The characteristic time for effecting changes in the absorbing column is then

$$t_{ion}10^{6\pm 1}\frac{r_{pc}^2}{L_{ion,44}}\text{ s}.$$
(9)
Different species have varying ionization timescales, depending on the details of their atomic structures. For example, if one examines the more abundant species in the relevant ionization conditions, SiXIV, FeXXV, and FeXXVI tend to have longer ionization times, while OVI and the ions around FeXX tend to have somewhat shorter timescales. The range is not terribly large, however: the $`\pm 1`$ in the exponent in equation 9 more or less comprises the range for the species of greatest interest (Krolik & Kriss 2001).

Timescales of days to years are, of course, very convenient timescales for human observations. Because we know which species should vary on which timescales (parameterized by the ionizing flux in the absorbing region), monitoring spectroscopy should be a very powerful tool for constraining the location and conditions of the absorbing gas.

### 4.5 How would the warm absorber look when seen from the side?

Finally, it is striking how much the warm absorber region resembles the reflection region that allows us to see the nuclei of type 2 Seyfert galaxies. Both regions appear to lie at distances $``$ pc; both are highly-ionized; both are flowing outward at speeds of several hundred km s-1; the total column densities in both are $`10^{22}`$$`10^{23}`$ cm-2. It would appear to make a great deal of sense to suppose that they are the same structure, simply viewed at different angles.


text:
# Dynamics of collapsing and exploding Bose-Einstein condensed vortex state

## I Introduction

Recent observation 1 ; ex2  of Bose-Einstein condensates (BEC) of dilute trapped bosonic atoms with repulsive interaction at ultra-low temperature has intensified theoretical and experimental studies on various aspects of the condensate 11 . One fascinating feature is the observation of quantized vortices exp  and vortex lattice exp1 ; expm  in the condensate as this is intrinsically related to the existence of superfluidity. Another interesting feature is the formation of a stable condensate composed of a finite number of attractive atoms less than a critical number $`N_{\text{cr}}`$ ex2 . The third noteworthy feature is the observation of Feshbach resonances in 23Na fbna , 85Rb ex3  and Cs fbcs  atoms, as in the presence of such resonances the effective atomic interaction can be varied in a controlled fashion by an external (background) magnetic field fbth .

For superfluid 4He(II) in a rotating container, no motion of the fluid is observed below a critical rotational frequency. Above this frequency quantized vortices appear in 4He(II) manifesting its superfluidity. However, because of the strong interaction, the theoretical description of this system is not easy. Quantized vortices have been observed exp ; exp1 ; expm  in trapped BEC and can be generated in theoretical mean-field models 2a ; 2c ; 2d ; 2f ; 2ff ; 2g ; fc  based on the Gross-Pitaevskii (GP) 8  equation. Different ways for generating vortices in a BEC have been suggested 2d , e.g., via spontaneous formation in evaporative cooling 2f , via controlled excitation to an excited state 2ff , by stirring a BEC using a laser with an angular frequency above a critical value 2c , or by the rotation of an axially symmetric trap with an angular frequency above a similar critical value 2g . In contrast to liquid 4He(II), a trapped BEC of small size is dilute and weakly interacting, which makes a mean-field analysis appropriate.

The observation of a condensate of attractive 7Li atoms and the subsequent measurement of the critical number $`N_{\text{cr}}`$ ex2  is in good agreement with the mean field analyses in a spherically symmetric trap skxx , although the agreement is not as good in the case of 85Rb unnum  atoms in an axially symmetric trap sk1 . If the number of atoms can somehow be increased beyond this critical number, due to interatomic attraction the condensate collapses emitting atoms until the number of atoms is reduced below $`N_{\text{cr}}`$. With a supply of atoms from an external source the condensate grows again beyond $`N_{\text{cr}}`$ and a sequence of collapses has been observed in 7Li by Gerton et al at Rice ex2 , where the number of atoms remain close to $`N_{\text{cr}}`$ and the collapse is driven by a stochastic process.

Recently, a more challenging experiment was performed by Donley et al. ex4  at JILA on a BEC of 85Rb atoms ex3  in an axially symmetric trap, where they varied the interatomic interaction by an external magnetic field near a Feshbach resonance fbth . Consequently, they were able to change the sign of the atomic scattering length, thus transforming a repulsive condensate of 85Rb atoms into a collapsing and highly explosive attractive condensate and studied the dynamics of the same ex4 . Immediately after the jump in the scattering length, one has a highly unstable BEC, where the number of atoms could be much larger than $`N_{\text{cr}}`$. Donley et al. ex4  have provided a quantitative estimate of the explosion by measuring the number of atoms remaining in the BEC as a function of time until an equilibrium is reached. Because this phenomenon of emission of a very large number of atoms in a small interval of time is reminiscent of an explosion and looks very much like a tiny supernova, or exploding star, the researchers dubbed it a “Bosenova”.

The essential aspects of the above experiments at Rice by Gerton et al ex2  and at JILA by Donley et al ex4  have been theoretically described by a variety of authors using the GP equation 8 . The theoretical analyses have not only produced the time-independent results, such as, the critical number $`N_{\text{cr}}`$ skxx ; sk1 , but also time-dependent results, such as, the variation of number of atoms of the BEC during collapse and explosion th1 ; th1a ; th2 ; th3 ; th4 ; th41 , both in reasonable agreement with experiment. This consolidates the use of the mean-field GP equation in describing the dynamics of collapsing and exploding BEC of small to medium size. These BEC’s composed of several thousand atoms can be considered dilute and weakly interacting and hence amenable to mean-field treatment. Motivated by the above success, using the GP equation we propose the numerical simulation of the dynamics of a rotating Bosenova with a single vortex composed of a small number (several thousands) of 85Rb atoms as in the experiment at JILA ex4 . We consider a single vortex state exp ; 2a ; sk1 , as appropriate for small condensates, opposed to vortex lattice for large condensates exp1 ; expm ; fc . A comparison of the present results with future experiment will provide a more stringent test for the mean-field GP equation.

In this paper we perform a mean-field analysis based on the time-dependent GP equation to understand the collapse and explosion of the attractive vortex state of 85Rb atoms in an axially symmetric trap. To account for the loss of atoms from the strongly attractive condensate we include an absorptive nonlinear three-body recombination term in the GP equation. The three-body recombination rate we use in numerical simulation is the same as used in a similar study with BEC’s with zero angular momentum th4  and is in agreement with previous experimental measurement k3  and theoretical calculation esry . In the present investigation we consider the complete numerical solution of the mean-field GP equation for an axially symmetric trap as in the experiment at JILA ex4 . As in the experiment at JILA with nonrotating condensates, we find that, also in the case of rotating vortex states, a large number of atoms could be emitted in a small interval of time and one could have a Bosenova-type explosion.

Throughout the present numerical simulation we make the assumption that the axial symmetry of the system is maintained. For small values of nonlinearity a dynamical quadrupole instability may cause an attractive BEC vortex state to split into two pieces that rotate around the axial direction ref1 . These pieces may unite to recover the original vortex and this split-merge cycle repeats. Similar instability is known to exist for BEC in a toroidal trap ref2 . Clearly, a full three-dimensional calculation of the collapsing phenomena taking into consideration the effect of the splitting of vortex states seems to be practically impossible at present and will be a welcome future work. In view of this here we present an axisymmetric model of the same which is expected to provide the essentials of the collapse dynamics of the vortex states.

In Sec. II we present the theoretical model and the numerical method for its solution. In Sec, III we present our results. Finally, in Sec. IV we present a brief discussion and concluding remarks.

## II The Gross-Pitaevskii Equation

### II.1 Theoretical Model

The time-dependent Bose-Einstein condensate wave function $`\mathrm{\Psi }(𝐫;\tau )`$ at position $`𝐫`$ and time $`\tau `$ allowing for atomic loss may be described by the following mean-field nonlinear GP equation 11 ; 8

$`[`$ $``$ $`i\mathrm{}{\displaystyle \frac{}{\tau }}{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V(𝐫)+𝒢N|\mathrm{\Psi }(𝐫;\tau )|^2{\displaystyle \frac{i\mathrm{}}{2}}(K_2N`$ (1)
$`\times `$ $`|\mathrm{\Psi }(𝐫;\tau )|^2+K_3N^2|\mathrm{\Psi }(𝐫;\tau )|^4)]\mathrm{\Psi }(𝐫;\tau )=0.`$
Here $`m`$ is the mass and $`N`$ the number of atoms in the condensate, $`𝒢=4\pi \mathrm{}^2a/m`$ the strength of interatomic interaction, with $`a`$ the atomic scattering length. The terms $`K_2`$ and $`K_3`$ denote two-body dipolar and three-body recombination loss-rate coefficients, respectively. There are many ways to account for the loss mechanism th1 ; th1a . Here we simulate the atom loss via the most important quintic three-body term $`K_3`$ th1 ; th2 ; th3 . The contribution of the cubic two-body loss term k3  is expected to be negligible th1 ; th3  compared to the three-body term in the present problem of the collapsed condensate with large density and will not be considered here.

The trap potential with cylindrical symmetry may be written as $`V(𝐫)=\frac{1}{2}m\omega ^2(r^2+\lambda ^2z^2)`$ where $`\omega `$ is the angular frequency in the radial direction $`r`$ and $`\lambda \omega `$ that in the axial direction $`z`$, with $`\lambda `$ the aspect ratio. We are using the cylindrical coordinate system $`𝐫(r,\theta ,z)`$ with $`\theta `$ the azimuthal angle. The normalization condition of the wave function is $`𝑑𝐫|\mathrm{\Psi }(𝐫;\tau )|^2=1.`$

The GP equation can easily accommodate quantized vortex states with rotational motion of the BEC around the $`z`$ axis. In such a vortex the atoms flow with tangential velocity $`L\mathrm{}/(mr)`$ such that each atom has quantized angular momentum $`L\mathrm{}`$ along the $`z`$ axis. This corresponds to an angular dependence of

$$\mathrm{\Psi }(𝐫,\tau )=\psi (r,z,\tau )\mathrm{exp}(iL\theta )$$
(2)
of the wave function, where $`\mathrm{exp}(iL\theta )`$ are the circular harmonics in two dimensions.

Now transforming to dimensionless variables defined by $`x=\sqrt{2}r/l`$, $`y=\sqrt{2}z/l`$, $`t=\tau \omega ,`$ $`l\sqrt{\mathrm{}/(m\omega )}`$, and

$$\frac{\phi (x,y;t)}{x}\sqrt{\frac{l^3}{\sqrt{8}}}\psi (r,z;\tau ),$$
(3)
we get from Eqs. (1) and (2)

$`[`$ $`i{\displaystyle \frac{}{t}}{\displaystyle \frac{^2}{x^2}}+{\displaystyle \frac{1}{x}}{\displaystyle \frac{}{x}}{\displaystyle \frac{^2}{y^2}}+{\displaystyle \frac{L^21}{x^2}}+{\displaystyle \frac{1}{4}}\left(x^2+\lambda ^2y^2\right)`$ (4)
$`+`$ $`\kappa n|{\displaystyle \frac{\phi (x,y;t)}{x}}|^2i\xi n^2|{\displaystyle \frac{\phi (x,y;t)}{x}}|^4]\phi (x,y;t)=0,`$
where $`n=Na/l`$, $`\kappa =8\sqrt{2}\pi `$, and $`\xi =4K_3/(a^2l^4\omega ).`$ From theoretical ver  and experimental k3  studies it has been found that for negative $`a,`$ $`K_3`$ increases rapidly as $`|a|^n`$, where the theoretical study favors $`n=2`$ and we represent this variation via this quadratic dependence. This makes the parameter $`\xi `$ above a constant th4  for an experimental set up with fixed $`l`$ and $`\omega `$ and in the present study we employ a constant $`\xi `$.

The normalization condition of the wave function for $`K_3=0`$ is

$$𝒩_{\text{norm}}2\pi _0^{\mathrm{}}\frac{dx}{x}_{\mathrm{}}^{\mathrm{}}𝑑y|\phi (x,y;t)|^2=1.$$
(5)
However, in the presence of loss $`K_3>0`$, $`𝒩_{\text{norm}}<1.`$ The number of remaining atoms $`N`$ in the condensate is given by $`N=N_0𝒩_{\text{norm}}`$, where $`N_0`$ is the initial number.

The root mean square (rms) sizes $`x_{\text{rms}}`$ and $`y_{\text{rms}}`$ are defined by

$`x_{\text{rms}}^2={\displaystyle \frac{2\pi }{𝒩_{\text{norm}}}}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑y|\phi (x,y;t)|^2x,`$ (6)
$`y_{\text{rms}}^2={\displaystyle \frac{2\pi }{𝒩_{\text{norm}}}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dx}{x}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑y|\phi (x,y;t)|^2y^2.`$ (7)
### II.2 Calculational Detail

We solve the GP equation (4) numerically using a split-step time-iteration method using the Crank-Nicholson discretization scheme sk1 ; th4 ; sk2 ; sk3 ; koo . We discretize the GP equation with time step 0.001 and space step 0.1 spanning $`x`$ from 0 to 15 and $`y`$ from $`35`$ to 35.

It is now appropriate to calculate the parameters of the present dimensionless GP equation (4) corresponding to the experiment at JILA for $`L=0`$ ex4 . As in that experiment we take the radial and axial trap frequencies to be $`\nu _{\text{radial}}=17.5`$ Hz and $`\nu _{\text{axial}}=6.8`$ Hz, respectively, leading to $`\lambda =0.389`$. The harmonic oscillator length $`l`$ of 85Rb atoms for $`\omega =2\pi \times 17.5`$ Hz and $`m79176`$ MeV is $`l=\sqrt{\mathrm{}/(m\omega )}=26070`$ Å. One unit of time $`t`$ of Eq. (4) is $`1/\omega `$ or 0.009095 s.

We consider a stable 85Rb condensate of $`N_0=16000`$ atoms with scattering length $`a_{\text{initial}}=7a_0`$, $`a_0=0.5292`$ Å. This wave function is obtained by time iteration of Eq. (4) employing the following normalized initial solution with a single central vortex sk1

$$\phi (x,y)=\left[\frac{\lambda }{2^{2L+3}\pi ^3(|L|!)^2}\right]^{1/4}x^{1+|L|}e^{(x^2+\lambda y^2)/4}$$
(8)
for $`n=0`$. In the course of above time iteration the nonlinear parameter $`n`$ was increased by steps of 0.0001 until its final value is attained. Then during an interval of time 0.1 ms the scattering length was ramped to $`a=a_{\text{collapse}}`$. The absorptive term $`\xi `$ was set equal to zero throughout above time iteration.

The final condensate is strongly attractive and unstable and undergoes a sequence of collapse and explosion. In our numerical simulation with $`L0`$ we consider a set of different values of $`a_{\text{collapse}}`$ $`(=263a_0,100a_0,30a_0,20a_0,`$ etc.) as well as $`N_0=6000`$.

For the simulation of collapse and explosion a nonzero value of $`\xi `$ $`(=2)`$ is chosen for different $`a_{\text{initial}}`$, $`a_{\text{collapse}}`$, and $`N_0`$ as in Ref. th4  and the time-evolution of the GP equation is continued. This value of $`\xi `$ reproduced the essentials of the experiment at JILA ex4  reasonably well for $`L=0`$ as well as produced th4  a $`K_3`$ in reasonable agreement with a previous experiment ($`K_3=4.24\times 10^{25}`$ cm6/s) k3  and theoretical calculation ($`K_3=6.7\times 10^{25}`$ cm6/s) esry  for $`a=370a_0`$. In particular we use $`K_3=9\times 10^{25}`$ cm6/s for $`a=370a_0`$. For smaller values of $`|a|`$, the $`K_3`$ values are scaled down using the relation $`K_3a^2`$.

## III Numerical Result

The numerical simulation using Eq. (4) with a nonzero $`\xi `$ as described above immediately yields the remaining number of atoms in the condensate after the jump in scattering length. The remaining number of atoms vs. time is plotted in Fig. 1 (a) for $`a_{\text{initial}}=7a_0`$, $`a_{\text{collapse}}=30a_0`$ and $`263a_0`$, $`N_0=16000`$, and $`L=0,1,2,4,`$ and 6. In Fig. 1 (b) the same results for $`N_0=6000`$ are plotted. In this figure we also plot some results of experiment at JILA for $`L=0`$ ex4 ; th4 . These experimental results are in agreement with the simulation for $`L=0`$. In Fig. 2 (a) we plot the particle loss curves for $`N_0=6000`$, $`a_{\text{initial}}=7a_0`$ and different $`a_{\text{collapse}}`$ for $`L=0`$. In Fig. 2 (b) we plot the same for $`L=1`$.

From Figs. 1 and 2 we find that for a fixed $`N_0`$, for a sufficiently small $`L`$ or a sufficiently large $`|a_{\text{collapse}}|`$, there could be collapse and explosion during a relatively short interval of time (called decay time) with the loss of a large fraction of the atoms. However, there is no collapse for a large enough $`L`$ or a small enough $`|a_{\text{collapse}}|`$. For example, for $`N_0=6000`$ in Fig. 1 (b) there is no collapse for $`L>2`$ for $`|a_{\text{collapse}}|=30a_0`$ and in Fig. 2 (a) there is no collapse for $`|a_{\text{collapse}}|<5a_0`$ for $`L=0`$.

In the experiment at JILA ex4  for $`L=0`$ it was observed that the strongly attractive condensate after preparation remains stable with a constant number of atoms for an interval of time $`t_{\text{collapse}}`$, called collapse time. This behavior is physically expected for medium to small values of $`|a_{\text{collapse}}|`$ $`(<50a_0)`$. Immediately after the jump in the scattering length to a negative value, the attractive condensate shrinks in size during $`t_{\text{collapse}}`$, until the central density increases to a maximum. Then the absorptive three-body term takes full control to initiate the explosion that last for few milliseconds. Consequently, the number of atoms remains constant for time $`t<t_{\text{collapse}}`$. The present results in Figs. 1(a) and (b) also show this behavior for $`|a_{\text{collapse}}|=30a_0`$. However, for larger $`|a_{\text{collapse}}|`$ $`(=263a_0)`$, the atomic attraction is very strong and the central density increases to a maximum quickly to start the explosion and $`t_{\text{collapse}}`$ is close to zero. In Figs. 2 (a) and (b) we see the dependence of particle loss and $`t_{\text{collapse}}`$ on $`|a_{\text{collapse}}|`$ for $`N_0=6000`$, $`a_{\text{initial}}=7a_0`$, and $`L=0`$ and $`L=1`$, respectively. From Figs. 1 and 2 we find that $`t_{\text{collapse}}`$ increases with $`L`$ for a fixed $`a_{\text{collapse}}`$ and with $`a_{\text{collapse}}`$ for a fixed $`L`$.

From Figs. 1 and 2 we find that after the collapse, the number of particles drop sharply during a small interval of time called decay time (few milliseconds), which means that the condensate emits a large number particles in an explosive fashion. This emission of particles is termed explosion.

After a sequence of collapse and explosion, for $`L=0`$ Donley et al. ex4  observed a “remnant” condensate of $`N_{\text{remnant}}`$ atoms at large times containing a fraction of the initial $`N_0`$ atoms. Figures 1 and 2 show such a behavior for different values of $`L`$ and $`a_{\text{collapse}}`$. In all cases the decay time during which the explosion takes place is small and of the order of few milliseconds. The decay time for vortex states ($`L0`$) is smaller than for nonrotating condensates $`(L=0)`$.

We studied the time evolution of the condensate for larger times. In Fig. 3 we plot the loss curves for $`N_0=16000`$, $`a_{\text{initial}}=7a_0`$, $`a_{\text{collapse}}=30a_0`$, and $`L=0,1,2,4,6`$ at larger times. The BEC continues to lose atoms if left for a long time but at a rate much slower than during the first explosion that we call primary. However, in the experiment at JILA Donley et al observed that a remnant condensate containing a fraction of the atoms survived with nearly constant number for more than one second. One possible reason for this discrepancy could be the following. In the actual experiment at JILA a major portion of the emitted atoms, called the burst atoms, remain trapped and oscillate around the central remnant. The presence of the burst atoms make the measurement of the number of atoms in the remnant a difficult task don . Some of these burst atoms may also rejoin the remnant to compensate for the three-body loss at large times. Such an effect is not included in the present model, which, hence, presents a larger loss for the remnant compared to experiment.

We also observe an interesting phenomenon in Fig. 3, e. g., the occurrence of smaller secondary and tertiary explosions after the primary one observed for small times. For $`L=2`$ after the primary collapse and explosion at $`t<10`$ ms with the loss of about 7000 atoms, there is another collapse and explosion with loss of about 2000 atoms at $`t150`$ ms. The primary and secondary explosions are separated by a large interval of time. We also see much weaker explosion(s) in the course of time evolution in Fig. 3, where the particle number varies in small steps. These explosions could be termed tertiary with the loss of few hundred atoms. It might be interesting to see if such secondary and tertiary explosion(s) could be observed experimentally.

Donley et al. ex4  provided a quantitative measurement of the variation of collapse time $`t_{\text{collapse}}`$ with the final scattering length $`a_{\text{collapse}}`$ for a given $`a_{\text{initial}}=0`$, $`N_0=6000`$, and $`L=0`$. We calculated this variation using our model for $`L=0,1`$ and 2. The $`t_{\text{collapse}}`$ vs. $`|a_{\text{collapse}}|/a_0`$ plots for $`L=0,1,2`$ are exhibited in Fig. 4 and compared with experimental data for $`L=0`$ ex4  as well as with another calculation using the mean-field GP equation in an axially symmetric trap for $`L=0`$ th3 . We see $`t_{\text{collapse}}`$ decreases with $`|a_{\text{collapse}}|/a_0`$ starting from an infinite value at $`|a_{\text{collapse}}|=a_{\text{cr}}`$, where $`a_{\text{cr}}`$ is the minimum value of $`|a_{\text{collapse}}|`$ that leads to collapse and explosion. The critical value $`a_{\text{cr}}`$ increases with $`L`$ and so does $`t_{\text{collapse}}`$ for a fixed $`|a_{\text{collapse}}|`$. For a given $`N_0`$, a critical value of $`nn_{\text{cr}}`$ can be defined via $`n_{\text{cr}}N_0a_{\text{cr}}/l`$. A necessary condition for collapse is $`N_0|a_{\text{collapse}}|/l>n_{\text{cr}}`$ sk1 . The value of $`n_{\text{cr}}`$ for a specific case is calculated as in Ref. sk1 . Consequently, $`a_{\text{cr}}/a_0`$ can be obtained. The value of $`a_{\text{cr}}/a_0`$ so evaluated for a specific $`L`$ is shown by an arrow near the curve for that particular $`L`$ in Fig. 4. The $`t_{\text{collapse}}`$ vs. $`|a_{\text{collapse}}|/a_0`$ curves should tend to infinity at the respective arrows and they do so in Fig. 4. There should not be any collapse for $`|a_{\text{collapse}}|<a_{\text{cr}}`$.

Donley et al. ex4  measured the number of remnant atoms for $`L=0`$, $`a_{\text{initial}}=7a_0`$ and different $`N_0`$ and $`a_{\text{collapse}}`$ th4 . We plot the same in Fig. 5 for $`L=0,1`$ and 2. The remnant number plotted in this figure is the number after the primary explosion(s) and not during or after possible secondary and tertiary explosions at larger times. In Figs. 1 and 2 the remnant number is obtained around $`t2030`$ ms and not at few hundred milliseconds. Our results in Fig. 5 for $`L=0`$ agree well th4  with the measurements of Donley et al. ex4 . In general the remnant number decreases with increasing $`L`$. However, there are some cases where the opposite trend has been observed in Figs. 1 and 2 and as well as in Fig. 5. For the smallest values of $`N_0`$ in Fig. 5, the condensate remains stable for $`L>0`$, and there is no collapse and explosion and hence no remnant numbers for $`L=1`$ and 2. For certain $`N_0`$ the results for all three $`L`$’s are not plotted as they coincide with other remnant numbers. The remnant number in some cases could be much larger than $`N_{\text{cr}}`$ for times on the order of tens of milliseconds.

Donley et al. ex4  observed that for $`L=0`$ the remnant condensate always oscillated in a highly excited collective state with approximate frequencies $`2\nu _{\text{axial}}`$ and $`2\nu _{\text{radial}}`$ being predominantly excited. This behavior emerges from the present simulation for all values of $`L`$. To illustrate this we plot in Fig. 6 sizes $`x_{\text{rms}}`$ and $`y_{\text{rms}}`$ vs. time for the condensate after the jump in the scattering length to $`30a_0`$ from $`7a_0`$ for $`N_0=16000`$ and $`L=2`$. We find a periodic oscillation in $`x_{\text{rms}}`$ and $`y_{\text{rms}}`$ with frequencies 13.6 Hz ($`2\nu _{\text{axial}}`$) and 35 Hz ($`2\nu _{\text{radial}}`$), respectively, as observed in experiment.

## IV Conclusion

In conclusion, we have employed a numerical simulation based on the solution sk1  of the mean-field Gross-Pitaevskii equation 8  with cylindrical symmetry to study the dynamics of collapse and explosion ex4 ; th4  of small attractive vortex states with $`L>0`$. The explosion is initiated by a sudden jump in the scattering length from a positive to negative value exploiting a Feshbach resonance fbna ; ex3 ; fbcs ; fbth . In the GP equation we include a quintic three-body nonlinear recombination loss term th1 ; th1a ; th2 ; th3  that accounts for the decay of the strongly attractive condensate. The results of the present simulation are to be considered as an extension of the experiment at JILA for $`L=0`$ ex4  to small vortex states.

We find the following features of this dynamics from the present numerical simulation: (1) The condensate undergoes collapse and explosion during a small interval of time of few milliseconds and finally stabilizes to a remnant condensate containing a fraction of initial number of atoms. The number in the remnant condensate for times on the order of tens of milliseconds can be much larger than the critical number for collapse $`N_{\text{cr}}`$ for the same atomic interaction. (2) In some cases after the primary explosion small secondary and tertiary explosions are observed. (3) The explosion takes place during a decay time of few milliseconds. This decay time for rotating Bosenova with vortex ($`L>0`$) is smaller than the same for nonrotating condensate with $`L=0`$. (4) The remnant condensate executes radial and axial oscillations in a highly excited collective state for a long time with frequencies $`2\nu _{\text{radial}}`$ and $`2\nu _{\text{axial}}.`$ (5) After the sudden change in the scattering length to a large negative value, the condensate needs an interval of time $`t_{\text{collapse}}`$ before it experiences loss via explosion. The interval $`t_{\text{collapse}}`$ increases with $`L`$ and $`a_{\text{collapse}}`$.

The simulation of the particle loss in strongly attractive rotating Bosenova, with a single axial vortex of small angular momentum per particle, may stimulate further theoretical and experimental studies. We have considered small vortex states as they can be well described by the mean-field GP equation. This will provide a test for the usefulness of this equation in handling particle loss. Otherwise, a similar study with a large Bose-Einstein condensed vortex lattice exp1 ; expm  is more challenging from both experimental and theoretical points.

###### Acknowledgements.

The work is supported in part by the CNPq and FAPESP of Brazil.


text:
# REFERENCES

SLAC-PUB-8813

April 19, 2002

Precision Measurement of the Proton and Deuteron Spin Structure Functions g2 and Asymmetries A2***Work supported by Department of Energy contract DE–AC03–76SF00515.

The E155 Collaboration P. L. Anthony,12 R. G. Arnold,1,†† T. Averett,15 H. R. Band,16 N. Benmouna,1 W. Boeglin,5 H. Borel,4 P. E. Bosted,1,†† S. L. B$`\ddot{\mathrm{u}}`$ltmann,14,§§ G. R. Court,6 D. Crabb,14 D. Day,14 P. Decowski,11 P. DePietro,1 H. Egiyan,15 R. Erbacher,12,∞ R. Erickson,12 R. Fatemi,14 E. Frlez,14 K. A. Griffioen,15 C. Harris,14 E. W. Hughes,2 C. Hyde-Wright,9 G. Igo,3 J. Johnson,12 P. King,15 K. Kramer,15 S. E. Kuhn,9 D. Lawrence,8 Y. Liang,1 R. Lindgren,14 R. M. Lombard-Nelsen,4 P. McKee,14 D.E. McNulty,14 W. Meyer,14,♭ G. S. Mitchell,16,§ J. Mitchell,13 M. Olson,10 S. Penttila,7‘ G. A. Peterson,8 R. Pitthan,12 D. Pocanic,14 R. Prepost,16 C. Prescott,12 B. A. Raue,5 D. Reyna,1,♡ P. Ryan,15 L. S. Rochester,12 S. Rock,1,†† O. Rondon-Aramayo,14 F. Sabatie,9,♠ T. Smith,7 L. Sorrell,1 S. St.Lorant,12 Z. Szalata,1,† Y. Terrien,4 A. Tobias,14 T. Toole,1,‡ S. Trentalange,3 F. R. Wesselmann,9 T. R. Wright,16 M. Zeier,14 H. Zhu,14 B. Zihlmann14

1American University, Washington, D.C. 20016 2California Institute of Technology, Pasadena, California 91125 3University of California, Los Angeles, California 90095 4DAPNIA-Service de Physique‘ Nucleaire, CEA-Saclay, F-91191 Gif/Yvette Cedex, France 5Florida International University, Miami, Florida 33199. 6University of Liverpool, Liverpool L69 3BX, United Kingdom 7Los Alamos National Laboratory, Los Alamos, New Mexico 87545 8University of Massachusetts, Amherst, Massachusetts 01003 9Old Dominion University, Norfolk, Virginia 23529 10St. Norbert College, De Pere, WI 54115 11Smith College, Northampton, Massachusetts 01063 12Stanford Linear Accelerator Center, Stanford, California 94309 13Thomas Jefferson National Accelerator Facility, Newport News, Virginia 23606 14University of Virginia, Charlottesville, Virginia 22901 15The College of William and Mary, Williamsburg, Virginia 23187 16University of Wisconsin, Madison, Wisconsin 53706 To be submitted to Physical Review Letters

## Abstract

We have measured the spin structure functions $`g_2^p`$ and $`g_2^d`$ and the virtual photon asymmetries $`A_2^p`$ and $`A_2^d`$ over the kinematic range $`0.02x0.8`$ and $`0.7Q^220`$ GeV2 by scattering 29.1 and 32.3 GeV longitudinally polarized electrons from transversely polarized NH3 and $`{}_{}{}^{6}\mathrm{LiD}`$ targets. Our measured $`g_2`$ approximately follows the twist-2 Wandzura-Wilczek calculation. The twist-3 reduced matrix elements $`d_2^p`$ and $`d_2^n`$ are less than two standard deviations from zero. The data are inconsistent with the Burkhardt-Cottingham sum rule if there is no pathological behavior as $`x0`$. The Efremov-Leader-Teryaev integral is consistent with zero within our measured kinematic range. The absolute value of A2 is significantly smaller than the $`A_2<\sqrt{R(1+A_1)/2}`$ limit.

The deep inelastic spin structure functions of the nucleons, $`g_1(x,Q^2)`$ and $`g_2(x,Q^2)`$, depend on the spin distribution of the partons and their correlations. The function $`g_1`$ can be primarily understood in terms of the quark parton model (QPM) and perturbative QCD with higher twist terms at low $`Q^2.`$ The function $`g_2`$ is of particular interest since it has contributions from quark-gluon correlations and other higher twist terms at leading order in $`Q^2`$ which cannot be described perturbatively. By interpreting $`g_2`$ using the operator product expansion (OPE) , it is possible to study contributions to the nucleon spin structure beyond the simple QPM.

The structure function $`g_2`$ can be written:

$$g_2(x,Q^2)=g_2^{WW}(x,Q^2)+\overline{g_2}(x,Q^2)$$
(1)
in which

$`g_2^{WW}(x,Q^2)`$ $`=`$ $`g_1(x,Q^2)+{\displaystyle _x^1}{\displaystyle \frac{g_1(y,Q^2)}{y}}𝑑y,`$ (2)
$`\overline{g_2}(x,Q^2)`$ $`=`$ $`{\displaystyle _x^1}{\displaystyle \frac{}{y}}\left({\displaystyle \frac{m}{M}}h_T(y,Q^2)+\xi (y,Q^2)\right){\displaystyle \frac{dy}{y}},`$ (3)
$`x`$ is the Bjorken scaling variable and $`Q^2`$ is the absolute value of the virtual photon four-momentum squared. The twist-2 term $`g_2^{WW}`$ was derived by Wandzura and Wilczek and depends only on $`g_1`$. The function $`h_T(x,Q^2)`$ is an additional twist-2 contribution that depends on the transverse polarization density in the nucleon. The $`h_T`$ contribution to $`\overline{g_2}`$ is suppressed by the ratio of the quark to nucleon masses $`m/M`$ and its effect is thus small for up and down quarks and is assumed negligible in the theoretical models referenced in this paper. The twist-3 part ($`\xi `$) comes from quark-gluon correlations and is the main focus of our study.

Low-precision measurements of $`g_2`$ and $`A_2`$ exist for the proton and deuteron , as well as for the neutron. In this Letter, we report new, precise measurements of $`g_2`$ and $`A_2`$ for the proton and deuteron.

Electron beams with energies of 29.1 and 32.3 GeV and longitudinal polarizations of $`P_b=(83.2\pm 3.0)`$% struck approximately transversely polarized NH3 (average polarization $`<P_t>=0.70)`$ or 6LiD $`(<P_t>=0.22)`$ targets. The beam helicity was randomly chosen pulse by pulse. Scattered electrons were detected in three independent spectrometers centered at 2.75, 5.5, and 10.5. The two small-angle spectrometers were the same as in SLAC E155, while the large-angle spectrometer had additional hodoscopes and a more efficient pre-radiator shower counter. Further information on the experimental apparatus can be found in references. The approximately equal amounts of data taken with the two beam energies and opposites signs of target polarization gave consistent results.

The measured asymmetry, $`\stackrel{~}{A_{}}`$, differs from $`A_{}`$ because the target polarizations were not exactly perpendicular to the beam line. For each kinematic bin $`\stackrel{~}{A_{}}`$ was formed using

$`\stackrel{~}{A_{}}={\displaystyle \frac{1}{f_{RC}}}\left[{\displaystyle \frac{C_1}{fP_t}}\left(\left({\displaystyle \frac{N_LN_R}{N_L+N_R}}\right){\displaystyle \frac{1}{P_b}}A_{EW}\right)+C_2\right]+A_{RC}`$ (4)
where $`N_L`$ and $`N_R`$ are the measured counting rates from the two beam helicities, including small corrections for pion and charge symmetric backgrounds, dead-time and tracking efficiency, and $`A_{EW}`$ is the electroweak asymmetry ($`8\times 10^5Q^2`$). The target dilution factor, $`f`$, is the fraction of free polarizable protons ($`0.13`$) or deuterons ($`0.18`$) for a given spectrometer acceptance. For the proton target, the nuclear correction $`C_10.98`$ is due to the polarization of the 15N and $`C_2`$=0. The deuteron data were extracted from the 6LiD results by applying a slightly $`x`$-dependent nuclear correction $`C_10.52`$ to account for the lithium and deuterium nuclear wave functions with 6Li$`\alpha +d`$. An additional correction $`C_2(x)\sigma _p/\sigma _d\stackrel{~}{A_{}^p}`$, ($`C_20.042`$) accounts for the $`4\%`$ polarized 7Li in the target. The quantities $`f_{RC}`$ and $`A_{RC}`$ are radiative corrections determined using a method similar to E143. The quantity $`1f_{RC}`$ was calculated as the proportion of events in a bin coming from elastic and quasi-elastic tails, and $`A_{RC}`$ included polarization-dependent elastic and quasi-elastic as well as inelastic and vertex corrections. The radiative dilution factor $`f_{RC}`$ has the effect of increasing the statistical errors at low $`x`$. Uncertainties in the radiative corrections were estimated by varying the input models over a range consistent with the measured data.

The detailed results for $`\stackrel{~}{A_{}}`$ are shown in Table I. Because $`\stackrel{~}{A_{}}`$ is close to zero, the relative statistical errors are always greater than 25%. The uncertainties due to target and beam polarization and dilution factor combined are 5.1% (proton) and 6.2% (deuteron). They are multiplicative and small compared to the statistical errors. The total systematic error is also shown in Table I.

We determined $`g_2(x,Q^2)`$ and $`A_2(x,Q^2)`$ from $`\stackrel{~}{A_{}}`$ (dominant contribution) and the previously measured $`g_1`$ (small contribution) using:

$`g_2`$ $`=`$ $`{\displaystyle \frac{yF_1}{2E^{}(\mathrm{cos}\mathrm{\Theta }\mathrm{cos}\alpha )}}\left[\stackrel{~}{A_{}}\nu {\displaystyle \frac{(1+ϵR)}{1ϵ}}{\displaystyle \frac{g_1}{F_1}}[E\mathrm{cos}\alpha +E^{}\mathrm{cos}\mathrm{\Theta }]\right]`$ (5)
$`A_2`$ $`=`$ $`\gamma (g_1+g_2)/F_1`$ (6)
where $`\mathrm{cos}\mathrm{\Theta }=\mathrm{sin}\alpha \mathrm{sin}\theta \mathrm{cos}\mathrm{\Phi }+\mathrm{cos}\alpha \mathrm{cos}\theta `$, $`\theta `$ is the spectrometer angle, $`\mathrm{\Phi }`$ is the angle between the spin plane and the scattering plane, $`\alpha =92.4^{}`$ is the angle of the target polarization with respect to the beam direction, $`y=\nu /E`$, $`\nu =EE^{}`$, $`E`$ and $`E^{}`$ are the incident and scattered electron energies, $`ϵ^1=1+2\left[1+1/\gamma ^2\right]\mathrm{tan}^2(\theta /2)`$, $`\gamma =\sqrt{Q^2/\nu ^2}`$ and $`F_1=F_2(1+4M^2x^2/Q^2)/[2x(1+R)].`$ We used a new $`Q^2`$-dependent parameterization of $`g_1`$  world data, the NMC fit to $`F_2(x,Q^2)`$  and the SLAC fit to $`R(x,Q^2)=\sigma _L/\sigma _T`$ . The structure functions for p, d, and n are related by $`g_2^d=(g_2^p+g_2^n)(11.5\omega _D)/2`$, where $`\omega _D=0.05`$, the fraction of D-wave in the deuteron wave function.

Results for $`A_2`$ and $`xg_2`$ for the three spectrometers and two energies are given in Table II with statistical errors. The systematic error on $`xg_2`$ is much smaller than the statistical error and is given approximately by $`a+bx`$ where $`a_p(a_d)`$=$`0.0016(0.0018)`$ and $`b_p(b_d)`$=$`0.0014(0.0021).`$ It includes the systematic errors on $`\stackrel{~}{A_{}}`$ as well as a 5% normalization uncertainty on $`g_1`$. The data cover the kinematic range $`0.02x0.8`$ and $`0.7Q^220`$ GeV2 with an average $`Q^2`$ of $`5`$ GeV2. Figure 1 shows the values of $`xg_2`$ as a function of $`Q^2`$ for several values of $`x`$ along with results from E143 and E155. The data approximately follow the $`Q^2`$ dependence of $`g_2^{WW}`$ (solid curve), although for the proton, the data points are lower than $`g_2^{WW}`$ at low and intermediate $`x`$ and higher at high $`x`$. The predictions of Stratmann are closer to the data.

To get average values at the average $`Q^2`$ for each $`x`$ bin we used the $`Q^2`$ dependence of $`g_2^{WW}`$: $`g_2(Q_{avg}^2)=g_2(Q_{exp}^2)g_2^{WW}(Q_{exp}^2)+g_2^{WW}(Q_{avg}^2).`$ These averaged results for $`A_2`$ and $`xg_2`$ are listed at the bottom of Table II. Figure 2 shows the averaged $`xg_2`$ of this experiment along with $`xg_2^{WW}`$ calculated using our parameterization of $`g_1.`$ The combined new data for p disagree with $`g_2^{WW}`$ with a $`\chi ^2/`$dof of 3.1 for 10 degrees of freedom. For d the new data agree with $`g_2^{WW}`$ with a $`\chi ^2/`$dof of 1.2 for 10 dof. The data for $`g_2^p`$ are also inconsistent with zero ($`\chi ^2/`$dof=15.5) while $`g_2^d`$ differs from zero only at $`x0.4.`$ Also shown in Fig. 2 is the Bag Model calculation of Stratmann which is in good agreement with the data, a chiral soliton model calculation which is too negative at $`x0.4`$ and the Bag Model calculation of Song which is in clear disagreement with the data.

The average values of $`A_2(x)`$, shown in Fig. 3, are consistent with zero at low $`x`$, increasing to about 0.1 at the highest $`x,`$ significantly different than zero. $`A_2^p`$ is many standard deviations lower than the Soffer limit of $`|A_2|<\sqrt{R(1+A_1)/2}0.4`$ for all values of $`x`$. The same is true for $`A_2^d`$, except at the highest $`x`$ value where the error is large.

The OPE allows us to write the hadronic matrix element in deep inelastic scattering in terms of a series of renormalized operators of increasing twist. The moments of $`g_1`$ and $`g_2`$ for even $`n2`$ at fixed $`Q^2`$ can be related to twist-3 reduced matrix element, $`d_n`$, and higher twist terms which are suppressed by powers of $`1/Q`$. Neglecting quark mass terms we find that:

$$d_n=2_0^1𝑑xx^n\left[\frac{n+1}{n}g_2(x,Q^2)+g_1(x,Q^2)\right]=2\frac{n+1}{n}_0^1𝑑xx^n\overline{g_2}(x,Q^2).$$
(7)
This relation is true to all orders of $`\alpha _s`$. The matrix element $`d_n`$ measures deviations of $`g_2`$ from the twist-2 $`g_2^{WW}`$ term. Note that some authors define $`d_n`$ with an additional factor of two. We calculated $`d_n`$ using $`d_2=3_0^1𝑑xx^2\overline{g_2}(x,Q^2)`$ (see eq.7) with the assumption that $`\overline{g_2}`$ is independent of $`Q^2`$ in the measured region. This is not unreasonable since $`d_n`$ depends only logarithmically on $`Q^2`$. The part of the integral for $`x`$ below the measured region was assumed to be zero because of the $`x^2`$ suppression. For $`x0.8`$ we used $`\overline{g_2}(1x)^m`$ where $`m`$=2 or 3, normalized to the data for $`x0.5`$. Because $`\overline{g_2}`$ is small at high $`x`$, the contribution was negligible for both cases. We obtained values of $`d_2^p=`$$`0.0025`$$`\pm `$0.0016 $`\pm `$0.0010 and $`d_2^d=`$$`0.0054`$$`\pm `$0.0023 $`\pm `$0.0005 at an average $`Q^2`$ of $`5`$ GeV2. We combined these results with those from SLAC experiments on the neutron (E142 and E154) and proton and deuteron (E143 and E155) to obtained average values $`d_2^p=`$$`0.0032`$$`\pm `$0.0017 and $`d_2^n=`$$`0.0079`$$`\pm `$0.0048 . These are consistent with zero (no twist-3) to within 2 standard deviations. The values of the 2nd moments alone are: $`_0^1𝑑xx^2g_2(x,Q^2)`$=$`0.0072`$$`\pm `$ 0.0005 $`\pm `$0.0003 (p) and $`0.0019`$$`\pm `$ 0.0007 $`\pm `$0.0001 (d).

Figure 4 shows the experimental values of $`d_2`$ for proton and neutron with their error, plotted along with theoretical models from left to right: Bag Models (Song, Stratmann, and Ji); sum rules (Stein, BBK, Ehrnsperger); chiral soliton models; and lattice QCD calculations ($`Q^2=5`$ GeV2, $`\beta =6.4`$). The lattice and chiral calculations are in good agreement with the proton data and two standard deviations below the neutron data. The sum rule calculations are significantly lower than the data. The Non Singlet $`3(d_2^pd_2^n)=`$ $`0.0141`$$`\pm `$0.0170 is consistent with an instanton vacuum calculation of $`0.001`$ .

The Burkhardt-Cottingham sum rule for $`g_2`$ at large $`Q^2`$, $`_0^1g_2(x)𝑑x=0,`$ was derived from virtual Compton scattering dispersion relations. It does not follow from the OPE since $`n=0`$. Its validity depends on the lack of singularities for $`g_2`$ at $`x=0,`$ and a dramatic rise of $`g_2`$ at low $`x`$ could invalidate the sum rule. We evaluated the Burkhardt-Cottingham integral in the measured region of $`0.02x0.8`$ at $`Q^2=5\mathrm{GeV}^2.`$ The results for the proton and deuteron are $`0.044`$$`\pm `$0.008 $`\pm `$0.003 and $`0.008`$$`\pm `$0.012 $`\pm `$0.002 respectively. Averaging with the E143 and E155 results which cover a slightly more restrictive $`x`$ range gives $`0.042`$$`\pm `$0.008 and $`0.006`$$`\pm `$0.011. This does not represent a conclusive test of the sum rule because the behavior of $`g_2`$ as $`x0`$ is not known. However, if we assume that $`g_2=g_2^{WW}`$ for $`x<0.02`$, and use the relation $`_0^xg_2^{WW}(y)𝑑y=x\left[g_2^{WW}(x)+g_1(x)\right]`$, there is an additional contribution of 0.020 (0.004) for the proton (deuteron).

The Efremov-Leader-Teryaev (ELT) sum rule involves the valence quark contributions to $`g_1`$ and $`g_2`$: $`_0^1x[g_1^V(x)+2g_2^V(x)]𝑑x=0.`$ Assuming that the sea quarks are the same in protons and neutrons, the sum rule takes a form $`_0^1x[g_1^p(x)+2g_2^p(x)g_1^n(x)2g_2^n(x)]𝑑x=0.`$ We evaluated this ELT integral in the measured region using our $`g_2`$ data and the fit to $`g_1`$. The result at $`Q^2=5\mathrm{GeV}^2`$ is $`0.013`$$`\pm `$ 0.008$`\pm `$0.002 , which is consistent with the expected value of zero. Including the data of E143  and E155 leads to $`0.011`$$`\pm `$ 0.008. The extrapolation to $`x`$=0 is not known, but is suppressed by a factor of $`x.`$

In summary, our results for $`g_2`$ follow approximately the twist-2 $`g_2^{WW}`$ shape, but deviate significantly at some values of $`x`$. The values obtained for the twist-3 matrix element $`d_2`$ from this measurement and the SLAC average are less than two standard deviations from zero. The data over the measured range are inconsistent with the Burkhardt-Cottingham sum rule if there is no pathological behavior as $`x0`$. The ELT integral is consistent with zero within our measured kinematic range.

We wish to thank the personnel of the SLAC accelerator and Experimental Facilities Departments for their efforts which resulted in the successful completion of the E155X experiment. This work was supported by the Department of Energy; the National Science Foundation; and the Centre National de la Recherche Scientifique and the Commissariat a l’Energie Atomique (French groups).


text:
# Northwestern University: N.U.H.E.P. Report No. 1000 April 4, 2002 Forward Compton Scattering, using Real Analytic Amplitudes

## 1 Introduction

The $`\rho `$-value is defined as the ratio of the real to the imaginary portion of the forward scattering amplitude. Block and Kadailov have shown that $`\rho _{\mathrm{nn}}=\rho _{\gamma p}`$ if one uses eikonals for the even portion of nucleon-nucleon scattering and $`\gamma p`$ scattering that have equal opacities, i.e., eikonals that have the same value at impact parameter $`b=0`$. This is the equivalent of the more physical statement that

$$\left(\frac{\sigma _{\mathrm{elastic}}(s)}{\sigma _{\mathrm{tot}}(s)}\right)_{\gamma p}=\left(\frac{\sigma _{\mathrm{elastic}}(s)}{\sigma _{\mathrm{tot}}(s)}\right)_{nn},\mathrm{for}\mathrm{all}s.$$
(1)
Block and Kaidalov have proved three factorization theorems:

1. $$\frac{\sigma _{\mathrm{nn}}(s)}{\sigma _{\gamma \mathrm{p}}(s)}=\frac{\sigma _{\gamma \mathrm{p}}(s)}{\sigma _{\gamma \gamma }(s)},$$

where the $`\sigma `$’s are the total cross sections for nucleon-nucleon, $`\gamma `$p and $`\gamma \gamma `$ scattering,
2. $$\frac{B_{\mathrm{nn}}(s)}{B_{\gamma \mathrm{p}}(s)}=\frac{B_{\gamma \mathrm{p}}(s)}{B_{\gamma \gamma }(s)},$$

where the $`B`$’s are the nuclear slope parameters for elastic scattering,
3. $$\rho _{\mathrm{nn}}(s)=\rho _{\gamma \mathrm{p}}(s)=\rho _{\gamma \gamma }(s),$$

where the $`\rho `$’s are the ratio of the real to imaginary portions of the forward scattering amplitudes,

with the first two factorization theorems having their own proportionality constant. These theorems are exact, for all  $`s`$ (where $`\sqrt{s}`$ is the c.m.s. energy), and survive exponentiation of the eikonal. The last theorem is valid independently of the model which takes one from $`nn`$ to $`\gamma p`$ to $`\gamma \gamma `$ reactions, as long as the respective eikonals have equal opacities, i.e., eq. (1) holds. We wish to demonstrate experimentally here the validity of the theorem that states that $`\rho _{\mathrm{nn}}(s)=\rho _{\gamma \mathrm{p}}(s)`$. However, no data are available in the hadronic sector for $`\rho _{\gamma \mathrm{p}}`$. The purpose of this note is to analyze forward Compton scattering at high energies in order to extract $`\rho _{\gamma p}`$ and then to compare it to $`\rho _{nn}`$.

Damashek and Gilman in 1970 have calculated $`\rho _{\gamma p}`$ using a singly-subtracted dispersion relation, up to a gamma ray laboratory energy $`\nu =20`$ GeV, i.e., to a c.m.s. energy of $`\sqrt{s}=6.2`$ GeV. Here we extend the $`\rho _{\gamma p}`$ evaluation up to $`\sqrt{s}=200`$ GeV and then compare the results to $`\rho _{nn}`$, the ratio of the even portions of the $`\overline{p}p`$ and $`pp`$ forward scattering amplitudes (for references, see ). In this paper we will calculate $`\rho _{\gamma p}`$ using real analytic amplitudes (see Section G. p.583 of ref. ). It is shown in ref.  that the numerical complexities of dispersion relations in analyzing $`\overline{p}p`$ and $`pp`$ scattering can be circumvented by direct use of analytic functions to fit forward $`\overline{p}p`$ and $`pp`$ scattering amplitudes, a technique first proposed by Bourrely and Fischer. We will introduce a variant of the Block and Cahn analysis appropriate for Compton scattering, $`\gamma p\gamma p`$.

## 2 Preliminaries

This work largely follows the procedures and conventions used by Block and Cahn. We use units where $`\mathrm{}=c=1`$. The variable $`s`$ is the square of the c.m. system energy, whereas $`\nu `$ is the laboratory system momentum. In terms of the even laboratory scattering amplitude $`f_+`$, where $`f_+(\nu )=f_+(\nu )`$, the total unpolarized Compton cross section $`\sigma _{\mathrm{tot}}`$ is given by

$`\sigma _{\mathrm{tot}}`$ $`=`$ $`{\displaystyle \frac{4\pi }{\nu }}\mathrm{Im}f_+(\theta =0),`$ (2)
where $`\theta `$ is the laboratory scattering angle. We will assume that our amplitudes are real analytic functions with a simple cut structure. For $`\gamma p`$ scattering, the assumed cut structure is a left-hand cut that begins at the gamma ray energy $`\nu _0`$ and a right-hand cut that begins at the gamma ray energy $`\nu _0`$, with a real amplitude on the real axis between $`\nu _0`$ and $`\nu _0`$. The threshold gamma ray energy for pion production is $`\nu _0=m_\pi +\frac{m_\pi ^2}{m}0.16`$ GeV, where $`m_\pi `$ and $`m`$ are the pion and proton masses, respectively. We will use an even amplitude for $`\gamma p`$ reactions in the high energy region $`\nu >>\nu _0`$, far above any cuts, (see ref., p. 587, eq. (5.5a), with $`a=0`$), where the even amplitude simplifies considerably and is given by

$$f_+=i\frac{\nu }{4\pi }\left\{A+\beta [\mathrm{ln}(s/s_0)i\pi /2]^2+cs^{\mu 1}e^{i\pi (1\mu )/2}\right\}+C_{\mathrm{subtraction}},$$
(3)
where $`A`$, $`\beta `$, $`c`$, $`s_0`$ and $`\mu `$ are real constants. The additional real constant $`C_{\mathrm{subtraction}}=f_+(0)`$ is the subtraction constant at $`\nu =0`$ needed in a singly-subtracted dispersion relation for the reaction $`\gamma +p\gamma +p`$ and is given by the Thompson scattering limit, i.e., $`f_+(0)=\alpha /m=3.03\mu \mathrm{b}\mathrm{GeV}`$. In eq. (3), we have assumed that the Compton cross section rises as $`\mathrm{ln}^2(s)`$ at ultra-high energies.

The real and imaginary parts of eq. (3) are given by

$`\mathrm{Re}{\displaystyle \frac{4\pi }{\nu }}f_+`$ $`=`$ $`\beta \pi \mathrm{ln}s/s_0c\mathrm{cos}(\pi \mu /2)s^{\mu 1}+{\displaystyle \frac{4\pi }{\nu }}f_+(0)`$ (4)
$`\mathrm{Im}{\displaystyle \frac{4\pi }{\nu }}f_+`$ $`=`$ $`A+\beta \left[\mathrm{ln}^2s/s_0{\displaystyle \frac{\pi ^2}{4}}\right]+c\mathrm{sin}(\pi \mu /2)s^{\mu 1}+,`$ (5)
where $`s=2m\nu +m^22m\nu `$ and with $`c`$ being a real constant. Using equations (2), (4) and (5), we find the total cross section for high energy Compton scattering is given by

$$\sigma _{\mathrm{tot}}=A+\beta \left[\mathrm{ln}^2s/s_0\frac{\pi ^2}{4}\right]+c\mathrm{sin}(\pi \mu /2)s^{\mu 1},$$
(6)
and that $`\rho `$, the ratio of the real to the imaginary portion of the forward scattering amplitude, is given by

$$\rho =\frac{\beta \pi \mathrm{ln}s/s_0c\mathrm{cos}(\pi \mu /2)s^{\mu 1}+\frac{4\pi }{\nu }f_+(0)}{\sigma _{\mathrm{tot}}},$$
(7)
with $`f_+(0)=3.03\mu \mathrm{b}\mathrm{GeV}.`$ We will use units of $`\nu `$ in GeV and $`s`$ in GeV2, and cross sections in $`\mu `$b. We have to fit the 5 real constants $`A`$, $`c`$, $`\beta `$, $`s_0`$ and $`\mu `$. If we assume that the term in $`c`$ is a Regge descending term, then $`\mu =1/2`$.

The high energy behavior that was assumed by Damashek and Gilman in 1970 when they calculated $`\rho _{\gamma p}`$ numerically using dispersion relations was that the cross section approached a constant value asymptotically, i.e., they assumed that $`\sigma _{\mathrm{tot}}`$, for $`\nu \mathrm{}`$, was given by

$$\sigma _{\mathrm{tot}}=A+\frac{c^{}}{\nu ^{1/2}},$$
(8)
with $`A=96.6\mu \mathrm{b}`$ and $`c^{}=70.2\mu \mathrm{b}\mathrm{GeV}^{1/2}`$, with $`\nu `$ measured in GeV. However, today we know experimentally that the cross section rises at high energies and that the rising term can be fit by a $`\mathrm{ln}^2s`$ term, as expressed in eq. (6).

As a check on our real analytic amplitude analysis, the highest energy $`\rho `$ values of ref. , where Damashek and Gilman used dispersion relations assuming the asymptotic cross section of eq. (8), can be simply reproduced from eq. (7) with $`\beta =0`$ and eq. (8) by the relation

$$\rho _{\mathrm{analytic}}=\frac{70.2/\nu ^{1/2}+38.07/\nu }{96.6+70.2/\nu ^{1/2}}.$$
(9)
The $`\rho `$-values calculated from eq. (9) and the results from the singly-subtracted dispersion relation of ref.  agree to better than 2% over the energy range $`10\mathrm{GeV}\nu 20`$ GeV. The numerical agreement is excellent—the simplicity and ease of calculation at high energies using real analytical amplitudes compared to a dispersion relation analysis is clear.

The $`\rho `$ values in ref.  were only calculated up to $`\nu =20`$ GeV, which corresponds to a c.m.s. energy of $`\sqrt{s}=6.2`$ GeV. We now make an amplitude analysis for $`\rho _{\gamma p}`$ using a rising cross section, asymptotically going as $`\mathrm{ln}^2s`$ , by fitting the experimental cross sections $`\sigma _{\mathrm{tot}}`$ in the energy interval $`5\mathrm{GeV}\sqrt{s}200`$ GeV to the parameters $`A`$, $`c`$, $`\beta `$ and $`s_0`$ of eq. (6), using a Regge descending trajectory with $`\mu =0.5`$.

## 3 Results and Conclusions

Since cross sections for $`\gamma p`$ scattering are now available for c.m.s. energies up to 200 GeV, we made a $`\chi ^2`$ fit to the experimental $`\sigma _{\mathrm{tot}}(\gamma p)`$ data in the c.m.s. energy interval $`5\mathrm{GeV}\sqrt{s}200`$ GeV. We find a reasonable representation of the data using eq. (6), with a $`\chi ^2`$ per degree of freedom of 0.98 for 40 degrees of freedom, with the coefficients:

$`A=115\pm 13`$, $`c=35.1\pm 97`$, $`\beta =1.07\pm 0.57`$, $`s_0=68.4\pm 80.2`$ GeV2,

using a fixed value of $`\mu =0.5`$ (cross sections in $`\mu `$b when $`s`$ is in GeV2). This fit, plotted as a function of c.m.s. energy, gives the dashed cross section curve $`\sigma _{\mathrm{tot}}(\gamma p)`$ in Fig. 1, as well as the dashed $`\rho _{\gamma p}`$ curve in Fig. 2, using eq. (6) and eq. (7), respectively. Several remarks are in order. The experimental data are taken from the Particle Data Group compilation and include the only three high energy points—in the neighborhood of 200 GeV—that are available. The very large errors of $`\beta `$ and $`s_0`$ clearly reflect the fact that these few points also have large systematic errors ranging from $`7`$ to 21%. The accuracy is sufficient to show that the cross section rises, but not much more. Obviously, much more precise data at high energies are required for pinning down the $`\mathrm{log}^2s`$ parameters.

In order to show visually the sensitivity of $`\sigma _{\mathrm{tot}}(\gamma p)`$ and $`\rho _{\gamma p}`$ to the parameters of the fit, we have also plotted in Fig. 1 the dotted curve (a slight variation of the parameters of $`A`$ and $`c`$ within their errors), where we have set $`A=107`$ and $`c=134`$. This curve has as its $`\rho _{\gamma p}`$ analog the dotted curve of Fig. 2.

Using eq. (7), Fig. 2 shows our result for $`\rho _{\gamma p}`$ compared to $`\rho _{nn}`$, the $`\rho `$-value for nucleon-nucleon scattering found in ref. , as a function of the c.m.s. energy $`\sqrt{s}`$, in GeV. The solid curve is $`\rho _{nn}`$; the dashed line is the $`\rho _{\gamma p}`$ curve which corresponds to the central values $`A=115`$, $`c=35.1`$, $`\beta =1.07`$, $`s_0=68.4`$ GeV2; the dotted line is the $`\rho _{\gamma p}`$ curve which uses the slightly varied parameters $`A=107`$ and $`c=134`$.

The agreement between the slightly modified $`\rho _{\gamma p}`$ and $`\rho _{nn}`$ over the energy interval $`5\mathrm{GeV}\sqrt{s}200`$ GeV lends experimental support, in a model independent way, for the three factorization theorems of Block and Kadailov. Of course, precision cross section data in the region $`20\mathrm{GeV}<\sqrt{s}<200`$ GeV would enable us to strengthen this conclusion.

## 4 Acknowledgments

The author would like to thank Professor J. D. Jackson for valuable critical comments and suggestions.