text:
# The Importance of Fourier Phases for the Morphology of Gravitational Clustering

## 1 Introduction

The standard theory of the origin of the large-scale structure of the universe involves the assumption that the structure we see today grew by gravitational instability from initial small fluctuations in the density field. In most popular variants of this model, particularly those involving cosmic inflation, the initial fluctuations are of a particularly simple form known as a Gaussian random field. Gaussian random fields are useful because many properties of Gaussian random density fields can be calculated analytically (e.g. Bardeen et al. 1986). Some direct motivation for such an assumption emerges from inflationary models, wherein the fluctuations are generated by quantum oscillations of the scalar field driving inflation (the “inflation”). Even if inflation turns out to be incorrect, however, the Central Limit Theorem tends to produce Gaussian fluctuations in any linear process, so that they are the most generic form for small initial conditions and a natural default assumption.

One particularly interesting property of Gaussian random fields is that the requirement for the density contrast $`\delta (𝐱)=[\rho (𝐱)\rho _0]/\rho _0`$ to be a Gaussian random field is equivalent to that the real and imaginary parts of its Fourier components $`\stackrel{~}{\delta }_k`$, where

$$\delta (𝐱)=\stackrel{~}{\delta }(𝐤)\mathrm{exp}(i𝐤𝐱),$$
(1)
are independently distributed. In other words, the Fourier modes $`\stackrel{~}{\delta }(𝐤)`$,

$$\stackrel{~}{\delta }(𝐤)=|\stackrel{~}{\delta }(𝐤)|\mathrm{exp}(i\varphi __𝐤),$$
(2)
possess phases $`\varphi _𝐤`$ which are independently distributed and uniformly random on the interval $`[0,2\pi ]`$. As the density field is simply a sum over a large number of Fourier modes and if the phases of each of Fourier modes are random, the Central Limit Theorem guarantees that the resulting superposition of the one-point probability distribution $`𝒫(\delta )`$ is close to Gaussian and that all of the field’s joint probability distributions are multivariate Gaussians. The statistical properties of an isotropic Gaussian random field are then completely specified by its second-order statistical quantity: the covariance function, or alternatively, its power spectrum $`P(k)=\stackrel{~}{\delta }^2(k)`$.

In the framework of gravitational instability, the growth of fluctuations can be understood analytically when the density fluctuation amplitude is small compared to the mean density; the linear perturbation theory tells us that each Fourier mode grows with the same rate independent of wavenumber and the statistical distribution of $`\delta `$ remains constant except its variance.

The linear theory breaks down when $`\delta ^2`$ is comparable with unity or beyond, and different Fourier modes start coupling. One way to look at the mode coupling is that $`\delta `$ is always constrained to value $`\delta 1`$. When the perturbation is small, the tail of $`𝒫(\delta )`$ in the part of $`\delta <1`$ assigned by Gaussian distribution is negligible, because probability of $`\delta <1`$ is small. When the density field evolves beyond the linear regime, i.e., $`\sigma ^2\delta ^21`$, a long tail at high $`\delta `$ is generated while lower bound is confined at $`\delta =1`$. Gaussian condition is therefore invalid, in that mode coupling effect causes the initial condition to skew. Terms in the evolution

equations for the Fourier modes that represent coupling between different modes are of second (or higher) order in $`\delta `$ and these are neglected when first-order perturbation theory is considered. Phases of Fourier modes therefore are therefore coupled together in a way which is yet to be fully elucidated but which has been recently investigated by Chiang & Coles (2000) and Coles & Chiang (2000).

One of the reasons for studying Fourier phases in depth is the question of possible primordial non-Gaussianity of the initial density distribution of the Universe. For example, there have been claims of non-Gaussianity from analysis results in the COBE DMR sky maps of cosmic microwave background using a number of different diagnostics: bispectrum analysis (Ferreira et al. 1998); the wavelet transform (Hobson et al. 1999); and Minkowski functionals (Schmalzing & Gorski 1998). The most direct practical approach, however, is through the distribution of Fourier phases. This diagnostic can also avoid the subtlety of more standard methods that depend on the distribution function $`𝒫(\delta )`$(Scherrer et al. 1991), in that a density field with Gaussian single-point density distribution $`𝒫(\delta )`$ is not necessarily a Gaussian field. Phase information is also important as a statistical diagnostic of non-linearity when the evolution of clustering in the non-linear regime which is generally intractable analytically. Indeed, , as we shall show, it is also closely related to the morphology of gravitational clustering and through this to the dynamical origin of structure.

The layout of this paper is as follows. In section 2 we present a visual demonstration of the link between Fourier phases and clustering morphology. Section 3 contains a brief discussion of some of the properties of phases and some pitfalls that must be avoided when using them as part of a clustering descriptor. In Section 4 we introduce a cross-correlation parameter $`S`$ and a rank-correlation $`\tau `$ as quantitative measures of the agreement between two distributions. In Section 5 we display the results of the correlation tests between sample distributions and phase-based reconstructions, and between simulations evolving from different initial power spectra but the same initial phase set. A brief discussion of the results follows in Section 6.

## 2 Visual Demonstration

To give a qualitative, visual description of the key ideas in this paper consider Fig. 1, in which we isolate the role of phases in determining clustering morphology. Plates (a) and (b) are two example realisations from two 2D N-body experiments, evolving from different initial power-law power spectra and different initial phase sets. Plate (a) is evolved from power spectral index $`n=1`$, and (b) from $`n=1`$. We perform a Fourier transform on both realisations, e.g., $`\delta ^a(𝐱)=\stackrel{~}{\delta }_𝐤^a\mathrm{exp}(i𝐤𝐱)`$, where $`\stackrel{~}{\delta }_𝐤^a=|\stackrel{~}{\delta }_𝐤^a|\mathrm{exp}(i\varphi _𝐤^a)`$. Plate (c) is obtained by taking the inverse Fourier transform from the combination of the phases of Fourier modes from (a), and the amplitudes from (b), i.e., $`^1[|\stackrel{~}{\delta }_𝐤^b|\mathrm{exp}(i\varphi _𝐤^a)]`$. Plate (d) is from the phases from (b) but amplitudes from (a). Therefore, plate (a) and (c) share the same phase configuration, so do (b) and (d). It is easy to see resemblance between (a) and (c), and between (b) and (d). Note that plate (a) and (d) have the same power spectrum, or equivalently, two-point correlation function, as do plate (b) and (c). Plate (e) is the inverse Fourier transform from only the phases from (a): $`^1[\mathrm{exp}(i\varphi _𝐤^a)]`$; it therefore retains only the phase information from (a). In plate (f), each mode keeps the same amplitude so its power spectrum is unchanged (i.e. the same as plate (a) and (d)) but the phases are redistributed randomly among the modes before inverse Fourier transform. Again, plate (e), the phase-only reconstruction (hereafter phase-only reconstruction), resembles the original distribution plate (a), and plate (c), the same-phase amplitude-swapped reconstruction (hereafter amplitude-swapped reconstruction), but (f) which has the same power spectrum as (a) but random phases (hereafter random-phase reconstruction), is featureless. This experiment suffices to show very clearly how phases determine morphology.

Another interesting property of Fourier phases in clustering morphology is that two realisations evolving from the same initial random phase set, though different power-law power spectra, will have their extrema at the same locations. The fact that Plate (a) and (b) are evolved from different initial phase sets is for demonstration purpose.

## 3 Quantifying Phase Shifts

Fourier phases reflect the locations of spatial ‘events’ more than Fourier amplitudes do. For a single hypothetical spike, represented by the Dirac $`\delta `$-function $`\delta _D(xx_0)`$, the amplitudes are constant and phases are $`kx_0`$. This has a white-noise spectrum, but very strong phase correlation. The phase configuration of a spike is very similar to that of a single density peak evolved from 1D Zel’dovich approximation; see Chiang & Coles (2000) for details. A translation $`x^{^{}}`$ of the Dirac $`\delta `$-function density field, $`\delta _D(xx_0x^{})`$, has no effect on the Fourier amplitudes, but the phases now become $`k(x_0+x^{})`$, which suffers a shift by a linear term proportional to wave number $`k`$. This dependency of phases on the choice of origin means that some care must be taken when trying to extract meaningful information. Some previous studies focused on the evolution of individual phases away from their initial values (Ryden & Gramann 1991; Soda & Suto 1992; Jain & Bertschinger 1998). The mean deviation from the initial phase can be defined as

$$\mathrm{\Delta }\varphi (k,a)=|\mathrm{\Delta }\varphi (𝐤,a)|=|\varphi (𝐤,a)\varphi (𝐤,a_i)|,$$
(3)
where the averages are performed over the different modes within a shell in k-space whose wave numbers lie in the range $`k0.5<|𝐤|<k+0.5`$. As long as there are enough modes, the maximal value of $`\mathrm{\Delta }\varphi (k,a)`$ is $`2\pi /3`$ (Jain & Bertschinger 1998), which can be understood as a variable obtained from the change from a random field to a hypothetical spike, i.e., from random phases to a very ordered state, $`\varphi (k,a)=kx_0`$( mod $`2\pi `$), where $`x_0`$ is the location of the spike and takes on any values except zero. This statistic will change following the translation of ‘events’, e.g., $`\mathrm{\Delta }\varphi (k,a)=\pi `$ if the hypothetical spike is shifted to the origin.

## 4 Cross and Rank Correlations

Although the examples displayed in Section 2 serve to demonstrate the importance of phases, it is necessary to measure objectively how well the phase-based reconstructions compare not only with the original distribution, but also with reconstruction from, say, random phases. The tool we use here is the cross-correlation statistic introduced by Coles, Melott & Shandarin (1993). A correlation coefficient is defined as

$$S_{\delta \delta ^r}=\frac{(\delta _{ij}\overline{\delta })(\delta _{ij}^r\overline{\delta }^r)}{(\delta _{ij}\overline{\delta })^2^{1/2}(\delta _{ij}^r\overline{\delta }^r)^2^{1/2}},$$
(4)
where $`\delta _{ij}`$ and $`\delta _{ij}^r`$ represent two density distribution and $`\overline{\delta }`$ and $`\overline{\delta }^r`$ are their mean densities, respectively. The indices $`i`$ and $`j`$ label the pixel positions in the two-dimensional simulations we use here for illustration. The parameter $`S`$ compares the density value of each grid point at $`(i,j)`$ in the original distribution $`\delta `$ with the corresponding grid point of the reconstruction $`\delta ^r`$; averages are taken over all grid points. $`S_{\delta \delta ^r}=1`$ denotes ‘completely positive correlation’, a perfect agreement between $`\delta `$ and $`\delta ^r`$, and a value of zero indicates the two distributions are uncorrelated. What is useful about this test is that it compares the morphology between two distributions, point by point, but does not take into account their variances. This is because $`S_{\delta \delta ^r}=1`$ when $`\delta =C\delta ^r`$, $`C`$ being a constant. If the structures of two distributions are different, $`|S|`$ is less than unity.

Cross-correlation tests the spatial correspondence of the locations of ‘events’, such as clusters of points or edges, between two distributions. Even the comparison by a small displacement between two identical periodic distributions will result in low value of $`S`$. It should be pointed out, therefore, that this grid-by-grid test is severe. As we have explained, phases are related to preservation of locations of ‘events’. Relative magnitudes of the ‘events’, however, are not preserved when the information of Fourier amplitudes is partially or totally lost. To relax the test, comparisons are made not only between raw density distributions $`\delta `$ and $`\delta ^r`$, but also between smoothed distributions $`\delta (𝐱,R)`$ and $`\delta ^r(𝐱,R)`$. For this purpose, a Gaussian window function is chosen to smooth the field,

