text:
# Proton-deuteron elastic scattering at low energies

## I Introduction

The three-nucleon continuum is an excellent testing ground for nuclear interactions, since both precise experiments and theoretical calculations are possible . Comparisons between measured three-nucleon observables and theoretical calculations provide stringent tests of the underlying nucleon-nucleon (NN) and three-nucleon (3N) interactions. Proton-deuteron elastic scattering at low energies offers two particular advantages: highly accurate experiments can be performed for a large number of observables, and theoretical calculations are simplified due to the energy being below the deuteron breakup threshold and the relatively low number of partial waves required.

Past experimental studies have focussed on the differential cross section and analyzing powers for p-d scattering at $`E_{c.m.}=2`$ MeV, just below the deuteron breakup threshold . Theoretical calculations of these p-d scattering observables using realistic NN and 3N interactions are now routinely possible . The calculations are found to reproduce the differential cross section and tensor analyzing powers (TAPs) $`T_{20}`$, $`T_{21}`$, and $`T_{22}`$ quite well, while the vector analyzing powers (VAPs) $`A_y`$ and $`iT_{11}`$ are underpredicted by $`30\%`$. The available 3N forces are found to have little influence, except for the $`J^\pi =\frac{1}{2}^+`$ partial wave . The situation in neutron-deuteron scattering is similar , although only the cross section and $`A_y`$ have been measured to date at low energies. By far the most glaring difference between theory and experiment is in the VAPs, a discrepancy which has come to be known as the “$`A_y(\theta )`$ puzzle” . This puzzle has been in existence for about ten years now, and the resolution is still unclear. It has been suggested that the discrepancy arises from inadequacies in the NN potential in the $`{}_{}{}^{3}P_{j}^{}`$ waves . Other workers have considered the role of 3N force effects . So far, no conclusive solution to the problem has been found, although most recent attention has been focussed on the role of the 3N interaction .

It is desirable to determine these observables at energies below $`E_{c.m.}=2`$ MeV, as the influence of higher partial waves is strongly reduced, allowing the dominant $`S`$\- and $`P`$-waves to be investigated with greater clarity. Also, when faced with disagreement between theory and experiment, it is interesting to study the energy dependence of the discrepancy, since this may provide clues as to its origin. New data will place strong constraints on any proposed solution. There is also considerable interest in obtaining low-energy data for the determination of the p-d scattering lengths .

In this paper we present measurements of the differential cross section and analyzing powers $`A_y`$, $`iT_{11}`$, $`T20`$, $`T_{21}`$ and $`T_{22}`$. for p-d scattering at $`E_{c.m.}=431.3`$ keV. These data are the first complete set of observables measured at such a low energy; previously only cross section  and $`A_y`$ data  have been obtained in this energy range. We also present an energy excitation function of $`iT_{11}(\theta _{c.m.}=87.8^{})`$ from $`E_{c.m.}=431.3`$ keV up to 2000 keV. Some of the data presented here have appeared previously in brief publications .

## II Experimental procedures

Measurements of the differential cross section and the analyzing powers $`A_y`$, $`iT_{11}`$, $`T_{20}`$, $`T_{21}`$, and $`T_{22}`$ were carried out at the Triangle Universities Nuclear Laboratory (TUNL). The analyzing power $`T_{22}`$ was deduced from measurements of $`T_{20}`$ and $`A_{yy}`$. Achieving a center-of-mass energy of 431.3 keV requires either a 650-keV proton beam incident on a deuteron target, or a 1300-keV deuteron beam incident on a proton target (the small corrections for energy loss in the target will be discussed later). The beams, targets, and detection methods are discussed below.

### A Low-energy proton beam

Beams of 650-keV protons were used for measurements of differential cross sections and $`A_y`$. Polarized or unpolarized beams of 72-keV $`{}_{}{}^{1}\mathrm{H}_{}^{}`$ ions were produced by the TUNL atomic beam polarized ion source . These ions were then accelerated to 450 keV by the minitandem accelerator . The final energy of 650 keV for the now positive ions was achieved by the $`200`$-kV bias on the 107-cm-diameter scattering chamber  where the measurements were performed.

Unpolarized proton beams were produced by bleeding hydrogen gas into the ionizer of the ion source. When polarized proton beams were used, a Wien filter downstream from the ion source was utilized to set the spin-quantization axis vertical in the laboratory. Two spin states were used, with polarizations $`p_Z\pm 0.7`$. The desired hyperfine states of atomic hydrogen were cycled approximately every second. This technique minimizes the effects of slow changes in beam position, target thickness, or amplifier gain on the measured analyzing powers. Proton beam polarization was determined using the $`{}_{}{}^{6}\mathrm{Li}(\stackrel{}{p},{}_{}{}^{3}\mathrm{He}){}_{}{}^{4}\mathrm{He}`$ reaction in a polarimeter  located at the rear of the scattering chamber. The polarization was measured several times throughout the measurements at an incident proton energy of 450 keV by lowering the chamber bias voltage. The proton polarization was found to be constant within $`\pm 3\%`$ throughout the measurements; the systematic error in the proton polarization is estimated to be $`\pm 4\%`$.

The energy calibration of the proton beam produced by the minitandem and high-voltage chamber system has been established to $`\pm 1`$ keV, using the 240.0- and 340.5-keV resonances in $`{}_{}{}^{19}\mathrm{F}(p,\alpha \gamma )`$ and the 405.4- and 445.8-keV resonances in $`{}_{}{}^{27}\mathrm{Al}(p,\gamma )`$. The resonance energies were taken from Ref. .

### B Deuteron beam

Polarized beams of 72-keV $`{}_{}{}^{2}\mathrm{H}_{}^{}`$ ions were produced using the same atomic-beam polarized ion source  used for protons. The beam was injected into the TUNL FN tandem accelerator, magnetically analyzed, and then delivered to a 62-cm diameter scattering chamber.

The same Wien filter mentioned previously was used to control the spin-quantization axis on target. The spin axis was longitudinal for the $`T_{20}`$ measurements, normal to the reaction plane for $`iT_{11}`$ and $`A_{yy}`$, and $`45^{}`$ offset from longitudinal in the reaction plane for $`T_{21}`$. Three deuteron beam polarization states were produced with the atomic beam polarized ion source: a maximum positive, a maximum negative, and an unpolarized state. The TAP data were obtained with $`p_{ZZ}\pm 0.70`$ and $`p_Z0.25`$, while beams with $`p_Z\pm 0.55`$ and $`|p_{ZZ}|0.05`$ were used for the $`iT_{11}`$ measurements. The spin states were also cycled approximately once every second. The beam polarization for the TAP data was determined using the $`{}_{}{}^{3}\mathrm{He}(\stackrel{}{d},p)`$ reaction in an online polarimeter located behind the scattering chamber . Deuteron beam vector polarization was determined online via the $`{}_{}{}^{12}\mathrm{C}(\stackrel{}{d},p)`$ reaction in a different polarimeter located behind the scattering chamber. The effective $`iT_{11}`$ for this reaction at $`E_d=1.3`$ MeV has been calibrated relative to the $`{}_{}{}^{3}\mathrm{He}(\stackrel{}{d},p)`$ polarimeter at $`E_d=12`$ MeV. More information on the use of the $`{}_{}{}^{12}\mathrm{C}(\stackrel{}{d},p)`$ reaction for deuteron vector polarimetry is available in Ref. . For the $`iT_{11}(87.8^{})`$ measurements the polarization was very stable over time, and was monitored periodically with the $`{}_{}{}^{3}\mathrm{He}(\stackrel{}{d},p)`$ polarimeter at $`E_d=12`$ MeV. The absolute uncertainty in all of the deuteron beam polarizations is estimated to $`\pm 3\%`$. The unpolarized deuteron and proton beams used for the $`{}_{}{}^{1}\mathrm{H}(d,d)`$ cross section measurements were produced by a direct extraction negative ion source.

The energy calibration of the beams produced by the FN tandem is determined within $`\pm 0.1\%`$ from the magnetic field in the $`52^{}`$ bending magnet, measured by an NMR magnetometer. The magnetometer was calibrated using the $`{}_{}{}^{7}\mathrm{Li}(p,n)`$ threshold ($`E_p=1880.443\pm 0.020`$ keV ) and the $`{}_{}{}^{19}\mathrm{F}(p,n)`$ threshold ($`E_p=4234.3\pm 0.8`$ keV ).

### C Targets

The measurements were performed by bombarding thin hydrogenated or deuterated carbon foils. The targets used for analyzing power measurements consisted of approximately $`1\times 10^{18}`$ and $`1.5\times 10^{18}`$ hydrogen isotope and carbon atoms/cm2, respectively. The targets used for the relative and absolute cross section measurements consisted of approximately $`5\times 10^{17}`$ and $`8\times 10^{17}`$ hydrogen isotope and carbon atoms/cm2, respectively. In order to calculate the energy loss in the targets it is necessary to know their carbon and hydrogen isotope thickness and the stopping power. The carbon content was determined by comparison with carbon foils of known thickness while the hydrogen thickness was determined either from the d-p elastic scattering cross section, or from the p-p elastic scattering cross section. For the analyzing power measurements the beam-energy loss in the targets is $`10`$ keV, while for the cross section measurements, the beam-energy loss was $`5`$ keV. In all cases the incident energy was adjusted so that the mean energy in the center of the target corresponded to $`E_{c.m.}=431.3\pm 0.8`$ keV. Here the quoted error includes the uncertainties in the incident energy and also the energy loss in the target. It should be noted that the use of thin targets is very important at low energies for minimizing energy loss and straggling effects. The differential cross section is particularly sensitive to the energy, a point which will be discussed quantitatively later in this paper.

### D Detection methods

The reaction products were detected using 100-$`\mu `$m-thick Si surface barrier detectors. Left-right symmetric detector configurations were utilized for all of the measurements. Count rates in the detectors were controlled by varying the solid angle (0.2-5 msr) and beam current (10-150 nA). In some cases, thin mylar foils (2-5 $`\mu `$m) were placed over the detectors to either (1) increase the energy separation between proton and deuteron peaks or (2) stop heavy recoils resulting from $`{}_{}{}^{12}\mathrm{C}(d,{}_{}{}^{12}\mathrm{C})d`$, $`{}_{}{}^{12}\mathrm{C}(d,{}_{}{}^{13}\mathrm{C})p`$, or $`{}_{}{}^{12}\mathrm{C}(d,{}_{}{}^{13}\mathrm{N})n`$ reactions. Measurements at far forward angles in the laboratory are limited by the high flux of elastic-scattered particles from carbon. For the proton-beam measurements in the high-voltage scattering chamber, detector signals were sent to ground potential via fiber-optic connections.

For the cross section, $`T_{20}`$, $`T_{21}`$, and $`A_{yy}`$ measurements the elastic scattering yields were determined from the energy spectra. Sample energy spectra are shown in Figs. 1-3. It should be noted that the low-energy heavy recoils seen in Figs. 2 and 3 can be eliminated if necessary by using the aforementioned mylar foils. The primary background under the p-d elastic scattering peaks at all angles is the low-energy tail from carbon elastic scattering. Pulsers were inserted into the detector electronics to facilitate dead-time corrections, which were typically less than 5% but always less than 10%. Peak yields were extracted using linear least-squares fits to the background on either side of the peak. For the analyzing power measurements, the background subtraction and dead-time corrections were performed separately for each spin state.

