T Muñoz-Britton

M Freer

I Ashwood

D Brown, W N Catford

N Curtis

S P Fox

B R Fulton

C W Harlin

M Laird

N I Ashwood

T A D Brown

W N Catford

A M Laird

P Mumby-Croft

A St. J Murphy

P Papka

D L Price

K Vaughan

D L Watson

D C Weisser

Search for the 2 + Excitation of the Hoyle State in 12 C Using the 12 C( 12 C,3α) 12 C Reaction

Keywords: numbers: 21.10.-k, 23.60.+e, 27.20.+n

published or not. The documents may come

Introduction

The structure of the nucleus 12 C is fascinating for a variety of reasons. For example, it lies at the present limit of the range of ab initio calculations such as the Quantum Monte Carlo approach [1,2], and provides an important test as it has a rather rich range of structural possibilities. The ground-state has an oblate-structure which is well represented in the shell-model. However, the next 0 + state at 7.65 MeV is not well described within this framework. Shell model calculations, for example those of Karataglidis et al. [3] reproduce, rather well, the energy of the first 2 + excitation, but in the region of the second state, the 0 + 2 (7.65 MeV), there is a void in the calculations; the energy of this state cannot be reproduced. This in itself points to the rather unusual structure. A similar conclusion is reached in no-core shell model calculations [4]. Indeed, in Ikeda's classification of the structure of light nuclei in 1968 it was assumed that this state should have a 3α cluster structure. Analysis of electron inelastic scattering data [5,6,7,8] indicates that the 7.65 MeV state has a radius which is ∼1.5 times that of the 12 C ground-state, i.e. a volume some 3-4 times larger. This larger volume reduces the overlap of the α-particles and may allow them to obtain their quasi-free characteristics in something approaching an α-particle gas, or perhaps a bosonic condensate [9]. It should be noted though that whilst FMD calculations predict a structure dominantly composed of three weakly bound α particles [10], they also suggest a 8 Be + α structure. This is in contradiction to a strict interpretation of a Bose-Einstein condensate and highlights the danger of taking the interpretation too literally.