$$\delta (𝐱,R)=d^2𝐱^{^{}}\delta (𝐱^{^{}})(\sqrt{2\pi }R)^2\mathrm{exp}(\frac{|𝐱𝐱^{^{}}|^2}{2R^2}).$$
(5)
As well as the linear cross-correlation $`S`$, we also use a non-parametric rank-correlation, Kendall’s $`\tau `$ parameter(Kendall and Gibbons 1990). This parameter uses the relative ordering of ranks to measure the degree of agreement between two compared distributions. For any two grid points $`(i,j)`$ and $`(i^{^{}},j^{^{}})`$, the given value is $`+1`$ if the relative ordering of the density values of $`\delta (i,j)`$ and $`\delta (i^{^{}},j^{^{}})`$ is the same as $`\delta ^r(i,j)`$ and $`\delta ^r(i^{^{}},j^{^{}})`$, i.e., either both decreasing or both increasing from $`(i,j)`$ to $`(i^{^{}},j^{^{}})`$; $`1`$ when one pair is increasing and the other is decreasing. The $`\tau `$ parameter compares $`N(N1)/2`$ pairs from total $`N^2`$ grid points and is normalized to $`[1,1]`$. Rank correlation measures such as this do not look specifically for linear association between the image and reconstruction, but for one-to-one ordering of the values in one relative to those in the other. If there is strong non-linear association, then $`\tau `$ can be close to unity even though $`S`$ may be small. Kendall’s $`\tau `$ parameter is used here as an auxiliary test, which is again “softened” with a Gaussian window function as in eq. (5).

## 5 Results

Fig. 2 shows cross-correlation coefficient $`S`$ and rank-correlation coefficient $`\tau `$ drawn against smoothing scale $`R`$ in computer grid units. Each panel includes correlations between the sample distribution and phase-only, amplitude-swapped reconstructions. Correlations of random-phase reconstructions are calculated only in $`S`$. On the top-left panel of Fig. 2, the sample distribution evolved from spectral index $`n=1`$(Fig.1a) is compared with phase-only reconstruction(Fig.1e), and with amplitude-swapped reconstruction (Fig.1c). The random-phase reconstruction (Fig.1f) is also compared for reference. The sample distribution of lower-left panel is evolved from $`n=1`$(Fig.1b). Even before smoothing, the higher degree of cross-correlation from phase-only and amplitude-swapped reconstructions than from random-phase reconstruction shows the ability of phases to retain the morphology of the sample distributions. Non-resemblance of random-phase reconstructions is shown through low value of $`|S|`$. That smoothing of the reconstructions increases (for amplitude-swapped it even decreases) certain level of $`S`$ indicates the limitation of phases on reconstruction with the loss of amplitude information.

The rank-correlations, on the other hand, show that phase-only reconstructions have higher agreement than amplitude-swapped ones on relative density magnitudes, which is due to the ‘twist’ on the latter’s power spectra. Surely, the loss of the Fourier amplitude information destroys part of the signal, as shown in the rank-correlation of the right panels. On the morphology of gravitational clustering, however, phases play a much more important role than amplitudes.

It emphasizes this point still further to test the correspondence between simulations evolving from different initial power spectra but the same initial phases. In Fig. (3) we compare the morphology for some relevant examples in terms of cross-correlation coefficient $`S`$ deployed above. In particular, we show realisations obtained by evolving $`N`$-body experiments from different 2D initial power-law power spectra $`n=1`$, $`0`$ $`1`$ and $`2`$. The left panel is the comparison of stage $`c`$ in which the scale of non-linearity is $`k_{NL}=64k_f`$, and the right panel is that of later stage $`f`$ with $`k_{NL}=8k_f`$(Chiang & Coles 2000, Beacom et al. 1991). We use the fundamental mode $`k_f=2\pi /L`$ as length unit, where $`L`$ is the length of the side of the simulation square.

Before smoothing, there is little or no correspondence between any of the realisations. This is due to the severe point-by-point test of the cross-correlation. After smoothing, however, there is a dramatic improvement between in the correspondence between $`n=1,0`$, $`n=0,1`$, and between $`n=1,2`$. We can examine the morphology evolving from the same phase set by the characteristic scale $`k_{NL}`$. For example, $`k_{NL}`$ for stage $`c`$ corresponds to 8 computer grid units, thus the smoothing on scales beyond 8 grid units erases non-linearities, while the larger scale structure remains in the linear regime. The linear growth of the density fluctuations set up by the same initial phases depends only on time. The comparison therefore indicates the intrinsic difference in clustering morphology arising from different initial power spectra. What we see on the scales beyond $`k_{NL}`$ is that the correspondence is low between $`n=1`$ and $`n=1,2`$.

The correspondence deteriorates when the difference between the spectral indices increases, i.e., the difference in intrinsic morphology is significant, or when the evolution goes into highly non-linear regime, where particles move away from their initial Lagrangian position and interact with each other non-linearly. Simulations evolving from the same initial phase configuration tend to develop nonlinear structures at or near the same spatial locations, but these structures appear with different contrast when the initial spectra are different. For example, filaments appear in the nonlinear regime in all cases, but for spectra with large $`n`$ these tend to be less well defined and broken up into clumps. Our statistic $`S`$ takes into account both the position and amplitude of structures that form so it indicates a deteriorating agreement for very different spectra. Nevertheless, it is clear that there is strong imprint upon the morphology of the initial phases resulting from the process of gravitational clustering. This can also be seen visually in the pictures shown by Beacom et al. (1991).

The cross-correlation coefficient being a single number, so one cannot infer from it precisely how the phase-based reconstructions perform relative to the morphology of the original distribution. One intriguing question is how the density regions from reconstructions can be compared to those from the original distributions. Is it $`\delta \delta ^r`$ that produces high value of $`S`$? From rank-correlation coefficient $`\tau `$, this is is not the case. In order to get more information between the two distributions, we also made a grid-by-grid comparison between them. A scatter diagram can be drawn in logarithmic scales for the comparison between each cell density from phase-based reconstructions against the corresponding one from the original distribution. If $`\delta \delta ^r`$, the points scatter along a straight line of slope equal to unity, and the relative magnitude of ‘events’ is preserved by a linear scale factor. Fig. 4 shows the scatter diagram between original distribution and phase-only, amplitude-swapped reconstructions for $`n=1`$ and $`n=1`$. Not all the points scatter along the straight line, as is suggested. Instead, most points scatter horizontally, particularly for phase-only cases. Also notice also that, in both reconstruction cases, the density values span only roughly one order of magnitude and the points that do not align well correspond to very small fluctuations. High values of $`S`$ result from phase configuration preserving the locations of high-density ‘events’ with large magnitudes in the original structure. The relatively flat phase-only reconstructed density distribution is caused by the flat power spectrum, in which the Fourier amplitude for each mode is squashed into unity. The lower two panels are scatter diagrams of the amplitude-swapped reconstructions against the originals, i.e., in Fig. 1 between (c) and (a) and between (d) and (b). The difference is amplitude-swapped reconstructions have the power spectra from the alternative sample distributions. The $`n=1`$ amplitude-swapped reconstruction has the power spectrum from the original distribution of $`n=1`$, which gives more power on small scales in the reconstruction, the $`n=1`$ reconstruction, on the other hand, has more power on large scales. This adverse effect causes the points of the low fluctuations to spread on the scatter diagram, which, however, doesn’t decrease much the cross-correlation coefficient $`S`$ between phase-only and amplitude-swapped reconstructions. This is in accord with with our visual impression that our eyes pick up the maxima between the distributions for comparison.

In Fig. 5 the scatter diagrams for the realisations evolving from the same phase set are also produced. Here only two cases are chosen, i.e., between $`n=1`$ and $`0`$, and between $`n=1`$ and $`2`$. The left two panels are comparison between $`n=1`$ and $`n=0`$, and the scattering is expected, which nonetheless follows the straight line. For the right panels, at early stage $`c`$, the deviation from the straight line is more isotropic, hence the correspondence is low. At late stage $`f`$, there is even more deviation, which comes from the mapping from the low-density regions of stage $`f`$, $`n=2`$. The reason is as follows. For the realisation of $`n=2`$ to reach the same level of non-linearity, the variance is higher than that of $`n=1`$, and there are a substantial number of voids. After smoothing to create continuous distribution, the void regions are smoothed as low-density regions, which is seen in the bottom-right panel.

## 6 Conclusion

We have shown the importance of Fourier phases on morphology of gravitational clustering via qualitative and quantitative demonstrations. It is interesting to remark upon the similarity of the results we have obtained here and those presented in Coles et al. (1993). The latter authors were interested in the ability of simple analytic methods to reproduce the clustering displayed by full N-body computations. They showed in particular that the Zel’dovich approximation (Zel’dovich 1970) could reproduce the full numerical results quite well, with cross-correlations similar to those we have found here. The Zel’dovich approximation works as well as it does in this respect because it places the caustic surfaces forming sheets and filaments near to the correct location in the N-body experiment. In other words, the Zel’dovich approximation has a high phase fidelity. It is less good at getting the amplitudes right. As we have shown, however, phases dominate the morphology.

With a new colour representation technique (Coles & Chiang 2000), phase information can be used to distinguish between non-Gaussianity induced by gravitational clustering and that by other mechanisms. For example, in the web page

```
http://www.nottingham.ac.uk/~ppzpc/phases/cmb.html

```

the phase configuration of the non-Gaussian temperature fluctuations on the CMB sky induced by cosmic strings can be seen to be intrinsically different from that of hierarchical clustering developed in N-body simulations. A new algorithm based on the analysis of the phase distribution of Fourier components is recently devised to extract noise from point sources from the CMB signal(Naselsky et al. 2000), which is another example practice of the close link between Fourier phases and morphology.

Future large-scale galaxy redshift surveys and microwave sky maps will reveal much morphological information about large-scale structure in the Universe. Existing statistical technology, however, is still dominated by second-order methods that are blind to phase information. Our success in extracting this information will therefore depend on the development of statistical methods sufficiently sensitive to the key ingredient: the distribution of Fourier phases.

## Acknowledgments

I thank Peter Coles for useful suggestions. This work was partly supported by Danmarks Grundforskningsfond through its support for TAC.


format: text

text:
# Calculations of three-body observables in 8B breakup

## I Introduction

Projectile breakup is an important reaction channel in the scattering of weakly bound nuclei. An accurate treatment of breakup is therefore a major ingredient in attempts to understand the properties of light exotic radioactive nuclei from reaction studies. The number of published experimental breakup studies, and also their accuracy, has increased rapidly. These include reactions in which both nuclear and Coulomb breakup effects are expected to be significant, e.g. . Until recently, the low intensities of available rare isotope beams has meant that many of the experiments were either designed to measure inclusive cross sections, with incomplete kinematics, or did not have adequate statistics to allow the extraction of exclusive observables. The cross sections extracted from the measurements were often integrated over fragment energies or angles, or both, and inevitably some details of the physical process were lost as a result. This is no longer the situation. Secondary beam intensities are becoming sufficiently high that coincidence experiments are now practical and, in the future, data will more routinely require a fully three- or more-body study, e.g. . The need for precise theoretical predictions of the breakup of two-body projectiles, and of their three-body observables, is the primary motivation for this work.

Theoretical reaction models, which treat breakup as an excitation of the projectile to a two-body continuum state, most naturally express their results as cross sections describing the center of mass (c.m.) and relative motions of the dissociated system, using 2-body kinematics. It has therefore been common for the experimental data to also be transformed to the c.m. frame, for ease of comparison, e.g. the theoretical calculations of  and the experimental data of . This process is ambiguous in the case of inclusive data. Much more important is that the three-body cross sections are explicitly coherent in contributions from different spin-parity excitations of the projectile and so have the potential to offer a far greater insight into the projectile structure and the reaction mechanism. An excellent example of this is the interference observed  in the cross section of the 7Be fragments, as a function of their component of momentum parallel to the beam direction, following 8B breakup on a heavy target at 44 MeV/nucleon.

In this paper we present calculations which are performed using full three-body kinematics. These calculations are carried out within the framework of the coupled discretised continuum channels (CDCC) methodology, e.g. , for breakup reactions of two-body projectiles. The interference between different excitation channels is shown to be significant for assessing the convergence of the calculations and those breakup excitations which contribute. The methods presented are applied to the breakup of 8B on a 58Ni target at $`E_{lab}=25.8`$ MeV, for which new measurements have been reported . We compare the results of the full CDCC analysis, and also distorted waves Born approximation (DWBA) and truncated coupled channels calculations, with these available data for the laboratory angle and energy distributions of the 7Be fragments. The calculations of Refs.  showed the importance of higher-order breakup couplings, the couplings between continuum states, upon the 8B center of mass cross section angular distribution. We will show in this work that these higher-order effects are manifest even more significantly in the energy distributions of the 7Be fragments following breakup.