A different technique was used for the VAPs $`A_y`$ and $`iT_{11}`$ since the magnitudes of these observables are roughly a factor of 10 smaller than for the TAPs. The scattered deuterons and protons were detected in coincidence using two pairs of silicon surface barrier detectors placed at symmetric angles on either side of the incident beam. The angles of the detectors were set to observe either protons or deuterons in the more forward detectors in coincidence with deuterons or protons detected in the more backward detector on the opposite side of the beam. A major benefit of this method is that data can be acquired at a much higher rate because the events resulting from carbon elastic scattering are essentially completely eliminated. Note that a factor of 10 reduction in magnitude of the analyzing power requires a factor of 100 increase in the number of counts if the same relative statistical accuracy is desired. An additional advantage of the technique is that the background under the peak of interest is much reduced, lowering the possibility of false asymmetries resulting from the asymmetry in the background. Histograms of the time difference between the fast timing signals from each coincident pair of detectors were stored for each spin state. The time resolution for the coincident proton-deuteron peaks was $`8`$ ns, with backgrounds $`<3\%`$. A sample time spectrum is shown in Fig. 4. Dead-time corrections ($`<3\%`$) were determined by sending test pulses to the detector preamplifiers with time delays adjusted to give distinct peaks in the time spectra. Peak yields for each spin state were extracted using linear least-squares fits to the background on either side of the peak, and were corrected for dead time. The possible influence of channel (time) dependent dead-time corrections was calculated and found to be negligible.

### E Measurements

The incident beam energy was adjusted to yield $`E_{c.m.}=431.3`$ keV in center of the target, except for the $`iT_{11}(87.8^{})`$ measurements where the incident deuteron energy varied over $`1.3E_d6.0`$ MeV. The beams were collimated to produce a 2 mm (horizontal) $`\times `$ 4 mm (vertical) beamspot on the center of the target. The $`0^{}`$ position for each detector was determined by measuring the elastic scattering from carbon on either side of the beam for angles near $`0^{}`$. The systematic uncertainty in the angular positioning is estimated to be $`\pm 0.1^{}`$. The number of incident particles was determined by beam current integration. For the analyzing power data and relative cross section data the targets were replaced approximately every 12 hours, or when the hydrogen isotope content decreased by $`30`$%.

Relative cross sections were determined with a proton beam in the high-voltage scattering chamber. The measurements were performed using a movable pair of left-right symmetric detectors, and a fixed pair located out of the plane of the other detectors at $`\theta _{lab}=37.5^{}`$. The yield of protons and deuterons from p-d elastic scattering in the fixed pair was used to normalize the yields in the movable pair.

Absolute p-d cross sections were determined relative to the proton-proton elastic scattering cross section at a nominal proton beam energy of 2600 keV. Detector pairs were placed at $`\theta _{lab}=25^{}`$ and $`35^{}`$ where both cross sections vary smoothly with angle. These angles give d-p elastic scattering measurements at $`\theta _{c.m.}=82.6^{}`$, $`110^{}`$, and $`130^{}`$. A direct extraction negative ion source was used to inject beam into the tandem accelerator for these measurements. The beam could be switched between 1300-keV deuterons to 2600-keV protons in about one minute by only adjusting two parameters: (1) the magnetic field in the bending magnet between the source and the accelerator, and (2) changing the accelerator terminal voltage. Note that the magnetic field in the $`52^{}`$ analyzing magnet, which defines the beam energy, is unchanged throughout this procedure. The deuteron and proton beams were switched several times for each target, using an integrated charge of $`1`$ $`\mu `$C per measurement. For these measurements accurate relative beam current integration is essential. Current was integrated from the target rod (biased to $`+100`$ V) and a plate at the rear of the chamber (suppressed by a shroud at $`100`$ V). As a test of the technique, we also measured using the same procedure $`{}_{}{}^{197}\mathrm{Au}(d,d)`$ and $`{}_{}{}^{197}\mathrm{Au}(p,p)`$ elastic scattering at $`\theta _{\mathrm{lab}}=140^{}`$, where the cross sections are described to a very high accuracy by the Rutherford formula. A target consisting of $`2\times 10^{17}`$ Au atoms/cm2 was utilized for these measurements. The results agreed within 0.8% with the ratio expected for Rutherford scattering ($`1\%`$ corrections due to electron screening  were taken into account).

## III Data Analysis and Results

### A Cross section measurements

The relative cross section data as a function of center-of-mass angle are shown in Fig. 5. These data have been normalized to the absolute cross section measurements.

Absolute p-d cross sections were determined at $`\theta _{lab}=25^{}`$ and $`35^{}`$ relative to the p-p scattering measurements at a nominal proton beam energy of 2600 keV. The p-p cross section for these energies and angles appears to be well understood. Our procedure is to utilize the Nijmegen energy-dependent partial-wave analysis  to generate cross sections for our energy and angles. A thorough discussion of the data upon which this analysis is based is given in Ref. . The p-p cross section data for $`E_p<10`$ MeV which are included in the Nijmegen analysis generally have absolute uncertainties $`<0.5\%`$ and are reproduced by the fit within this error. In addition there are the p-p cross section data of Knecht, Dahl, and Messelt  which include measurements at $`E_p=2425`$ keV, very close to the energy of the present experiment. The quoted experimental uncertainties vary in the range 0.1-0.3%. Although not used in the Nijmegen analysis, the Nijmegen energy-dependent fit agrees with the results of Ref.  at 2425 keV within 0.4% for the angles of our interest. The systematic uncertainty in the absolute p-d cross sections is estimated to be $`\pm 1.1\%`$; the contributions to this uncertainty are summarized in Table I. The assumed p-p cross sections and the p-d cross sections derived from them are given in Table II.

### B Analyzing power measurements

Analyzing powers were determined from the particle yields measured for each spin state. The yields for each spin state were corrected for background and dead time, and were normalized by the number of incident particles. Tests for false asymmetries using tensor-polarized deuteron beams were carried out by measuring 197Au(d,d) scattering at $`\theta _{lab}=40^{}`$ under identical conditions as the p-d TAP measurements. For $`E_d=1.3`$ MeV, all of analyzing powers for 197Au(d,d) are expected to be $`<10^4`$ . The results were consistent with zero at the level of $`5\times 10^4`$, the statistical uncertainty of the measurement.

For each pair of left-right symmetric detectors, the analyzing powers were determined as follows. We first define the following ratios

$`L`$ $`=`$ $`Y_L^{(1)}/Y_L^{(2)}\text{ and}`$ (1)
$`R`$ $`=`$ $`Y_R^{(1)}/Y_R^{(2)},`$ (2)
where $`Y_L^{(i)}`$ and $`Y_R^{(i)}`$ are the particle yields in the left and right detector, respectively, while $`i=1,2`$ denotes the spin state. Then

$`T_{20}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[{\displaystyle \frac{L+R2}{p_{ZZ}^{(1)}\frac{1}{2}(L+R)p_{ZZ}^{(2)}}}\right],`$ (3)
$`T_{21}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}\left[{\displaystyle \frac{1L}{p_{ZZ}^{(1)}Lp_{ZZ}^{(2)}}}{\displaystyle \frac{1R}{p_{ZZ}^{(1)}Rp_{ZZ}^{(2)}}}\right],`$ (4)
$`A_{yy}`$ $`=`$ $`{\displaystyle \frac{L+R2}{p_{ZZ}^{(1)}\frac{1}{2}(L+R)p_{ZZ}^{(2)}}},`$ (5)
$`iT_{11}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}\left[{\displaystyle \frac{LR}{2p_Z^{(1)}(L+R)p_Z^{(2)}}}\right],\text{ and}`$ (6)
$`A_y`$ $`=`$ $`{\displaystyle \frac{2(L1)(R1)}{(p_Z^{(2)}p_Z^{(1)})(LR)}},`$ (7)
where $`p_Z^{(i)}`$ and $`p_{ZZ}^{(i)}`$ are the vector and tensor polarizations for spin state $`i`$. It should be noted that Eqs. (5) and (6) are approximate; for the analyzing powers and polarizations encountered in this experiment the corrections are negligible. Values for the TAP $`T_{22}`$ were derived from $`A_{yy}`$ and $`T_{20}`$ using $`T_{22}=A_{yy}/\sqrt{3}T_{20}/\sqrt{6}`$. The analyzing power data are shown in Figs. 6 and 7. The measured $`iT_{11}(87.8^{})`$ excitation function is shown in Fig. 8.

## IV Discussion

These experimental results are compared to calculations utilizing the Pair-Correlated Hyperspherical Harmonic basis  to construct the scattering wave function, and the Kohn variational principle to determine the scattering matrix elements . The calculations were performed using the AV18 NN potential  and with AV18 plus the 3N interaction of Urbana (UR-IX) , and are shown in Figs. 67, and 9. All of the cross sections in Fig. 9 have been divided by the cross section calculated using the AV18+UR-IX potential, so that the small differences are more apparent. The cross section measurements are seen to be in excellent agreement with the calculations using the AV18+UR-IX potential. The calculation using the AV18 potential alone differs from the measurements by 3-4%, a small but statistically-significant amount. As shown in Fig. 6, the TAP data are in reasonably good agreement with calculations using either AV18 or AV18+UR-IX. It is, however, significant to note that the AV18+UR-IX potential gives slightly better agreement for all of the TAPs.

The VAP results are shown Fig. 7. As we have previously reported , the calculations using both the AV18 and AV18+UR-IX potentials underpredict the data by $`40`$%. This discrepancy has now been observed for a wide range of energies for both $`nd`$ and $`pd`$ scattering. The origin of this “$`A_y`$ puzzle” is not clear. Specific efforts to study this problem have been undertaken recently by employing new forms of 3N forces . Moreover in Ref.  a new NN interaction constructed from chiral perturbation theory has been shown to give a better description of the VAPs at low energy.

The measured energy dependence of $`iT_{11}(87.8^{})`$, shown in Fig. 8, will have to be explained by any proposed solution to the $`A_y`$ puzzle. It is seen that the AV18+UR-IX potential consistently underpredicts $`iT_{11}`$ for all energies below the deuteron breakup threshold. Although not shown, the AV18 potential calculation lies $`10`$% lower than for AV18+UR-IX at these energies as well. We also note the the present measurements at $`E_{c.m.}=1.67`$ and 2.00 MeV are in excellent agreement with previous data .

The comparison between the present data and theoretical calculations at $`E_{c.m.}=431.3`$ keV can be made more quantitative through the use of the $`\chi ^2`$ parameter. For the $`N=94`$ data points the AV18 calculation yields $`\chi ^2=844`$, while the AV18+UR-IX calculation yields $`\chi ^2=268`$. The result that $`\chi ^2/N`$ is much greater than one even for the AV18+UR-IX calculation is expected since there are large discrepancies for the VAPs as well as smaller but statistically-significant differences in the other observables. A more detailed discussion of $`\chi ^2`$ comparisons and also phase-shift analysis by $`\chi ^2`$ minimization will be the subject of another publication. Here we comment on two systematic uncertainties in the experimental data and how they affect the calculated $`\chi ^2`$. Considerable attention was given to the determination of the energy of the experiment, $`E_{c.m.}=431.3\pm 0.8`$ keV. In Fig. 10 we show the sizable effect on $`\chi ^2`$ of varying the energy of the theoretical calculations with the AV18+UR-IX interaction. It is clearly important for the experimental energy to be determined as precisely as possible. We have also investigated the importance of the absolute normalization of the experimental cross sections. In Fig. 11 we show the effect on $`\chi ^2`$ of varying the normalization of the cross section. Here $`\chi ^2`$ is calculated for the 22 cross section points. Again very significant effects are observed.