The 7.65 MeV 0 + state has a rather central role in stellar nucleosynthesis. It acts as the gateway through which 12 C is formed via the triple-α process. Its existence was first predicted by Hoyle [11] in order to explain the observed abundance of 12 C and was later experimentally confirmed by Fowler and co-workers [12]. It has subsequently become called the Hoyle-state. In order to arrive at a detailed understanding of the structure associated with this state it is important to determine the nature of possible collective excitations. Determining of location of the associated 2 + will allow the moment of inertia, and hence a measure of the deformation, to be found. To date no firm evidence for such an excitation has unambiguously been found. However, Itoh et al. have reported evidence for a 2 + state close to 9.7 MeV with a width which is a few hundred keV [13]. These measurements were performed using the 12 C(α,α') inelastic scattering reaction, and multipole decomposition of the angular distributions was used to extract the possible 2 + strength. Similar measurements using the 12 C(p,p ) reaction also indicate a broad (600 keV) 2 + state close to 9.6 MeV [14]. Both FMD [10] and AMD [15] calculations also suggest the possible existence of a 2 + 2 excitation in the region of 9-12 MeV. Here we present measurements of the 12 C( 12 C, 12 C * ) reaction in which the angular correlation technique is used to search for a 2 + signature close to 9.7 MeV.

Experimental Details

The present measurement was performed using the 14UD tandem Van der Graaff Pelletron accelerator facility at the the Australian National University (ANU). A 101.5 MeV 12 C beam was used to bombard a 50 µgcm -2 12 C foil target; the beam current was typically 10 enA and the charge state of the beam was 6 + .

The reaction studied was 12 C( 12 C,3α) 12 C and the three α-particles were detected in coincidence in an array of four detector telescopes (two primary and two secondary). These provided an angular coverage in the range of θ * ≈ 5 • -30 • (in the laboratory frame); this suitably covered the range of excitation energies being investigated. The array allowed coverage of the complete azimuthal range. The primary detector telescopes consisted of three elements: a (50 × 50) mm 2 , 70 µm thick, double-sided Si strip detector (DSSD) placed in front of a (50 × 50) mm 2 , 500 µm thick, position-sensitive Si resistive strip detector (RS) which was in turn located in front of a 10 mm thick CsI scintillator. The DSSDs were divided into two sets of 16 independent 3 mm wide strips, the 16 strips on the front face having a horizontal orientation and the 16 on the back face vertical. The RS detectors were separated into 16 position sensitive strips; these provided position information with a (FWHM) resolution of ∼ 0.3 mm (in-plane) and ± 1.5 mm (out-ofplane). The CsI scintillators were used to stop any highly energetic particles that passed completely through both the Si detectors. These latter detectors provided a software veto, removing any high energy events in which there was not full energy deposition in the Si elements (due to the CsI scintillators' inability to resolve multiple hits into a telescope). The two secondary telescopes consisted of a (50 × 50) mm 2 , 70 µm thick, Si quad detector in front of a 500 µm thick DSSD (with the same strip arrangement as in the primary detectors), which were located in front of a 10 mm thick CsI scintillator. The Si quadrant detectors were separated into 4 (25 × 25) mm 2 segments (quadrants). The elements in both telescopes provided E -E particle identification information to distinguish α-particles as required. The detectors and amplifiers were calibrated using a combination of two techniques: a 3-line α-source ( 239 Pu, 241 Am and 244 Cm), and 45 MeV 12 C ions elastically scattered from a 197 Au target. The energy resolution was ∼ 150 keV (FWHM) for the silicon detectors and ∼ 1.5 % (FWHM) for the CsI scintillators.

Analysis and Results

The 12 C( 12 C,3α) 12 C reaction was used to investigate the possible existence of a 2 + excitation of the Hoyle state in the 9.7 MeV region. This was performed by considering the reverse of the 'triple-alpha' process by which 12 C is formed in stellar nucleosynthesis. Techniques were then employed which allowed the suppression of a known 3 -resonance which also exists at this energy. Analysis was primarily undertaken in decays through the 8 Be g.s. + α reaction channel as this is a pre-requisite for the angular correlation technique used.

In order to improve statistics, events were selected in which only 3α-particles were detected. i.e. the 12 C recoil particle was undetected. This was performed using E -E particle identification spectra. Following this selection, the total energy (E tot ) was then calculated from

E tot = E 1 + E 2 + E 3 + E rec = E beam + Q.
Here E 1-3 correspond to the energies of the three detected α-particles, E beam the beam energy, Q the Q-value of the break-up reaction (Q=-7.272 MeV) and E rec the recoil energy of the undetected 12 C nucleus. E rec was determined by applying momentum conservation between the beam and the final-state particles, whilst assuming the mass of the recoil particle, m rec = 12 amu.

The E tot spectrum is shown in Figure 1 and shows a series of peaks. The highest energy peak corresponds to a Q-value of ∼ -7.5 MeV which is consistent with the Qvalue of α-break-up of 12 C (-7.272 MeV). The remaining peaks correspond to excited states of the 12 C recoil nucleus. The peak at ∼ 90 MeV is consistent with the recoil 12 C being emitted in its first excited state at 4.4 MeV, the second member of a rotational band built on the ground state configuration. The peak at ∼ 84 MeV corresponds to the recoil nucleus being emitted in the 9.64 MeV 3 -state. By gating on the highest- energy peak, only events in which the recoil nucleus was formed in its (0 + ) ground state were selected; this is a requirement of the angular correlation method used later. The background events that extend to E tot values greater than the beam energy are believed to most likely arise from pile-up events and from events in which the detected particles are produced from different beam-target interactions.

The break-up of 12 C into 3α-particles occurs sequentially in two steps:

where the lifetime of the 8 Be nucleus is sufficiently short so as not to affect the overall kinematics. The decay of the intermediate 8 Be nucleus was reconstructed and compared with the Q-value for the 8 Be break-up (Q Be =92 keV). The decay energy was calculated using

E decay = 2 i=1 E α i - P Be 2 2M Be , (1)
where E α i are the energies of two of the three α-particles, P Be is the momentum of the 8 Be nucleus, and M Be is the mass of the 8 Be nucleus. All possible permutations of 2α-pairs were checked and events filtered to include those for which an appropriate Q-value existed for only one pairing. A 12 C excitation energy spectrum was constructed by a similar method to that used to reconstruct the 8 Be break-up: The dashed lines correspond to the energy at which the 2 + excitation was found in Ref. [13,14].

E ex = 3 i=1 E α i - P C 2 2M C + E th , (2)
where the terms are analogous to those in equation 1 with the addition of E th , which accounts for the missing break-up energy (E th = 7.272 MeV). Figure 2 shows the excitation energy spectrum for the break-up for 12 C for decays via 8 Be gs ; a number of known energy levels are labeled, along with the proposed location of the 2 + resonance associated with the Hoyle state identified in the inelastic scattering measurements [13,14]. The widths of the peaks is dominated by the experimental resolution which is 90 keV for the 7.65 MeV peak and 175 keV (FWMH) at 9.64 MeV.

In an attempt to 'look beneath' the 3 -resonance, which coincides with the proposed location of the 2 + state, angular correlation analysis has been used. The correlations depend on the angles which describe the two different stages of the reaction process (see Figure 3). The angle θ * describes the centre-of-mass emission angle of the excited 12 C nucleus, measured with respect to the beam axis. The second angle, ψ, corresponds to the emission angle of the first α-particle in the centre-of-mass frame of the 12 C nucleus, again measured with respect to the beam axis. As all of the initial and final state particles are spin zero, then for a scattering angle θ * =0 the 12 C excited state is constrained to be populated in the m=0 magnetic substate. The subsequent α-particle emission thus follows that of a squared Legendre-Polynomial of order J (the spin of the state in 12 C); |P J (cos(ψ))| 2 . Away from zero degrees other m-substates may contribute which gives rise to a phase shift in the Legendre polynomial; |P J (cos(ψ + ∆ψ))| 2 , where ∆ψ = ∆θ * l i -J J .

Here l i is the entrance channel grazing angular momentum [START_REF] Silviera | Proc. 14 th Winter Meeting on Nucl. Phys. Borneo[END_REF]. This results in a ridge structure in the θ * -ψ plane. Thus, the angular correlation spectrum is characterised by the spin of the state both in terms of the frequency of the oscillations as a function of ψ and also the gradient of the loci in the θ * -ψ plane [START_REF] Freer | [END_REF]. A more complete discussion of the correlation technique can be found in Ref. [START_REF] Freer | [END_REF]. This implies there will be regions in which the contribution from a 3 -state will be minimal and a 2 + state maximal. The simulated angular correlations for a 2 + and 3 -state are shown in the top part of Fig. 4. The Monte Carlo simulations performed simulated both the initial 12 C + 12 C interaction and also the resulting two stage sequential 12 C -→ 3α break-up.

The simulations considered energy lost by the beam in the target, as well as the beam spot size and beam divergence; for the decay fragments energy loss, energy straggling and angular straggling in the target and also energy loss and angular straggling in the detectors was considered. Finally, positional resolution and energy resolution of the detectors was included. These simulations reproduce the 12 C( 12 C,3α) reaction and the decay via the 8 Be ground state and reproduce the observed width of the 3 -peak (175 keV). It can be seen that at the location of the minimum in the 3 -correlation that the ridge in the 2 + distribution is close to being maximal. Thus, gating in the "valley" of the 3 -correlation pattern should enhance any 2 + contribution. The experimental angular correlations for the 3 -region are also shown in Fig. 4. The ridge pattern is similar to that observed in the simulations. To gain an impression of the factor by which the 3 -contribution is suppressed by selecting data in the correlation minimum the variation in intensity has been measured. Between the peak of most intense ridge and the valley the yield is reduced by a factor of 2.5, whilst for the less intense ridge at lower values of θ * the reduction is by a factor of 1.3. Thus, on average the 3 - strength is reduced by a factor of ∼ 2 (calculated from the simulations). A software window (window (ii), Figure 4) was placed over the region where the ratio of the 2 + to 3 - strength was expected to be largest. An excitation energy spectrum for the 12 C break-up was created for events located within this window. This spectrum was normalised by dividing the channel contents by the integrated number of counts in the 3 -peak. The resulting spectrum was compared with a similarly normalised spectrum obtained from a region containing most of the data at an energy range across the 3 -peak (window (i), Figure 4). The two spectra in Figure 5 appear to indicate a difference between the peak shapes with the data in the 3 -minimum showing a very slight enhancement at higher values of E ex . When the normalised peaks are subtracted from one another then there is a small remainder in the candidate region for a 2 + state.

In order to characterise this difference a fitting function was devised which was composed of two Gaussian functions (equation 4). This was then fitted to the experimental difference data. The fitting function used was

f (E x ) = exp -Ex-c 1 √ 2σ 1 2 + α exp -Ex-c 2 √ 2σ 2 2 √ 2π (σ 1 + ασ 2 ) - exp -Ex-c 1 √ 2σ 1 2 + β exp -Ex-c 2 √ 2σ 2 2 √ 2π (σ 1 + βσ 2 ) (4)
where c 1 and c 2 are the centroid positions of the peaks (c 1 = 9.64 MeV), σ 1 and σ 2 are controlled by the widths of the peaks (the width of the 3 -state, 150 keV (FWHM), was limited by the experimental resolution), and α and β are the 2 + contribution in each software window. Whilst c 1 , σ 1 were known, c 2 , σ 2 and β-α were optimised by the fitting routine. In order to simplify Equation 4it can be re-written as

f (E x ) = α -β (σ 1 + ασ 2 ) (σ 1 + βσ 2 )     - exp -Ex-c 1 √ 2σ 1 2 √ 2π σ 1 + exp -Ex-c 2 √ 2σ 2 2 √ 2π σ 2     . (5)
If ασ 2 and βσ 2 are both smaller than σ 1 then equation 5 may be written as

f (E x ) (α -β)     - exp -Ex-c 1 √ 2σ 1 2 √ 2πσ 1 + σ 2 2 σ 2 1 exp -Ex-c 2 √ 2σ 2 2 √ 2πσ 2     . ( 6
)
This shows that the fitting function should have a bipolar form, with a negative going component described by a Gaussian whose width is represented by that of the 3 -peak and amplitude by the difference α -β, whilst the positive component is associated with a second Gaussian with a width determined by σ 2 and amplitude (α -β)σ 2 2 /σ 2 1 . The above function was fitted to the difference data shown in Fig. 5 and produced values of c 2 = 9.67-9.71 MeV, σ 2 = 32-83 keV (FWHM = 75-195 keV) and β -α = -0.3 to +0.2 %. The quantity β -α evaluates the change in relative 2 + and 3 -contributions for the two selected regions. A value of zero indicates no additional contribution. It should be noted that the range of widths explored in the present analysis exceed the experimental resolution and so the Gaussian line-shape that has been used is only an approximation to the resonance profile. This approximation does not substantially affect the conclusions.

Simulations

In order to understand the above result more precisely Monte Carlo simulations were performed which replicated the 3 -resonance and background contributions as seen in the experimental data. These simulations were used to explore the sensitivity of the current technique and also any effect of the precision of the alignment of the detection system.

In the first stage of this analysis separate simulations for a 3 -and 2 + decay were performed using the correlation function described by Eqn. 5, where a value of 9 . Result of subtracting a spectrum from the dominant 3 -region from a spectrum of a region at which a J π =2 + resonance would be maximal.

l i = 22 was used. The results of these simulations are shown in Fig. 4. The events from the simulations were analysed using the same analysis procedures as used for the experimental data. In this way different contributions from a 2 + state could be added to the simulated 3 -decay data. In addition, a smooth background was added at the level that was observed below the 9.64 MeV peak in Fig. 5 and also a contribution from the 10.84 MeV, 1 -, state was included at the level observed in the experimental data. The simulations and background contributions were combined such that the number of counts were the same as in the data such that the same statistical sensitivity was achieved in the simulations. A number of simulations were performed for the generated 2 + resonance; varying centroid positions (9.70 MeV, 9.85 MeV and 10.0 MeV) and of varying widths (FWHM = 250 keV, 500 keV and 750 keV). Each 2 + resonance was then combined with the 3 -and background simulation in a ratio of 15:85 (i.e. a 15 % level). Similar calculations were also performed with 10, 5 and 2 % contributions from a 2 + state.

The resulting angular correlations were processed in the same way as described for the experimental data; i.e. gating in windows (i) and (ii), normalising the two corresponding excitation energy spectra and computing the difference. A number of initial 'guess' parameters were tried with the fitting routine in order to probe the sensitivity of the fit (equation 4), and by minimising the χ 2 value the optimal fit was obtained. At the 15% level the fitting routine was successfully able to distinguish the

Centroid (MeV) Width (keV) 15% 10% 5% 2% 250 × × × 9.70 500 × × × 750 × × × 250 × × 9.85 500 × × × 750 × × × 250 10.0 500 × 750 × × Table 1.
Levels at which a resonance can be detected for various centroids and widths.

2 + from the background under all conditions. The ratio was then altered (10, 5, 2 %), reducing the size of the 2 + resonance, until the fitting routine was no longer able to verify its existence (defined as when the fit parameters no longer matched the known characteristics of the peak entered into the simulation). This allowed a relationship to be found between the position/centroid and the strength of a resonance in order that it might be detected by the techniques employed in this analysis (table 1). Examples of some of the simulations and optimal fits are given in Fig. 6. Table 1 indicates that for a 2 + state well separated from the 9.64 MeV, 3 -, state (e.g. at 10.0 MeV and 250 keV width) it is possible to unambiguously pick out this contribution down to very low levels (2 %). Larger contributions (15 %) can be distinguished for all widths at all peak separations down to 60 keV. Unambiguous identification of the peak parameters becomes more difficult the smaller the separation from the 3 -state and the larger the width -as might be expected.

A number of Monte Carlo simulations were performed in order to investigate the sensitivity of the results with regard to the position of the detector telescopes; simulations were performed which introduced displacements in the position of both a detector-telescope (Fig. 7) and also the position of the beam with respect to the detectorarray (Fig. 8). Figure 7 shows the unaltered Monte Carlo simulation (i), the effect of one telescope being offset from the intended position (ii and iii), and the effect of one telescope being rotated by a small amount (iv and v). The effect of one telescope being offset is independent of the direction in which the offset occurs and causes the single 3 - peak to separate in two components. The effect is however more pronounced if the offset is in a direction parallel to the direction of the detector strips (ii) than in a direction perpendicular to the detector strips (iii). The effect of rotating one detector telescope by a small amount is dependent on the point about which the rotation occurs but is consistent in that it causes a broadening of the 3 -peak, and to a greater extent in the suppressed 3 -subset than in the entire data set. Simulations to investigate moving the position of the beam suggest a similar effect to those seen when a single telescope is offset, but also indicate the sensitivity is much greater. It also appears that the effects (i)

(ii) (iii) (iv) (v)
Figure 7. The effect that shifting a detector telescope can have on the excitation energy spectrum as calculated by using Monte Carlo simulations; (i) the unaltered 3 - peak, (ii) shifting one telescope 2 mm away from the beam in the direction of the detector strips, (iii) shifting one telescope 2 mm in a direction perpendicular to the detector strips, (iv) rotating a telescope 2 • clockwise about the beam axis, (v) rotating a telescope 2 • clockwise about its centre.

Discussion and Conclusions

The experimental difference data (Fig. 5) show a small amplitude deviation from zero (± 0.01). This is best described by an additional contribution to the 3 -peak at an energy of 9.67-9.71 MeV which corresponds to a very small change in intensity between the two windows used in the analysis (β-α = -0.3 to 0.2 %). This would correspond to a very small 2 + contribution. Note the negative value would indicate a reduction in the 2 + strength. The width found from the fit was in the range 75 to 195 keV, close to the width of the 3 -peak. These values are inconsistent with those extracted in the inelastic scattering data [13,14]. It is thus likely that the deviation observed in the experimental difference data (Fig. 5) is actually related to an experimental artifact, e.g. a slight misalignment of the detectors -as observed in the simulations. If a 600 keV wide 2 + state exists at 9.6 MeV, then it is populated in the present reaction at a strength which lies below the present experimental sensitivity. The summary of the Monte Carlo simulations shown in Table 1 indicates that for a state of the above parameters then this strength is less than ∼ 10 % of the strength of the 3 -state. If the width and centroid of the additional peak in equation 4 are set to 600 keV and 9.6 MeV respectively (rather than allowed to be free parameters) then the values of β-α are smaller than 0.1 %, i.e. representing virtually no change in the line-shape between the two windows used in the analysis of the angular correlations. In other words a broad 2 + state with a strength above 15 % is not found in the present data. It can be concluded from this analysis that if an additional 2 + resonance exists in the present data which has a width of less than 250 keV, then it must have a population strength of less than ∼ 10 % of the 3 -peak and an excitation energy of less than 9.85 MeV. If the width lies in the range of 500 to 750 keV then to be clearly identified it should have an excitation energy of ∼10 MeV and a population strength of 5 %. If it lies at lower energy much stronger strengths are required for unambiguous identification.

In conclusion, a search has been made for a 2 + state close to 9.6 MeV in 12 C using the 12 C( 12 C, 12 C[3α]) reaction. Angular correlation techniques have been used to enhance any 2 + contribution to the excitation energy spectrum, compared with a dominant 3 -contribution. These measurements fail to find any evidence for a new state at this energy. It is clear that in order to provide much more stringent limits increased experimental sensitivity is required. Improved excitation energy resolution would help define the excitation energy spectrum and improved energy and angular resolution would help better emphasise the minima in the angular correlation pattern.

Recent R-matrix analysis of the decay of 12 N and 12 B suggest that the 2 + excitation of the Hoyle state may lie at a higher energy (∼11 MeV) [15]. We have performed an analysis of this region using the angular correlation technique described here and can find no strong evidence for a 2 + state.

Inelastic scattering experiments using either the 12 C(p,p') 12 C or 12 C(α,α') 12 C reaction are both candidate methods for better exploring the 3 -region. A recent 12 C(p,p') 12 C experiment has found evidence suggesting a broad 2 + state lying slightly lower in energy than the 9.64 MeV 3 -state [14]. Future measurements are also planned which would combine the α-particle decay measurements with a magnetic spectrometer. The enhanced excitation energy resolution afforded by this technique would permit a more sensitive correlation analysis to be performed.

Figure 1 .

 1 Figure1. Total energy spectrum. The marked peak corresponds to events in which the recoil nucleus was formed in its ground state. Peaks at lower energy correspond to the 4.4 and 9.6 MeV excitations of the recoil nucleus.

Figure 2 .

 2 Figure 2. Excitation energy spectrum for the break-up of 12 C into 3α-particles.The dashed lines correspond to the energy at which the 2 + excitation was found in Ref.[13,14].

Figure 3 .

 3 Figure 3. The definition of the correlation angles θ * and ψ used in the present analysis (from Ref. [17]).

Figure 4 .

 4 Figure 4. Angular correlations. Top left Simulation of the angular correlation for the decay of a J π = 2 + state at 9.7 MeV. Top right Simulation of the angular correlation for the decay of a J π = 3 -state at 9.64 MeV. Bottom The angular correlation for the experimental data in the region of the 9.64 MeV 3 -resonance. The windows indicate positions of software gates used (see text for further details).

Figure 5 .

 5 Figure5. Result of subtracting a spectrum from the dominant 3 -region from a spectrum of a region at which a J π =2 + resonance would be maximal.

C -→ 8 Be + α ,8 Be -→ α + α ,

of shifting the beam in different directions are the same as if a telescope was offset in the same direction. These simulations show that for the most part the types of structures misalignments of the detectors would cause have a different type of structure to that which would occur if an additional state would exist below the 9.64 MeV, 3 -, peak.

Figure 8. The effect of the beam not passing through the central axis of the detector array: (i) Monte Carlo simulation of a 3 -peak with the beam aligned with the central axis, (ii) the same simulation as (i) but with the beam shifted 1 mm towards one detector, (iii) the same simulation as (i) but with the beam shifted 1 mm in a direction perpendicular to the detector-strips.


P Salamon
email: salamon@atomki.hu

A T Kruppa

Curvature Correction in the Strutinsky's Method

Keywords: numbers: 21.10.Dr, 21.10Ma, 21.10.Pc

come

Introduction

Masses, or, equivalently ground state energies of nuclei are very important basic quantities. For example, shell structure changes around drip lines or existence of new magic numbers beyond the heaviest known nuclei can be studied with the help of mass measurements. Important applications are the determination of possible decay modes of a system and the calculation of energy release in a given reaction. Mass values are important not only for nuclear physics but for astrophysics and particle physics as well.

Nowadays not only the precision of mass measurements has improved dramatically but it is possible to measure masses of short lived nuclei [1,2]. Very accurate experimental mass values are now available for almost every part of the periodic table and this is a great challenge for theoretical description.

Global mass calculations have been carried out using the so called macroscopicmicroscopic method [3,4]. In this approach the binding energy is the sum of two terms E(N, Z) = E mac (N, Z) + E mic (N, Z). The first term E mac (N, Z) is a smooth function of the neutron (N ) and proton number (Z). It can be determined using a macroscopic model. The simplest macroscopic expression is derived from the liquid drop model. This form can be generalized in several respects. E.g. a well known macroscopic term is the formula of the finite range droplet model [4].

Beside the smooth trend the binding energy formula contains an oscillatory part with pure quantum origin. The microscopic part E mic (N, Z) consists of two terms shell correction and pairing energy E mic (N, Z) = E shell (N, Z)+E pair (N, Z). The pairing term is usually calculated using BCS or Lipkin-Nogami scheme and the shell correction by the Strutinsky's method [5].

The shell correction method have been generalized in order to describe special excited states. High-spin states of a rotating nuclei can be described by the cranked shell-model. The modification of this model was developed [6,7] in order to include the shell effects. This spin dependent shell correction approach is used in many investigations [8][9][10][11]. An another important application of the shell correction method is the calculation of fission barriers [12][13][14][15][16]. This quantity is a key element for understanding of many phenomena e.g. stable super-heavy elements. The importance of the shell effects at an outer saddle-point barrier is recently emphasized [17]. The shell correction method have been used not only in nuclear physics but for description of mesoscopic systems. Metallic clusters have pronounced shell structures and these systems are an ideal ground for the application of the shell correction method [18][19][20].

The first step in the calculation of shell correction is to take a single particle Hamiltonian. The energy E of the ground state is simply sum of the energies of occupied single particle states. This energy can be decomposed into a sum of a smooth Ē and a fluctuating term Ẽ, E = Ē + Ẽ. It can be argued using the Strutinsky's energy theorem [21] that the quantities E shell (N, Z) and Ẽ are identical.

Recently new facts have emerged concerning the fluctuating part of the energy. The shell correction is a rapidly changing function of the particle number. Earlier, the average value of Ẽ was expected to be zero. Recently it has been shown that the average value of the fluctuating part of the energy is non zero and it contributes to the smooth part of the energy [22]. The extra smooth trend of the average value of Ẽ can be analytically calculated [22,23] in case of simple model potentials. It would be desirable to minimize this extra smooth trend of the shell correction. A study of the nuclear mass error distribution [24], using results of the Strutinsky type macroscopicmicroscopic calculation of [4], has been revealed a regular correlation and it has been suggested that refinement of the Strutinsky's method is needed. The focus of the present paper is how to calculate the smooth part of the energy Ē or equivalently the shell correction. The basic concepts are the quantum single particle level density (SPLD) and the smooth SPLD. The quantum SPLD is simply sum of Dirac-delta functions centered at the single particle energies. The smooth SPLD is free from singularities and it is a smooth function. It may be determined using semiclassical methods [START_REF] Brack | Semiclassical Physics[END_REF]. Shell correction calculations [3,[START_REF] Nazarewicz | [END_REF] were carried out using the smooth semi-classical SPLD based on Wigner-Kirkwood expansion. Smooth SPLD can be derived using only the single particle energies and a special smoothing method called Strutinsky's smoothing or standard averaging method [START_REF] Nazarewicz | [END_REF]. The Strutinsky type [4] and semi-classical [3] macroscopic-microscopic models differ on the way how the smooth SPLD is calculated. Beside the energy averaging particle number space averagings were developed [27][28][29] in order to calculate Ē and even the combination of averaging spaces was considered [30]. The Strutinsky's smoothing of quantum SPLD is carried out with a well defined procedure. First, a normalized, positive and even function is selected. The smoothing means that the convolution of the quantum SPLD and the kernel function is calculated [21,31]. This simple procedure however is insufficient. We have to demand selfconsistency of the smoothing process. If the smoothing is applied to a sufficiently smooth function (polynomial) then the result has to be identical with the original function. This property is referred to as curvature corrected smoothing.

In the most standard form of the Strutinsky's method the curvature correction is taken into account by multiplying the starting kernel function with an appropriately chosen polynomial [21]. Brack and Pauli have introduced [31] an alternative method to achieve the effect of the curvature correction. In this paper we introduce two new curvature correction methods.

The organization of the paper is the following. In section 2 we describe the Strutinsky's shell correction method and than in section 3 we give precise meaning of the curvature correction. Two standard methods and our new methods for curvature correction are also presented here. Finally at section 4 we compare the methods using harmonic oscillator and Saxon-Woods potentials. Summary is given at section 5.

Strutinsky's shell correction method

In the Strutinsky's method first an appropriate single particle Hamiltonian is taken. The eigenvalues of bound states and the corresponding degeneracies are denoted by i and d i , respectively. For simplicity the { i } series is arranged in increasing order. An important notion of the method is the single particle level density (SPLD). The so called exact or quantum mechanical SPLD is of the form

g( ) = N b i=1 d i δ( -i ), (1)
where δ(i ) is the Dirac-delta function and the number of bound states is N b . The central notion of the Strutinsky's method is the smooth SPLD. It is denoted by ḡ( ). The smoothed Fermi-level λ is obtained from the particle number equation

N τ = λ -∞ ḡ( )d , (2)
where N τ is the neutron or proton number of the considered nucleus since the shell correction is separately calculated for neutrons and protons. The smoothed total energy Ē can be calculated by an integral

Ē = λ -∞ ḡ( )d . (3)
The shell correction is defined by the equation

Ẽ = E sp -Ē, (4)
where the total single particle energy E sp is taken in the spirit of the independent particle shell model

E sp = Nτ i=1 i . (5)
The new energy series { i } is derived from the series { i } by repeating the eigenvalues i d i times.

Smoothing with convolution

Convolution is an ideal tool to smooth an oscillatory function y(x). We take a bell shaped kernel function K(x) and the smoothed ȳ(x) function is calculated by the convolution of K(x) and y(x)

y(x) = 1 γ ∞ -∞ K x -x γ y(x )dx , (6)
where the positive real number γ defines the smoothing width. The smooth SPLD can be obtained from ( 6) by identifying y with the quantum SPLD g. The result is

ḡ( ) = 1 γ N b i=1 d i K -i γ . (7)
The particle number equation takes the form

N τ = N b i=1 d i ν i -∞ K(x)dx (8)
and the smoothed total energy is

Ē = N b i=1 d i γ ν i -∞ xK(x)dx + i ν i -∞ K(x)dx , (9)
where

ν i = ( λ -i )/γ.

Curvature correction

We expect that a proper definition of smoothing has the following property: the smoothing does not change a smooth function. Here we identify smooth functions with polynomials [5,21]. We require that the mapping y(x) → y(x) is the identity map on the set of polynomials of order at most (2m + 1) i.e.

p m (x) = p m (x) = 1 γ ∞ -∞ K x -x γ p m (x )dx , (10)
where p m (x) is an arbitrary polynomial of order at most (2m + 1). We will say that a kernel K satisfying (10) has curvature correction of order m. If we want to emphasize this property the kernel will be denoted by K m instead of K. Other papers use the notation p = 2m for the curvature correction. A kernel with curvature correction of order zero will be called starting kernel and this function is assumed to be an even one. If we want to stress actual values of the smoothing range γ and the order of the curvature correction m we will use the notation ḡγ,m ( ) for the smooth SPLD.

Standard curvature correction methods

For completeness we review two standard curvature correction methods. In the most widely used method a kernel with curvature correction of order m is searched in the form [21]

K m (t) = P m (t)K 0 (t), (11)
where P m (t) is a polynomial of order 2m

P m (t) = m i=0 c i t 2i . (12)
The coefficients c i are determined by the condition (10). The following set of equations can be gained [21]

m i=0 M 2(i+j) (K 0 )c i = δ j,0 0 ≤ j ≤ m, (13)
where the moments of an arbitrary kernel K is defined by

M i (K) = ∞ -∞ t i K(t)dt, i = 0, 1, 2, . . . . ( 14
)
This approach is called polynomial curvature correction (PCC). Brack and Pauli [31] have taken the kernel with curvature correction of order m in the form

K m (t) = m i=0 a 2i d 2i dt 2i K 0 (t). ( 15
)
We will refer to this method as Brack and Pauli curvature correction (BPCC) method.

The condition (10) uniquely determines the coefficients a 2i . The value of a 0 is 1. The other expansion coefficients satisfy the following set of equations

m j=1 B i,j a 2j = -b 2i (i = 1, . . . m), (16)
where

B i,i = 1, B i,j = 0 (i < j), B i,j = b 2(i-j) (1 < i < j) and b 2i is related to the moments of K 0 b 2i = M 2i (K 0 ) (2i)! . ( 17
)

Curvature correction with many widths

The error term of the smoothing

δy(x) = y(x) -y(x) ( 1 8 )
can be cast into the form

δy(x) = ∞ n=1 (-1) n γ n 1 n! M n (K)y (n) (x), (19)
where the n-th derivative of y(x)is denoted by y (n) (x) and it is assumed M 0 (K) = 1.

If the moments of the kernel function K satisfy the following set of equations M 0 (K) = 1, M i (K) = 0, i = 1, 2, . . . , (2m + 1) then from (19) it is clear that such a kernel function has curvature correction of order m.

Let us assume that we have a kernel K m and suppose that it is an even function. We claim that a kernel with curvature correction of order (m + 1) can be written in the form

K m+1 (t) = K m (t) -c -2m-3 K m t c 1 -c -2m-2 , m ≥ 0. ( 20
)
Here c is an arbitrary real number but c = 0 and c = 1.

In order to prove our claim it is enough to show the validity of the equations

∞ -∞ K m+1 (t)dt = 1 (21) and ∞ -∞ t i K m+1 (t)dt = 0, i = 1, 2, . . . , (2m + 3). ( 22
)
The proofs are very simple. We have to substitute (20) into ( 21) and ( 22), use the linearity of the integration and in the second term of the left hand sides change the integration variable from t to u = t/c. If the same c value is used at each iteration step (starting from m = 0) then the following closed expression for the kernel K m (m ≥ 1) can be derived

K m (t) = A m m k=0 (-1) k c (m-k)(m-k+2) m k c 2 K 0 t c k , (23)
where the notation

m r q =      1 i fr = 0 (1-q m )(1-q m-1 )...(1-q m-r+1 ) (1-q)(1-q 2 )...(1-q r ) if 1 ≤ r ≤ m 0 i fr > m (24)
stands for the Gaussian coefficient (q-binomial coefficient) and

A m = m k=0 c 2k 1 c (m+2) 2 -4 m k=1 1 1 -c -(2k+2) . (25)
Kernels derived using (20) or ( 23) are called many width kernels and the smoothing method is named as many width curvature correction (MWCC). The name comes from the form of the smooth SPLD of the MWCC method

ḡγ,m ( ) = A m m k=0 (-1) k c (m-k)(m-k+2)+k m k c 2 ḡγc k ,0 ( ). ( 26
)
The smooth SPLD corresponding to a many width kernel with curvature correction of order m is a linear combination of smooth SPLD's derived from the starting kernel having smoothing width γc k (k = 0, . . . m). The individual terms ḡγc k ,0 ( ) behave differently, some terms may show the shell structure, other terms may over smooth but the proper linear combination can make a balance.

Derivative curvature correction

The recursion (20) is not valid at c = 1. If the limit c → 1 is taken using the L'Hospital rule then a new kernel with curvature correction of order m emerges

K m+1 (t) = 2m + 3 2m + 2 K m (t) + 1 2m + 2 t d dt K m (t). ( 27
)
Kernels obtained in this way are called derivative kernels and the method is called derivative curvature correction (DCC).

With the aid of the recursion (27) a closed expression can be derived for the derivative kernels (m ≥ 0)

K m (t) = m k=0 a(m + 1, k + 1) m!2 m t k d k dt k K 0 (t), (28)
where the coefficients a(n, m) obey the recursion

a(n, m) = (2n + m -2)a(n -1, m) + a(n -1, m -1) (29)
with initial conditions a(1, 1) = 1, a(n, 0) = 0 (n is arbitrary) and a(n, m) = 0 for n < m. Apart from the multiplicative factor t k the kernel ( 28) is similar to the kernel (15). In contrast to the BPCC method the linear combination coefficients in (28) are independent from the starting kernel. In the BPCC method the linear combination coefficients are calculated using (16).

The recursion (27) simplifies the calculation of auxiliary quantities of the shell correction method. For the calculation of the particle number we need the integral (see (8))

N m (x) = x -∞ K m (t)dt, m ≥ 0. ( 30
)
We can easily get the following recursion for the function N m (x)

N m+1 (x) = N m (x) + 1 2m + 2 xK m (x). ( 31
)
This recursion is valid provided lim t→-∞ tK m (t) = 0, m ≥ 0. For the calculation of the smoothed energy we need the integral (see ( 9))

E m (x) = x -∞ tK m (t)dt, m ≥ 0. ( 32
)
We can derive the following recursion

E m+1 (x) = 2m + 1 2m + 2 E m (x) + 1 2m + 2 x 2 K m (x). ( 33
)
Here we assumed lim t→-∞ t 2 K m (t) = 0, m ≥ 0.

Results

At the beginning of an application a starting kernel has to be chosen. The most widely used form corresponds to a Gaussian K 0 (t) = exp(-t 2 )/ √ π. The kernel K 0 (t) = 1/(2 Cosh 2 (t)) was considered in [31] using the BPCC method. We will use both functions and call them Gaussian and Cosh-type kernels respectively. An interesting observation have to be mentioned. In the case of Gaussian starting kernel the methods PCC, BPCC and DCC lead to the same kernel K m . The MWCC method however furnishes new kernel function. In the case of Cosh-type starting kernel each presented method is different.

The MWCC method in addition to the smoothing width γ and the order of the curvature correction m contains a new dimensionless technical parameter c. For a light nucleus (N τ =20) the shell correction as the function of the parameter c is displayed in Figure 1 using three different value of the smoothing width. According to this figure there is an optimal region for c (the shell correction is stable with respect to c). This region is roughly between 0.8 and 1 in the presented case if each γ value is considered.

In all calculations we have done we have found an optimal region for c and in the majority of cases this region was below one. Since the shell correction method already contains two technical parameters (γ and m) it would be unwise to increase them so we determined a procedure how to fix the value of c.

At each fixed values of γ and m the value of c is determined with the following procedure. An averaging region around the investigated particle number N τ is selected. Shell corrections of all nuclei belonging to the averaging window are calculated. The c parameter is used to get the average value of the shell corrections as close to zero as possible. This procedure is inspired by recent works [22,23]. The shell correction is a rapidly changing function of the particle number. The average value of Ẽ was expected to be zero but it turned out that the average of the fluctuating part of the energy is non zero and it contributes to the smooth part of the energy [22]. The value of c is determined in such a way to minimize this extra smooth trend of the shell correction.

In case of a light and a medium heavy nuclei we studied the dependence on the averaging window. The optimal value of c is determined the way described above. The order of the curvature correction was three. The upper bounds of the averaging windows correspond to shell closers and the lower bounds are a start of a new shell. The results are presented in Figure 2. The shell correction hardly depends on the averaging window (notice the scale on the vertical axis). We have found that in addition to the own shell of the considered particle number one shell below and one shell above is enough for the determination of the optimal value of c.

The shell correction method contains two technical parameters the smoothing width γ and the order of the curvature correction m. The result is supposed to be independent from the values of these parameters if they are selected in a reasonable way. In the case of harmonic oscillator or similar potentials (infinite number of bound states) it was demonstrated [31] that there is a large region of the smoothing width where the shell correction is practically constant (plateau condition). Therefore the new curvature correction methods have to be checked in the case of harmonic oscillator potential to see if the plateau is remained or not.

The parameter of the harmonic oscillator potential ( ω =6 MeV), the considered particle number (N τ =70) and the order of the curvature correction (m = 3) were chosen to be the same as in [31]. In this work the PCC (Gaussian type kernel) and the BPCC (Cosh-type kernel) methods were used. The shell correction as a function of the smoothing width is displayed in Figure 3 using Gaussian starting kernel. Since the shell spacings are ω the shell corrections are calculated only for γ > ω otherwise the effect of shells can be noticed in the smooth SPLD [37].

The curve signed by PCC m = 3 in Figure 3 agrees with the result of [31]. The large absolute value of the shell correction is due to fact that the particle number 70 is a magic number in the case of a harmonic oscillator potential. In order to see the dependence on the order of the curvature correction the calculation is repeated for m = 5 and it is displayed also in Figure 3. On the considered γ range the plateau region of the MWCC method is roughly two times larger than the plateau of the PCC method. Furthermore the increase of the order of the curvature correction from 3 to 5 does not modify the result of the MWCC method but there is large change in the case of the PCC method.

In the case of Cosh-type starting kernel the shell correction is shown in Figure 4. Using this kernel the quality of the DCC, MWCC and BPCC curvature correction methods can be judged. In this case also the MWCC method is proved to the best, it gives the largest plateau region. The results of the BPCC and MWCC methods differ only at large γ values. In this large γ region only the MWCC method keeps the plateau property. The importance of the large plateau region derives from the fact that in global mass calculation for a given nucleus only one γ value is considered. It would be desirable to have this fixed γ value in the plateau region.

It remains to study the properties of the curvature correction methods in the case of realistic potential. We restrict ourselves for spherically symmetric mean field. We will consider Saxon-Woods potential with spin-orbit term. The geometrical and strength parameters of the Saxon-Woods potential are taken from the parametrization given in [32]. There is a great difference between harmonic oscillator and realistic potentials. The first one has only bound states whereas the Saxon-Woods potential supports both bound and scattering (continuum) states.

It is known that the presence of continuum states and the finiteness of the number of bound states cause difficulties in the shell correction method [START_REF] Nazarewicz | [END_REF]33]. The formula (1) for the quantum SPLD contains only the contribution of bound states. The exact form of the continuum SPLD is known [33] and it was used in [33,34]. The continuum SPLD is expressed by the help of the scattering phase shift. In shell correction calculations the bound states usually are obtained by diagonalizing the Hamiltonian on a given square integrable basis. On such a basis however it is difficult to calculate the continuum SPLD. An elegant way to solve this problem was developed in [35] and was used in shell correction calculation [36]. In this method two diagonalizations have to be carried out. The first diagonalization is the standard one (we take the Hamiltonian as it is). In the second diagonalization the interaction is switched off and only the kinetic energy is diagonalized using the same basis as in the first diagonalization. The results of the second diagonalization are denoted by

( ) = 1 γ i d i K -i γ -d (0) i K - (0) i γ , (34)
where the summation over i extends to all states gained from the diagonalizations. The value of the parameter c is determined in the same way as in the harmonic oscillator case. Since now we have only finite number of bound states the upper bound of the averaging window has a natural value N b i=1 d i . It corresponds to a nucleus where all bound states are occupied. The lower bound of the averaging window was N τ = 50. We have checked that the result is hardly depends on this value.

The neutron shell correction for 146 Gd is displayed in Figures 5 and6. In the case of Gaussian starting kernel (Figure 5) the order of the curvature correction was three and five. Results using Cosh-type starting kernel are given in Figure 6 here the order of the curvature correction was seven. From Figures 5 and6 the same conclusion can be drawn. In the case of the method MWCC a very large plateau region is present (the shell correction is practically independent from the value of γ). The other shell correction methods have no plateau, they satisfy only the local plateau condition (the derivative of the shell correction with respect to γ is zero at some γ point). The plateau property can be made quantitative if the total variation of the shell correction is determined in a given γ region. Considering Figure 5 the total variation are 0.3 MeV and 3 MeV for the MWCC and the PCC methods. In the case of the harmonic oscillator the total variation of the shell correction is only 0.01 MeV for the MWCC method. The plateau property in the realistic potential case is not as good as in the harmonic oscillator case.

The dependence on the order of the curvature corrections is the same as in the harmonic oscillator case. The results of the MWCC method is practically independent from the parameter m. In the case of the methods PCC, DCC and BPCC the shell correction is sensitive to m.

In the shell correction method the order of the curvature correction and the smoothing range are determined by the plateau condition. In some cases for example from some light nuclei, the plateau condition could not be found [4]. As a light nuclei the 40 Ca is considered now. The neutron shell correction is displayed in Figure 7. The results are calculated using the methods PCC and MWCC. The order of the curvature correction was 3, 4 and 5. The number of particles now is only twenty but we can observe the same pattern as in the earlier calculations when the particle number was eighty two. The shell correction as the function of γ is more stable in the case of the method MWCC than in the method PCC. The dependence on the order of the curvature correction is not as good as in the case of the heavier nucleus 146 Gd. However it is clear from this calculation that the method MWCC performs better than the method PCC even in the light mass region.

Summary

In the Strutinsky's shell correction method the smooth single particle level density is calculated by a convolution. The aim of the curvature correction is to satisfy a self consistency condition, the smoothing (convolution) must leave a smooth function unchanged. We are aware only two curvature correction methods but in global mass calculations only one of them is used. We call the standard method polynomial curvature correction since the curvature correction is achieved by multiplying the kernel of the convolution by a polynomial. Two new curvature correction methods are introduced. One of the new methods (we call it many width curvature correction (MWCC)) is proved to be the best among the investigated curvature correction methods.

For potentials without scattering states (harmonic oscillator and similar potentials) the shell correction is practically constant (plateau condition) on a very large region over the smoothing width in each curvature correction methods. However the length of the plateau in the MWCC method is much larger than in the other methods.

For potentials supporting scattering states (realistic potentials) the advantage of the MWCC method is more pronounced. In the case of each methods except the MWCC method only the local plateau condition is satisfied. Due to the lack of convergence in the order of the curvature correction it is difficult to define a unique value for the shell correction. If the calculation is carried out by the MWCC method a large plateau region exists even in the case of realistic potentials and the shell correction is quite stable with respect to the increase of the order of the curvature correction. It seems to us that the MWCC method is better suited for shell correction calculation than the standard procedures.

  smooth SPLD which contains the effect of scattering states reads ḡ

Figure 1 .

 1 Figure captions

Figure 2 .

 2 Figure 2.Shell correction for the harmonic oscillator potential ( ω = 6 MeV) as a function of the smoothing width using different averaging windows for the optimal determination of the value of the parameter c. The order of the curvature correction is three. The kernels are derived from Gaussian starting kernel using the MWCC method. The particle number is N τ =20 for part (a) and N τ =70 for part (b).

Figure 3 .

 3 Figure 3. Shell correction for the harmonic oscillator potential ( ω = 6 MeV, N τ =70) as a function of the smoothing width using different curvature correction methods.The order of the curvature correction is three and five. The kernels are derived from Gaussian starting kernel. In this case the methods PCC, BPCC and DCC give the same kernel for the smoothing.

Figure 4 .

 4 Figure 4. Same as Figure 3 but the kernels are derived from Cosh-type starting kernel. The order of the curvature correction is three.

Figure 5 .

 5 Figure 5.Neutron shell correction for 146 Gd as a function of the smoothing width using different curvature correction methods. The single particle mean field was a Saxon-Woods potential with spin-orbit term. The order of the curvature correction is three and five. The kernels are derived from Gaussian starting kernel. In this case the methods PCC, BPCC and DCC give the same kernel for the smoothing.

Figure 6 .

 6 Figure 6. Same as Figure 5 but the kernels are derived from Cosh-type starting kernel. The order of the curvature correction is seven.

Figure 7 .

 7 Figure 7.Neutron shell correction for 40 Ca as a function of the smoothing width using the curvature correction method PCC and MWCC. The single particle mean field was a Saxon-Woods potential with spin-orbit term. The order of the curvature correction is three, four and five. The kernels are derived from Gaussian starting kernel.

Acknowledgments

This work has been supported by the Hungarian OTKA fund No. K72357.


Shung-Ichi Ando

Michael C Birse

Effective Field Theory of 3 He

3 He and the triton are studied as three-body bound states in the effective field theory without pions. We study 3 He using the set of integral equations developed by Kok et al. which includes the full off-shell T -matrix for the Coulomb interaction between the protons. To leading order, the theory contains: two-body contact interactions whose renormalized strengths are set by the NN scattering lengths, the Coulomb potential, and a three-body contact interaction. We solve the three coupled integral equations with a sharp momentum cutoff, Λ, and find that a three-body interaction is required in 3 He at leading order, as in the triton. It also exhibits the same limit-cycle behavior as a function of Λ, showing that the Efimov effect remains in the presence of the Coulomb interaction. We also obtain the difference between the strengths of the three-body forces in 3 He and the triton.

I. INTRODUCTION

Since Weinberg first proposed applying the ideas of effective field theory (EFT) to nuclear forces [1], much effort has gone into this approach. (For reviews, see Refs. [START_REF] Beane | At the frontier of particle physics: handbook of QCD[END_REF][START_REF] Bedaque | [END_REF][4].) Although there is still some debate about how best to implement it at energies where pion-exchange forces are resolved, the picture is clearer at lower energies. Here few-nucleon systems can be described by a "pionless" EFT based on two-and three-body contact interactions [5][6][7][8].

The resulting expansion of the two-body force is just that of the effective-range expansion [9], but the EFT framework makes it possible to extend this consistently to other effective operators and three-body forces. This theory has been applied extensively to two-body systems, where it has been extended to include the effects of the Coulomb interaction on proton-proton scattering [10][11][12][13][14]. In that system, it corresponds to a distorted-wave or "modified" version of the effective-range expansion [9].

Confidential: not for distribution. Submitted to IOP Publishing for peer review 28 June 2010 Three-body systems have also been studied in the pionless EFT [5,8]. In these, exchange of one particle between an interacting pair and the third particle leads to a long-range force that can be either attractive or repulsive depending on the overall symmetry of the wave function. Attractive channels, such as the triton, display the Efimov effect [15] in the limit of infinite two-body scattering length. This consists of an infinite tower of bound states with energies in a constant ratio. It results from a discrete scale invariance in these systems and corresponds to a limit cycle of the renormalisation group (RG) for the contact three-body force [16].

The 3 He nucleus is of particular interest since it can potentially provide access to properties of the neutron that require a polarised target. It should be possible to describe its low-energy properties within the pionless EFT, extended to include the Coulomb interaction. However the only previous application of the theory to the pd system is in the work of Rupak and Kong, who treated Coulomb effects perturbatively and only considered the spin-3/2 channel [17].

A nonperturbative treatment may not be essential for the 3 He ground state, given the typical momenta in the wave function. However it will be necessary if the EFT is to be extended to describe pd scattering near threshold. Also, it is needed to answer the question of whether the Efimov effect survives in the presence of the Coulomb interaction. Arguments based on the degrees of the singularities in the potentials suggest it should and this has been checked by Hammer and Higa for a simpler two-body model [18]. However it has not previously been confirmed in a three-body system.

Like most other EFT studies of the triton and other three-body systems [START_REF] Bedaque | [END_REF], our treatment is formulated in terms of a set of integral equations. However a quite different approach, based on the resonating-group method, was recently applied to three-and four-nucleon systems by Kirscher et al. [19]. While this also treats the Coulomb interaction to all orders, limitations on the range of cutoffs mean that it is not able to test whether the Efimov effect occurs.

Here we report the results of a study of 3 He in the pionless EFT with the Coulomb interaction treated nonperturbatively. Our approach is based on that used by Kok et al. [20] to study 3 He with separable potentials and the full off-shell Coulomb T -matrix. The resulting set of integral equations is very similar to that developed by Skornyakov and ter-Martirosian [21] which have more recently been used in EFT studies of the triton and other three-body systems [START_REF] Bedaque | [END_REF].

Our results show that the Efimov effect does indeed occur in the presence of the Coulomb interaction. We also determine the strength of the three-body force needed to reproduce observed 3 He binding energy. The difference between this and the corresponding force needed for the triton is of the expected size for an electromagnetic effect.

II. LAGRANGIAN

In the pionless EFT, the strong forces between nucleons are described by two-and threebody contact interactions. In the present context it is convenient to represent the two-body interactions in terms of dibaryon fields. The resulting effective Lagrangian for the threenucleon system can be written as [8,22]

L = L N + L s + L t + L 3 , (1)
where L N is the standard one-nucleon Lagrangian in the heavy-baryon formalism,

L N = N † iv • D + 1 2m N (v • D) 2 -D 2 N . (2)
Here v µ is a velocity vector satisfying a condition v 2 = 1, D µ is the covariant derivative, and m N is the nucleon mass. The terms L s and L t are dibaryon effective Lagrangian for spin singlet and triplet parts, respectively, and these read

L s = s † a iv • D + ∆ s(a) s a -y s s † a N T P ( 1 S 0 ) a N + h.c. , (3)
L t = t † i {iv • D + ∆ t } t i -y t t † i N T P ( 3 S 1 ) i N + h.c. , (4)
where s a and t i are the corresponding dibaryon fields. The strengths of the corresponding two-body interactions depend on ∆ s(a) and ∆ t , the mass differences between the dibaryons and two nucleons, and y s and y t , the coupling constants for the dibaryon-nucleon-nucleon vertices. The projection operators for the two-nucleon 1 S 0 and 3 S 1 states are

P ( 1 S 0 ) a = 1 √ 8 τ 2 τ a σ 2 , P ( 3 S 1 ) i = 1 √ 8 τ 2 σ 2 σ i , (5)
where τ a and σ i are Pauli matrices for isospin and spin, respectively.

The present work is based on the leading-order (LO) terms in the expansion of this EFT of Q. Here Q stands for any of the low-energy scales in the problem: the relative momenta, 1/a where a is any of the NN scattering lengths, and the inverse of the Bohr radius κ = α EM m N /2. More precisely, we keep all terms of order Q -1 , noting that the strong attraction in the NN S waves enhances the wave functions near the origin and hence promotes contact interactions proportional to 1/a or κ to this order [6,7,23,[START_REF] Birse | Proc. of the 6th International Workshop on Chiral Dynamics[END_REF].

The Coulomb interaction leads to a well-known logarithmic divergence in the NN loop integrals [10][11][12]. This is renormalised by a LO counterterm that contributes to the pp scattering length. Given that we need to include isospin breaking in this channel, we have also taken account of charge-independence breaking between the 1 S 0 np and nn scattering lengths, even though this is formally of higher order. This allows us to present the general structures of the three-body equations, which will be needed for extensions of this work to higher orders. We have not included the neutron-proton mass difference which is not only of higher order but is smaller than might be expected as a result of cancellations between electromagnetic and strong contributions [START_REF] Gasser | [END_REF]. To the order we work, observables depend only on the combinations ∆ s(a) /y 2 s and so we have chosen to leave the couplings y s isospininvariant.

The three-body force is expressed as a dibaryon-nucleon contact interaction. It is given by the effective Lagrangian,

L 3 = m N H(Λ) 3Λ 2 y 2 t N † ( σ • t) † ( σ • t)N -y s y t N † ( σ • t) † ( τ • s)N + h.c. + y 2 s N † ( τ • s) † ( τ • s)N , (6)
where H(Λ) is the coupling constant, which runs with the scale Λ of the cutoff we impose on the coupled integral equations. Note that this interaction contributes only to three-body channels with total spin 1/2 and total isospin 1/2.

The building blocks needed to construct the three-body integral equations from the Lagrangian are as follows. In the one-body sector, we have the non-relativistic nucleon propagator,

iS N (p) = i p 0 -p 2 2m N + i . (7)
In the two-body sector, we need the dressed dibaryon propagators with and without the Coulomb interaction. The propagators for the np and nn channels, which have no Coulomb = + + + ... interaction, are represented by the diagrams in Fig. 1 and are given by

iD s,t E - 3q 2 4m N = 4π m N y 2 s,t -i -µ -4π∆s,t m N y 2 s,t + 3 4 q 2 -m N E -i + i . ( 8
)
The two-nucleon loop diagrams here have been dimensionally regularised using the PDS scheme with the subtraction scale µ. The dependence on this can be removed by renormalising the constants 4π∆ s,t /(m N y 2 s,t ) using

µ + 4π∆ s(np) m N y 2 s = 1 a np , µ + 4π∆ s(nn) m N y 2 s = 1 a nn , µ + 4π∆ t m N y 2 t = γ , (9)
where a np and a nn are the scattering lengths for the spin-singlet np and nn channels, respectively, and γ is related to the deuteron binding energy

B 2 through γ = √ m N B 2 . Working
to leading order, we neglect corrections from the effective ranges.

= + + + ... In the pp channel, the Coulomb interaction dresses the two-nucleon Green's function in the bubble diagram for the dressed dibaryon propagator, as shown in Fig. 2. The resulting propagator is given by

iD s(pp) E - 3q 2 4m N = 4π m N y 2 s -i - 4π∆ R s(pp) my 2 s -2κH(κ/p ) , (10)
where

H(η) = ψ(iη) + 1 2iη -ln(iη) , (11)
ψ is the logarithmic derivative of the Γ function, and

-ip = 3 4 q 2 -m N E -i . The renormalised constant ∆ R s(pp)
for the pp channel cancels both the linear and logarithmic divergences of the loop diagram. It is related to the Coulomb scattering length a C by

1 a C = 4π∆ R s(pp) m N y 2 s = 4π∆ s(pp) m N y 2 s + µ -2κ 1 -C E + ln µ 4κ , (12)
where C E = 0.577215 • • • is Euler's constant. Note that the logarithmic divergence means that it is impossible to make a model-independent decomposition of 1/a C into strong and electromagnetic contributions [10][11][12]. This implies that within this EFT there is no unambiguous way to separate Coulomb from other isospin-breaking effects.

The final building block for the integral equations is the off-shell Coulomb T-matrix. A convenient form for it is the integral representation [26] for negative energies, k 2 /m N < 0:

p |T C (k 2 /m N )|p = e 2 (p -p) 2 1 -4iη 1 0 dt t iη 4t -(1 -t) 2 (x 2 -1) , (13)
where η = κ/k and

x 2 = 1 + (p 2 -k 2 )(p 2 -k 2 ) k 2 (p -p) 2 . ( 14
)

III. INTEGRAL EQUATIONS

Now we construct the integral equations for scattering of a third nucleon off a deuteron, concentrating on the channels where the third nucleon is in an S wave, and the total spin and isospin of the three particles are both 1/2. The negative-energy solutions of these describe the bound triton and 3 He states.

We work the center of mass frame for the Nd scattering and use the notation of Bedaque et al. [8]. In the isospin-symmetric case, this process can be described by an amplitude a(p, k) for scattering into states with a 3 S 1 dibaryon, and b(p, k) for scattering into those with a 1 S 0 dibaryon. Here k denotes the initial (on-shell) relative momentum of the deuteron and the third nucleon, and p the final (off-shell) momentum.

The amplitudes satisfy coupled integral equations which can be represented diagrammatically as in Fig. 3. The kernels of these equations consist of two terms: one-nucleon-exchange, which provides a long-range force between the dibaryon and the third nucleon, and the threebody contact interaction.

In the limit of infinite two-body scattering lengths (and zero deuteron binding energy) the long-range force is scale-independent. It is this feature that leads to the Efimov effect in channels where the force is attractive. When we include the Coulomb interaction, isospin symmetry is broken and we must introduce separate amplitudes b + (p, k) and b 0 (p, k) for pd scattering into states with 1 S 0 pp and np dibaryons, respectively. These and a(p, k) satisfy a set of three coupled integral equations. For completeness, we also present the corresponding set of three equations that arise when we allow for isospin-breaking effects in the nd case, involving the amplitude b -(p, k) for states with nn dibaryons. As already mentioned, these are not essential for the present LO calculations but will be needed for extensions to higher orders.

A. Integral equations for the nd channel

Before presenting the equations for pd scattering, we first look at the simpler equations for nd scattering. Allowing for isospin breaking in the scattering lengths, the three integral equations for the S = 1/2, T = 1/2 channel are

a(p, k) = K (a) (p, k; E) + 2 H(Λ) Λ 2 + 1 π Λ 0 dq q 2 K (a) (p, q; E) + 2 H(Λ) Λ 2
a(q, k)

-γ + 3 4 q 2 -mE + 1 π Λ 0 dq q 2 K (a) (p, q; E) + 2 3
H(Λ) Λ 2 b 0 (q, k) -1 anp + 3 4 q 2 -mE + 2 π Λ 0 dq q 2 K (a) (p, q; E) + 2 3
H(Λ) Λ 2 b -(q, k) -1 ann + 3 4 q 2 -mE , ( 15
) b 0 (p, k) = 3K (a) (p, k; E) + 2 H(Λ) Λ 2 + 1 π Λ 0 dq q 2 3K (a) (p, q; E) + 2 H(Λ) Λ 2
a(q, k)

-γ + 3 4 q 2 -mE + 1 π Λ 0 dq q 2 -K (a) (p, q; E) + 2 3
H(Λ) Λ 2 b 0 (q, k) -1 anp + 3 4 q 2 -mE + 2 π Λ 0 dq q 2 K (a) (p, q; E) + 2 3
H(Λ) Λ 2 b -(q, k) -1 ann + 3 4 q 2 -mE , (16) b
-(p, k) = 3K (a) (p, k; E) + 2 H(Λ) Λ 2 + 1 π Λ 0 dq q 2 3K (a) (p, q; E) + 2 H(Λ) Λ 2
a(q, k)

-γ + 3 4 q 2 -mE + 1 π Λ 0 dq q 2 K (a) (p, q; E) + 2 3
H(Λ) Λ 2 b 0 (q, k) -1 anp + 3 4 q 2 -mE + 2 π Λ 0 dq q 2 2 3 H(Λ) Λ 2 b -(q, k) -1 ann + 3 4 q 2 -mE , (17)
where

K (a) (p, q; E) = 1 2pq ln p 2 + q 2 + 2pq -m N E p 2 + q 2 -2pq -m N E . ( 18
)
If we set a nn = a np and b 0 = b -= b, these reduce to the isospin-symmetric forms of the equations in Ref. [8].

The integrals over the relative momentum q are all cut off at q = Λ. The resulting dependence on Λ can be cancelled by the three-body force, whose strength is H(Λ)/Λ 2 .

B. Integral equations for the pd channel

We now turn to the corresponding equations for pd scattering. These differ from the ones in the previous subsection by having b + instead of b -and, more importantly, adding the Coulomb interaction between the two protons. They can be obtained from the equations developed by the Groningen group [20], by omitting the separable form factors and instead imposing a sharp cutoff on the relative momenta.

The Coulomb interaction leads to additional long-range terms in the kernels of these The coupled integral equations can be written

a(p, k) = K (a) (p, k; E) + K (c) (p, k; E) + 2K (d) (p, k; E) + 2 H(Λ) Λ 2 + 1 π Λ 0 dq q 2 K (a) (p, q; E) + K (c) (p, k; E) + 2K (d) (p, k; E) + 2 H(Λ) Λ 2 × a(q, k) -γ + 3 4 q 2 -mE + 1 π Λ 0 dq q 2 K (a) (p, q; E) + K (c) (p, k; E) + 2 3
H(Λ) Λ 2 b 0 (q, k) -1 anp + 3 4 q 2 -mE + 2 π Λ 0 dq q 2 K (a) (p, k; E) + K (b) 13 (p, k; E) + 2 3
H(Λ) Λ 2 b + (q, k) -1 a C -2κH(κ/p ) , ( 19
) b 0 (p, k) = 3K (a) (p, k; E) + 3K (c) (p, k; E) + 2 H(Λ) Λ 2 + 1 π Λ 0 dq q 2 3K (a) (p, q; E) + 3K (c) (p, q; E) + 2 H(Λ) Λ 2 a(q, k) -γ + 3 4 q 2 -mE - 1 π Λ 0 dq q 2 K (a) (p, q; E) + K (c) (p, q; E) + 2K (d) (p, q; E) - 2 3
H(Λ) Λ 2 × b 0 (q, k) -1 anp + 3 4 q 2 -mE + 2 π Λ 0 dq q 2 K (a) (p, q; E) + K (b)
13 (p, q; E) + 2 3

H(Λ) Λ 2 b + (q, k) -1 a C -2κH(κ/p ) , (20) b
+ (p, k) = 3K (a) (p, k; E) + 3K (b) 31 (p, k; E) + 2 H(Λ) Λ 2 + 1 π Λ 0 dq q 2 3K (a) (p, q; E) + 3K (b) 31 (p, q; E) + 2 H(Λ) Λ 2 a(q, k) -γ + 3 4 q 2 -mE + 1 π Λ 0 dq q 2 K (a) (p, k; E) + K (b) 31 (p, q; E) + 2 3
H(Λ) Λ 2 b 0 (q, k) -1 anp + 3 4 q 2 -mE + 2 π Λ 0 dq q 2 2 3 H(Λ) Λ 2 b + (q, k) -1 a C -2κH(κ/p ) , (21)
where K (a) is the one-nucleon exchange kernel defined in Eq. ( 18) above, and K 13 , the subscripts 1 and 3 refer to initial and final dibaryons, 1 denoting np and 3 pp.)

In diagram (b), the dibaryon-nucleon-nucleon vertex projects out only the S-wave part of the Coulomb T-matrix. This considerably simplifies the corresponding terms in the equations by reducing the number variables that need to be integrated from six to two. The corresponding terms in the kernels are

K (b) 13 (p, q; E) = 1 2 1 -1 dy φ (b)
13 (p, q, y; E) p 2 + q 2 -mE + pqy ,

K (b) 31 (p, q; E) = 1 2 1 -1 dy φ (b) 31 (p, q, y; E) p 2 + q 2 -mE + pqy , (22)
where y = p • q = cos θ and

φ (b)
13 (p, q, y; E) = -

4κβ p 2 + β 2 1 0 dt t κ/β 1 1 + t 2 -2 p 2 -β 2 p 2 +β 2 t , (24) φ (b)
31 (p, q, y; E) = -

4κβ p 2 + β 2 1 0 dt t κ/β 1 1 + t 2 -2 p 2 -β 2 p 2 +β 2 t . ( 25
)
Here we have defined the momentum variables

p = p + 1 2 q , p = q + 1 2 p , β = 3 4 q 2 -mE , β = 3 4 p 2 -mE . (26)
The kernels for diagrams (c) and (d) involve six dimensional integrals but we can analytically evaluate two of these (over the azimuthal angles). In the case of K (c) , it is convenient to introduce the vectors,

a = p + q , b = p -q , (27)
Choosing axes where a lies along the z-direction, these can be written

a = a (0, 0, 1) , b = b (sin θ , 0, cos θ ) , l = l (sin θ cos φ , sin θ sin φ , cos θ ) , ( 28
)
in terms of their magnitudes

a = p 2 + q 2 + 2pq cos θ , b = p 2 + q 2 -2pq cos θ , (29)
and the angles defined by sin θ = 2pq sin θ

(p 2 -q 2 ) 2 + 4p 2 q 2 sin 2 θ , cos θ = p 2 -q 2 (p 2 -q 2 ) 2 + 4p 2 q 2 sin 2 θ . ( 30
)
The kernel then becomes

K (c) (p, q; E) = - κ π ∞ 0 dl 1 -1 d(cos θ) 1 -1 d(cos θ ) 1 0 dt (κ/β)t κ/β-1 bl d 2 + 4t(1 -t) -2 l 2 β 2 × 1 cos θ cos θ + 2 bl d 2 + 4t(1 -t) -2 l 2 β 2 1 1 -y 2 1 - 1 cos θ cos θ -2 bl d 2 + 4t(1 -t) -2 l 2 β 2 1 1 -y 2 2 , (31)
where

y 1,2 = sin θ sin θ cos θ cos θ ± 2 bl d 2 + 4t(1 -t) -2 l 2 β 2 , d = l 2 + p 2 + q 2 + pq cos θ - 3 2 al cos θ -m N E , β 2 = 3 4 (l 2 + a 2 + 2al cos θ ) -m N E . ( 32
)
The last term in the kernel, arising from diagram (d), requires the most care in its numerical evaluation. This is because the long-range photon exchange gives rise to a logarithmic IR singularity. In the present work, where we are interested only in bound states, the finite extent of the wave functions helps to regulate this singularity. Calculations of pd scattering amplitudes will, however, require different numerical methods. Following Ref. [20], we split K (d) up into three pieces,

K (d) (p, q; E) = K (d,V ) (p, q; E) + K (d,V ) (p, q; E) + K (d,T -V ) (p, q; E) , (33)
where K (d,V ) and K (d,V ) denote the contributions the from singular and nonsingular parts of one-photon exchange, and K (d,T -V ) the remaining (also nonsingular) terms from iterating the Coulomb potential. These have pieces have the forms

K (d,V ) (p, q; E) = - κ 4pqβ ln p 2 + q 2 + 2pq p 2 + q 2 -2pq , (34)
K (d,V ) (p, q; E) = -κ 1 -1 d(cos θ) 1 p 2 + q 2 -2pq cos θ 1 p 2 + q 2 -2pq cos θ × arctan -p 2 + 2q 2 -pq cot θ 2β p 2 + q 2 -2pq cos θ + arctan 2p 2 -q 2 -pq cos θ 2β p 2 + q 2 -2pq cos θ - 1 2β , (35)
K (d,T -V ) (p, q; E) = -κ 1 -1 d(cos θ) 1 p 2 + q 2 -2pq cos θ × 1 π ∞ 0 dl l 2 1 -1 d(cos θ ) 1 0 dt (κ/β l )t κ/β l -1 D √ 1 -c 2 - 1 p 2 + q 2 -2pq cos θ arctan -p 2 + 2q 2 pq cos θ 2β p 2 + q 2 -2pq cos θ + arctan 2p 2 -q 2 -pq cos θ 2β p 2 + q 2 -2pq cos θ , (36)
where

β = 3 4 p 2 -m N E , β = 3 4 q 2 -m N E , β l = 3 4 l 2 -m N E , D = (l 2 + q 2 + ql cos θ -m N E)(l 2 + p 2 + pl cos θ -m N E) +4t(1 -t) -2 β 2 l (p 2 + q 2 -2pq cos θ) , c = (l 2 + q 2 + ql cos θ -m N E)pl sin θ sin θ /D . (37)

IV. NUMERICAL RESULTS FOR BOUND STATES

We solve for the three-body bound states by taking the homogeneous parts of the coupled integral equations and treating them as nonlinear eigenvalue problems,

K(E) u = u , (38)
where u T = (a, b 0 , b ± ) and the 3 × 3 matrices of integral operators K(E) can be extracted from the expressions in the previous section.

To convert K(E) into a matrix we use Gauss-Legendre quadrature to replace the momenta p and q by discrete variables. We have examined three versions of this discretisation, in terms of variables x defined by: (1) p = m tan(x) where m = 140 MeV, (2) p = x (no change of variable), and (3) p = m 0 exp(x) where m 0 = 1 MeV. The variables p and q run from zero to the cutoff Λ, except in the third of these where we also need to impose a cutoff at the very small momentum m 0 .

We the construct the determinant of the matrix K(E) -1 and search on E to find the energy eigenvalues at which det[K(E) -1] = 0. For one of these values of E, we then solve the eigenvalue problem Eq. ( 38) and determine the eigenvector corresponding to the unit eigenvalue of the matrix K(E).

For given values of the two-body input parameters we first set the three-body force H(Λ) to zero and adjust the value of the cutoff until the shallowest bound-state eigenvalue reproduces the observed three-body binding energy. Starting from this cutoff, Λ 0 , we then determine the values of H(Λ) needed to reproduce this binding energy for other cutoffs.

A. Triton

In this case, we work with the determinant constructed from the homogeneous part of the coupled integral equations in Eqs. ( 15), ( 16), and ( 17). The two-body input is provided by γ = 45.703 MeV , a np = -23.749 ± 0.008 fm [27] , a nn = -18.5 ± 0.4 fm [28] .(39)

We adjust the three-body force to fit the triton binding energy,

B 3 ( 3 H) = 8.48182 MeV , (40)
which corresponds to a momentum scale α = √ m N B 3 = 89.239 MeV.

For 16 Gaussian integration points (corresponding to a 48×48 matrix) we find that our results have converged and, for the three choices of variable discussed above, they agree to 6 significant figures.

For H(Λ) = 0, we reproduce the experimental binding energy with the cutoff Λ 0 = 380.689 MeV. The values of the three-body force needed for other cutoffs are shown in Fig. 5.

As we take Λ → ∞, the cutoff becomes much larger than the physical scales in the system (1/a) and H(Λ) displays the limit-cycle behaviour found in Ref. [8]. This reflects the presence of the Efimov effect in this three-body system, a tower of deeply bound states with energies in a constant ratio [15]. In three-nucleon systems only the shallowest of these states lies within the domain of the EFT and so corresponds to a physical state.

The asymptotic form of H(Λ) in this limit is [8]

H(Λ) = - sin[s 0 ln(Λ/Λ * ) -arctan(1/s 0 )] sin[s 0 ln(Λ/Λ * ) + arctan(1/s 0 )] . (41)
where s 0 = 1.00624 • • • and the scale parameter Λ * is defined by

Λ 0 = Λ * exp[(1/s 0 ) arctan(1/s 0 )] . ( 42
)
This form is shown by the dashed line in Fig. 5. The numerical results deviate from it only for small values of Λ where the finite scattering lengths can influence the renormalisation of the three-body force.

B. 3 He channel

For the 3 He channel, we need to the pp Coulomb scattering length as a two-body input parameter, a C = -7.8063 ± 0.0026 fm [29].

The experimental binding energy of 3 He is

B 3 ( 3 He) = 7.71804 MeV , (44)
corresponding to a scale α = 85.127 MeV

As already noted, the K (d,V ) term contains a logarithmic singularity that needs to be dealt with carefully. We handle it with the method outlined in Appendix B in Ref. [20]. In the case of a bound state, the other factors in the integral are regular and so the logarithmic singularity has a finite integral. We can take the interval around the singular point and, treating the other factors as constant, explicitly integrate the logarithm over it. The result can be expressed in a similar form to a Gaussian quadrature, b a dq ln(p j -q) 2 y(q) = N i =j

w i ln(p j -q i ) 2 y(q i ) + w j y(p j ) , (45)
where y(q) is a smooth function and the w i are the usual weights. The modified weight for the singular point is given by w j = q j> q j< dq ln(p j -q) 2 = -2w j + 2(q j> -p j ) ln(q j> -p j ) + 2(p j -q j< ) ln(p j -q j< ) ,

where q j< = a + j-1 k=1 w k and q j> = a + j k=1 w k . To examine the size of Coulomb effects in the three-nucleon system, we first consider an isospin-symmetric three-body force and take H(Λ 0 ) for the same value Λ 0 = 380.689 MeV as in the triton. We use N 4 = 32 Gaussian points for the four integrals needed to evaluate the kernels K (c) and K (d) . The integral equations are discretised using up to N = 32 points (leading to 96 × 96 matrices). We fit the results for N = 16 and 32 to the form α(N ) A + B/N in order to extrapolate to N = ∞.

Using the first choice for the discretisation variable (p = m tan(x)), we get α = 84.8084 MeV, or a binding energy of B 3 = 7.66038 MeV. This differs from the triton energy for the same three-body force by 0.82 MeV. It is within 1% of the observed 3 He energy, indicating that the isospin-violating three-body force is indeed a higher-order contribution. This is as expected, given the absence of any modification of the renormalisation, since the isospin-violating term is suppressed by one power of the inverse Bohr radius κ relative to the LO force. Our difference between the triton and 3 He binding energies is significantly larger than the 0.66 ± 0.03 MeV found by Kirscher et al. [19]. However, although working to NLO, those authors took a charge-independent value for the 1 S 0 scattering length. Their results are consistent with other estimates of the pure Coulomb contribution [30,31].

Because we take different values for the nn and np 1 S 0 scattering lengths, our calculations do include some higher-order isospin-breaking effects. Given that we have not done a consistent calculation to N 2 LO (or even NLO), the best we can say is that the ∼160 keV difference between our result and that of Ref. [19] is comparable in size to the contributions of 70-110 keV that have been estimated using other approaches [30][31][32][33]. One caveat that should be made about these comparisons is that, as noted above, there is no model-independent way to separate Coulomb and other isospin-breaking effects in this EFT.

By adjusting the strength of the three-body force, we can fit the experimental binding energy. As in the triton case, we have done this for a range of cutoffs. The results are plotted in Fig. 6. They display the same limit-cycle behaviour as in the triton channel. Moreover the differences between the forces needed in the two channels are very small, reflecting the fact that a symmetric force can give a very good account of 3 He.

To quantify the size of the isospin breaking in the three-body forces, it is better not to consider the values of H(Λ), since these cycle between ±∞. The scale Λ 0 or, equivalently Λ * , is a more appropriate measure. Its value of Λ 0 = 382.46 MeV differs from that for the triton by 0.46%. The contribution of this force to the triton-3 He splitting is about 50 keV. This is much larger than the 5 keV estimate of the three-body contribution in Ref. [33], which may reflect the fact that other isospin-breaking effects are missing from our calculatation. Also, the results of Ref. [33] were obtained using an EFT that includes pion fields.

V. SUMMARY

We have studied 3 He within the framework of the pionless EFT, treating the Coulomb interaction nonperturbatively. We solve the set of integral equations developed by Kok et al. [20], which can be thought of as extending the equation of Skornyakov and ter-Martirosian [21] to include the full off-shell Coulomb T -matrix.

We find that a three-body interaction is required at leading order in 3 He, as in the triton.

This force also exhibits the same limit-cycle behavior as a function of the cutoff Λ as found by Bedaque et al. [8] for systems with purely short-range interactions. These results show that the Coulomb interaction is not singular enough to alter the deeply bound "Efimov" states in the three-nucleon system.

We find that an isospin symmetric three-body force, fit to the triton binding energy, can give the 3 He binding to better than 1%. The scale parameter of the force that fits 3 He differs from the corresponding value for the triton by about 0.5%. This is of the expected size for an order-α electromagnetic effect.

Our results demonstrate that the Coulomb interaction has no nonperturbative effect on the renormalisation of the three-body force in 3 He. Moreover the overall contributions of the Coulomb interaction to the binding energy of 3 He are small, perhaps unsurprisingly given the typical momenta involved, and so could have been calculated perturbatively. This will not be the case if our approach is extended to describe low-energy pd scattering. However, as noted above, such an extension will require improved numerical techniques to handle the singularity of the Coulomb interaction.

Our approach could also be used to determine the 3 He wave function from the corresponding eigenvector of the set of integral equations. This could then be used in calculations of electromagnetic properties of 3 He, based on the methods in Ref. [34].

Our present results do include certain other isospin-breaking effects beyond LO in this EFT, through the nn and np 1 S 0 scattering lengths. A proper treatment of these in the threenucleon system will require a consistent higher-order calculation. We have provided here the basic forms of the integral equations allowing for isospin-breaking in the NN scattering lengths. The relevant power counting for the terms needed in a full calculation can be determined using the techniques in Refs. [35,36].

FIG. 1 :

 1 FIG. 1: Diagrams for the dressed dibaryon propagator, denoted by the double line with filled circle. Single lines represent nucleon propagators; double lines undressed dibaryon propagators.

FIG. 2 :

 2 FIG. 2: Diagrams for dressed dibaryon propagator for the pp channel, denoted by the double line with shaded circle. The shaded ovals here represent the two-nucleon Green's function dressed with the Coulomb interaction.

FIG. 3 :

 3 FIG.3: Diagrams corresponding to the coupled integral equations for nucleon-dibaryon scattering.Thin lines denote nucleons, double and thick lines 3 S 1 and 1 S 0 dibaryons. Those with filled circles are dressed dibaryon propagators, as in Fig.1. Ovals with stripes denote the off-shell amplitude a(p, k) with initial and final spin-triplet dibaryons, and those with cross stripes the amplitude b(p, k) with initial spin-triplet and final spin-singlet dibaryon fields.

  equations. The diagrams for the full long-range kernel are shown in Fig. 4. The one-nucleon exchange term of diagram (a) is supplemented by diagram (b) if a proton is exchanged, and by (c) for neutron exchange. Finally, diagram (d) represents the Coulomb interaction between an np dibaryon and a proton.

FIG. 4 :

 4 FIG. 4: Diagrams corresponding to the long-range terms in the kernel for pd scattering. Here, the shaded oval denotes the off-shell Coulomb T -matrix (since contributions where the two protons do not interact are already contained in (a)). Diagram (b) has been shown for the case of an initial np dibaryon and a final pp one. The corresponding diagram for the inverse process is the mirror image of this.

  ) , and K(d) are the pieces corresponding to the Coulomb diagrams (b), (c), (d) in Fig. 4. (In K (b) 31 and K (b)

FIG. 5 :

 5 FIG. 5: Values of the three-body force H as a function of the cutoff Λ in MeV. For comparison, the dashed curve shows the asymptotic, scale-free form of H(Λ).

FIG. 6 :

 6 FIG.6: Evolution with Λ of the three-body forces H(Λ) for 3 He (solid line) and the triton (dashed).

Acknowledgments

This work was supported by STFC grants PP/F000448/1 and ST/F012047/1. We are grateful to L. Kok and his colleagues in Groningen for providing copies of the reports [20].


Tadeusz Wibig
email: wibig@zpk.u.lodz.pl

The non-extensivity Parameter of a Thermodynamical Model of Hadronic Interactions at LHC energies

Keywords: numbers: 12.40.Ee, 13.85.Hd, 13.85.Rm, 89.75.-k

The LHC measurements above SPS and Tevatron energies give an opportunity to test predictions of the non-extensive thermodynamical picture of hadronic interaction to examine the measured transverse momentum distributions for new a interaction energy range. We determined the Tsallis model non-extensivity parameter for the hadronization process before short-lived particles decay and distort the initial p ⊥ distribution. We have shown that it follows exactly the smooth rise determined at lower energies below the present LHC record. The energy dependence of the q parameter is consistent with expectations and the evidence of the asymptotic limit may be seen.

The thermodynamical concept of the hadronization process following the fifty-year old Hagedorn idea was recently successfully used in [1] to describe multiplicities in pp, e + e -, and N N high energy interactions. However, the main problem of such a treatment: the nonexponential tails of transverse momentum distributions remains. One explanation is proposed in [2] introducing the pre-hadronization dynamics of fireballs. A different, more general, solution is to extend the standard, Boltzmann statistics following the Tsallis idea [3]. The generalised, non-extensive statistics technique is often used to describe different kinds of phenomena: from interstellar dynamics to the cosmic ray source puzzle [4]. Among them high energy interaction phenomenology is one of the fields where the new statistics is very successful.

It was shown that such non-extensive statistics works well for p ⊥ spectra in e + e - [5], and could be used for pp and pp interactions too [6,7]. Here we show that the results of the Tsallis generalisation of the Becattini canonical thermodynamics works well when applied also to the LHC energies: 2.36 and 7 TeV.

In Refs. [9] and [10] the function of the form inspired by the Tsallis formula

E d 3 N ch dp 3 ∼ dN ch dη 1 + (q -1) E ⊥ T -1 q-1 (1)
was used to fit the p ⊥ data. The shape of the invariant transverse momentum distribution is well described by the above formula with the values of the q parameter of 1.13 and 1.15 (with a systematic uncertainty of order 0.005) and T : 130 and 145 MeV for energies 2.36 and 7 TeV, respectively. We would like to use the same transverse momentum data, but not only to adjust the simple and convenient formula parameters. We try to say something about the mechanism of the hadronization process. The fits of the form of Eq.1 show the direction one can follow, but the physical conclusions should be based on rather more detailed analysis of the thermodynamical picture. We would like to remind that the word 'thermodynamical' here should be taken with care. As it was pointed out already by Hagedorn, (e.g., [8]) the meaning of the temperature T in the "hadronic matter near the boiling point" is different than the ordinary temperature of the ordinary boiling soup.

The formalism we used in the present work is described in detail in [7]. The canonical partition function is generalised by replacing the exponential Boltzmann factor by the respective Tsallis formula:

exp - E T -→ 1 + (q -1) T p 2 + m 2 j -1 q-1 , (2)
where q is called the non-extensivity parameter. The transverse momentum distribution is given by the partition function

Z f p 2 T dp 2 T = j ∞ n=1 (-1) n+1 γ sj n (2J j + 1) V 2π 2 × × dp L x j Z(Q 0 -nq j ) Z(Q 0 ) dp 2 T , (3)
where x j is the statistical weight (Eq.( 2) of the particular state (with classical Boltzmann or modified Tsallis factor) γ sj is the strangeness suppression factor, J j is the spin of the j-th type hadron (of quantum numbers q j ).

V denotes the hadronization volume.

Respective numerical integrations for the partition function determination have been performed. Then summation and integration of Eq.( 3) have been performed taking into account almost one hundred species of daughter particles. All unstable particles were then decayed at the end and the final p ⊥ distribution was compared with the distributions measured by CMS at two energies: 2.36 [9] and 7 TeV [10]. The T parameter is sensitive to small p ⊥ 's while q describes the non-exponential tails. No significant and systematic energy dependence of T was found, and eventually we found that this parameter can be fixed at 130 MeV as in our previous paper [7]. Results of the fits are shown in Fig. 1. Values of the q-parameter found for the √ s of 2.36 and 7 TeV are 1.109±0.001 and 1.222±0.001, respectively. The quality of both fits (χ 2 /ndf) is quite good: 19/16 and 10/26, respectively.

The old data at 900 GeV from the SPS experiment [11] are shown together with the fit from [7]. As it is seen the p -4 ⊥ steepness expectations at high p ⊥ are still far and asymptotic possibility.

The non-extensive thermodynamics works well at lower energies, as it has been shown e.g. in [7], not only describing each individual energy transverse momentum spectrum, but the resulting parameter q behaviour exhibits the regular and expected smooth energy dependence. The low energy Boltzmann limit of q = 1 (observed for 100 GeV/c data) together with the upper (asymptotic, for infinite energy) limit of 1.25 [3] were combined in the functional form

q(E) = 1.25 -A s -B >1 . (4)
Parameters A and B were adjusted to the data up to Tevatron energies in [7] and the result is shown in Fig. 2. The values of A and B are 0.33 and 0.054, respectively. In Fig. 2 The open symbol represents the values obtained in [7] for energies up to Tevatron. Two solid circles show values adjusted to CMS data, their uncertainties are much smaller than the symbols diameters.

ter could remain equal to 130 MeV and a smooth energy dependence of the non-extensivity parameter q found for lower energies is confirmed by the new data.

FIG. 1 :

 1 FIG.1: Transverse momentum distributions measured by CMS at 2.36 and 7 TeV and SPS at 900 GeV and respective non-extensive statistics fits (described in the text).

FIG. 2 :

 2 FIG. 2: Energy dependence of the non-extensivity parameter.The open symbol represents the values obtained in[7] for energies up to Tevatron. Two solid circles show values adjusted to CMS data, their uncertainties are much smaller than the symbols diameters.