## II Theoretical Considerations

We consider the breakup reaction $`pc+v`$ in which the projectile nucleus $`p`$ is a bound state of a core particle $`c`$, of spin $`I`$ and projection $`\mu `$, and a valence particle $`v`$, of spin $`s`$ and projection $`\sigma `$. These particles are, presently, assumed structureless and so their internal wave functions are represented by the spinors $`𝒳_I`$ and $`𝒳_s`$. The total angular momentum of the ground state of $`p`$ is $`J_p`$, with projection $`M`$, the relative orbital angular momentum of the two constituents is $`\mathrm{}_0`$, and their separation energy is $`_0(>0)`$. The incident wave number of the projectile in the center of mass (c.m.) frame of the projectile and target is $`\stackrel{}{K}_0`$ and the co-ordinate $`z`$–axis is chosen in the incident beam direction. The target $`t`$ is assumed to have spin zero and no explicit target excitation is included. Target excitation is therefore present only through the complex effective interactions of $`c`$ and $`v`$ with the target. Our three-body solution of the Schrödinger equation calculates an approximate description of the projection of the full many-body $`p+t`$ wave function onto the ground states of the target, core and valence nuclei. This three-body wave function is denoted $`\mathrm{\Psi }_{\stackrel{}{K}_0M}(\stackrel{}{r},\stackrel{}{R})`$ where $`\stackrel{}{R}`$ is the position of the c.m. of $`p`$ relative to the target and $`\stackrel{}{r}`$ is the position of $`v`$ relative to the core $`c`$. The particle masses are $`m_p=m_c+m_v`$ and $`m_t`$.

### A Construction of continuum bin states

In the coupled discretised continuum channels (CDCC) method  the breakup of $`p`$ is assumed to populate a finite set of selected $`c+v`$ excited configurations, with quantum numbers $`J_p^{},\mathrm{},j`$, where $`\stackrel{}{j}=\stackrel{}{\mathrm{}}+\stackrel{}{s}`$ and $`\stackrel{}{J}_p^{}=\stackrel{}{j}+\stackrel{}{I}`$. Here, each of these spin-parity excitations will be assumed diagonal in all of these angular momentum labels. The excitations are also assumed to extend to some maximum relative energy $`_{max}(J_p^{})`$ or wave number $`k_{max}`$. This momentum range is then divided into a number $`𝒩(J_p^{})`$ of intervals or bins, each with a width $`\mathrm{\Delta }k_i=[k_ik_{i1}]`$. We label each such momentum bin by $`\alpha (i,J_p^{},\mathrm{},j,s,I)`$.

In each of these relative motion bins a single representative wave function is constructed from those $`c+v`$ scattering states $`f_\alpha (k,r)`$ internal to the bin, with assumed angular momentum coupling

$`\widehat{\varphi }_\alpha ^M^{}(\stackrel{}{r})=\left[\left[Y_{\mathrm{}}(\widehat{\stackrel{}{r}})𝒳_s\right]_j𝒳_I\right]_{J_p^{}M^{}}u_\alpha (r)/r.`$ (1)
The radial functions $`u_\alpha `$ are square integrable and are a superposition

$`u_\alpha (r)=\sqrt{{\displaystyle \frac{2}{\pi N_\alpha }}}{\displaystyle _{k_{i1}}^{k_i}}g_\alpha (k)f_\alpha (k,r)𝑑k`$ (2)
of the scattering states, eigenstates of the $`c+v`$ internal Hamiltonian $`H_p`$, with weight function $`g_\alpha (k)`$. $`N_\alpha =_{k_{i1}}^{k_i}|g_\alpha (k)|^2𝑑k`$ is a normalisation constant. The $`f_\alpha `$ are defined here such that, for $`r\mathrm{}`$,

$$f_\alpha (k,r)\left[\mathrm{cos}\delta _\alpha (k)F_{\mathrm{}}(kr)+\mathrm{sin}\delta _\alpha (k)G_{\mathrm{}}(kr)\right],$$
(3)
where $`k\alpha `$ and $`F_{\mathrm{}}`$ and $`G_{\mathrm{}}`$ are the regular and irregular partial wave Coulomb functions. So the $`f_\alpha `$ are real when using a real $`c+v`$ two-body interaction. An optimal discretisation of the continuum requires a consideration of the number, the boundaries $`k_i`$, the widths $`\mathrm{\Delta }k_i`$ and the weights $`g_\alpha `$ in the bins, which may depend on the $`J_p^{}`$ configuration. Energy conservation relates the $`c+v`$ c.m. wave numbers $`K_\alpha `$ and corresponding bin state energies $`\widehat{}_\alpha `$, as

$`{\displaystyle \frac{\mathrm{}^2K_\alpha ^2}{2\mu _{pt}}}+\widehat{}_\alpha ={\displaystyle \frac{\mathrm{}^2K_0^2}{2\mu _{pt}}}_0,`$ (4)
where we define each bin energy by $`\widehat{}_\alpha =\widehat{\varphi }_\alpha |H_p|\widehat{\varphi }_\alpha `$ and where $`\mu _{pt}`$ is the projectile-target reduced mass.

For non-$`s`$-wave bins typically one uses $`g_\alpha (k)=1`$ for a non-resonant continuum in which case $`N_i=\mathrm{\Delta }k_i`$ and $`\widehat{}_i=\mathrm{}^2\widehat{k}_i^2/(2\mu _{cv})`$ with $`\widehat{k}_i^2=[k_i^3k_{i1}^3]/(3\mathrm{\Delta }k_i)`$. For $`s`$-wave bins it is an advantage to use $`g_\alpha (k)=k`$. This stabilises the extraction of the three-body transition amplitude at low relative breakup energies, discussed later in Eq. (15). In this case $`N_i=\widehat{k}_i^2\mathrm{\Delta }k_i`$ and the bin energies are $`\widehat{}_i=\mathrm{}^2[k_i^5k_{i1}^5]/(10\mu _{cv}\mathrm{\Delta }k_i\widehat{k}_i^2)`$.

### B Coupled channels amplitudes

These bin states $`\widehat{\varphi }_\alpha `$ provide an orthonormal relative motion basis for the coupled channels solution of the three-body $`c+v+t`$ wave function. The bins and the coupling potentials $`\widehat{\varphi }_\alpha |U(\stackrel{}{r},\stackrel{}{R})|\widehat{\varphi }_\beta `$ are constructed, and the coupled equations are solved, using the coupled channels code fresco . Here $`U(\stackrel{}{r},\stackrel{}{R})`$ is the sum of the interactions of $`c`$ and $`v`$ with the target, which is expanded to a maximum specified multipole order $`q`$. The coupled equations solution generates the (two-body) scattering amplitudes, summed over partial waves, for populating each bin state $`J_p^{},M^{}`$ from initial state $`J_p,M`$, as a function of the angle $`\theta _K`$ of the c.m. of the excited projectile in the c.m. frame

$`\widehat{}_{M^{}M}(\stackrel{}{K}_\alpha )`$ $`=`$ $`{\displaystyle \frac{4\pi }{K_0}}\sqrt{{\displaystyle \frac{K_\alpha }{K_0}}}{\displaystyle \underset{LL^{}J}{}}(L0J_pM|JM)(L^{}MM^{}J_p^{}M^{}|JM)`$ (5)
$`\times `$ $`\mathrm{exp}(i[\sigma _L+\sigma _L^{}]){\displaystyle \frac{1}{2i}}\widehat{𝒮}_{LJ_p:L^{}J_p^{}}^J(K_\alpha )Y_L^0(\widehat{\stackrel{}{K}}_0)Y_L^{}^{MM^{}}(\widehat{\stackrel{}{K}}_\alpha ).`$ (6)
Here $`\sigma _L`$ and $`\sigma _L^{}`$ are the Coulomb phases appropriate to the initial and final state c.m. energies and the $`\widehat{𝒮}_{LJ_p:L^{}J_p^{}}(K_\alpha )`$ are the partial wave $`S`$–matrices for exciting bin state $`\alpha `$ with c.m. wave number $`K_\alpha `$. When calculated using fresco , these amplitudes are expressed in a coordinate system with $`x`$–axis in the plane of $`\stackrel{}{K}_0`$ and $`\stackrel{}{K}_\alpha `$, such that the azimuthal angle $`\varphi _{K_\alpha }`$ of $`\stackrel{}{K}_\alpha `$ is zero. When discussing three-body observables, it is more convenient to define the coordinate system with respect to the fixed positions of the detectors in the laboratory. For such a general $`x`$–coordinate axis the coupled channels amplitudes must subsequently be multiplied by $`\mathrm{exp}(i[MM^{}]\varphi _K)`$, with $`\varphi _K`$ referred to the chosen $`x`$-axis.

For use in the following, the two-body inelastic amplitudes of Eq. (6) are re-normalised to that of the $`T`$-matrix by removal of their two-body phase space factors, so that

$`\widehat{𝒯}_{M^{}M}^\alpha (\stackrel{}{K}_\alpha )={\displaystyle \frac{2\pi \mathrm{}^2}{\mu _{pt}}}\sqrt{{\displaystyle \frac{K_0}{K_\alpha }}}\widehat{}_{M^{}M}(\stackrel{}{K}_\alpha ).`$ (7)
Throughout, we adopt scattering state and $`T`$-matrix normalisations such that, asymptotically, the plane wave states $`\mathrm{exp}(i\stackrel{}{k}\stackrel{}{r})`$ that enter are multiplied by unity.

It follows that the inelastic differential cross section angular distribution, in the center of mass frame, for excitation of a given bin state is

$`{\displaystyle \frac{d\sigma (\alpha )}{d\mathrm{\Omega }_K}}`$ $`=`$ $`{\displaystyle \frac{1}{2J_p+1}}\left[{\displaystyle \frac{\mu _{pt}}{2\pi \mathrm{}^2}}\right]^2{\displaystyle \frac{K_\alpha }{K_0}}{\displaystyle \underset{MM^{}}{}}\left|\widehat{𝒯}_{M^{}M}^\alpha (\stackrel{}{K}_\alpha )\right|^2`$ (8)
$`=`$ $`{\displaystyle \frac{1}{2J_p+1}}{\displaystyle \underset{MM^{}}{}}\left|\widehat{}_{M^{}M}(\stackrel{}{K}_\alpha )\right|^2.`$ (9)
### C Three-body breakup amplitudes

Less obvious is the relationship of the CDCC two-body inelastic amplitudes $`\widehat{𝒯}_{M^{}M}^\alpha (\stackrel{}{K}_\alpha )`$ to the breakup transition amplitudes $`T_{\mu \sigma :M}(\stackrel{}{k},\stackrel{}{K})`$ from an initial state $`J_p,M`$ to a general physical three-body final state of the constituents . This is needed to make predictions for experiments with general detection geometries, since each detector configuration and detected fragment energy involves a distinct final state c.m. wave vector $`\stackrel{}{K}`$, breakup energy $`_k`$, and relative motion wave vector $`\stackrel{}{k}`$.

To clarify this connection we make the CDCC approximation to the exact (prior form) breakup transition amplitude, by replacing the exact $`c+v+t`$ three-body wave function, $`\mathrm{\Psi }_{\stackrel{}{K}_0M}(\stackrel{}{r},\stackrel{}{R})`$, by its CDCC approximation $`\mathrm{\Psi }^{CD}`$, as

$`T_{\mu \sigma :M}(\stackrel{}{k},\stackrel{}{K})=\varphi _{\stackrel{}{k}\mu \sigma }^{()}(\stackrel{}{r})e^{i\stackrel{}{K}\stackrel{}{R}}|U(\stackrel{}{r},\stackrel{}{R})|\mathrm{\Psi }_{\stackrel{}{K}_0M}^{CD}(\stackrel{}{r},\stackrel{}{R}).`$ (10)
Here $`\varphi _{\stackrel{}{k}\mu \sigma }`$ is the $`c+v`$ final state. Upon inserting the set of all included bin-states, which are assumed complete within the model space used, then