## V Conclusions

We have measured differential cross sections and the analyzing powers $`A_y`$, $`iT_{11}`$, $`T_{20}`$, $`T_{21}`$, and $`T_{22}`$ at $`E_{c.m.}=431.3`$ keV, as well as an excitation function of $`iT_{11}(\theta _{c.m.}=87.8^{})`$ for $`431.3E_{c.m.}2000`$ keV. The 431.3-keV data comprise the lowest-energy complete set of cross section and analyzing power observables for nucleon-deuteron scattering. The cross section and TAPs are found to be in good agreement with calculations using the AV18 NN interaction and the Urbana-IX 3N interaction. The VAPs $`A_y`$ and $`iT_{11}`$ are found to be underpredicted by $`40`$%. At this time it is not clear whether these discrepancies result from inadequacies in the assumed NN interaction, the 3N interaction, or some other source. In this context, the very low energy data presented here represent a very stringent constraint for any new theoretical model of the nuclear interaction.

## Acknowledgments

The authors would like to thank B. J. Crowe and B. M. Fisher for their assistance in the data collection process, and M.C.M. Rentmeester for supplying the Nijmegen p-p scattering predictions. This work was supported in part by the U.S. Department of Energy, Office of High Energy and Nuclear Physics, under grant No. DE-FG02-97ER41041.


format: text

text:
# Ab initio description of tunnel junctions

## 1 Introduction

The discovery of Giant MagnetoResistance  in 1988 initiated a renewed interest in the phenomenon of spin dependent tunneling. The first time Julliere in 1975  succeded to measure the spin polarized current in a tunnel junction composed of two ferromagnetic layers separated by an insulating barrier. Using improved experimental techniques it is now possible to produce high quality tunnel junctions showing a remarkable tunnel magnetoresistance . A widely used model to interpret the results in terms of the polarization of the ferromagnetic leads was proposed by Julliere . It was based on earlier studies of Meservey and Tedrow on tunneling between superconductors and ferromagnets . The model assumed that the tunnel current is proportional to the product of the effective tunneling density of states of the considered spin direction for both ferromagnetic leads. In this paper we want to focus on the ab initio description of TMR of a Fe/Vacuum/Fe junction using a supercell approach. We will discuss the properties of the local and total density of states and the influence of the finite thickness of the electrode on current polarization and TMR ratio.

## 2 Transport

To evaluate the tunnel current in an Fe/Vacuum/Fe system we use a supercell description. We consider unit cells with a fixed number of empty sites to simulate the vacuum barrier with a thickness of 0.55 nm (4 ML) and a variing number of Fe sites describing the electrodes. Based on density functional theory the electronic structure of the system was calculated selfconsistently using a Screened Korringa-Kohn-Rostoker (KKR) method . This method is especially advantageous for prolonged supercells and we treated supercells with up to 108 atoms for this investigation.

We consider ballistic transport perpendicular to the layers in the limit of zero bias. Starting from the Landauer formula the conductance through a barrier between two reservoirs can be expressed by the following Fermi surface integral

$$G(𝐧)=e^2\frac{A}{2}\underset{\sigma ,k}{}\delta (E_k^\sigma E_F)|𝐧𝐯_k^\sigma |,$$
(1)
if we assume ideal transmission for all states at the Fermi energy $`E_F`$. $`A`$ denotes the cross section of the junction. $`𝐧𝐯_k^\sigma `$ is the component of the Fermi velocity in current direction. The sum in Eq.(1) goes over all states $`k`$ and over all spin directions $`\sigma `$, that means, we assume that the electron spin is conserved during the transmission process. The finger-print of magnetoelectronics is that the conductance changes as a function of the magnetic configuration from antiparrallel (AP) into parallel (P) under the influence of an external magnetic field. For this reason, we consider the electronic structure and the conductance of the superlattices in both magnetic configurations (see Fig. 1).

The TMR ratio is finally defined by

$$TMR=1\frac{G^{AP}}{G^P}.$$
(2)
## 3 Density of states and polarization

Julliere  assumed that the TMR ratio is proportional to the polarization of an effective density of states of the ferromagnetic leads. The question is which density is appropriate to satisfy this model. The polarization of any density of states is defined to be

$$P_n=\frac{n^{}n^{}}{n^{}+n^{}},$$
(3)
with $`n^\sigma `$ being a spin-projected local density of states depending on the position in the unit cell or the total density of states for the whole system. Results for the polarization of the local densities of states of an Fe/Vacuum/Fe system with a vacuum barrier of 1.1 nm (8 ML) are shown in Fig. 2.

In Fig. 2 the local density of states at the Fermi level at the different sites in the supercell is shown. The Fe polarization is positive and bulk-like inside the Fe layers and changes sign at the Fe interface due to the existence of surface resonances. Inside the vacuum barrier the polarization remains negative and tends to zero for larger distances from the surface.

## 4 TMR

Fig. 3a shows the corresponding conductance for P and AP configuration in dependence on the Fe layer thickness with a barrier consiting of 4 vacuum layers (0.55 nm). In the asymptotic region (for thick Fe layers) the conductance is a superposition of various oscillatory contributions around a constant value of 0.256 for P and 0.077 for AP configuration, respectively. All values of G are given in atomic units of the Rydberg system ($`G_0=1.75210^{14}\mathrm{\Omega }^1m^2`$). Characteristic oscillation periods are 2.3, 3.5 and about 7.5 ML. These periods are related to stationary points on the bulk Fe Fermi surface and describe wave vectors at which an remarkable amount of states at $`E_F`$ is quantized .

The Polarization PG of the conductance in P configuration (Fig. 3b) behaves non-monotonic. The mean value in the asymptotic region is about -0.2. This is consistent with the negative polarization of the local density of states at the Fe surface and in the vacuum barrier region. MacLaren et al.  obtained a positive polarization near 1 for a semi-infinite Fe/Vacuum/Fe junction. The difference should first be related to the structure of the barrier consisting of ASA potentials of constant height and secondly to the different geometry of the leads. In our calculation the potentials are calculated self-consistently which leads to different shapes of the barrier potential. The barrier profile, however, can influence the transmission coefficients drastically as it is known from free electron models.

The TMR ratio in dependence on the Fe layer thickness is shown in Fig. 3c. The horizontal line for thicknesses larger than 25 indicates the average of about 0.7. To elucidate the origin of the oscillations of conductance and TMR we investigated the total density of states of the periodic system at the Fermi level as a function of the Fe layer thickness (Fig. 4).

Obviously, their is no difference for the density of states in P and AP configuration of the junction, that is, TMR is not caused by changes in the density of states. The enhancement of the density of states due to surface resonances dominates the minority band for thin Fe layers. For thicker Fe layers the spin-projected density of states tends towards the Fe bulk values indicated by the horizontal lines. The variations of the total density for large thicknesses are very small and can not account for the oscillations of the conductances. That illustrates the fact that the conductance perpendicular to the layers is not caused by the huge number of states at the Fermi energy but by a few states with a large velocity in current direction. The relative position of these states with respect to the Fermi energy can change drastically with the Fe layer thickness due to quantum confinement and consequently causes the conductance and TMR oscillations.

## 5 Summary

The ballistic conductance was calculated for an Fe/Vacuum/Fe superlattice using the Landauer formula. Although we considered very large Fe layer thicknesses up to 50 ML (about 7 nm) the conductance and the TMR ratio still oscillates. The origin of the oscillations is the finite size of the Fe electrodes. The calculated mean value of the conductance polarization $`P_G=0.2`$ agrees with the polarization of the local density of states at the Fe surface. The calculated average of the TMR ratio is 0.7.


text:
# Moments of Inertia of Nuclei in the Rare Earth Region: A Relativistic versus Non-Relativistic Investigation

## I Introduction

One of the oldest problems in our understanding of collective motion of nuclei is the moments of inertia of ground state rotational bands in well deformed nuclei. They depend in a very sensitive way on the collective properties such as deformations and on pairing correlations of these many-body systems. Since rotational bands have been detected in nuclei nearly fifty years ago and since the first microscopic calculations of the moments of inertia by Inglis , these quantities have been used as a testing ground for nearly all microscopic theories of collective motion. They describe the response of the strongly interacting nuclear many-body system to an external Coriolis field breaking time reversal symmetry. They are, therefore, in some sense comparable to the static magnetic susceptibility in condensed matter physics.

The earliest microscopic calculations were based on a mean field of a deformed harmonic oscillator . In these calculations, residual interactions were neglected. In this way one founds the values of the moments of inertia identical to those of a rigid body with the same shape, in strong disagreement with the experimentally observed values, which where considerably smaller. It has been pointed out already very early  that residual two-body interactions would lower the values of the moment of inertia obtained in the Inglis model. The most important correlations causing such a reduction are pairing correlations . In fact, Belyaev  showed that a simple extension of the Inglis formula in the framework of the BCS theory is able to reduce the theoretical moments of inertia dramatically because of the large energy gap in the spectrum of quasiparticle excitations occurring in the denominator of the Belyaev formula. Therefore, the small moments of inertia of the rotational bands provided one of the most important experimental hints for a superfluid behaviour of these heavy open shell nuclei. Extended calculations using the theory of Belyaev have been carried out by Nilsson and Prior  using the BCS model based on the single particle spectrum of the Nilsson potential.

Apart from the fact that the results of these calculations were relatively successful, there are, as we know today, a number of open problems, namely:

i) Belyaev’s formula is based on generalized mean field theory violating essential symmetries. It was pointed out already by Migdal , that Galileian invariance is broken. He therefore modified the Belyaev formula by taking into account more complicated correlations to correct the violation of this symmetry. The question of the restoration of the broken Galileian invariance in the particle-hole and particle-particle channels has been later discussed in a number of articles, see for example Refs.  and references therein.

ii) Since Belyaev’s formula describes only quasiparticles moving independently, higher order correlations have to be taken into account. This has been done by Thouless and Valatin  who considered all orders of the interaction in a theory describing the linear response of the system to the external Coriolis field. Marshalek and Weneser  showed that the method of Thouless and Valatin preserves all the symmetries violated in the mean field approximation in linear order. In that sense Migdal’s formula was just a special case to deal with Galilean invariance. Marshalek showed in a series of papers (see for example Ref. ) that this is just the linear approximation of a more general theory based on Boson expansion techniques treating the symmetries appropriately in all orders .

iii) Much more elaborated versions of the cranked Nilsson model  showed that the $`𝒍^2`$-term in this model, which corrects in an elegant way the fact that realistic potentials for heavy nuclei are much flatter than an oscillator in the nuclear interior, introduces a strong spurious momentum dependence. This leads to the values for the moments of inertia deviating considerably from the experimental values. However, this problem is to large extent cured either by Strutinsky renormalization of the moments of inertia  or by an additional term to the cranked Nilsson potential that restores the local Galilean invariance .