$`T_{\mu \sigma :M}(\stackrel{}{k},\stackrel{}{K})={\displaystyle \underset{\alpha ,M^{}}{}}\varphi _{\stackrel{}{k}\mu \sigma }^{()}|\widehat{\varphi }_\alpha ^M^{}\widehat{\varphi }_\alpha ^M^{}e^{i\stackrel{}{K}\stackrel{}{R}}|U(\stackrel{}{r},\stackrel{}{R})|\mathrm{\Psi }_{\stackrel{}{K}_0M}^{CD}(\stackrel{}{r},\stackrel{}{R}),`$ (11)
where the sum is over all bins $`\alpha `$ which contain wave number $`k`$. We should now recognise that the matrix elements $`\widehat{𝒯}_{M^{}M}^\alpha (\stackrel{}{K}_\alpha )`$ of Eq. (7), obtained from the coupled channels solution, are precisely the transition matrix elements appearing in Eq. (11), i.e.

$`\widehat{𝒯}_{M^{}M}^\alpha (\stackrel{}{K}_\alpha )=\widehat{\varphi }_\alpha ^M^{}e^{i\stackrel{}{K}_\alpha \stackrel{}{R}}|U(\stackrel{}{r},\stackrel{}{R})|\mathrm{\Psi }_{\stackrel{}{K}_0M}^{CD}(\stackrel{}{r},\stackrel{}{R}),`$ (12)
but calculated on the grid of $`\theta _\alpha `$ and $`K_\alpha `$ values. For the first term in Eq. (11) one obtains

$`\varphi _{\stackrel{}{k}\mu \sigma }^{()}|\widehat{\varphi }_\alpha ^M^{}={\displaystyle \frac{(2\pi )^{3/2}}{k\sqrt{N_\alpha }}}{\displaystyle \underset{\nu }{}}(i)^{\mathrm{}}(\mathrm{}\nu s\sigma |jm)(jmI\mu |J_p^{}M^{})\mathrm{exp}[i\overline{\delta }_\alpha (k)]g_\alpha (k)Y_{\mathrm{}}^\nu (\widehat{\stackrel{}{k}}),`$ (13)
where $`\overline{\delta }_\alpha (k)=\delta _\alpha (k)+\sigma _\alpha (k)`$ is the sum of the nuclear and Coulomb phase shifts for $`c+v`$ scattering at relative wave number $`k`$. It follows that the three-body breakup $`T`$–matrix can be written

$`T_{\mu \sigma :M}(\stackrel{}{k},\stackrel{}{K})`$ $`=`$ $`{\displaystyle \frac{(2\pi )^{3/2}}{k}}{\displaystyle \underset{\alpha \nu }{}}(i)^{\mathrm{}}(\mathrm{}\nu s\sigma |jm)(jmI\mu |J_p^{}M^{})\mathrm{exp}[i\overline{\delta }_\alpha (k)]`$ (14)
$`\times `$ $`Y_{\mathrm{}}^\nu (\widehat{\stackrel{}{k}})g_\alpha (k)T_{M^{}M}(\alpha ,\stackrel{}{K}).`$ (15)
Here the $`T_{M^{}M}(\alpha ,\stackrel{}{K})`$ will be interpolated from the matrices $`\widehat{𝒯}_{M^{}M}^\alpha (\stackrel{}{K}_\alpha )`$, available on the calculated $`K_\alpha `$ and $`\theta _{K_\alpha }`$ grid. Specifically, for each value of $`\stackrel{}{K}`$, we evaluate

$`T_{M^{}M}(\alpha ,\stackrel{}{K})=\mathrm{exp}(i[MM^{}]\varphi _K)\left[\widehat{𝒯}_{M^{}M}^\alpha (\stackrel{}{K})/\sqrt{N_\alpha }\right]`$ (16)
where the value of the bracketed term on the right hand side is interpolated from the coupled channels solution.

In practice this interpolation is carried out as a function of the deviation of $`K`$ from the threshold center of mass wave number. For non-$`s`$-wave breakup, the amplitudes are constrained to vanish at the breakup threshold $`K_{thr}`$, i.e.

$`\widehat{𝒯}_{M^{}M}^\mathrm{}0(\stackrel{}{K}_{thr})=0,{\displaystyle \frac{\mathrm{}^2K_{thr}^2}{2\mu _{pt}}}={\displaystyle \frac{\mathrm{}^2K_0^2}{2\mu _{pt}}}_0.`$ (17)
We note that in Eqs. (15) and (16) the functional dependence of the $`T`$-matrix on the angles of $`\stackrel{}{k}`$, the phase shifts $`\overline{\delta }_\alpha (k)`$, and the azimuthal angle $`\varphi _K`$ are all treated exactly. The grid of $`\theta _K`$ values can also be very fine without computing cost. The most important requirement is therefore that the number of bin states used to describe each $`[0k_{max}]`$ $`J_p^{}`$ excitation must be sufficient to allow an accurate interpolation of the amplitudes in the value of $`\mathrm{\Delta }K=|KK_{thr}|`$, or alternatively in $`k`$.

### D Three-body observables

The three-body amplitudes of Eq. (15) are used to compute the triple differential cross sections for breakup. If the energy of the core particle is measured then

$$\frac{d^3\sigma }{d\mathrm{\Omega }_cd\mathrm{\Omega }_vdE_c}=\frac{2\pi \mu _{pt}}{\mathrm{}^2K_0}\frac{1}{(2J_p+1)}\underset{\mu \sigma M}{}\left|T_{\mu \sigma :M}(\stackrel{}{k},\stackrel{}{K})\right|^2\rho (E_c,\mathrm{\Omega }_c,\mathrm{\Omega }_v).$$
(18)
With our $`T`$-matrix normalisations, and non-relativistic kinematics, the necessary three-body phase space factor $`\rho (E_c,\mathrm{\Omega }_c,\mathrm{\Omega }_v)`$, the density of states per unit core particle energy interval for detection at solid angles $`\mathrm{\Omega }_v`$ and $`\mathrm{\Omega }_c`$, is

$$\rho (E_c,\mathrm{\Omega }_c,\mathrm{\Omega }_v)=\frac{m_cm_v\mathrm{}k_c\mathrm{}k_v}{(2\pi \mathrm{})^6}\left[\frac{m_t}{m_v+m_t+m_v(\stackrel{}{k}_c\stackrel{}{K}_{tot})\stackrel{}{k}_v/k_v^2}\right].$$
(19)
Here $`\mathrm{}\stackrel{}{k}_c`$ and $`\mathrm{}\stackrel{}{k}_v`$ are the core and valence particle momenta in the final state and $`\mathrm{}\stackrel{}{K}_{tot}`$ the total momentum of the system, all evaluated in the frame, c.m. or laboratory, of interest. The association with the appropriate $`T`$-matrix elements in Eq. (18) is made through

$$\stackrel{}{K}=\stackrel{}{k}_c+\stackrel{}{k}_v\frac{m_p}{m_p+m_t}\stackrel{}{K}_{tot},\stackrel{}{k}=\frac{m_c}{m_p}\stackrel{}{k}_v\frac{m_v}{m_p}\stackrel{}{k}_c.$$
(20)
As the data under discussion here are inclusive with respect to the valence particle (proton) angles, the calculated triple differential cross sections must be integrated numerically over $`\mathrm{\Omega }_v`$. The presented observables are also integrated and averaged over the solid angles $`\mathrm{\Delta }\mathrm{\Omega }_c`$ subtended by the core particle detectors, with the stated detector efficiency profiles $`\epsilon (\mathrm{\Omega }_c)`$ , i.e.

$$\frac{d^2\sigma }{d\mathrm{\Omega }_cdE_c}=\frac{1}{\mathrm{\Delta }\mathrm{\Omega }_c}_{\mathrm{\Delta }\mathrm{\Omega }_c}𝑑\mathrm{\Omega }_c\left\{\epsilon (\mathrm{\Omega }_c)𝑑\mathrm{\Omega }_v\frac{d^3\sigma }{d\mathrm{\Omega }_cd\mathrm{\Omega }_vdE_c}\right\}.$$
(21)
It is most convenient to choose the $`xz`$ plane to be that defined by the beam and the core particle detector.

## III Application to sub-Coulomb breakup

The method detailed above is applied to the breakup of 8B on 58Ni at energy $`E_{lab}=25.8`$ MeV, for which new data are available . A first experiment was performed in 1996 at the Nuclear Structure Laboratory of the University of Notre Dame (ND) , one motivation being to clarify the importance of the E2 contribution to the Coulomb dissociation process, an issue which is still not completely resolved . In that first experiment, the measured 7Be fragments were detected at only one laboratory angle ($`40^{}`$), assumed to be free from the influence of strong interaction contributions. However, as a result of theoretical predictions  of a strong nuclear peak beyond $`40^{}`$, and claims also of Coulomb-nuclear interference at around $`40^{}`$, a more complete experiment was recently carried out using the now upgraded ND facility. Measurements were obtained of an angular distribution of the 7Be fragments  and also of their energy distributions  for the range of measured laboratory angles. Although the removed proton is not observed, since the heavy fragment energies are identified, the presented 7Be fragment distributions are known to contain no contribution from proton transfer reactions to bound states of 59Cu. There may nevertheless be contributions from knockout or stripping processes in which the proton excites the target and is absorbed. Such contributions are not calculated in this work. Proton transfer reactions to near-threshold (unbound) states of 59Cu, if present, could also contribute. We comment briefly below on the latter.

### A The CDCC model space

Model space parameters common to all the CDCC calculations are as follows. Partial waves up to $`L_{max}=1000`$ and radii $`R`$ up to $`R_{coup}=500`$ fm were used for the computation of the projectile-target relative motion wave functions. The continuum bins were calculated using radii $`rr_{bin}=60`$ fm. The 7Be intrinsic spin was neglected, the core being assumed to behave as a spectator. Thus we set $`I=0`$. The proton spin, $`s=1/2`$, was included and hence $`J_p^{}=j`$.

In the final calculations presented all $`J_p^{}`$ states consistent with relative orbital angular momenta $`\mathrm{}3`$, i.e. $`J_p^{}`$ up to $`f_{7/2}`$, were included. We show that the effects of the $`g`$-wave continuum are small. The bin state discretisation was carried out up to maximum relative energy $`_{max}=10`$ MeV for each state. The number of bins in the $`s_{1/2}`$-continuum was 32. For each of the other $`J_p^{}`$, 16 bins were used. These had equally spaced $`k_i`$ from $`k=0`$ to $`k_{max}`$. In the case of the DWBA calculations shown the model space is the same, however, the bin states are coupled to the ground state in first-order only. Calculations using potential multipoles $`q4`$ in constructing the coupling potentials will be shown but the final calculations require $`q3`$.

For the 7Be-58Ni system, the interaction of Moroz et al.  was used, as in the earlier analysis of Ref. . The proton-7Be binding potential was taken from Esbensen and Bertsch (EB) . The model of Kim et al.  is also considered. The potential used to construct the bin states was the same (real) potential as was used to bind the 8B ground state, assumed a pure $`p_{3/2}`$ proton single-particle state. The proton-58Ni potential is first taken from the global parameterization of Becchetti and Greenlees (BG) , but is also discussed below.

### B Results of calculations

It is important to note from the outset that the total breakup cross section angular distribution of the c.m. of the excited projectile, the sum of the two-body inelastic differential cross sections of Eq. (9), is incoherent in the different bin components. This is not the case for the three-body amplitude of Eq. (15) and the triple differential cross sections, Eq. (18). The practical convergence of the calculation, i.e. the dependence of the observables on the assumed model space, is therefore much more subtle in this case.

The three-body calculations are found to require a more extended set of bins, excitation energies, and potential multipoles. Whereas the use of energy bins up to only 3 MeV of relative energy, and multipoles $`q2`$, e.g. in Ref. , gives stable (converged) c.m. differential cross sections, in the sense of Eq. (9), this is not the case for the calculations of the triple differential cross sections and the energy and angle integrated distributions. We need the enlarged coupled channels model space, as detailed above, with bins extending beyond $`_{max}=8`$ MeV to obtain a converged result for these three-body observables. Furthermore, even when the extended range of continuum energies is included, the bin discretisation may itself not be fine enough so that the basis of bin states is sufficiently complete. We have therefore verified the stability of our results, with regard to the bin size, by doubling the number of bins and confirming that the same results are produced.

#### 1 Angular distributions