Realistic applications of the Thouless-Valatin theory are by no means trivial. They should be based on self-consistent solutions of the mean field equations, because only for those solutions the RPA theory preserves the symmetries . In addition, they require the inversion of the RPA-matrix. Meyer-ter-Vehn et al.  have carried out such calculations in a restricted configuration space replacing the self-consistent mean field in an approximate way by the Woods-Saxon potential. As residual interaction they used density dependent Migdal forces $`F^\omega `$ in the $`ph`$-channel and $`F^\xi `$ in the $`pp`$-channel. These interactions have been carefully adjusted to experimental data for the underlying configuration space. The results of these calculations showed, that there are indeed effects originating from both channels, each of them modifying the Belyaev values, but cancelling themselves to a large extent. Therefore one could understand why older calculations  based on the generalized mean field model gave reasonable results as compared to the experiment.

Nowadays there are theories available where the Hartree-(Fock)-Bogoliubov equations can be solved in a fully self-consistent way in the rotating frame for finite angular velocity $`\Omega `$. Using the resulting wavefunctions $`|\Phi _\Omega `$ the Thouless-Valatin moment of inertia can be found as

$`J={\displaystyle \frac{d}{d\Omega }}\Phi _\Omega |\widehat{J}_x|\Phi _\Omega |_{\Omega =0}.`$ (1)
In this way one avoids the inversion of the full RPA-matrix, as task which is so far technically impossible for realistic forces in a full configuration space. Among these theories the properties of rotating nuclei are described in a way free from adjustable parameters only in the Cranked Relativistic Hartree-Bogoliubov (CRHB) theory  and non-relativistic density-dependent Cranked Hartree-Fock-Bogoliubov (CHFB) theory with finite range Gogny forces . Several realistic investigations of the moments of inertia in normal- and in super-deformed bands have been carried out in the literature in the framework of non-relativistic CHFB theory with Gogny forces . Similar investigations in the relativistic framework have been performed without pairing in the $`A60`$ , 80  and $`140150`$  regions of superdeformation where the pairing correlations are expected to be considerably quenched at high spin. The recently developed formalism of the CRHB theory has been applied so far only for the description of the moments of inertia in the $`A190`$ mass region of superdeformation . A very successful description of the moments of inertia has been obtained in the framework of these two theories. The aim of the present investigation is to find the similarities and differences between these two theories using in a systematic way the moments of inertia of rare-earth nuclei as a testing ground.

## II Theoretical tools

The CRHB theory  is an extension of cranked relativistic mean field theory  to the description of pairing correlations in rotating nuclei. It describes the nucleus as a system of Dirac nucleons which interact in a relativistic covariant manner through the exchange of virtual mesons : the isoscalar scalar $`\sigma `$ meson, the isoscalar vector $`\omega `$ meson, and the isovector vector $`\rho `$ meson. The phonon field $`(A)`$ accounts for the electromagnetic interaction. The CRHB equations for the fermions in the rotating frame are given in one-dimensional cranking approximation by

$`\left(\begin{array}{cc}h\mathrm{\Omega }_x\widehat{J}_x& \widehat{\mathrm{\Delta }}\\ \widehat{\mathrm{\Delta }}^{}& h^{}+\mathrm{\Omega }_x\widehat{J}_x^{}\end{array}\right)\left(\begin{array}{c}U_k\\ V_k\end{array}\right)=E_k\left(\begin{array}{c}U_k\\ V_k\end{array}\right)`$ (2)
where $`h=h_D\lambda `$ is the Dirac Hamiltonian $`h_D`$ for the nucleon with mass $`m`$

$`h_D=𝜶(i\mathbf{}𝑽(𝒓))+V_0(𝒓)+\beta (m+S(𝒓))`$ (3)
minus the chemical potential $`\lambda `$ defined from the average particle number constraint

$`\Phi _\Omega |\widehat{N}|\Phi _\Omega =N.`$ (4)
The Dirac Hamiltonian contains a repulsive vector potential $`V_0(𝒓)`$, an attractive scalar potential $`S(𝒓)`$ and the magnetic potential $`𝑽(𝒓)`$ which leads to non-vanishing currents in the systems with broken time-reversal symmetries . These currents play an extremely important role in the description of the moments of inertia  and thus they are taken into account fully self-consistently in the calculations. In Eq. (2), $`U_k`$ and $`V_k`$ are quasiparticle Dirac spinors and $`E_k`$ denotes the quasiparticle energies. Furthermore, $`\widehat{J}_x`$ and $`\mathrm{\Omega }_x`$ are the projection of total angular momentum on the rotation axis and the rotational frequency.

The time-independent inhomogeneous Klein-Gordon equations for the mesonic fields are given by

$`\left\{\mathrm{\Delta }(\Omega _x\widehat{L}_x)^2+m_\sigma ^2\right\}\sigma (𝒓)`$ $`=`$ $`g_\sigma \left[\rho _s^p(𝒓)+\rho _s^n(𝒓)\right]g_2\sigma ^2(𝒓)g_3\sigma ^3(𝒓),`$ (5)
$`\left\{\mathrm{\Delta }(\Omega _x\widehat{L}_x)^2+m_\omega ^2\right\}\omega _0(𝒓)`$ $`=`$ $`g_\omega \left[\rho _v^p(𝒓)+\rho _v^n(𝒓)\right],`$ (6)
$`\left\{\mathrm{\Delta }(\Omega _x(\widehat{L}_x+\widehat{S}_x))^2+m_\omega ^2\right\}𝝎(𝒓)`$ $`=`$ $`g_\omega \left[𝒋^p(𝒓)+𝒋^n(𝒓)\right],`$ (7)
with source terms involving the various nucleonic densities and currents

$`\rho _s^i(𝒓)={\displaystyle \underset{k>0}{}}[V_k^i(𝒓)]^{}\widehat{\beta }V_k^i(𝒓),`$ $`\rho _v^i(𝒓)={\displaystyle \underset{k>0}{}}[V_k^i(𝒓)]^{}V_k^i(𝒓)`$ (8)
$`𝒋^i(𝒓)`$ $`=`$ $`{\displaystyle \underset{k>0}{}}[V_k^i(𝒓)]^{}\widehat{𝜶}V_k^i(𝒓).`$ (9)
The sums over $`k>0`$ run over all quasiparticle states corresponding to positive energy single-particle states (no-sea approximation) and the indexes $`i`$ could be either $`n`$ (neutrons) or $`p`$ (protons). For simplicity, the equations for the $`\rho `$ meson and the Coulomb fields are omitted in Eqs. (7) since they have the structure similar to the equations for the $`\omega `$ meson, see Refs.  for details. Since the coupling constant of the electromagnetic interaction is small compared with the coupling constants of the meson fields, the Coriolis term for the Coulomb potential $`A_0(𝒓)`$ and the spatial components of the vector potential $`𝑨(𝒓)`$ are neglected in the calculations.

The pairing potential $`\mathrm{\Delta }`$ in Eq. (2) is given by

$`\mathrm{\Delta }\mathrm{\Delta }_{ab}={\displaystyle \frac{1}{2}}{\displaystyle \underset{cd}{}}V_{abcd}^{pp}\kappa _{cd}`$ (10)
where the indices $`a,b,\mathrm{}`$ contain the space coordinates $`𝒓`$ as well as the Dirac and isospin indices $`s`$ and $`t`$. It contains the pairing tensor $`\kappa `$

$`\kappa \kappa (𝒓,s,t,𝒓^{},s^{},t)={\displaystyle \underset{E_k>0}{}}V_k^{}(𝒓,s,t)U_k(𝒓^{},s^{},t)`$ (11)
and the matrix elements $`V_{abcd}^{pp}`$ of the effective interaction in the $`pp`$-channel. In the present version of CRHB theory, pairing correlations are only considered between the baryons, because pairing is a genuine non-relativistic effect, which plays a role only in the vicinity of the Fermi surface. The central part of the Gogny interaction containing two different finite range terms (see Eq. (14)) is employed in the $`pp`$ (pairing) channel.

The CRHB calculations have been performed with the NL1 parametrization  of the RMF Lagrangian. The D1S set of parameters  is used for the Gogny force in the pairing channel. The CRHB-equations are solved in the basis of an anisotropic three-dimensional harmonic oscillator in Cartesian coordinates. A basis deformation of $`\beta _0=0.3`$ has been used. All fermionic and bosonic states belonging to the shells up to $`N_F=13`$ and $`N_B=16`$ are taken into account in the diagonalisation and the matrix inversion, respectively. This truncation scheme provides reasonable numerical accuracy for the physical observables which as estimated in the calculations with larger fermionic basis is on the level of $`1.5\%`$ or better for kinematic moment of inertia $`J^{(1)}`$ and charge quadrupole moments $`Q_0`$. In order to calculate the derivative with respect to $`\mathrm{\Omega }`$ in Eq. (1), all CRHB calculations have been performed at rotational frequency $`\mathrm{\Omega }_x=0.05`$ MeV.

The starting point of the non-relativistic CHFB theory based on the Gogny force is the phenomenological finite range two-body interaction of the form

$`V^{pp}(1,2)`$ $`=`$ $`{\displaystyle \underset{i=1,2}{}}e^{[(𝒓_1𝒓_2)/\mu _i]^2}(W_i+B_iP^\sigma H_iP^\tau M_iP^\sigma P^\tau )+`$ (12)
$`+`$ $`iW_{LS}(\mathbf{}_{12}\delta (𝒓_1𝒓_2)\mathbf{}_{12})(𝝈_1+𝝈_2)`$ (13)
$`+`$ $`t_3(1+P^\sigma x_0)\delta (𝒓_1𝒓_2)[\rho (𝑹)]^{1/3}`$ (14)
which is used simultaneously both in $`pp`$ and $`ph`$ channels. In Eq. (14), $`𝑹=(𝒓_1+𝒓_2)/2`$. The transformation to the rotating frame  leads to equations similar to (2,10,11). The only difference is that the Dirac Hamiltonian of Eq. (3) is replaced by the non-relativistic Hartree-Fock Hamiltonian $`h_{ij}`$ containing the density-dependent Gogny force and the rearrangement term $`\mathrm{\Gamma }_{ij}`$, stemming from the density dependence of the force

$`h_{ij}h_{ij}+\mathrm{\Gamma }_{ij}=t_{ij}+{\displaystyle \underset{qq^{}}{}}\stackrel{~}{\upsilon }_{iqjq^{}}\rho _{q^{}q}+{\displaystyle \frac{\delta H^{}}{\delta \rho }}f_{ij}(𝑹)`$ (15)
In the above expression $`f_{ij}(𝒓)`$ is the quantity appearing in the second quantization form of the density operator $`\rho (𝒓)=_{ij}f_{ij}(𝒓)c_i^+c_j`$. The parameter set D1S  has been used in the present calculations. The CHFB-equations are again solved in the basis of an anisotropic three-dimensional harmonic oscillator in Cartesian coordinates with the oscillator length $`b_0=1.98`$ fm and the deformation of basis $`\beta _0=0.3`$. Only single-particle states satisfying the condition

$`\mathrm{}\omega _xn_x+\mathrm{}\omega _yn_y+\mathrm{}\omega _zn_zN_{max}\mathrm{}\omega _0\mathrm{with}N_{max}=11.1`$ (16)
have been included in the basis. The HFB equation has been solved with the Conjugated Gradient Method .