The convergence of the three-body calculations with $`_{max}`$ is clearly illustrated in Fig. 1. Here we show the 7Be laboratory differential cross section angular distributions from calculations that include continuum bins up to $`_{max}=3,4,6,8`$, and 10 MeV. The calculations for this convergence test use multipoles $`q2`$ and $`\mathrm{}3`$. The calculations use the BG proton-target potential and the EB proton-7Be potential. For the larger $`_{max}`$ the bins have been constructed so as not to alter their low energy discretisation. The calculation of the three-body cross sections thus provides a different interpretation of the reaction mechanism, and evidence for significantly higher energy excitations than would be deduced from the earlier calculations and their comparison with the 8B c.m. cross section. We will show that these high relative motion excitations are reflected in the calculated breakup energy distributions for 7Be and the proton.

Figs. 2 and 3 present the calculated 7Be laboratory differential cross section angular distribution, integrated over energy and proton angles and averaged over the core detector solid angles, and compare this with the data . The 7Be detectors were circular, subtending a solid angle $`\mathrm{\Delta }\mathrm{\Omega }_c`$ comprising a circle of radius $`6^{}`$ about the nominal laboratory angle $`\theta _{lab}`$. They have a stated Gaussian efficiency profile $`\epsilon (\theta )`$ with full width half maximum of 10.9 . Here $`\theta `$ is measured from the nominal $`\theta _{lab}`$ setting.

The convergence of the calculations with multipole order, and also with the included continuum partial waves, is shown in Fig. 2. Here the long dashed curve is the result shown in Fig. 1, converged with respect to excitation energy, with $`q2`$ and $`\mathrm{}3`$. The solid curve includes also the effects of the $`q=3`$ multipole couplings for $`\mathrm{}3`$. The dash-dot curve is a calculation where $`q=4`$ multipole couplings and the $`\mathrm{}=4`$ breakup partial wave are included. The additional effects are small and the remaining calculations therefore use the truncated model space with $`q3`$ and $`\mathrm{}3`$.

The solid curve in Fig. 3 uses the BG proton-target potential and the EB proton-7Be potential. In Ref.  it was shown that different 7Be-58Ni potential models give essentially the same shape for the 8B c.m. angular distribution, while the cross section normalisation depends on the size of the 8B g.s. wave function. The long-dashed curve in Fig. 3 shows the results of using the proton-7Be interaction of Kim et al. . Consistent with earlier work, the cross section is enhanced due to the larger predicted 8B root mean squared (r.m.s.) radius in this model.

The Becchetti-Greenlees  proton-58Ni potential, used above and previously, has surface imaginary strength and geometry parameters $`W`$=12 MeV, $`r_W`$=1.32 fm, and $`a_W`$=0.534 fm when computed at 3 MeV proton energy. Experience tells us  that the BG parameters give reasonable fits to data only down to approximately 10 MeV. An alternative global parameterization, tailored for use below 20 MeV, has a similar imaginary strength but somewhat smaller radius and diffuseness parameters $`r_W`$=1.25 fm, and $`a_W`$=0.47 fm  and leads to very similar results. There are, however, also potential parameters fitted to elastic scattering data at 5.45 MeV . This analysis uses a Gaussian surface term and obtains a much reduced absorptive strength, $`W`$=3.5 MeV, $`r_W`$=1.23 fm, and $`a_W`$=1.2 fm. We will refer to this as the VG potential. There is therefore some uncertainty in this potential input. The dot-dashed curve in Fig. 3 shows the calculated 7Be angular distribution from the VG potential. The cross section is changed only slightly at smaller angles. At the larger angles the calculated cross section is enhanced and is consistent with the experimental angular distribution data.

Our calculations show that the 8B structure (size) and proton-target potential uncertainties affect the calculations in characteristically different ways. The former produces an overall scaling while the latter produces, principally, a large angle enhancement. The data, currently, do not allow these effects to be discriminated further. In the final event, the overall agreement between the calculations and the data in Fig. 3 is qualitatively similar to the comparisons made in Ref. . There the calculated 8B c.m. cross sections  are compared with an approximate transformation of the measured 7Be data of Fig. 3 to the c.m. frame. Such approximate comparisons, however, are not necessary.

We observe that the results of our calculations are qualitatively quite different to those presented in Ref. , where an isotropic approximation was assumed in calculating the 7Be fragment laboratory cross sections. Those calculations show a radical change of shape of the angular distribution at forward angles which is not present in the calculations of Figs. 1, 2, and 3 in which the angular dependences are treated exactly.

#### 2 Energy distributions

In Fig. 4 we show the calculated breakup energy distributions of the 7Be fragments, together with the data from Ref. , for four measured laboratory configurations. For the smallest angle, $`20^{}`$, the calculations and the data are the average of the distributions at $`\theta _{lab}=19^{}`$ and $`\theta _{lab}=21^{}`$. For the largest angles, $`50/60^{}`$, the curves and data are similarly the average of the distributions obtained at $`\theta _{lab}=50^{}`$ and $`\theta _{lab}=60^{}`$. The measured cross sections are zero outside of the range of the data points shown. The solid curves use the BG proton distortion and the EB proton-7Be potential. The general features of the data, their magnitude, centroids, and widths, are well described by the calculations. The long-dashed curves are the results using the Kim proton-7Be potential. They show an enhanced cross section, discussed earlier, but a very similar shape. The dot-dashed curves are calculated using the VG proton distortion and the EB proton-7Be potential. The small arrows on the energy axis in Fig. 4 (and Fig. 5) indicate 7/8 of the 8B energy for elastic scattering at each laboratory angle. An overall reduction in the mean energy of the heavy fragments within the breakup reaction is evident.

Further insight is gained by looking at the results of DWBA calculations, and also calculations in which a subset of the continuum couplings are switched off, shown in Figs. 5(a–d). The long-dashed lines show the DWBA calculations. The dot-dashed lines are the results of CDCC coupled channels calculations but in which all continuum-continuum (CC) couplings between bin states are removed. The solid lines are the full calculations, as were shown in Fig. 4. We see that the calculations in the absence of CC couplings, both DWBA and truncated coupled channels, show energy distributions which are strongly asymmetric and have an enhanced high energy peak. This is very similar to what is observed in the 7Be fragment parallel momentum distributions from 8B breakup observed at higher-energy . As in that case, we show in Fig. 6 that this asymmetry has its origin in the interference between the E1 transitions to even breakup partial waves, and the E2 transitions to odd breakup partial waves. These $`E\lambda `$ amplitudes, which individually give approximately symmetric energy distributions, interfere to give strongly asymmetric responses. The very nearby kinematic cutoff in our case distorts the symmetry somewhat. The E2/E1 amplitude ratio in this lower energy case is also greater and so the asymmetry is enhanced compared to higher energies.

In the full CDCC calculations these asymmetries are essentially removed as a result of the higher-order couplings. This higher-order coupling induced suppression of the E1/E2 interference asymmetry was also a feature of the (higher energy) dynamical calculations in Ref. . The suppression is more complete at the lower energy discussed here. Figure 7 shows the analogue of Fig. 6(a), the calculated cross sections to odd and even breakup partial waves, from the full CDCC calculations using EB and BG potentials. Evident is the interference, both within and between the odd and even partial wave excitations. We note that the analogue of the $`E2`$ cross section, the $`p+f`$ wave contribution, is not itself suppressed, and is in fact large. The interference between the two contributions in Figs. 7 and 6(a) is however very different in the two cases.

Also evident in these two figures is that the odd breakup partial waves contribution in the CDCC calculation is significantly narrower than that calculated using DWBA. This narrowing is already manifest in $`s+p`$ wave two-step ($`q2`$ Coulomb) calculations and arises there from interference between the first-order E2 and second-order E1 amplitudes for populating the $`p`$ wave continuum. The importance of these particular interfering paths was also noted in Ref. , there in connection with a reduction in the calculated 8B decay energy spectrum at higher energy, when going beyond first-order Coulomb excitation theory. The calculated energy distributions reveal even more clearly than those for the angular distribution the importance of a full treatment of the dynamical couplings within the continuum.

#### 3 Additional calculations and comments

Since the proton separation energy from the 59Cu(g.s.) is $`S_p=3.42`$ MeV, proton transfer to the 59Cu(g.s.) would produce 7Be fragments with $`26`$ MeV of kinetic energy in the c.m. frame, and so such events are not part of the energy distributions measured. Those transfers that might contribute to the energy spectra of Fig. 4 would therefore be to excited (resonant) proton levels in 59Cu at around 9 MeV of excitation. If the proton-58Ni interaction supported one or more potential resonances, then the CDCC reaction mechanism would include their dynamical effects since breakup, by projectile excitation and by transfer to unbound states, are not distinguishable mechanisms in the three-body reaction model used. Clearly, however, the ability of the proton-58Ni interaction to support such resonance strength, and its absorptive content, are closely related questions. As was noted earlier, Fig. 3, use of the VG proton-target potential calculates an enhanced large angle cross section. Clarifying this sensitivity, and the possible role of such final-state resonances, requires further study and fine tuning of the proton-target potential. A full discussion of this topic is beyond the scope and motivation of the present article.

With this sensitivity to the proton-target potential in mind, however, in Fig. 8, we show the calculated proton laboratory angular distributions from the EB and Kim 8B wave functions, and the BG and VG proton distorting potentials. We note that the magnitude, but not the shape, of the proton cross section angular distribution shows a significant sensitivity to the assumed absorption in the proton-target system. Precise data could therefore verify and constrain this element of the calculations.

The shape of the calculated proton energy distribution, like that for the 7Be fragments, shows little sensitivity to the absorptive content of the proton distortion or to the choice of 8B binding potential. The calculations in Fig. 9 use the EB (solid) and Kim (long dashed) models for the proton-7Be interaction and the BG proton-target interaction. The dot-dashed curve uses the EB proton-7Be interaction and the VG proton-target interaction. The calculated proton energy distributions, integrated over all 7Be fragment angles, peak for $`E_p3.8`$ MeV and have a width $`\mathrm{\Gamma }4`$ MeV. The tail of the energy distribution is seen to extend to high energy, reflecting the high relative energy excitations of the 8B discussed earlier in connection with the convergence of the CDCC calculations. Figure 10 shows the energy distributions predicted when the 7Be fragments emerge at laboratory angles of 20, 30 and 40. In this case the arrows on the different curves indicate 1/8 of the 8B energy for elastic scattering at each laboratory angle. The calculations show an increased average energy (acceleration) of the removed protons from the dynamics of the breakup process.

## IV Summary and conclusions

In this paper we have calculated the most exclusive three-body breakup observables of a two-body projectile using the coupled channels CDCC methodology. The formalism is applied to investigate the angular and energy distributions of the 7Be fragments resulting from the sub-Coulomb breakup of 8B on a 58Ni target, the subject of recent experiments. We show that the convergence of the CDCC calculations of these observables is more subtle than that for the cross section of the c.m. motion of the 8B and requires a significantly more extended space of 8B excitation energies. The required excitation energy range is clarified.

Our calculations show that the 8B structure and the absorptive content of the proton-target potentials affect the calculated 7Be fragment angular distributions differently, the former producing an overall scaling, and the latter a large angle enhancement. Reducing the strength of the imaginary part of the proton potential in line with a phenomenological study , provides agreement with the larger angle data. The full CDCC calculations are shown to provide a good description of the measured 7Be fragment energy distributions. The widths and positions of these distributions are found to be rather insensitive to the details of the potentials used within the calculations. The presence of coupling between the continuum states is shown to be crucial to understand both the magnitudes of these energy distributions and their measured energy centroids. The absorptive content of the proton-target potentials affect the magnitudes of the calculated proton angular and energy distributions significantly, although their shapes are little affected. The calculated proton (7Be) fragment energy distribution reveals an overall increased (reduced) average energy of the fragment from the dynamics of the breakup process.

The application of these techniques to calculate the parallel momentum distribution of the heavy breakup fragments following the nuclear dissociation of the two-body system 11Be will be reported elsewhere . Further applications to systems with significant Coulomb dissociation strength, such as for 8B breakup at energies of 40 MeV/nucleon and greater, are also in progress.

###### Acknowledgements.