We also consider in this investigation the fluctuations in the pairing field by using the technique of an approximate particle number projection before the variation introduced by Lipkin and Nogami (further APNP(LN)) and discussed in detail in non-relativistic case in Refs. . In the relativistic case, the same approximate particle number projection is used but only the $`pp`$-part of the interaction is taken into account for the Lipkin-Nogami procedure, see Ref.  for details.

## III Results and discussion

In order to see the dependence of the results on proton and neutron numbers several nuclei in the Gd, Dy, Er and Yb isotope chains, the ground state rotational bands of which are close to rotational limit ($`E(4^+)/E(2^+)3.3`$), have been selected for the present study. The results of relativistic and non-relativistic calculations with and without APNP(LN) are presented in Tables 1-3 and Figs. 1-2 and compared with the experiment. Such quantities as charge quadrupole moments \[deformations\], pairing energies and moments of inertia are discussed below.

The calculated and experimental charge quadrupole moments $`Q`$ and quadrupole deformation parameters $`\beta `$ derived from $`Q`$ by

$`Q=\sqrt{{\displaystyle \frac{16\pi }{5}}}{\displaystyle \frac{3}{4\pi }}ZR_0^2\beta \mathrm{where}R_0=1.2A^{1/3}`$ (17)
are shown in Table 1 and Fig. 1. The general feature is that the charge quadrupole moments $`Q`$ obtained in relativistic calculations are larger than the ones of the non-relativistic calculations. In the non-relativistic case, the $`Q`$ values calculated with APNP(LN) are slightly smaller than the ones obtained without APNP(LN) because of the larger pairing correlations (see Table 2), which in general favours more spherical configurations. This trend also persists in the relativistic case, but there are the cases (164,166Er, 164,166Yb) in which APNP(LN) leads to larger charge quadrupole moments as compared with unprojected calculations. As shown in Fig. 1, the non-relativistic results are somewhat closer to the experiment than the relativistic ones. However, considering that the experimental values of $`Q`$ are subject of considerable experimental errors , one can conclude that both theories describe experimental charge quadrupole moments reasonable well which allows us to proceed further with the study of more sensitive quantities such as moments of inertia.

In Hartree-(Fock)-Bogoliubov calculations the size of the pairing correlations is usually measured in terms of the pairing energy defined as

$`E_{pair}={\displaystyle \frac{1}{2}}\text{Tr}(\mathrm{\Delta }\kappa ).`$ (18)
This is not an experimentally accessible quantity, but it is a measure for the size of the pairing correlations in the theoretical calculations. These quantities are shown in Table 2 and Fig. 2 for protons and neutrons separately. Both in relativistic and non-relativistic calculations, we observe that APNP(LN) leads to an increase of the pairing energies. This increase shows large variations as a function of the proton and neutron numbers. In general, this increase is larger in relativistic calculations. For example, proton pairing energies increase on average by factor 2.16 (with minimal and maximal increase being equal to 1.53 (156Dy) and 3.36 (160Gd)). In non-relativistic calculations, the average increase of proton pairing energies is only 1.55 (with minimal and maximal increase being equal to 1.11 (156Dy) and 1.97 (170Yb)). Neutron pairing energies behave in a similar way but there the difference between relativistic and non-relativistic calculations is smaller: the average increase of neutron pairing energies due to APNP(LN) is 1.81 in the relativistic and 1.73 in the non-relativistic calculations. In some cases, as for example in the Gd isotopes, the increase of neutron pairing energies due to APNP(LN) is larger in non-relativistic calculations. This increase of pairing energies due to APNP(LN) will lead to an increase of the pairing gaps, as it is well known from many phenomenological calculations using the monopole pairing force . We also see that with few exceptions the pairing energies are smaller in relativistic calculations. An additional effect of APNP(LN) is the increase of absolute values of binding energies. In relativistic calculations, APNP(LN) provides an additional binding by $`2.5`$ MeV.

Calculated moments of inertia are given in Table 3 and Fig. 1. Comparing the results of calculations without APNP(LN), it is clear that the moments of inertia are systematically larger in the relativistic case. Although one cannot completely exclude that this feature is to some extent connected with a different angular momentum content of single-particle orbitals in relativistic and non-relativistic calculations, a detailed analysis of pairing energies and moments of inertia suggests that this fact can be explained in a more realistic way by the different effective masses of the two theories: $`m^{}/m0.6`$ in RMF theory and $`0.7`$ in the non-relativistic theory. Thus the corresponding level density in the vicinity of the Fermi level is smaller in the relativistic theory which in general leads to weaker pairing correlations (see Table 2 and discussion above) as compared with non-relativistic calculations and as a result to larger moments of inertia. APNP(LN) restores to large extent the correct size of pairing correlations and thus its effect is larger in relativistic calculations. The average decrease of the moments of inertia due to APNP(LN) over the considered set of nuclei is 1.35 and 1.15 in relativistic and non-relativistic calculations, respectively. It is also clearly seen that APNP(LN) improves in average and especially in relativistic case the agreement between experimental and calculated moments of inertia. The level of agreement between calculations with APNP(LN) and experiment is similar in both theories, however, some discrepancies still remain.

## IV Conclusions

In conclusion, the moments of inertia, charge quadrupole moments and pairing energies of well-deformed nuclei in the rare-earth region have been investigated within relativistic and non-relativistic mean field theories with and without approximate particle number projection by means of the Lipkin-Nogami method. With no adjustable parameters it was possible to obtain good description of experimental data. It was found that the particle number projection plays a more important role in the relativistic calculations most likely reflecting the lower effective mass. In addition, it has larger impact on the moments of inertia and the pairing energies as compared with the charge quadrupole moments. Remaining deviations from experimental data could be related either to the parametrization of the mean field or to the interaction in the pairing channel or to the approximate character of the particle number projection. Further and more systematic investigations are needed for clarification of the main source of discrepancies between theory and experiment.

A.V.A. acknowledges support from the Alexander von Humboldt Foundation. This work has been supported in part by DGICyT, Spain under project PB97-0023 and the Bundesministerium für Bildung und Forschung under the project 06 TM 875. P.R. wishes to express his gratitude to the Spanish Ministry for Education for support of his work at the Universidad Autónoma de Madrid in the framework of the A. von Humboldt - J.C. Mutis program for the Scientific Cooperation between Spain and Germany.

## V Figure caption

Fig. 1. Experimental and calculated charge quadrupole moments $`Q`$ (top panels) and moments of inertia $`2J^{(1)}`$ (bottom panels) of well deformed rare-earth nuclei. The experimental data are shown by solid unlinked circles. The results of calculations with APNP(LN) are shown by the lines without symbols. The lines with open symbols are used to indicate the results of calculations without APNP(LN). Solid and dashed lines are used for relativistic and non-relativistic results, respectively.
Fig. 2. Calculated neutron (top panels) and proton (bottom panels) pairing energies. The results of calculations with (without) APNP(LN) are shown by the lines without (with) symbols. Solid and dashed lines are used for relativistic and non-relativistic results, respectively.


text:
# The Effect of Weak Interactions on the Ultra-Relativistic Bose-Einstein Condensation TemperatureImperial/TP/00-1/6, hep-ph/0011286, 11th September 2001

## 1 Introduction

The inclusion of finite densities of conserved charges in thermal field theory poses well known problems in the study of phase transitions. On one hand, perturbative analyses of charged scalar fields  give information about the phase structure, but cannot probe the critical point. Standard perturbative calculations are plagued with infrared divergences and after dealing with these, the asymptotic expansion breaks down. On the other hand, finite charges cannot be easily represented in lattice Monte-Carlo simulations. Encoding a non-zero charge in the grand canonical ensemble renders the action complex and therefore useless as a statistical weighting . These obstacles invoke a demand for alternative methods.

Several authors have recently considered the effect of repulsive interactions on the condensation temperature of a dilute Bose gas, a non-relativistic problem. It was pointed out by  that the leading correction could be isolated in the static Matsubara mode and a mean field calculation could be performed on this mode alone. Since then the static mode has been studied using the $`1/N`$ expansion , with the linear delta expansion , and now numerically .

In this paper we consider the effect of interactions on the transition temperature for scalars in the ultra-relativistic limit. By this we mean that we are at the high temperature limit, but that the typical self-energy corrections, $`\mathrm{\Sigma }`$, the chemical potential, $`\mu `$, or the cube root of the charge density, $`\rho `$, can be of the same order $`T\mathrm{\Sigma }\mu \rho ^{1/3}`$. This limit is appropriate for the study of high temperature symmetry breaking where Bose-Einstein condensation and spontaneous symmetry breaking have interesting similarities. The relevant finite temperature and density 4D Feynman diagrams are difficult to handle and only Jones and Parkin  include the setting-sun diagram for the self-energy. To avoid these problems, we use dimensional reduction  and this is the main focus of our paper. In our regime, we will show that dimensional reduction for a relativistic model is much more complicated than in the non-relativistic limit of our model which studied in , but we will show it is still manageable. As in the non-relativistic case, we take advantage of the good infrared behaviour of the non-static Matsubara modes. This allows us to integrate these modes perturbatively in an attempt to simplify our calculation.

Following dimensional reduction, the problem of relativistic Bose-Einstein condensation has been reduced to the study of a phase transition in an effective three-dimensional theory at zero temperature and zero density. Standard non-perturbative methods can then be employed to study this model but we do not investigate the many alternatives here, merely choosing one for exemplary purposes. Here we use the linear delta expansion for its simplicity and resemblance to standard perturbation theory.

By way of contrast, the only other study of a relativistic finite density system using dimensional reduction we know of is a study of QCD , rather than the Higgs sector studied here. The resulting 3D action is complex, making subsequent numerical analysis of  more complicated.

## 2 Dimensional Reduction

We begin by considering a relativistic system of bosons described by a complex scalar field theory. We encode a conserved charge by working in the grand canonical ensemble with chemical potential $`\mu `$. The partition function is then

$$Z=[d\mathrm{\Psi }^{}][d\mathrm{\Psi }]\mathrm{exp}\{S\}$$
(1)
where field integrations are periodic over imaginary time $`\beta =1/T`$ and the action is given by

$`S`$ $`=`$ $`{\displaystyle _0^\beta }𝑑\tau {\displaystyle d^3x\left\{[(_\tau +\mu )\mathrm{\Psi }^{}][(_\tau \mu )\mathrm{\Psi }]+\mathrm{\Psi }^{}\mathrm{\Psi }+m^2\mathrm{\Psi }^{}\mathrm{\Psi }+\lambda (\mathrm{\Psi }^{}\mathrm{\Psi })^2\right\}}.`$ (2)
The charge density is obtained from $`Z`$ as follows:

$$\rho =\frac{T}{V}\frac{1}{Z}\frac{Z}{\mu }$$
(3)
where we let the volume $`V`$ tend to infinity.

The periodicity of the fields is made explicit by a mode expansion and the non-static modes are integrated perturbatively111Alternatively, one can calculate static quantities in the 3D effective theory and in the full 4D theory. By matching these results one can then relate the coefficients of the 3D effective theory to those of the 4D one . However there is a lack of results in the literature for our model at high densities.. This not only gives an overall factor to $`Z`$ but also renormalizes the parameters of the static mode, $`\mathrm{\Phi }`$. The result is an effective three-dimensional theory whose dependence on the temperature and chemical potential is contained within its mass and coupling:

$$Z=\mathrm{exp}\{\beta VF(\mu )\}[d\mathrm{\Phi }^{}][d\mathrm{\Phi }]\mathrm{exp}\{S_{3\mathrm{D}}\}.$$
(4)
Writing $`F=F^{(0)}+\lambda F^{(1)}+𝒪(\lambda ^2)`$, the quadratic part of the action can be integrated to give

$$F^{(0)}(\mu )=T\underset{n0}{}\frac{d^3p}{(2\pi )^3}\mathrm{ln}\mathrm{\Delta }(\omega _n,p).$$
(5)
where $`\mathrm{\Delta }(\omega _n,p)=[(\omega _ni\mu )^2+p^2+m^2]^1`$ and, since we are dealing with bosons, $`\omega _n=2\pi nT`$. The leading perturbative correction is given by the figure-of-eight diagram:

$$F^{(1)}(\mu )=2T^2\left[\underset{n0}{}\frac{d^3p}{(2\pi )^3}\mathrm{\Delta }(\omega _n,p)\right]^2.$$
(6)
We use dimensional regularization in the $`\overline{\mathrm{MS}}`$ scheme, making the replacement

$$\frac{d^3p}{(2\pi )^3}_\mathrm{p}=\left(\frac{e^\gamma M^2}{4\pi }\right)^ϵ\frac{d^{32ϵ}p}{(2\pi )^{32ϵ}}$$
(7)
and subtracting only the terms which are divergent as $`ϵ0`$. $`M`$ is an arbitrary renormalization scale and $`\gamma `$ is the Euler-Mascheroni constant.

The ultra-relativistic limit is defined by $`\rho m^3`$ or, equivalently, $`Tm`$. In order to avoid confusion over the different expansion parameters we will set $`m=0`$. Since we shall handle the infrared region non-perturbatively, this does not create additional problems. We keep corrections of $`𝒪(\mu ^2/T^2)`$ as we will find these are in fact $`𝒪(\lambda )`$ for the critical theory. Ignoring any $`\mu `$-independent terms, the factor $`F`$ is then given by (see  for details)

$$F(\mu )=\frac{\mu ^2T^2}{6}\left[1\frac{\mu ^2}{4\pi ^2T^2}+\frac{\lambda }{4\pi ^2}+𝒪\left(\lambda ^2\right)\right],$$
(8)
where by $`𝒪(\lambda ^2)`$ we mean $`𝒪(\mu ^4/T^4,\lambda \mu ^2/T^2,\lambda ^2)`$.

The non-static modes also give corrections to the static mode action. We write the effective 3D action as

$$S_{3\mathrm{D}}=d^3x\left[\mathrm{\Phi }^{}\mathrm{\Phi }+r\mathrm{\Phi }^{}\mathrm{\Phi }+u(\mathrm{\Phi }^{}\mathrm{\Phi })^2\right],$$
(9)
where we neglect the higher dimensional operators $`(\mathrm{\Phi }^{}\mathrm{\Phi })^n`$ for $`n3`$. This is firstly on the grounds that their coefficients are high order in $`\lambda `$ and further that all non-renormalizable interactions are suppressed by factors of $`1/T^2`$ (for further discussion of this point see ). It can be checked that in the calculation which follows, these neglected terms contribute at higher order in the $`\lambda `$-expansion than is considered.

The fields should be scaled by $`\sqrt{T}`$ in order to remove an overall factor of $`1/T`$ from the action and to the order we work they undergo no field renormalization from the non-static mode integration.

Referring to figure 1, the 3D mass is given by

$$r=\mu ^2+\mathrm{\Sigma }_1+\mathrm{\Sigma }_2+\mathrm{\Sigma }_3+\mathrm{\Sigma }_4+\mathrm{\Sigma }_5+𝒪(\lambda ^3),$$
(10)
where

$`\mathrm{\Sigma }_1`$ $`=`$ $`4\lambda T{\displaystyle \underset{n0}{}}{\displaystyle _\mathrm{p}}\mathrm{\Delta }(\omega _n,p)`$ (11)
$`=`$ $`{\displaystyle \frac{\lambda T^2}{3}}\left[1{\displaystyle \frac{3\mu ^2}{2\pi ^2T^2}}+𝒪\left(\lambda ^2\right)\right];`$ (12)
$`\mathrm{\Sigma }_2`$ $`=`$ $`16\lambda ^2T^2{\displaystyle \underset{n_10}{}}{\displaystyle _{\mathrm{p}_1}}\mathrm{\Delta }(\omega _{n_1},p_1){\displaystyle \underset{n_20}{}}{\displaystyle _{\mathrm{p}_2}}\mathrm{\Delta }^2(\omega _{n_2},p_2)`$ (13)
$`=`$ $`{\displaystyle \frac{\lambda ^2T^2}{6\pi ^2}}\left[{\displaystyle \frac{1}{2ϵ}}+2\mathrm{ln}{\displaystyle \frac{Me^\gamma }{4\pi T}}+1\gamma +{\displaystyle \frac{\zeta ^/(1)}{\zeta (1)}}+𝒪\left(\lambda \right)\right].`$ (14)
Though $`\mathrm{\Sigma }_4`$ involves a static internal line it should be included since it will not arise from the effective static 3D theory222Note that at this order, integrating out the heavy modes actually leaves one with a non-local action with terms such as $`d^3xd^3x^{}|\mathrm{\Phi }(x)|^2B(x,x^{})|\mathrm{\Phi }(x^{})|^2`$ where $`B`$ is an $`O(\lambda ^2)`$ bubble diagram. It is when approximating the 3D theory by the purely local one (9) that contributions, such as $`\mathrm{\Sigma }_4`$, coming purely from non-local terms in the exact non-local effective theory must not be forgotten. For instance, the same type of non-local $`B`$ term should also lead to a diagram similar to $`\mathrm{\Sigma }_2`$ but with a light petal on the top and one has to check that this is of lower order than required. Such problems suggest that the matching of Green functions approach to dimensional reduction  might be simpler.. This contribution is in fact well behaved in the infrared and to leading order we find (see appendix)

$$\mathrm{\Sigma }_3+\mathrm{\Sigma }_4=0+𝒪(\lambda ^3T^2).$$
(15)
Finally we include the diagram with the one-loop vertex counterterm which is given by $`(5\lambda ^2T/8\pi ^2ϵ)(\mathrm{\Phi }^{}\mathrm{\Phi })^2`$. This is

$`\mathrm{\Sigma }_5`$ $`=`$ $`4{\displaystyle \frac{5\lambda ^2T}{8\pi ^2}}{\displaystyle \frac{1}{ϵ}}{\displaystyle \underset{n0}{}}{\displaystyle _\mathrm{p}}\mathrm{\Delta }(\omega _n,p)`$ (16)
$`=`$ $`{\displaystyle \frac{5\lambda ^2T^2}{12\pi ^2}}\left[{\displaystyle \frac{1}{2ϵ}}+\mathrm{ln}{\displaystyle \frac{Me^\gamma }{4\pi T}}+1\gamma +{\displaystyle \frac{\zeta ^/(1)}{\zeta (1)}}+𝒪\left(\lambda \right)\right].`$ (17)
We may now sum all the contributions to equation (10) giving

$$r=\mu ^2+\frac{\lambda T^2}{3}\left[1\frac{3\mu ^2}{2\pi ^2T^2}+\frac{3\lambda }{4\pi ^2}\left(\frac{1}{2ϵ}+\frac{1}{3}\mathrm{ln}\frac{Me^\gamma }{4\pi T}+1\gamma +\frac{\zeta ^/(1)}{\zeta (1)}\right)+𝒪\left(\lambda ^2\right)\right].$$
(18)
We shall find the coupling to be given with sufficient accuracy by

$$u=\lambda T+𝒪(\lambda ^2).$$
(19)
The coupling is now dimensionful due to the scaling of the static fields by $`\sqrt{T}`$.

Use of equations (3) and (4) with (8), (9), (18) and (19) gives

$`\rho `$ $`=`$ $`{\displaystyle \frac{F}{\mu }}+2\mu T\left[1+{\displaystyle \frac{\lambda }{2\pi }}+𝒪\left(\lambda ^2\right)\right]\mathrm{\Phi }^{}\mathrm{\Phi }`$ (20)
$`=`$ $`{\displaystyle \frac{\mu T^2}{3}}\left[1{\displaystyle \frac{\mu ^2}{2\pi ^2T^2}}+{\displaystyle \frac{\lambda }{4\pi ^2}}+𝒪\left(\lambda ^2\right)\right]+2\mu T\left[1+{\displaystyle \frac{\lambda }{2\pi ^2}}+𝒪\left(\lambda ^2\right)\right]\mathrm{\Phi }^{}\mathrm{\Phi },`$ (21)
where $`\mathrm{\Phi }^{}\mathrm{\Phi }`$ denotes the Green function evaluated in the effective 3D theory. This we cannot calculate perturbatively since the expansion will break down when probing large length scales greater than $`1/u`$.

## 3 Linear Delta Expansion

To evaluate the quantity $`\mathrm{\Phi }^{}\mathrm{\Phi }`$ we need a non-perturbative method as the 3D sector retains all the infrared divergences of the 4D theory. The effective 3D theory is studied at zero temperature and density so the problem is greatly simplified. Any standard non-perturbative method can be used at this point.

We will use LDE, the Linear Delta Expansion, though the method is also known by several other names (see  for a brief summary). LDE has been used successfully in many situations, including studies of scalar theories at non-zero density such as . In toy models, where exact results are achievable, LDE is known to produce convergent results and to do so much faster than alternatives, for instance see  and references therein. In full QFT, LDE has also often proved to be better than other methods .

We begin by defining $`S_\delta `$ which interpolates between our 3D action (from which we now drop the subscript) and some soluble action $`S_0`$ as $`\delta `$ varies from 1 to 0:

$$SS_\delta =\delta S+(1\delta )S_0.$$
(22)
We are free to choose

$$S_0=d^3x\left[\mathrm{\Phi }^{}\mathrm{\Phi }+\mathrm{\Omega }^2\mathrm{\Phi }^{}\mathrm{\Phi }\right]$$
(23)
such that

$$S_\delta =d^3x\left[\mathrm{\Phi }^{}\mathrm{\Phi }+(\mathrm{\Omega }^2\delta \mathrm{\Omega }^2+\delta r)\mathrm{\Phi }^{}\mathrm{\Phi }+\delta u(\mathrm{\Phi }^{}\mathrm{\Phi })^2\right].$$
(24)
Any physical quantity $`𝒫`$ is evaluated as a power series in $`\delta `$ to some finite order. This quantity will generally have some dependence on $`\mathrm{\Omega }`$ which we fix by some specified criterion which we take to be the principal of minimal sensitivity (PMS) :

$$\frac{d𝒫}{d\mathrm{\Omega }}|_{\delta =1,\mathrm{\Omega }=\overline{\mathrm{\Omega }}}=0.$$
(25)
This variational procedure allows for the emergence of non-perturbative behaviour.

The Green function we require can now be written as

$$\mathrm{\Phi }^{}\mathrm{\Phi }=\frac{d^3p}{(2\pi )^3}G_\delta (p),$$
(26)
where

$$G_\delta (p)=\left[p^2+\mathrm{\Omega }^2\delta \mathrm{\Omega }^2+\delta r+\mathrm{\Sigma }_\delta (p)\right]^1.$$
(27)
The criterion defining the transition temperature is that the correlation length in the original theory be infinite. This can be expressed as $`G_\delta ^1(0)|_{\delta =1}=0`$ which is satisfied by imposing

$$\delta r+\mathrm{\Sigma }_\delta (0)=0.$$
(28)
Use of this relation in equation (26) gives

$$\mathrm{\Phi }^{}\mathrm{\Phi }=\frac{d^3p}{(2\pi )^3}\left[p^2+\mathrm{\Omega }^2\delta \mathrm{\Omega }^2+\mathrm{\Sigma }_\delta (p)\mathrm{\Sigma }_\delta (0)\right]^1$$
(29)
and expanding to second order in $`\delta `$ we have

$$\mathrm{\Phi }^{}\mathrm{\Phi }=\frac{d^3p}{(2\pi )^3}\frac{1}{p^2+\mathrm{\Omega }^2}\left[1+\frac{\delta \mathrm{\Omega }^2}{p^2+\mathrm{\Omega }^2}+\frac{\delta ^2\mathrm{\Omega }^4}{(p^2+\mathrm{\Omega }^2)^2}\frac{\mathrm{\Sigma }_\delta (p)\mathrm{\Sigma }_\delta (0)}{p^2+\mathrm{\Omega }^2}\right]+𝒪(\delta ^3).$$
(30)
The only momentum dependent contribution to the self energy results from the setting sun diagram (see figure 2):

$$\mathrm{\Sigma }_3^{3\mathrm{D}}(p)=8\delta ^2u^2\frac{d^3k}{(2\pi )^3}\frac{d^3q}{(2\pi )^3}\frac{1}{(k^2+\mathrm{\Omega }^2)}\frac{1}{(q^2+\mathrm{\Omega }^2)}\frac{1}{((k+q+p)^2+\mathrm{\Omega }^2)}.$$
(31)
We continue our use of dimensional regularization in the $`\overline{\mathrm{MS}}`$ scheme, using the same scale $`M`$ as in our 4D heavy mode calculations. The required integrals (see, for example,  and ) are

$`{\displaystyle _\mathrm{p}}{\displaystyle \frac{1}{p^2+\mathrm{\Omega }^2}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{4\pi }}\left[1+2ϵ\left(\mathrm{ln}{\displaystyle \frac{M}{2\mathrm{\Omega }}}+1\right)+𝒪(ϵ^2)\right],`$ (32)
$`\mathrm{\Sigma }_3^{3\mathrm{D}}(0)`$ $`=`$ $`{\displaystyle \frac{8\delta ^2u^2}{(4\pi )^2}}\left[{\displaystyle \frac{1}{4ϵ}}+\mathrm{ln}{\displaystyle \frac{M}{3\mathrm{\Omega }}}+{\displaystyle \frac{1}{2}}+𝒪(ϵ)\right],`$ (33)
$`{\displaystyle _\mathrm{p}}{\displaystyle \frac{1}{p^2+\mathrm{\Omega }^2}}\mathrm{\Sigma }_3^{3\mathrm{D}}(p)`$ $`=`$ $`{\displaystyle \frac{2\delta ^2u^2}{(4\pi )^3\mathrm{\Omega }}}\left[{\displaystyle \frac{1}{2ϵ}}+3\mathrm{ln}{\displaystyle \frac{M}{2\mathrm{\Omega }}}+12\mathrm{ln}2+𝒪(ϵ)\right].`$ (34)
Further integrals may be derived from (32) by successive differentiations with respect to $`\mathrm{\Omega }`$.

Summing the contributions to the Green function, we find

$$\mathrm{\Phi }^{}\mathrm{\Phi }=\frac{\mathrm{\Omega }}{4\pi }+\frac{\delta \mathrm{\Omega }}{2(4\pi )}+\frac{\delta ^2\mathrm{\Omega }}{8(4\pi )}\frac{4\delta ^2u^2}{(4\pi )^3\mathrm{\Omega }}\mathrm{ln}\frac{4}{3}.$$
(35)
The divergences cancel and there is no need to invoke any counter terms. We now apply the PMS condition to $`\mathrm{\Phi }^{}\mathrm{\Phi }`$ and arrive at

$$\overline{\mathrm{\Omega }}=\pm \frac{u}{\pi }\left(\frac{2}{3}\mathrm{ln}\frac{4}{3}\right)^{\frac{1}{2}}.$$
(36)
The question of which solution to choose can be answered by comparison with the solution in the large-$`N`$ case for $`N2`$. Referring to  and  we find that the positive solution is appropriate, giving

$`\mathrm{\Phi }^{}\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{u\overline{f}}{(4\pi )^2}},`$ (37)
$`\overline{f}`$ $`=`$ $`\left(6\mathrm{ln}{\displaystyle \frac{4}{3}}\right)^{\frac{1}{2}}1.314.`$ (38)
Inserting this into equation (21) we have for the critical density

$$\rho =\frac{\mu T^2}{3}\left[1\frac{\mu ^2}{2\pi ^2T^2}+\frac{\lambda }{(4\pi )^2}\left(46\overline{f}\right)+𝒪\left(\lambda ^2\right)\right].$$
(39)
The chemical potential is an unwanted free variable in this expression and we must use relation (28) to constrain $`\mu `$ at the transition temperature. Along with the setting sun contribution to the self energy, we also have to order $`\delta ^2`$ (see figure 2)

$`\mathrm{\Sigma }_1^{3\mathrm{D}}`$ $`=`$ $`4\delta u{\displaystyle _\mathrm{p}}{\displaystyle \frac{1}{p^2+\mathrm{\Omega }^2}},`$ (40)
$`\mathrm{\Sigma }_2^{3\mathrm{D}}`$ $`=`$ $`16\delta ^2u^2{\displaystyle _\mathrm{p}}{\displaystyle \frac{1}{(p^2+\mathrm{\Omega }^2)^2}}{\displaystyle _\mathrm{q}}{\displaystyle \frac{1}{q^2+\mathrm{\Omega }^2}},`$ (41)
$`\mathrm{\Sigma }_4^{3\mathrm{D}}`$ $`=`$ $`4\delta ^2u(\mathrm{\Omega }^2r){\displaystyle _\mathrm{p}}{\displaystyle \frac{1}{(p^2+\mathrm{\Omega }^2)^2}}.`$ (42)
Evaluating and summing these contributions gives

$$\mathrm{\Sigma }_\delta (p=0)=\frac{\delta u\mathrm{\Omega }}{\pi }\frac{8\delta ^2u^2}{(4\pi )^2}\left[\frac{1}{4ϵ}+\mathrm{ln}\frac{M}{3\mathrm{\Omega }}\frac{1}{2}\right]+\frac{\delta ^2u\mathrm{\Omega }}{2\pi }\frac{\delta ^2ur}{2\pi \mathrm{\Omega }}+𝒪(\delta ^3),$$
(43)
which upon substituting into (28) gives

$$\delta r=\frac{\delta u\mathrm{\Omega }}{\pi }\frac{\delta ^2u\mathrm{\Omega }}{2\pi }+\frac{\delta ^2u^2}{2\pi ^2}\left[\frac{1}{4ϵ}+\mathrm{ln}\frac{M}{3\mathrm{\Omega }}+\frac{1}{2}\right]+𝒪(\delta ^3).$$
(44)
We apply the PMS condition to $`r`$, giving $`\overline{\mathrm{\Omega }}=u/\pi `$, and insert this value into equation (44):