We thank Dr Valdir Guimarães and Prof. Jim Kolata for providing the data presented in tabular form and for detailed discussions of the experimental arrangement. The financial support of the United Kingdom Engineering and Physical Sciences Research Council (EPSRC) in the form of Grants Nos. GR/J95867 and GR/M82141 and Portuguese support from Grant FCT PRAXIS/PCEX/P/FIS/4/96 are gratefully acknowledged.


text:
# Complex Scaling of the Faddeev Equations

## 1 Introduction

The complex scaling method  invented in early 70-s remains one of the most effective approaches to calculation of resonances in few-body systems. This method is applicable to an $`N`$-body problem in the case where inter-particle interaction potentials are analytic functions of coordinates. The complex scaling gives a possibility to rotate the continuous spectrum of an $`N`$-body Hamiltonian in such a way that certain sectors of unphysical sheets neighboring the physical one turn into a part of the physical sheet for the resulting non-selfadjoint operator. Resonances appear to be complex eigenvalues of this operator  while the binding energies stay fixed during the scaling transformations. Therefore, when searching for the resonances within the complex scaling approach one may apply the methods which are usually employed to locate the binding energies. Some reviews of the literature on the complex scaling and its many applications can be found, in particular, in . Here we only mention that there is a rigorous mathematical proof  that for a rather wide class of interaction potentials the resonances given by the complex scaling method coincide with the “true scattering resonances”, i. e. the poles of the analytically continued scattering matrix in the unphysical sheets.

Along with the complex scaling, various different methods are also used for calculations of the resonances. Among the methods developed to calculate directly the scattering-matrix resonances we, first, mention the approach based on the momentum space Faddeev integral equations  (see, e. g., Ref.  and references cited therein). In this approach one numerically solves the equations continued into an unphysical sheet and, thus, the three-body resonances arise as the poles of the continued T-matrix. Another approach to calculation of the scattering-matrix resonances is based on the explicit representations  for the analytically continued T- and S-matrices in terms of the physical sheet. From these representations one infers that the three-body resonances can be found as zeros of certain truncations of the scattering matrix only taken in the physical sheet. Such an approach can be employed even in the coordinate space .

To the best of our knowledge there are no published works applying the complex scaling to the Faddeev equations. Therefore, we consider the present investigation as a first attempt undertaken in this direction. However, the purpose of our work is rather two-fold. On the one hand, we make the complex scaling of the Faddeev differential equations. On the other hand we compare the complex scaling method with the scattering-matrix approach suggested in . We do this making use of both the approaches to examine resonances in a model system of three bosons having the nucleon masses and in the three-nucleon ($`nnp`$) system itself.

## 2 Formalism

First, we recall that, after the scaling transformation, the three-body Schrödinger operator reads as follows

$$H(\vartheta )=\mathrm{e}^{2\vartheta }\mathrm{\Delta }_X+\underset{\alpha =1}{\overset{3}{}}V_\alpha (\mathrm{e}^\vartheta x_\alpha )$$
(1)
where $`\vartheta =\mathrm{i}\theta `$ is the scaling parameter with $`\theta `$. By $`\mathrm{\Delta }_X`$ we understand the six–dimensional Laplacian in $`X(𝐱_\alpha ,𝐲_\alpha )`$ where $`𝐱_\alpha ,𝐲_\alpha `$ are the standard Jacobi variables, $`\alpha =1,2,3`$. Notation $`V_\alpha `$ is used for the two-body potentials which are assumed to depend on $`x_\alpha =|𝐱_\alpha |`$ but not on $`\widehat{𝐱}_\alpha =𝐱_\alpha /x_\alpha `$.

The corresponding scaled Faddeev equations which we solve read

$`[\mathrm{e}^{2\vartheta }\mathrm{\Delta }_X+v_\alpha (\mathrm{e}^\vartheta x_\alpha )z]\mathrm{\Phi }^{(\alpha )}(z;X)+V_\alpha (\mathrm{e}^\vartheta x_\alpha ){\displaystyle \underset{\beta \alpha }{}}\mathrm{\Phi }^{(\beta )}(z;X)=f_\alpha (X),`$
$`\alpha =1,2,3.`$ (2)
Here $`f=(f_1,f_2,f_3)`$ is an arbitrary three-component vector with components $`f_\alpha `$ belonging to the three-body Hilbert space $`L_2(^6)`$.

The partial-wave version of the equations (2) for a system of three identical bosons at the zero total angular momentum $`L=0`$ reads

$`\mathrm{e}^{2\mathrm{i}\theta }H_0^{(l)}\mathrm{\Phi }_l(z;x,y)z\mathrm{\Phi }_l(z;x,y)+V(\mathrm{e}^{\mathrm{i}\theta }x)\mathrm{\Psi }_l(z;x,y)=f^{(l)}(x,y)`$ (3)
where $`x>0`$, $`y>0`$ and $`H_0^{(l)}`$ denotes the partial-wave kinetic energy operator,

$$H_0^{(l)}=\frac{^2}{x^2}\frac{^2}{y^2}+l(l+1)\left(\frac{1}{x^2}+\frac{1}{y^2}\right),l=0,2,4,\mathrm{},$$

while $`\mathrm{\Psi }_l`$ stands for the partial-wave component of the total wave function,

$$\mathrm{\Psi }_l(z;x,y)=\mathrm{\Phi }_l(z;x,y)+\underset{l^{}}{}_1^{+1}𝑑\eta h_{ll^{}}(x,y,\eta )\mathrm{\Phi }_l^{}(z;x^{},y^{}).$$
(4)
Here, $`x^{}=\sqrt{\frac{1}{4}x^2+\frac{3}{4}y^2\frac{\sqrt{3}}{2}xy\eta }`$ and $`y^{}=\sqrt{\frac{3}{4}x^2+\frac{1}{4}y^2+\frac{\sqrt{3}}{2}xy\eta }`$. Explicit expression for the geometric function $`h_{ll^{}}(x,y,\eta )`$ can be found, e. g., in .

The partial-wave equations (3) are supplied with the boundary conditions

$$\mathrm{\Phi }_l(z;x,y)|_{x=0}=0\text{and}\mathrm{\Phi }_l(z;x,y)|_{y=0}=0.$$
(5)
For compactly supported inhomogeneous terms $`f^{(l)}(x,y)`$ the partial-wave Faddeev component $`\mathrm{\Phi }_l(z;x,y)`$ also satisfies the asymptotic condition

$$\begin{array}{ccc}\mathrm{\Phi }_l(z;x,y)\hfill & =\hfill & \delta _{l0}\psi _d(\mathrm{e}^{\mathrm{i}\theta }x)\mathrm{exp}(\mathrm{i}\sqrt{zϵ_d}\mathrm{e}^{\mathrm{i}\theta }y)\left[\mathrm{a}_0(z)+o\left(y^{1/2}\right)\right]\hfill \\ & +\hfill & \frac{\mathrm{exp}(\mathrm{i}\sqrt{z}\mathrm{e}^{\mathrm{i}\theta }\rho )}{\sqrt{\rho }}\left[A_l(z;y/x)+o\left(\rho ^{1/2}\right)\right],\hfill \end{array}$$
(6)
For simplicity it is assumed in this formula that the two-boson subsystem has only one bound state with the energy $`ϵ_d`$, and $`\psi _d(x)`$ represents its wave function. The values of $`\mathrm{a}_0`$ and $`A_l(y/x)`$ are the main asymptotical coefficients effectively describing the contributions to $`\mathrm{\Phi }_l`$ from the elastic $`(2+12+1)`$ and breakup $`(2+11+1+1)`$ channels, respectively. Hereafter, by $`\sqrt{\zeta }`$, $`\zeta `$, we understand the main (arithmetic) branch of the function $`\zeta ^{1/2}`$.

In the scaling method a resonance is looked for as the energy $`z`$ which produces a pole to the quadratic form

$$Q(\theta ,z)=\left[H_F(\theta )z\right]^1f,f$$

where $`H_F(\theta )`$ is the non-selfadjoint operator resulting from the complex-scaling transformation of the Faddeev operator. The latter operator is just the operator constituted by the l. h. s. parts of Eqs. (2). The resonance energies should not, of course, depend on the scaling parameter $`\theta `$ and on the choice of the terms $`f^{(l)}(x,y)`$.

In the scattering-matrix approach we solve the same partial-wave Faddeev equations (3) with the same boundary conditions (5) and (6) but for $`\theta =0`$ and

$$f^{(l)}(x,y)=V(x)_1^{+1}𝑑\eta h_{l0}(x,y,\eta )\psi _d(x^{})\mathrm{sin}(\sqrt{zϵ_d}y^{}).$$

The resonances are looked for as zeroes of the truncated scattering-matrix (see  for details) $`\mathrm{s}_0(z)=1+2\mathrm{i}\mathrm{a}_0(z)`$, where the $`(1+11+1)`$ elastic scattering amplitude $`\mathrm{a}_0(z)`$ for complex energies $`z`$ in the physical sheet is extracted from the asymptotics (6).

For numerical solution of the boundary-value problem (3 – 6) we employ its finite-difference approximation in the hyperradius-hyperangle coordinates. A detail description of the finite-difference algorithm used can be found in Ref. .

## 3 Results

In the table we present our results obtained for a complex-scaling resonance in the model three-body system which consists of identical bosons having the nucleon mass. To describe interaction between them we employ a Gauss-type potential of Ref.

$$V(x)=V_0\mathrm{exp}[\mu _0x^2]+V_b\mathrm{exp}[\mu _b(xx_b)^2]$$

with $`V_0=55`$ MeV, $`\mu _0=0.2`$ fm-2, $`x_b=5`$ fm, $`\mu _b=0.01`$ fm-2 and $`V_b=1.5`$. The figures in the table correspond to the roots of the inverse function $`[Q(\theta ,z)]^1`$ for $`L=0`$ and $`l=0`$ only taken into account. In the present calculation we have taken up to 400 knots in both hyperradius and hyperangle variables while for the cut-off hyperradius we take 40 fm. One observes from the table that the position of the resonance depends very weakly on the scaling parameter $`\theta `$ which confirms a good numerical quality of our results. We compare the resonance values of the table to the resonance value $`z_{\mathrm{res}}=5.9520.403\mathrm{i}`$ MeV obtained for the same three-boson system with exactly the same potentials but in the completely different scattering-matrix approach of Ref. . We see that, indeed, both the complex scaling and the scattering matrix approaches give the same result.

| $`\theta `$ | $`z_{\mathrm{res}}`$ (MeV) | $`\theta `$ | $`z_{\mathrm{res}}`$ (MeV) |
| --- | --- | --- | --- |
| 0.25 | $`5.95250.4034\mathrm{i}`$ | 0.50 | $`5.95260.4032\mathrm{i}`$ |
| 0.30 | $`5.95260.4033\mathrm{i}`$ | 0.60 | $`5.95260.4033\mathrm{i}`$ |
| 0.40 | $`5.95260.4032\mathrm{i}`$ | 0.70 | $`5.95260.4034\mathrm{i}`$ |

We also watched the trajectory of the above resonance when the barrier amplitude $`V_b`$ varied (see. Fig. 1). While the complex scaling method was applicable it gave practically the same positions for the resonance. For the barrier amplitudes $`V_b`$ smaller than 1.0 only the scattering-matrix approach allows to locate the resonance (which finally, for $`V_b<0.85`$, turns into a virtual level).

As to the $`nnp`$ system in the $`S`$ – state where we employed the MT I–III  potential model, both the methods applied give no resonances on the two-body unphysical sheet (see ). Moreover, we have found no resonances in the part of the three-body sheet accessible via the complex scaling method. Thus, at least in the framework of the MT I–III model we can not confirm the experimental result of Ref.  in which the point $`1.5\pm 0.3\mathrm{i}(0.3\pm 0.15)`$ MeV was interpreted as a resonance corresponding to an exited state of the triton 3H.

The triton virtual state can be only calculated within the scattering-matrix method but not in the scaling approach. Our present improved scattering-matrix result for the triton virtual state is $`2.690`$ MeV (i. e. the virtual level lies 0.47 MeV below the two-body threshold). This result has been obtained with the MT I-III potential on a grid having 1000 knots in both hyperradial and hyperradial variables and with the value of cut-off hyperradius equal to 120 fm. Notice that some values for the virtual-state energy obtained by different authors can be found in  and all of these values are about 0.5 MeV below the two-body threshold.