$`r`$ $`=`$ $`{\displaystyle \frac{u^2}{2\pi ^2}}\left[{\displaystyle \frac{1}{4ϵ}}+\mathrm{ln}{\displaystyle \frac{M}{u}}+c_r\right],`$ (45)
$`c_r`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{\pi }{3}}+{\displaystyle \frac{3}{2}}1.546.`$ (46)
The dependence of $`r`$ on the scale $`M`$ is exact because of superrenormalisability and agrees with that noted elsewhere (e.g. ) even though it appears here in the context of a particular non-perturbative calculation. Comparing with equation (18) we have

$`\mu ^2+{\displaystyle \frac{\lambda T^2}{3}}\left[1{\displaystyle \frac{3\mu ^2}{2\pi ^2T^2}}+{\displaystyle \frac{3\lambda }{4\pi ^2}}\left({\displaystyle \frac{1}{2ϵ}}+{\displaystyle \frac{1}{3}}\mathrm{ln}{\displaystyle \frac{Me^\gamma }{4\pi T}}+1\gamma +{\displaystyle \frac{\zeta ^/(1)}{\zeta (1)}}\right)+𝒪(\lambda ^2)\right]`$
$`={\displaystyle \frac{\lambda ^2T^2}{2\pi ^2}}\left[{\displaystyle \frac{1}{4ϵ}}+\mathrm{ln}{\displaystyle \frac{M}{\lambda T}}+c_r\right].`$ (47)
The divergences cancel and we can obtain $`\mu `$ in terms of the critical temperature:

$`\mu `$ $`=`$ $`{\displaystyle \frac{\sqrt{\lambda }T}{\sqrt{3}}}\left[1+{\displaystyle \frac{\lambda }{8\pi ^2}}\left(\alpha _{\mathrm{mu}}\mathrm{ln}{\displaystyle \frac{T}{M}}+\eta _{\mathrm{mu}}\mathrm{ln}\lambda +a_{\mathrm{mu}}\right)+𝒪(\lambda ^2)\right],`$ (50)
$`\alpha _{\mathrm{mu}}=5,\eta _{\mathrm{mu}}=6,`$
$`a_{\mathrm{mu}}=\mathrm{ln}\left({\displaystyle \frac{e^\gamma }{4\pi }}\right)+1+3\left(\gamma +{\displaystyle \frac{\zeta (1,1)}{\zeta (1)}}\right)6c_r3.2706c_r6.007.`$
It is now clear that for the critical theory, $`\mu ^2/T^2𝒪(\lambda )`$. Substituting into the equation for the charge density gives

$`\rho `$ $`=`$ $`{\displaystyle \frac{\sqrt{\lambda }T^3}{3\sqrt{3}}}\left[1+{\displaystyle \frac{\lambda }{8\pi ^2}}\left(\alpha _{\mathrm{r}ho}\mathrm{ln}\left({\displaystyle \frac{T}{M}}\right)+\eta _{\mathrm{r}ho}\mathrm{ln}(\lambda )+a_{\mathrm{r}ho}\right)+𝒪(\lambda ^2)\right],`$ (53)
$`\alpha _{\mathrm{r}ho}=5,\eta _{\mathrm{r}ho}=6,`$
$`a_{\mathrm{r}ho}=a_{\mathrm{mu}}+{\displaystyle \frac{2}{3}}3\overline{f}3.92366c_r3\overline{f}21.759.`$
Equation (53) relates the critical density to the critical temperature for the Bose-Einstein condensation of a complex scalar field in the ultra-relativistic limit. Though (53) looks like an expansion in $`\lambda `$ we stress that the result is non-perturbative because of the severe IR problems in calculating $`\mathrm{\Phi }^{}\mathrm{\Phi }`$ and self-energies at the critical point in the three-dimensional theory. We should expect $`\rho `$ to have some renormalization scale dependence since we have included a one-loop vertex counterterm in the calculation.

## 4 Discussion

The leading term in our expressions for the critical chemical potential (50) or equivalently for the critical density (53), are the usual leading high temperature results, $`\mu _0=\sqrt{\lambda /3}T`$ and $`\rho _0=\sqrt{\lambda /27}T^3`$.

The $`\lambda ^2\mathrm{ln}(T/M)`$ term is exactly that expected from the running of the 4D coupling $`\lambda `$ in the leading term using perturbation theory where one finds

$$\lambda (M_2)=\lambda (M_1)+\frac{5\lambda ^2}{8\pi ^2}\mathrm{ln}\left(\frac{M_2^2}{M_1^2}\right).$$
(54)
The perturbative result is appropriate as the leading behaviour comes only from the heavy modes, and these are dealt with perturbatively in this calculation. Thus we find that once the implicit scale dependence of $`\lambda `$ is accounted for, our results for $`\mu `$ and $`\rho `$ are actually independent of the scale $`M`$ as all exact physics results should be.

The $`\lambda ^2\mathrm{ln}(\lambda )`$ term comes directly from our expression (45) for $`r`$ and in particular comes from dependence on the scale $`M`$ which is exact for for the super renormalisable theory. In the context of a Green function matching approach, as discussed for instance in , this term might be described as running the 3D mass $`r`$ from the scale $`MT`$, used when dealing with the heavy modes, down to $`u=\lambda T`$ appropriate for the 3D static theory. Overall then our expressions for the critical chemical potential and density agree with our expectations from other calculations.

In terms of actual numbers, for $`\lambda =1/8`$ the fractional corrections to this leading term are of the order of a few percent, e.g. in units of the scale $`M`$, for $`T=10.0`$, the critical chemical potential is $`\mu =\mu _0\times (10.010)`$, $`\mu _0=2.04`$ and the critical density is $`\rho =\rho _0\times (10.058)`$, $`\rho _0=68.0`$.

A good test of our central result, the dimensional reduction, is to compare our formulae for the critical chemical potential with that extracted from the results of Jones and Parkin . They also use the linear delta expansion but they apply it directly to the full four-dimensional theory anywhere in the symmetric phase. They can in principle study any temperature and density, though they limit their analysis to high temperatures. We use dimensional reduction so our method is always limited to high temperatures. However, this brings several benefits to us as further fields, including fermions, can be added with great ease and our non-perturbative effort is much less. The effective three-dimensional theory derived here can in principle be studied using any non-perturbative technique, including Monte-Carlo since the effective action is real. While we have only looked at the critical point, we have also given the critical density which is the actual measured quantity. Finally, our results are completely analytic, while Jones and Parkin can only find numerical solutions to their equations, though the numerics are relatively straight forward.

Turning to the details of the results, we find that the Jones-Parkin method gives the same qualitative behaviour as our results for $`\lambda =1/8`$, $`T/M=0.5\mathrm{}10.0`$. The fractional correction to the leading $`\mu `$ behaviour, $`\mathrm{\Delta }\mu :=\mu /\mu _01`$ is shown in figure 3. However their results for these parameters  are best fitted with slightly different coefficients333A shift from MS to $`\overline{\mathrm{MS}}`$ scales is needed., namely $`\alpha _{\mathrm{r}ho,JP}/\alpha _{\mathrm{r}ho}=1.2(1)`$, and choosing $`\eta _{\mathrm{r}ho,JP}=\eta _{\mathrm{r}ho}`$ we find $`a_{\mathrm{r}ho,JP}+\eta _{\mathrm{r}ho,JP}\mathrm{ln}\lambda =17.8(2)`$ compared with our value $`18.5`$. In our calculation the $`\alpha _{\mathrm{m}u}`$ value was set by perturbation theory and can not be exact. Since Jones and Parkin use an entirely non-perturbative method, even though both methods ought to be valid in the parameter range considered, there are likely to be small differences between the two results. These ought to be $`O(\lambda )`$ fractional corrections so the results seem to be consistent as far as they go.

Finally, we can illustrate the power of our work by noting that any non-perturbative calculation can be used with our formulae444Note that the normalisation of fields and the definition of $`u`$ or $`\lambda `$ often differ from those used here by simple constant factors. for dimensional reduction in the presence of a large density. For instance, since this work was completed, two lattice Monte Carlo studies of the 3D action relevant here have appeared . While these numerical results were produced for the non-relativistic 4D Bose-Einstein condensation problem, we can just as easily apply them to our relativistic case. The two studies are completely consistent but for definiteness we use the values in  which gives $`c_{\mathrm{r},lat}=0.671(1)`$ and $`\overline{f}_{\mathrm{l}at}=0.57(1)`$. In this case, the $`\alpha `$ and $`\eta `$ terms of both $`\mu `$ and $`\rho `$ expressions are fixed by dimensional reduction so the only difference is in the $`a`$ parameters, and we find that the lattice data leads to $`a_{\mathrm{m}u,lat}=0.756`$ and $`a_{\mathrm{r}ho,lat}=13.2`$ compared to our values of $`6.01`$ and $`21.8`$ respectively. While these appear to be large differences, in physical quantities such as the critical density $`\rho `$ the constant coming from $`a`$ is in fact overwhelmed by the contribution from the $`\mathrm{ln}\lambda `$ term for $`\lambda 1`$ which is where the dimensional reduction process is valid. The difference is more of a comment on the efficacy of different non-perturbative approximations in 3D calculations (see  for a good comparison) than particularly important to our results.

In conclusion, we have shown that by organising the modes into those which can be pertubatively integrated and those which cannot, we minimise the non-perturbative effort needed to study Bose-Einstein condensation at relativistic temperatures and densities. The method is economical and provides a reliable estimate of the critical temperature and critical density.

## Acknowledgements

We would like to thank H.F.Jones and P.Parkin for providing us with some of their data. DJB was supported by PPARC.

## Appendix

In this appendix we consider how the sunset diagram can be decomposed into contributions from different modes.

We begin by considering the $`\mu =0`$ case and state the results (see, for example, )

$`{\displaystyle \underset{n_1,n_2=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle _{\mathrm{p}_1,\mathrm{p}_2}}{\displaystyle \frac{1}{(\omega _{n_1}^2+p_1^2+m^2)}}{\displaystyle \frac{1}{(\omega _{n_2}^2+p_2^2+m^2)}}{\displaystyle \frac{1}{((\omega _{n_1}+\omega _{n_2})^2+(p_1+p_2)^2+m^2)}}`$
$`={\displaystyle \frac{1}{(4\pi )^2}}\left[{\displaystyle \frac{1}{4ϵ}}+\mathrm{ln}{\displaystyle \frac{M}{3m}}+{\displaystyle \frac{1}{2}}\right]+𝒪\left({\displaystyle \frac{m}{T}}\right)`$ (A.1)
and

$$_{\mathrm{p}_1,\mathrm{p}_2}\frac{1}{(p_1^2+m^2)}\frac{1}{(p_2^2+m^2)}\frac{1}{((p_1+p_2)^2+m^2)}=\frac{1}{(4\pi )^2}\left[\frac{1}{4ϵ}+\mathrm{ln}\frac{M}{3m}+\frac{1}{2}\right].$$
(A.2)
We shall split (A.1) into purely non-static internal lines ($`n_10;n_20;n_1n_2`$), one static internal line, ($`n_10;n_20;n_1=n_2`$ and $`n_10;n_2=0`$ and $`n_1=0;n_20`$), and purely static internal lines, ($`n_1=0;n_2=0`$). This covers all the possibilities and we may write

$$\underset{n_1,n_2=\mathrm{}}{\overset{\mathrm{}}{}}=\underset{\stackrel{n_1,n_20}{n_1n_2}}{}+6\delta _{n_1}\underset{n_2>0}{}+\delta _{n_1}\delta _{n_2}$$
(A.3)
Denoting these contributions to (A.1) as $`I_{\mathrm{nonstatic}}`$, $`I_{\mathrm{mixed}}`$ and $`I_{\mathrm{static}}`$, we immediately see that $`I_{\mathrm{static}}`$ is given by equation (A.2).

The contribution with one static line can be calculated by taking the static line to be massless. This does not cause infrared divergences and is appropriate for the critical theory. We may also set $`m=0`$ in the other propagators since $`Tm`$. We thus have

$`I_{\mathrm{mixed}}`$ $`=`$ $`6{\displaystyle \underset{n>0}{}}{\displaystyle _{\mathrm{p}_1,\mathrm{p}_2}}{\displaystyle \frac{1}{p_1^2}}{\displaystyle \frac{1}{p_2^2+\omega _n^2}}{\displaystyle \frac{1}{(p_1+p_2)^2+\omega _n^2}}+𝒪\left({\displaystyle \frac{m}{T}}\right)`$ (A.4)
$`=`$ $`6{\displaystyle \underset{n>0}{}}\left({\displaystyle \frac{e^\gamma M^2}{4\pi }}\right)^ϵ{\displaystyle \frac{\mathrm{\Gamma }^3(1/2+ϵ)}{(4\pi )^{3/2ϵ}\mathrm{\Gamma }(1+2ϵ)}}{\displaystyle _{\mathrm{p}_1}}{\displaystyle \frac{1}{p_1^2}}{\displaystyle \frac{1}{(p_1^2+\omega _n^2)^{1/2+ϵ}}}+𝒪\left({\displaystyle \frac{m}{T}}\right)`$ (A.5)
$`=`$ $`6{\displaystyle \underset{n>0}{}}\left({\displaystyle \frac{e^\gamma M^2}{4\pi \omega _n^2}}\right)^{2ϵ}{\displaystyle \frac{\mathrm{\Gamma }^2(1/2+ϵ)}{2ϵ(ϵ1/2)(4\pi )^{32ϵ}}}+𝒪\left({\displaystyle \frac{m}{T}}\right)`$ (A.6)
$`=`$ $`\left({\displaystyle \frac{e^\gamma M^2}{16\pi ^3T^2}}\right)^{2ϵ}{\displaystyle \frac{3\zeta (4ϵ)\mathrm{\Gamma }^2(1/2+ϵ)}{ϵ(ϵ1/2)(4\pi )^{32ϵ}}}+𝒪\left({\displaystyle \frac{m}{T}}\right)`$ (A.7)
$`=`$ $`{\displaystyle \frac{3}{(4\pi )^2}}\left[{\displaystyle \frac{1}{4ϵ}}+\mathrm{ln}{\displaystyle \frac{M}{2T}}+{\displaystyle \frac{1}{2}}\right]+𝒪\left({\displaystyle \frac{m}{T}}\right).`$ (A.8)
Finally, the purely non-static contribution can be found by subtracting the other contributions from (A.1):

$$I_{\mathrm{nonstatic}}=\frac{3}{(4\pi )^2}\left[\frac{1}{4ϵ}+\mathrm{ln}\frac{M}{2T}+\frac{1}{2}\right]+𝒪\left(\frac{m}{T}\right).$$
(A.9)
Turning to the case where $`\mu 0`$, $`I_{\mathrm{nonstatic}}`$ is unchanged at leading order in the high temperature expansion since $`\mu T`$ and we may clearly take the $`\mu 0`$ limit without causing infrared divergences. $`I_{\mathrm{mixed}}`$ is also unchanged upon choosing the static line to be critical. Up to corrections of $`𝒪(\mu ^2/T^2)`$ and $`𝒪(m/T)`$, the non-static and mixed contributions to the sunset diagram sum to zero.