text:
# SCALING VIOLATIONS IN THE Q² LOGARITHMIC DERIVATIVE OF F₂

## 1 Introduction

The physics of small $`Q^2`$ and small $`x`$ is associated with the search for the scale where gluon saturation , implied by s-channel unitarity , sets in. Gluon saturation signals the transition from the perturbative (hard) to the non perturbative (soft) QCD regime. We expect this transition to be preceded by signatures of screening corrections (SC) which should be experimentally visible even though the relevant scattering amplitude has not yet reached the unitarity (black disk) limit. Although, the general theoretical framework for saturation is reasonably well understood , the specific kinematic domain where we expect to see evidence of the role it plays, is as yet not determined. Consequently, this QCD component depends on the relevant DIS experimental data and the associated phenomenology. While the global analysis of the proton structure function $`F_2`$ (or $`\sigma _{tot}^{\gamma ^{}p}(W,Q^2)`$) data shows no significant deviations from DGLAP , there are dedicated HERA measurements of the $`Q^2`$ logarithmic derivative of the proton structure function $`F_2`$ , which add to our knowledge as they provide detailed information on the rate of change of the logarithmic $`Q^2`$ dependence of $`F_2`$, and, hence, they magnify possible deviations from the expected DGLAP behaviour at small x and small $`Q^2`$.

The interest in ($`\frac{F_2}{\mathrm{ln}Q^2}`$) stems from the observation that:

1. In leading order DGLAP evolution the $`Q^2`$ logarithmic slope of $`F_2`$ at low $`x`$, is directly proportional to the gluon structure function , since

$$\frac{F_2(x,Q^2)}{\mathrm{ln}Q^2}=\frac{2\alpha _S}{9\pi }xG^{DGLAP}(x,Q^2),$$
(1)
where $`xG^{DGLAP}(x,Q^2)`$ denotes the DGLAP gluon distribution of the proton.
2. $`\frac{F_2}{\mathrm{ln}Q^2}`$ is also related to the dipole cross section

$$\frac{F_2(x,Q^2)}{\mathrm{ln}Q^2}=\frac{Q^2}{3\pi ^3}\sigma _{dipole}(r_{}^2,x)=\frac{Q^2}{3\pi ^2}𝑑b^2Ima_{el}(r_{}^2,x;b);$$
(2)
where $`a_{el}(r_{}^2,x)`$ is the elastic scattering amplitude at fixed impact parameter of a dipole of size $`r_{}=2/Q`$ with energy $`x`$. For the amplitude $`a_{el}(r_{}^2,x)`$ we have a non-linear equation in the region of high density QCD (see Refs.) as well as the unitarity constraint $`a_{el}(r_{}^2,x;b)\mathrm{\hspace{0.17em}1}`$. Using Eq. (2) one can extract the dipole - target amplitude directly from the experimental data.

In this letter we present a detailed study of $`\frac{F_2(x,Q^2)}{\mathrm{ln}Q^2}`$. Our goal is to utilize this quantity, which is relatively well measured, to examine whether the system of partons has reached a QCD regime of gluon saturation. The complexity and difficulty of this investigation is illustrated by Fig. 1 which shows two different limits of $`\frac{F_2}{\mathrm{ln}Q^2}`$: the DGLAP limit at low $`x`$ for different parton distribution functions (pdf)  solutions of the evolution equations; and the unitarity limit at very long distances (low $`Q^2`$). From Fig. 1 one can see that (i) both the DGLAP equation and gluon saturation which leads to the unitarity bound displayed in Fig. 1 have the same qualitative behaviour decreasing at low $`Q^2`$; and (ii) the saturation scale ($`Q_s(x)`$) at which the DGLAP solution crosses the unitarity bound, differs for different parameterizations. $`Q_s(x)45\mathrm{GeV}^2`$ for GRV98 and $`Q_s(x)12\mathrm{GeV}^2`$ for CTEQ5 and MRS99 parameterizations. Such a low value for the saturation scale indicates that we shall have difficulty attempting to distinguish between models with gluon saturation effects, and models that include soft and hard interactions, since the matching of the components occurs at similar values of $`Q^2`$ ($`1GeV^2`$).

In our study we will incorporate the following initial observations and constraints:

1. Recent measurements of $`F_2`$ , extending to very low $`Q^2`$, find that it has a monotonic behaviour down to the photo production limit. These data suggest a smooth transition from the perturbative (DGLAP, high $`Q^2`$) to the non perturbative (Regge or classical gluon field, exceedingly small $`Q^2`$) regime.
2. Preliminary HERA data  on $`\frac{F_2}{\mathrm{ln}Q^2}`$ at small $`Q^2`$ and small $`x`$ do not show any turnover when plotted at fixed $`x`$ or fixed $`Q^2`$.
3. An approach which provides a way of successfully describing $`\frac{F_2}{\mathrm{ln}Q^2}`$ should also be able to reproduce data in other channels which are sensitive to $`xG(x,Q^2)`$ such as $`F_2^c`$, $`F_L`$ and the photo and DIS production of heavy vector mesons.

The main goal of this letter is to investigate and compare:

1. The ability of conventional DGLAP, using the pdfs , to reproduce the available data on $`\frac{F_2}{\mathrm{ln}Q^2}`$.
2. The role of screening corrections (SC). In this approach , DGLAP is modified by SC due to the $`q\overline{q}`$ percolation through the target, as well as by modifications due to screening in the gluon distribution $`xG(x,Q^2)`$.
3. Within the framework of pQCD no allowance is made for explicit soft contributions, both in the non-screened and screened formulations of DGLAP. The latter provides a smooth interpolation between the high and low $`Q^2`$ regimes. In the present investigation we have checked the validity of our methods at low values of $`Q^2`$, i.e. just above the DGLAP $`Q_0^2`$ threshold (defined for each pdf input). For the case with SC, we have also considered an extrapolation for $`Q^2Q_0^2`$. Our aim is to determine the values of $`Q^2`$ that one can successfully describe without the help a soft component. To this end we compare with the Donnachie and Landshoff (DL) two Pomeron model , which explicitly uses the sum of a hard and a soft Pomeron.
4. We conclude with a critical discussion and caution regarding attempts that have been made to determine the gluon saturation scale using data with highly correlated variables such as in the Caldwell and fixed W plots . We also comment on the general problems associated with gluon saturation and screening corrections relating to other models which include SC .

## 2 The small $`Q^2`$ and small $`x`$ behavior of $`\frac{F_2}{\mathrm{ln}Q^2}`$ in various models

As is well known, a global DGLAP analysis of the data with the recent pdfs  is in good agreement with the experimental data. A study  comparing the screened and non screened DGLAP calculations of $`F_2(x,Q^2)`$ showed only a small difference due to SC even in the small $`Q^2`$ and $`x`$ region attained by present HERA measurements. A significant deviation of the data from Eq. (1), where $`xG^{DGLAP}`$ is obtained from the global $`F_2`$ analysis, may serve as an experimental signature indicating the growing importance of unitarity corrections. This was recently suggested by Caldwell , showing a rather complicated plot of $`\frac{F_2}{\mathrm{ln}Q^2}`$ in which each measured point had different correlated values of $`Q^2`$ and $`x`$. The Caldwell plot suggested a dramatic turnover of $`\frac{F_2}{\mathrm{ln}Q^2}`$ corresponding to $`Q^2`$ of about $`3GeV^2`$ and $`x<510^3`$, in contrast to the behavior expected from GRV94  at sufficiently small $`Q^2`$ and $`x`$. The problem with this presentation is that, as suggestive as it may seem, it does not discriminate between different dynamical interpretations . It is, actually compatible with an overall data generator  as well as the latest pdfs which were adjusted to account for this observation.

Better discrimination is achieved if we carefully study the small $`Q^2`$ and $`x`$ dependences of $`\frac{F_2}{\mathrm{ln}Q^2}`$ at either fixed $`Q^2`$ or fixed $`x`$ values, so as to be free of the kinematic correlation between $`Q^2`$ and $`x`$ that is endemic in the data displayed in the Caldwell plot. Such preliminary HERA data have recently became available . Note that the method that has been used to specify the logarithmic derivative of the structure function $`F_2`$ in  is by fitting the $`F_2`$ data at fixed $`x`$, to the expression

$$F_2(x,Q^2)=a(x)+b(x)\mathrm{ln}Q^2+c(x)[\mathrm{ln}Q^2]^2$$
(3)
and then make the identification

$$F_2/\mathrm{ln}Q^2=b(x)+2c(x)\mathrm{ln}Q^2$$
(4)
H1  compared the results using Eq. (4) with the results obtained by taking the derivative of their QCD parameterization (H1 QCD fit). They find that for $`x\mathrm{\hspace{0.17em}0.0002}`$, the L.H. and R.H. sides of Eq. (4) match very well, however, for $`x<2\times 10^4`$ there is a discrepancy. We note that the experimental data of  are correlated, such that the lower values of $`x`$ are associated with the lower values of $`Q^2`$.

At the initial stage of our investigation we examined whether DGLAP (using the pfds  in the NLO approximation) reproduce the data. The results are shown in Fig. 2. We note that:

1. GRV98 overestimates the low $`Q^2\mathrm{\hspace{0.17em}5}GeV^2`$ data. For higher $`Q^2`$, the reproduction of the data is good.
2. MRS99 underestimates the data at $`Q^2`$ = 1.9 and 40 $`GeV^2`$. It is above the data in the $`3Q^210GeV^2`$ range and does well for $`Q^2`$ = 12 and 30 $`GeV^2`$.
3. CTEQ5 (we have used CTEQ5HQ as it reproduces the energy dependence of $`J/\psi `$ photo production ) provides a good reproduction of the experimental data over the entire range $`1.9Q^240GeV^2`$. A few words of explanation are called for regarding the numerical NLO calculation of $`F_2`$ from the CTEQ pdf (we discuss only CTEQ in this context since both GRV98 and MRS99 supply a code for calculating the structure function, from which the parton densities have been parameterized). When calculating $`F_2`$ in NLO one should take care when inserting the threshold for charm production, so that $`F_2`$ remains a smooth function of $`Q^2`$. The effect of the threshold is small if one looks at the structure function, but it can be rather large when examining it’s $`Q^2`$ derivative. We found that by following the consistent treatment of charm in , we obtain a smooth behaviour of both $`F_2`$ and $`\frac{F_2}{\mathrm{ln}Q^2}`$.

.

In order to study the role of SC in our calculations, we follow the eikonal SC formalism presented in Ref.111For theoretical and phenomenological status of this approach see Ref. ., where screening is calculated in DLA in both the quark sector, to account for the percolation of a $`q\overline{q}`$ through the target, and the gluon sector, to account for the screening of $`xG(x,Q^2)`$. The factorizable result that we obtain is

$$\frac{F_2^{SC}(x,Q^2)}{\mathrm{ln}Q^2}=D_q(x,Q^2)D_g(x,Q^2)\frac{F_2^{DGLAP}(x,Q^2)}{\mathrm{ln}Q^2}.$$
(5)
SC in the quark sector are given by

$$D_q(x,Q^2)\frac{F_2^{DGLAP}(x,Q^2)}{\mathrm{ln}Q^2}=\frac{Q^2}{3\pi ^2}𝑑b^2\left(1e^{\kappa _q(x,Q^2;b^2)}\right),$$
(6)
$$\kappa _q=\frac{2\pi \alpha _S}{3Q^2}xG^{DGLAP}(x,Q^2)\mathrm{\Gamma }(b^2).$$
(7)
The calculation is simplified if we assume a Gaussian parameterization for the two gluon non perturbative form factor,

$$\mathrm{\Gamma }(b^2)=\frac{1}{R^2}e^{b^2/R^2}.$$
(8)
SC in the gluon sector are given by

$$xG^{SC}(x,Q^2)=D_g(x,Q^2)xG^{DGLAP}(x,Q^2),$$
(9)
where

$$xG^{SC}(x,Q^2)=\frac{2}{\pi ^2}_x^1\frac{dx^{}}{x^{}}_0^{Q^2}𝑑Q^{\mathrm{\hspace{0.17em}2}}𝑑b^2\left(1e^{\kappa _g(x^{},Q^{\mathrm{\hspace{0.17em}2}};b^2)}\right),$$
(10)
and $`\kappa _g=\frac{9}{4}\kappa _q`$.

An obvious difficulty in the above calculation of $`xG^{SC}`$ stems from the fact that the $`Q^{\mathrm{\hspace{0.17em}2}}`$ integration spans not only the short (pQCD), but also the long (npQCD) distances. To overcome this difficulty we assume that

$$xG\left(x,Q^2<\mu ^2\right)=\frac{Q^2}{\mu ^2}xG(x,\mu ^2).$$
(11)
Our choice of the above extrapolation is motivated by the gauge invariance requirement that $`xGQ^2`$ when $`Q^2\mathrm{\hspace{0.17em}0}`$, and is supported by recent ZEUS measurements of $`F_2`$ at exceedingly small $`Q^2`$ .

The SC calculation described above can be applied to any given input pdf where the only adjusted parameters are $`R^2`$ and $`\mu ^2`$. For the radius $`R^2`$, see Eq. (8), we use the value 8.5 $`GeV^2`$ which is determined directly from the forward slope of $`J/\mathrm{\Psi }`$ photo production . $`\mu ^2`$ is conveniently fixed at $`Q_0^2`$, the lowest $`Q^2`$ value accessible for the input pdf we use. Once we have chosen our pdf, our SC calculation is essentially parameter free. We have checked that our output results are not sensitive to the fine tuning of these fixed parameters.

Before presenting our results we would like to recall two important features of SC:

1. They can only dampen (reduce) the results obtained with the unscreened input dpfs.
2. They are negligible for large enough $`Q^2`$ and/or $`x`$. This can be easily deduced from the functional dependence of $`\kappa _q`$ and $`\kappa _g`$ on $`Q^2`$ and $`x`$ .

Throughout the SC calculation we have used as input the $`\overline{MS}`$ version of NLO GRV98. The results obtained using GRV98 DIS are very similar.

Our fixed $`Q^2`$, and fixed $`x`$ results are displayed in Fig. 3 and Fig. 4, respectively. As can be seen in the limit of small $`Q^2`$ and $`x`$ there is a significant difference between the screened and non screened values of $`\frac{F_2}{\mathrm{ln}Q^2}`$. As expected the SC results are smaller and softer than the non screened input. All in all, our reproduction of the experimental data for $`Q^2>\mathrm{\hspace{0.17em}2}GeV^2`$ is very good. The ZEUS fixed $`Q^2`$ = 1.9 $`GeV^2`$ data are somewhat softer than our predictions. This feature is conspicuous when we compare with the very small $`Q^2`$ = 0.3 and 0.75 $`GeV^2`$ data, where we have made use of the extrapolation given in Eq.(11). In Fig. 4 we see that, using this extrapolation, screened GRV’98 is also successful in describing the data at values of $`Q^2\mathrm{\hspace{0.17em}1}GeV^2`$, where in the DL parameterization (see Fig. 6) the soft component dominates.

To further examine the role of the soft (non perturbative) component in $`\frac{F_2}{\mathrm{ln}Q^2}`$, we compare the above DGLAP results with the DL two Pomeron model . The DL parameterization is based on the Regge formalism, and consists of the coherent sum of contributions from a hard and soft Pomeron, a normal Reggeon and an additional contribution from the charmed sector which is proportional to the hard Pomeron. Each of these fixed j-poles are multiplied by a fitted $`Q^2`$ form factor. The hard Pomeron with a fixed trajectory has an intercept of 1.43 at t=0, while the soft Pomeron has an intercept of 1.08 at t=0. As can be seen in Fig. 3, the DL model reproduces the data at all values of $`Q^2`$.

We conclude that one can obtain a good description of $`\frac{F_2}{\mathrm{ln}Q^2}`$ data with $`Q^2\mathrm{\hspace{0.17em}1.9}`$ $`GeV^2`$ using any of the following options:

1. Conventional DGLAP evolution using CTEQ5 pdfs as input.
2. DGLAP evolution using GRV98 pdfs as input, provided SC are included.
3. The two Pomeron (DL) model combining hard and soft components.

For smaller $`Q^2(<Q_0^2)`$, a soft contribution is called for, for which there are several parameterizations available, but not yet a precise theory.

Recently, Golec-Biernat and Wuesthoff  have suggested studying the $`Q^2`$ and $`x`$ behaviour of $`\frac{F_2}{\mathrm{ln}Q^2}`$ at fixed W as a method to determine the gluon saturation scale from the anticipated turnover in these plots. The new ZEUS low $`Q^2`$ presentations  of these plots show, indeed, the anticipated turnover structure in these figures, suggesting that gluon saturation is attained at $`Q^2`$ a few Ge$`V^2`$. We advocate using the invariant variables $`x`$ and $`Q^2`$ when describing $`F_2`$. Introducing other variables such as W to study the structure function gives rise to spurious effects which are predominantly kinematic. In the specific procedure suggested by Golec-Biernat and Wuesthoff, the combination of the kinematic relation between $`x,Q^2`$ and W with the very general behaviour of $`xG(x,Q^2)`$ is sufficient to produce a turnover. Its exact location depends on the details of the numerical input. Consequently, the suggested fixed W plots do not convey any new dynamical information on the saturation problem even if such information is contained in the analyzed data. We note that a fixed $`W`$ plot is natural and informative for Regge type models. It appears that any parameterization for $`F_2`$ (or $`xG`$) which has factorizable $`Q^2`$ and $`x`$ dependences, and provides a reasonable description of the data such as the Buchmueller- Haidt (BH) model , gives rise to the fixed W turnover effect.

To illustrate this point we consider the fixed W behaviour in three parameterizations which have completely different dynamics:

1. DGLAP with CTEQ5 pdf. As noted, the results of this parameterization are close to the experimental data for $`\frac{F_2}{\mathrm{ln}Q^2}`$ at fixed $`Q^2`$ and at fixed $`x`$ values over the kinematic range of interest. This parameterization does not include a specific soft component.
2. The GLMN model , which is a pure pQCD dipole model with SC. As such, the model relates indirectly to gluon saturation, even though it is constructed so as to include unitarity corrections below actual saturation. This parameterization has no explicit soft component.
3. The DL two Pomeron parameterization .

All three models, as well as the BH parameterization, follow the experimental behaviour of $`\frac{F_2}{\mathrm{ln}Q^2}`$ at fixed W rather well, including the observed turnover. We demonstrate this in Fig. 5 by comparing the ZEUS data  with the results of our SC model. In Fig. 6 we plot the DL model predictions versus the data at fixed W. We also display in this figure the contribution of the hard DL Pomeron.

We conclude with three general comments:

1. Gluon saturation is not unique in producing a turnover at fixed $`W`$.
2. The gluon saturation scale may be estimated theoretically from the contours produced at the boundary of $`\kappa _g=\mathrm{\hspace{0.17em}1}`$, as discussed in our papers . This theoretical analysis suggests that gluon saturation occurs at $`Q^2`$ 1 $`GeV^2`$. Note that $`\kappa _g`$ depends on $`xG(x,Q^2)`$, which is determined from a DGLAP pdf fit to the global $`F_2`$ data.
3. We note that in the BH and DL parameterizations the soft contributions are concentrated at a rather large typical scale $`\mathrm{\hspace{0.17em}\hspace{0.17em}2}GeV^2`$. This observation supports the idea that the soft Pomeron originates from rather short distances .

## 3 Discussion and Conclusions

The main experimental observation regarding $`\frac{F_2}{\mathrm{ln}Q^2}`$ is that at fixed $`x`$ it is a monotonic decreasing function of $`Q^2`$, and for fixed $`Q^2`$ it increases as $`x`$ becomes smaller. This is the observed pattern even for the lowest measured value of $`Q^2`$ = 0.3 $`GeV^2`$ . The main objective of our investigation was to check whether the preliminary H1  and ZEUS  data contains a decisive experimental signature indicating signs of gluon saturation. This effect is associated with the transition from the relatively well understood pQCD hard domain to the more complex, and less understood soft domain. We have shown in this paper that the gluon saturation is capable of describing all HERA data on $`\frac{F_2}{\mathrm{ln}Q^2}`$ and, therefore, the data do not contradict the idea that HERA has reached a new QCD regime: high parton density QCD.

However, the final conclusion of our analysis is that the present data can be described by a number of models, which are very different in their basic construction. Specifically, we show that all available data with $`Q^2>Q_0^2`$, of $`\frac{F_2}{\mathrm{ln}Q^2}`$ can be described by DGLAP using the CTEQ5 pdf as input, as well as GRV98 pdf which have been corrected for screening. The DL model which is Regge motivated and includes a significant soft component, also provides a good description of the data.

We have shown that the turnover seen in $`\frac{F_2}{\mathrm{ln}Q^2}`$ at fixed $`W`$, as a function of either fixed $`Q^2`$ or $`x`$, is not necessarily connected with saturation effects and cannot be used to discriminate between models. There has also been a successful reproduction of the Caldwell plot by Kaidalov et al.  within the framework of the CKMT model , and by Forshaw et al. using a colour dipole model formalism. We are of the opinion that fitting the Caldwell as well as the fixed W plots is a necessary first step, but a comparison with all the detailed measurements is essential to test the adequacy of a given model.

In Fig. 7 we attempt to clarify the similarities and differences between the non screened CTEQ5 and the screened GRV98 DGLAP calculations whose results are very similar, (see Figs. 2 and 3). To this end we display the calculated distributions of $`\frac{F_2}{\mathrm{ln}Q^2}`$ and $`\alpha _sxG(x,Q^2)`$ at fixed values $`x=\mathrm{\hspace{0.17em}10}^3,10^4,10^5`$ and $`10^6`$. As can be seen the results using CTEQ5 (non screened) and GRV98 (screened) in the limit of very small $`x`$ are very similar. It is therefore suggestive, that CTEQ5 contains significant screening effects, that are absent in the boundary conditions used in GRV98. Note that CTEQ5 is not defined for $`x<\mathrm{\hspace{0.17em}10}^5`$. Our estimates of SC for the CTEQ5 parameterization (based on Eq. (5) - Eq. (10)) is about 10%, which should be considered as the uncertainty in this parameterization. We also show the DL predictions for $`\frac{F_2}{\mathrm{ln}Q^2}`$, which differ from the hard partonic DGLAP at small enough $`x`$.

The ZEUS data at exceedingly low $`Q^2`$  are of particular interest when investigating the transition from the DGLAP dominated region to the non perturbative (low $`Q^2`$) region. This transition was expected to be observed experimentally since $`F_2Q^2`$ as $`Q^2\mathrm{\hspace{0.17em}0}`$ due to gauge invariance requirements, and for large $`Q^2`$, $`F_2(Q^2)^\gamma `$ from DGLAP evolution, where $`\gamma `$ is the anomalous dimension. This transition has been seen by ZEUS  in their measurements of $`F_2(x,Q^2)`$ at small values of $`x`$ and $`Q^2`$, where the transition appears to be at $`Q^2\mathrm{\hspace{0.17em}1}GeV^2`$, and is compatible with our early theoretical estimates . This provides an indirect indication that this is the scale for the onset of gluon saturation, but we still lack a decisive signature for this effect.

Based on our present investigation we conclude that the behaviour of $`\frac{F_2}{\mathrm{ln}Q^2}`$ as measured in the kinematic region presently accessible at HERA, does not display unambiguous signs of saturation. This is compatible with the information displayed in Fig. 1 where we see that the anticipated gluon saturation scale is close to $`Q_0^2`$, and to the scale at which the soft contribution becomes significant.

We are extending our search for possible effects of gluon saturation to the channel of photo and DIS production of $`J/\mathrm{\Psi }`$ , where the cross section should be very sensitive to this signal, as it is proportional to $`[xG(Q^2,x)]^2`$.

Acknowledgments:

UM wishes to thank UFRJ and FAPERJ (Brazil) for their support. EG and EL would like to acknowledge the support and hospitality of the Theory Group at DESY. EN would like to thank H. Abramowicz for a helpful discussion. This research was supported by in part by the Israel Academy of Science and Humanities and by BSF grant # 98000276.