format: text

text:
# Radiative alpha-capture cross sections from realistic nucleon-nucleon interactions and variational Monte Carlo wave functions

## I Introduction

Electromagnetic captures of alpha particles on $`{}_{}{}^{3}\mathrm{H}`$ and $`{}_{}{}^{3}\mathrm{He}`$ are important processes in astrophysics. Together, they are responsible for all 7Li production in the standard big bang nucleosynthesis (BBN) calculation. Because their cross sections are also difficult to measure in the laboratory at the relevant energies (20–500 keV), they are the major sources of uncertainty in the calculated primordial 7Li abundance . $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$ is also important for predicting the production rate of $`{}_{}{}^{8}\mathrm{B}`$ and 7Be neutrinos in the sun. Accurate knowledge of its rate at solar energies ($`20`$ keV) is therefore important for studies of the solar neutrino problems .

Recent theoretical studies of these cross sections have used two closely-related approaches. Potential models  treat alpha particles and tritons (or $`{}_{}{}^{3}\mathrm{He}`$ nuclei – the isospin symmetry of these systems makes the two processes almost identical) as point particles and the final states as bound states of these point particles. Authors of such models generate wave functions from potentials that fit experimental scattering phase shifts and bound state properties (binding energies, electromagnetic moments, etc.), and then compute the cross section as a direct capture. The resonating-group method (RGM)  is fully microscopic, in that it solves an explicitly seven-body problem with a nucleon-nucleon potential, the parameters of which are adjusted to reproduce bound-state and resonance properties for the particular problem being solved. The name is derived from the choice of basis states for solving the Schrödinger equation, which consists of one or more partitions of the particles into clusters with internal harmonic oscillator structure. The potential models now have a well-founded justification in the resonating-group work, in the form of the microscopic potential model . Although the RGM models have apparently been very successful in describing these reactions (the calculation of Kajino  correctly predicted both energy dependence and normalization for the experiment of Brune et al. ), it is not clear that this is the final word, because the agreement with the data seems to be spoiled when the model space is expanded .

It was shown in the early work of Christy and Duck  that low-energy radiative captures on light nuclei may be treated to good approximation as external direct captures, that is, as one-step processes in which most of the matrix element arises outside the nuclear interaction radius. It remains true in more detailed models that the largest contributions to the matrix elements arise in regions well outside the range of the internuclear forces. In principle, then, the cross section energy dependence is given by convolutions of (positive- and negative-energy) Coulomb wave functions with the current operator, while the normalization of the cross section is determined by the asymptotic normalization of the bound state in the appropriate clusterization channel.

However, fine details are obviously missing from the external direct capture model, motivating application of more detailed models that can predict both energy dependence (which is relatively easy to measure) and normalization (which is not easy to measure, and which must be inserted by hand in an external direct capture calculation after determination by some other means). Small effects involving the short-range ($`<5`$ fm) behavior of the nuclear wave functions, both in the alpha-trinucleon channel and in other channels, can affect the cross sections by several percent. Such effects are probably behind the differences between models pointed out in Ref. . In principle, there are also small corrections to the leading-order current operators. Because of the astrophysical importance of determining these cross sections, and especially the need for low-energy extrapolation of $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$ for solar physics, it is important to apply new approaches to this problem as they become available, and compare the results with past efforts.

In this context, recent developments in the physics of light nuclei are particularly interesting. There now exist “realistic” nucleon-nucleon potentials which describe the $`np`$ and $`pp`$ scattering data, as well as the deuteron, with high precision (e.g. Ref. ). Further interactions not describable by two-body potentials are described by three-nucleon potentials, which have been adjusted to reproduce residual effects in the energy spectrum of light nuclei and properties of nuclear matter . Wave functions have been developed for these potentials in systems with up to eight nucleons . This provides an opportunity to approach the problem of radiative captures on light nuclei using realistic potentials and the computational techniques that have been developed to utilize them.

Conversely, astrophysical interest in these processes has resulted in relatively precise measurements, which make cross section calculations useful tests of the wave functions in addition to the reproduction of the static moments, electron scattering properties, and energy spectra to which they have been compared in the past. We have already reported the application of these wave functions to a radiative capture calculation in a paper on the process $`d(\alpha ,\gamma )^6\mathrm{Li}`$ .

The remainder of this paper describes cross section calculations for the reactions $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ and $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$, using bound-state wave functions derived from realistic potentials by the variational Monte Carlo method. It is organized as follows: In Section II, we describe the wave functions used to compute the cross sections. In Section III, we describe the electromagnetic current operators and the methods used to compute their matrix elements. In Section IV, we describe the results for cross sections and branching ratios. In Section V, we examine the implications of our results.

## II Wave functions

### A Bound states

We used ground states of $`{}_{}{}^{3}\mathrm{H}`$, $`{}_{}{}^{3}\mathrm{He}`$, $`{}_{}{}^{4}\mathrm{He}`$, 7Li, and 7Be which were found by the variational Monte Carlo (VMC) technique for the Argonne $`v_{18}`$ two-nucleon potential (hereafter AV18)  and the Urbana IX three-nucleon potential (UIX) . The radiative captures can go to either the ground state or the (bound) first excited state in both 7Li and 7Be, so the first excited states of these nuclei were also needed. These wave functions were generated by the same VMC method as the ground states. The bound-state wave functions have been reported in Refs.  (triton and 4He) and  (modified here as in Ref.  to obtain 7Li and 7Be bound states with desired asymptotic properties).

The VMC method proceeds by constructing wave functions as products of pair and triplet correlations between nucleons, and adjusting free parameters in these correlations to minimize energy expectation values which are computed by a Monte Carlo integration. The bound state wave functions are built from central and operator correlations between nucleons, acting on a Jastrow wave function,

$$|\mathrm{\Psi }_T=\left[1+\underset{i<j<k}{}\stackrel{~}{U}_{ijk}^{TNI}\right]\left[𝒮\underset{i<j}{}(1+U_{ij})\right]|\mathrm{\Psi }_J,$$
(1)
where $`U_{ij}`$ and $`\stackrel{~}{U}_{ijk}^{TNI}`$ are two- and three-body correlation operators that include spin and isospin dependence and $`𝒮`$ is a symmetrization operator, needed because the $`U_{ij}`$ do not commute. The sums and products throughout this Section are over all nucleons. For 4He, the Jastrow part takes a relatively simple form:

$$|\mathrm{\Psi }_J=\underset{i<j<k4}{}f_{ijk}\underset{i<j4}{}f(r_{ij})|\mathrm{\Phi }_\alpha (0000)_{1234},$$
(2)
where $`f(r_{ij})`$ and $`f_{ijk}`$ are pair and triplet functions of relative position only, and $`\mathrm{\Phi }_\alpha (0000)`$ is a determinant in the spin-isospin space of the four particles. Jastrow wave functions for $`{}_{}{}^{3}\mathrm{He}`$ and for the triton are constructed analogously, but the parameters of the $`f(r_{ij})`$ for these nuclei have been chosen to minimize energy expectation values of the three-body nuclei rather than of the alpha particle. The triton and 3He are identical in our calculation, except for their isospin vectors. In cases where the distinction is unimportant, we refer to both nuclei as “the trinucleon,” and denote them both by $`\tau `$ in subscripts that label clusters.

For larger nuclei, spatial dependences must be introduced to place some particles in the p-shell. The $`A=7`$ Jastrow wave function is constructed from scalar correlations multiplying a shell model wave function,

$`|\mathrm{\Psi }_J`$ $`=`$ $`𝒜\{{\displaystyle \underset{i<j<k4}{}}f_{ijk}^{sss}{\displaystyle \underset{n4}{}}{\displaystyle \underset{5m<7}{}}{\displaystyle \underset{m<p7}{}}f_{nmp}^{spp}{\displaystyle \underset{i<j4}{}}f_{ss}(r_{ij})`$ (5)
$`{\displaystyle \underset{k4}{}}{\displaystyle \underset{5n7}{}}f_{sp}(r_{kn})f_{pp}(r_{56})f_{pp}(r_{57})f_{pp}(r_{67})`$
$`{\displaystyle \underset{LS[n]}{}}(\beta _{LS[n]}|\mathrm{\Phi }_7(LS[n]JMTT_3)_{1234:567})\},`$
where $`𝒜`$ is an antisymmetrization operator over all partitions of the seven particles into groups of four and three. For the central pair and triplet correlations, $`f_{xy}(r_{ij})`$ and $`f_{ijk}^{xyz}`$, the $`xyz`$ denote whether the particles are in the s- or p-shell. The shell model wave function $`|\mathrm{\Phi }_7(LS[n]JMTT_3)`$ has orbital angular momentum $`L`$, spin $`S`$, and spatial symmetry $`[n]`$ coupled to total angular momentum $`J`$, projection $`M`$, isospin $`T`$, and charge state $`T_3`$, and is explicitly written as

$`|\mathrm{\Phi }_7(LSJM[n]TT_3)_{1234:567}`$ $`=`$ $`|\mathrm{\Phi }_\alpha (0000)_{1234}\varphi _p^{LS[n]}(R_{\alpha 5})\varphi _p^{LS[n]}(R_{\alpha 6})\varphi _p^{LS[n]}(R_{\alpha 7})`$ (9)
$`\{[Y_{1m_l}(\mathrm{\Omega }_{\alpha 5})Y_{1m_l^{}}(\mathrm{\Omega }_{\alpha 6})Y_{1m_l^{\prime \prime }}(\mathrm{\Omega }_{\alpha 7})]_{LM_L}`$
$`\times [\chi _5({\displaystyle \frac{1}{2}}m_s)\chi _6({\displaystyle \frac{1}{2}}m_s^{})\chi _7({\displaystyle \frac{1}{2}}m_s^{\prime \prime })]_{SM_S}\}_{JM}`$
$`\times [\nu _5({\displaystyle \frac{1}{2}}t_3)\nu _6({\displaystyle \frac{1}{2}}t_3^{})\nu _7({\displaystyle \frac{1}{2}}t_3^{\prime \prime })]_{TT_3}.`$
The $`Y_{LM}(\mathrm{\Omega })`$ are spherical harmonics, $`\chi (s,m_s)`$ are spinors, and $`\nu (t,t_3)`$ are spinors in isospin, while brackets with subscripts denote angular momentum and isospin coupling.

Particles 1–4 are placed in the s-shell core with only spin-isospin degrees of freedom, while particles 5–7 are placed in p-wave orbitals $`\varphi _p^{LS[n]}(R_{\alpha k})`$ that are functions of the distance between the center of mass of the core and particle $`k`$. Different amplitudes $`\beta _{LS[n]}`$ in Eq. (5) are mixed to obtain an optimal wave function; for the $`J^\pi ,T=3/2^{},1/2`$ ground state of 7Be, the p-shell can have $`\beta _{1\frac{1}{2}[3]}`$, $`\beta _{1\frac{1}{2}[21]}`$, $`\beta _{1\frac{3}{2}[21]}`$, $`\beta _{2\frac{3}{2}[21]}`$, and $`\beta _{2\frac{1}{2}[21]}`$ terms. By far the largest contribution from these terms, as expected and as derived by diagonalization of the variational wave functions, is $`\beta _{1\frac{1}{2}[3]}`$. This is true of all the mass-seven, $`T=1/2`$ bound states, which are the final states of the radiative captures in question.

The two-body correlation operator $`U_{ij}`$ is defined as:

$$U_{ij}=\underset{p=2,6}{}\left[\underset{ki,j}{}f_{ijk}^p(𝐫_{ik},𝐫_{jk})\right]u_p(r_{ij})O_{ij}^p,$$
(10)
where the $`O_{ij}^{p=2,6}`$ = $`𝝉_i𝝉_j`$, $`𝝈_i𝝈_j`$, $`𝝈_i𝝈_j𝝉_i𝝉_j`$, $`S_{ij}`$, and $`S_{ij}𝝉_i𝝉_j`$, and the $`f_{ijk}^p`$ is an operator-independent three-body correlation. The six radial functions $`f_{ss}(r)`$ and $`u_{p=2,6}(r)`$ are obtained from two-body Euler-Lagrange equations with variational parameters as discussed in detail in Ref. . They are taken to be the same in the p-shell nuclei as in 4He, except that the $`u_{p=2,6}(r)`$ are forced to go to zero at large distance by multiplying in a cutoff factor, $`\left[1+\mathrm{exp}[R_u/a_u]\right]/\left[1+\mathrm{exp}[(rR_u)/a_u]\right]`$, with $`R_u`$ and $`a_u`$ as variational parameters. The $`f_{sp}`$ correlation is constructed to be similar to $`f_{ss}`$ for small separations, but goes smoothly to a constant of order unity at large distances ($`r>5`$ fm):

$$f_{sp}(r)=\left[a_{sp}+\frac{b_{sp}}{1+\mathrm{exp}[(rR_{sp})/a_{sp}]}\right]f_{ss}(r)+c_{sp}(1\mathrm{exp}[(r/d_{sp})^2]),$$
(11)
where $`a_{sp}`$, $`b_{sp}`$, etc. are additional variational parameters. The $`f_{pp}(r)`$ correlation in the mass-seven nuclei is the same as the $`f(r_{ij})`$ correlation in the trinucleon, so that when the three p-shell nucleons are all far from the s-shell core, they look very much like a trinucleon.

These choices for $`f_{ss}`$, $`f_{sp}`$, $`f_{pp}`$, and $`u_{p=2,6}`$ guarantee that when the three p-shell particles are all far from the s-shell core, the overall wave function factorizes as:

$$\mathrm{\Psi }_T[f_{sp}(r_{\alpha \tau })]^{12}[\varphi _p^{LS[n]}(r_{\alpha \tau })]^3\psi _\alpha \psi _\tau ,$$
(12)
where $`\psi _\alpha `$ is the variational 4He wave function, $`\psi _\tau `$ is the variational trinucleon wave function, and $`𝐫_{\alpha \tau }`$ denotes the separation between the centers of mass of the $`\alpha `$ and trinucleon clusters. Provided that $`[f_{sp}(r_{\alpha \tau })]^{12}`$ goes to a constant quickly enough and smoothly enough, the long-range correlation between clusters is proportional to $`[\varphi _p^{LS[n]}(r_{\alpha \tau })]^3`$.

The single-particle functions $`\varphi _p^{LS[n]}(R_{\alpha k})`$ describe correlations between the s-shell core and the p-shell nucleons, and have been taken in previous work  to be solutions of a radial Schrödinger equation for a Woods-Saxon potential and unit angular momentum, with energy and Woods-Saxon parameters determined variationally. It is important for low-energy radiative captures that these functions reproduce faithfully the large-separation behavior of the wave function, because the matrix elements receive large contributions at cluster separations greater than 10 fm. In fact, at 20 keV, more than 10% of the cross section for $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$ comes from cluster separations beyond 20 fm. At these distances, well outside the nuclear interaction distance, the clusterization with the lowest cluster-separation energy should be the most important. We have therefore modified the bound-state wave functions for the capture calculation to enforce cluster-like behavior, matching laboratory cluster separation energies, when the three p-shell nucleons are all far from the s-shell core.

In general, for light p-shell nuclei with an asymptotic two-cluster structure, such as $`\alpha d`$ in 6Li or $`\alpha t`$ in 7Li, we want the large separation behavior to be

$$[\varphi _p^{LS[n]}(r\mathrm{})]^nW_{km}(2\gamma r)/r,$$
(13)
where $`W_{km}(2\gamma r)`$ is the Whittaker function for bound-state wave functions in a Coulomb potential (see below) and $`n`$ is the number of p-shell nucleons. We achieve this by solving the equation

$$\left[\frac{\mathrm{}^2}{2\mu _{41}}\left(\frac{d^2}{dr^2}\frac{\mathrm{}(\mathrm{}+1)}{r^2}\right)+V(r)+\mathrm{\Lambda }(r)\right]r\varphi _p^{LS[n]}(r)=0,$$
(14)
with $`\mathrm{}=1`$, $`\mu _{41}`$ the reduced mass of one nucleon against four, and $`V(r)`$ a parametrized Woods-Saxon potential plus Coulomb term:

$$V(r)=\frac{V_0}{1+\mathrm{exp}[(rR_0)/a_0]}+\frac{2(Z2)}{n}\frac{e^2}{r}F(r).$$
(15)
Here $`V_0`$, $`R_0`$, and $`a_0`$ are variational parameters, $`(Z2)/n`$ is the average charge of a p-shell nucleon, and $`F(r)`$ is a form factor obtained by folding $`\alpha `$ and proton charge distributions together. The $`\mathrm{\Lambda }(r)`$ is a Lagrange multiplier that enforces the asymptotic behavior at large $`r`$, but is cut off at small $`r`$ by means of a variational parameter $`c_0`$:

$$\mathrm{\Lambda }(r)=\lambda (r)\left[1\mathrm{exp}\left((r/c_0)^2\right)\right].$$
(16)
The $`\lambda (r)`$ is given by

$$\lambda (r)=\frac{\mathrm{}^2}{2\mu _{41}}\left[\frac{1}{u_L}\frac{d^2u_L}{dr^2}\frac{2}{r^2}\right]\frac{2(Z2)}{n}\frac{e^2}{r},$$
(17)
where $`u_L`$ is directly related to the Whittaker function (solution to the radial Schrödinger equation for negative-energy states in a purely Coulomb potential):

$$u_L/r=(W_{km}(2\gamma r)/r)^{1/n}.$$
(18)
Here $`\gamma ^2=2\mu _{4n}B_{4n}/\mathrm{}^2`$, with $`\mu _{4n}`$ and $`B_{4n}`$ the appropriate two-cluster reduced mass and binding energy, $`k=2(Z2)e^2\mu _{4n}/\mathrm{}^2\gamma `$, and $`m=L+\frac{1}{2}`$.

For the 7Li ground state, the largest contribution has $`B_{43}=2.47`$ MeV (binding energy of 7Li relative to $`\alpha `$ and $`t`$ clusters) and $`L=1`$ corresponding to the asymptotic $`P`$-wave of the 7Li ground state, or amplitude $`\beta _{1\frac{1}{2}[3]}`$ in Eq. (5). None of the other possible amplitudes $`\beta _{LS[n]}`$ correspond to asymptotic $`\alpha t`$ clusterizations. However, there is no reason for them not to be present in compact configurations of the nucleons. Including such components in the wave functions improves the binding energies of the mass-seven bound states by about $`0.2`$ MeV. The asymptotic forms of $`\varphi _p^{LS[n]}(r)`$ in the lower-symmetry channels are set to match the threshold for $`{}_{}{}^{7}\mathrm{Li}^6\mathrm{Li}+n`$. Analogous descriptions hold for the other bound states (the 7Li excited state and the two 7Be bound states), with the appropriate thresholds substituted for $`B_{43}`$.

The 7Be ground and first excited state Jastrow functions have been treated in previous development of the variational Monte Carlo wave functions as the isospin rotations ($`T_3=+1/2`$ instead of $`1/2`$) of the corresponding 7Li shell-model-like wave functions. In this work, the 7Li and 7Be bound state Jastrow functions also differ by the choices of $`B_{4n}`$ for the asymptotic cluster behavior of the wave functions, which match the cluster breakup thresholds as described above in each case. This choice of $`B_{4n}`$ has the result that in configurations where the p-shell nucleons are far from the s-shell core, the energy is the sum of the Coulomb potential, the kinetic energy contributed by the $`\varphi _p^{LS[n]}`$, and the cluster binding energies (well-reproduced because the core resembles an alpha particle and the p-shell has been constructed to resemble a trinucleon). Because the $`\varphi _p^{LS[n]}`$ matches the laboratory binding energy for the known Coulomb potential, the local energies at large particle separations match the known binding energies in the mass-seven nuclei. This agreement has been confirmed numerically.

In Refs.  the $`f_{ijk}^{sss}`$ three-body correlation of Eq.(5) was a valuable and computationally inexpensive improvement to the trial function, but no $`f_{ijk}^{ssp}`$ or $`f_{ijk}^{spp}`$ correlations could be found that were of any benefit. However, for the types of wave functions used here, it is found that the correlation

$$f_{nij}^{spp}=1+q_1[f_{ss}(r_{ij})/f_{pp}(r_{ij})1]\mathrm{exp}[q_2(r_{ni}+r_{nj})],$$
(19)
with $`i,j`$ labels of p-shell nucleons and $`q_{1,2}`$ as variational parameters, is very useful . It effectively alters the central pair correlations between pairs of p-shell nucleons from their trinucleon-like forms to be more like the pair correlations within the s-shell when the two particles are close to the core. This correlation improves the binding energy by $`0.25`$ MeV in 7Li.

The authors of Refs.  reported energies both for the trial function $`\mathrm{\Psi }_T`$ of Eq.(1), and for more sophisticated variational wave functions, $`\mathrm{\Psi }_V`$, which add two-body spin-orbit and three-body spin- and isospin-dependent correlation operators. The $`\mathrm{\Psi }_V`$ gives improved binding compared to $`\mathrm{\Psi }_T`$ in both the mass-four and the mass-seven wave functions considered, but is significantly more expensive to construct because of the numerical derivatives required for the spin-orbit correlations. In the case of an energy calculation, the derivatives are also needed for the evaluation of $`L`$-dependent terms in AV18, so the cost is only a factor of two in computation. However, for the evaluation of other expectation values the relative cost increase is $`6A`$. Thus in the present work we choose to use $`\mathrm{\Psi }_T`$ for non-energy evaluations; this proved quite adequate in studies of 6Li form factors  and of the six-body radiative capture $`d(\alpha ,\gamma )^6\mathrm{Li}`$.

The variational Monte Carlo (VMC) energies and point proton RMS radii obtained with $`\mathrm{\Psi }_T`$ are shown in Table I along with the results of essentially exact Green’s function Monte Carlo (GFMC) calculations  and the experimental values. We note that the underbinding of the $`A=7`$ nuclei in the GFMC calculation arises from the AV18/UIX model and not the many-body method; it can be improved by the introduction of more sophisticated three-nucleon potentials .

Although the present variational trial functions with the imposed Coulomb asymptotic correlations produced a variational improvement in the case of 6Li, they give approximately the same energies as the older shell-model-like correlations of Refs.  for 7Li and 7Be. Unfortunately, because the variational energies of the $`A=7`$ bound states are not below those of separated alpha and trinucleon clusters, it is possible to lower the energy significantly by making the wave functions more diffuse. Therefore, the variational parameters were constrained to give RMS charge radii which agree reasonably with experiment, as seen in Table I. Given $`\varphi _p^{LS[n]}`$, the easiest way to constrain the RMS radius in the variational procedure is to choose the parameters of the $`f_{sp}`$ correlations between s- and p-shell nucleons to adjust the probability of finding the p-shell nucleons far from the s-shell core. There was considerable freedom in the specific form of these correlations as long as asymptotic properties of the wave functions were not being tested, because of the insensitivity of energy expectation values to the tails of the wave functions. However, the form of the correlation in the $`\alpha \tau `$ channel depends on $`f_{sp}^{12}`$, as seen in Eq. (12), so the large-cluster-separation parts of the wave function are very sensitive to the choice of $`f_{sp}`$. (The twelfth power arises because there are four particles in the s-shell and three in the p-shell, and thus $`4\times 3=12`$ sp pairs.) Prior to this work, VMC wave functions had cluster distributions that dropped by an extra factor of two beyond 5 fm, relative to the drop expected on the basis of the clusterization arguments given above. The extra drop had two closely-related effects: the asymptotic normalization coefficients (see below) for two-cluster breakup were too small, and so were the cross sections which we initially computed from them – by a factor of about two. The present wave functions perform more poorly than previous VMC wave functions on the ordering of the mass-seven bound states (a perennial difficulty for both VMC and GFMC because of the close spacing of the states), but they give larger asymptotic normalizations at a reasonable cost in binding energy. The relationship between nuclear size (as measured by quadrupole moments) and cross sections for direct radiative captures has been noticed before, and applied usefully both to these  and other reactions, most notably $`{}_{}{}^{7}\mathrm{Be}(p,\gamma )^8\mathrm{B}`$ . We note that our $`S(0)`$ and quadrupole moments for both mass-seven systems fit the general trends shown in Ref. .

The asymptotic two-cluster behavior of the seven-body wave functions can be studied by computing the two-cluster $`\alpha \tau `$ distribution functions, $`𝒜\psi _\alpha \psi _t^{m_t},𝐫_{\alpha t}\psi _{\mathrm{Li}}^{m_7}`$ for 7Li, and its analogs for the other seven-body bound states we consider. These functions are described in Ref. . They can be expressed in terms of Clebsch-Gordan factors, spherical harmonics, and the radial functions $`R(r_{\alpha \tau })`$ plotted in Figs. 1 and 2. At large $`r_{\alpha \tau }`$, $`r_{\alpha \tau }R(r_{\alpha \tau })`$ should be proportional to a Whittaker function, as described above. The proportionality constant is the asymptotic normalization constant $`C_1`$. We have extracted $`C_1`$ from the overlap functions by a least-squares procedure, matching them to the appropriate Whittaker functions. We find that for all the mass-seven wave functions considered here, $`C_1`$ becomes asymptotic at $`r_{\alpha \tau }`$ 7–9 fm. In Table II we present asymptotic normalization coefficients for the $`A=7`$ bound states, based on fitting overlaps in the region 7–10 fm as shown in Figs. 1 and 2.

Our values of $`C_1`$ for the mass-seven nuclei tend to be somewhat smaller than the values found in the literature . Only in the case of the 7Li ground state are experimental determinations of the $`C_1`$ of reasonable quality. Brune et al.  find a “world average” of $`3.55\pm 0.27\mathrm{fm}^{1/2}`$, relying mainly on theoretical models, and giving less weight to the partially experimentally-based evaluations of Igamov. For the other states, we rely on the Igamov et al.  extraction of ANCs from Kajino’s RGM calculations with the MHN potential. These calculations match the radiative capture data for $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ very well in both energy dependence and normalization. However, Igamov et al. report “nuclear vertex constants,” which differ from ANCs by prefactors whose definitions are ambiguous in the literature; the numbers presented in the second row of Table II should be used with caution. The present results for the ANCs are also subject to correlated uncertainties characteristic of the Monte Carlo integration method used to compute the overlaps, and the uncertainties are therefore difficult to estimate reliably (as discussed in more detail with regard to matrix element densities in Sec. III A below).

As noted in the Introduction, low-energy direct captures can be treated to good approximation by considering only cluster separations beyond a few fm, and using only the longer-range parts of the bound- and initial-state wave functions to compute matrix elements. This approach requires the provision of “spectroscopic factors,” or more precisely, the ANCs discussed above. Because we have imposed the condition that the large-cluster-separation part of the ground state match its expected form, it is true that for our cross section calculations, the VMC method provides $`C_1`$. However, it also provides the inner few fm of the wave function, which requires some model of what is going on inside the nuclear interaction radius, and may be important for understanding differences in the logarithmic derivatives of the $`S`$-factor found in various theoretical studies.

### B Scattering states

The initial-state wave functions are taken to be elastic-scattering states of the form

$$|\psi _{\alpha \tau };LSJM=𝒜\left\{\varphi _{\alpha \tau }^{JL}(r_{\alpha \tau })Y_{LM_L}(\widehat{𝐫}_{\alpha \tau })\underset{ij}{}G_{ij}|\psi _\alpha \psi _\tau ^{m_S}\right\}_{LSJM},$$
(20)
where curly braces indicate angular momentum coupling, $`𝒜`$ antisymmetrizes between clusters, $`\psi _\alpha `$ is the 4He ground state, and $`\psi _\tau ^{m_S}`$ is the trinucleon ground state in spin orientation $`m_S`$.

The $`G_{ij}`$ are identity operators if the nucleons $`i`$ and $`j`$ are in the same cluster. Otherwise, they are a set of both central and non-central pair correlation operators which introduce distortions in each cluster, under the influence of individual nucleons from the other cluster. They are derived from solutions for nucleon-nucleon correlations in nuclear matter , and become the identity operator at pair separations beyond about 2 fm. (These correlations have been included in the definition of the overlap functions shown in Figs. 1 and 2.)

The correlations $`\varphi _{\alpha \tau }^{JL}`$ are derived phenomenologically. The variational seven-body bound-state wave functions do not give the correct energies with respect to cluster breakup, so we do not expect to be able to use the variational technique to solve for these correlations. Instead, we generate the $`\varphi _{\alpha \tau }^{JL}`$ as solutions to Schrödinger equations, from cluster-cluster potentials that describe phase shifts of $`\alpha `$-$`{}_{}{}^{3}\mathrm{H}`$ and $`\alpha `$-$`{}_{}{}^{3}\mathrm{He}`$ scattering as scattering of point particles. Because of the small amount of available laboratory data and the large amount of work that has already been put into generating potentials that reproduce them, we take cluster-cluster potentials from the literature . We now point out the main features of these potentials.

We treat all of these models as (and many have been explicitly constructed as) descriptions of both the $`\alpha t`$ and $`\alpha ^3\mathrm{He}`$ systems, with appropriate Coulomb potentials and laboratory masses. Each of the potentials we use to generate the $`\varphi _{\alpha \tau }^{JL}`$ has a deep, attractive central term, and a spin-orbit term. The spin-orbit terms are constrained mainly by the spacing between the $`P_{3/2}`$ and $`P_{1/2}`$ bound states, and between the resonances at 2.16 and 4.21 MeV in $`F_{7/2}`$ and $`F_{5/2}`$ $`\alpha t`$ scattering, respectively (as well as their analogs in $`\alpha ^3\mathrm{He}`$ scattering). The other interesting feature in the scattering of the odd partial waves is the apparent hard-sphere behavior of the phase shifts in $`P`$-wave scattering. This comes about because the wave functions between the clusters must respect the Pauli principle by allowing antisymmetrization of nucleon wave functions between clusters. In the $`P`$ waves, this takes the form of wave functions that have a single node whose location is almost independent of scattering energy, and it therefore gives rise to phase shifts which look like scattering from a hard sphere . (See Figs. 1 and 2, where the steep dips in absolute values of the cluster distributions at 2 fm correspond to nodes). In the cluster-cluster potentials, this requires that the central term be large and negative, so that the ground state wave function has one node. Of course, a more tightly-bound “forbidden” state with zero nodes exists for such a potential, but it does not allow antisymmetrization of the nucleon wave functions. We note that the requirement of a particular nodal structure only models approximately the effects of the Pauli principle on the inter-cluster correlations. We did not use potentials that enforced the hard-sphere behavior with repulsive short-range terms .

Because of the requirements of the Pauli principle, one also expects different potentials to describe the odd- and even-parity scattering. The even-parity phase shifts are unfortunately lacking in details which models must match (see Fig. 3), beyond the apparent hard-sphere scattering in the $`S`$-wave data, corresponding to the Pauli-required minimum of two nodes in the wave function . The measured $`D`$-wave phase shifts are even worse, being consistent with zero (or $`180^{}`$, with the knowledge that there must be at least one node in the wave function) throughout the region below 8 MeV for both $`\alpha `$-t and $`\alpha `$-$`{}_{}{}^{3}\mathrm{He}`$ systems . This lack of features is particularly unfortunate in that by far the most important reaction mechanisms for radiative capture in this system are $`E1`$ captures from $`S`$\- and $`D`$-wave scattering states. The most stringent tests of any of the $`D`$-wave potentials have been comparison of their phase shifts with the results of more elaborate (RGM) theoretical models.

Finally, we note the difficulty of reproducing the published models, which is due to omitted descriptions of details of the potentials, particularly handling of the Coulomb potential at short range. Rather than guess how to fix up each potential, we restrict attention to potentials that allowed good descriptions of the low-energy scattering on the first try (with no short-range cutoff in the $`1/r`$ Coulomb term unless explicitly described by the potential’s authors). Relatively small differences in the cluster-cluster potentials are directly connected to the size of the scattering wave function at energies less than 1 MeV and cluster separations less than 20 fm, resulting directly in differences in the normalization of our computed capture cross sections from one potential to the next. However, it is possible to eliminate the worst potentials on the basis of the low-energy phase shifts; potentials which underpredict the phase shifts also produce radiative capture cross sections that are too low by as much as a factor of two, relative to those generated from other cluster-cluster potentials. Potentials which were created for use in orthogonal cluster models (as opposed to simple cluster models) were eliminated from application to our problem on this criterion.

## III Operators

The cross sections were computed by a multipole expansion of the electromagnetic current operator ,

$$\sigma (E_{\mathrm{c}.\mathrm{m}.})=\underset{LSJ_iJ_f\mathrm{}}{}\frac{8\pi }{2J_f+1}\frac{\alpha }{v_{\mathrm{rel}}}\frac{q}{1+q/m_7}\left[\left|E_{\mathrm{}}^{LSJ_iJ_f}(q)\right|^2+\left|M_{\mathrm{}}^{LSJ_iJ_f}(q)\right|^2\right],$$
(21)
where $`\alpha `$ is the fine structure constant ($`\alpha =e^2/\mathrm{}c`$), $`v_{\mathrm{rel}}`$ is the $`\alpha `$$`\tau `$ relative velocity, $`m_7`$ is the mass of the final state, and $`E_{\mathrm{}}^{LSJ_iJ_f}(q)`$ and $`M_{\mathrm{}}^{LSJ_iJ_f}(q)`$ are the reduced matrix elements (RMEs) of the electric and magnetic multipole operators with multipolarity $`\mathrm{}`$ connecting the scattering states in channel $`LSJ_i`$ to bound states of 7Li or 7Be with angular momentum $`J_f`$. The center-of-mass energy of the emitted photon is given by

$`q`$ $`=`$ $`m_7\left[1+\sqrt{1+{\displaystyle \frac{2}{m_7}}(m_\tau +m_\alpha m_7+E_{\mathrm{c}.\mathrm{m}.})}\right],`$ (22)
$``$ $`m_\alpha +m_\tau m_7+E_{\mathrm{c}.\mathrm{m}.}`$ (23)
where $`m_\tau `$, $`m_\alpha `$, and $`m_7`$ are the rest masses of the trinucleon, 4He, and the appropriate 7Li or 7Be state, respectively. The astrophysical $`S`$-factor is then related to the cross section via

$$S(E_{\mathrm{c}.\mathrm{m}.})=E_{\mathrm{c}.\mathrm{m}.}\sigma (E_{\mathrm{c}.\mathrm{m}.})\mathrm{exp}(2Z_1Z_2\pi \alpha /\mathrm{}v_{\mathrm{rel}}),$$
(24)
where $`Z_1`$ and $`Z_2`$ are the charges of the two initial-state nuclei.

The dominant reaction mechanism, for creation of both the excited state and the ground state, is an $`E1`$ (electric dipole) transition from an $`S`$-wave scattering state to a bound state. Capture into the excited state is followed immediately by electromagnetic decay to the ground state, so the total cross section of interest for astrophysics is the sum of the cross sections for captures into the two states. At energies above about 500 keV, $`E1`$ capture from $`D`$ waves becomes important. We computed transitions originating from scattering states with orbital angular momentum $`L=0,1,2`$ and 3, via $`M1,M2,E1,E2,`$ and $`E3`$ transitions, and found that up to the 0.1% level, only $`E1`$ captures originating from $`S`$\- and $`D`$-wave scattering states matter at energies below 1.5 MeV.

With the exception of the $`E1`$ term, all RMEs were computed in the standard long-wavelength approximation (LWA), keeping only the lowest-order term in photon wavenumber of the modified spherical Bessel functions appearing in the RME integrals. The LWA is valid to reasonable accuracy because at the low energies under consideration ($`E<1`$ MeV), the ratio of system size to photon wavelength is less than $`10\mathrm{fm}/200\mathrm{fm}=0.05`$.

In a previous study , we developed code to examine the isospin-forbidden $`E1`$ transition in $`d(\alpha ,\gamma )^6\mathrm{Li}`$. We have applied this code to compute corrections to the LWA for the $`E1`$ transitions under consideration here. An examination of these corrections for the mass-seven system is in principle of interest for the problem of extrapolating cross sections to low energies. However, we find that all but one of them provide contributions of less than 0.05% of the total cross section. (See Ref.  for a list and detailed discussion of the corrections we applied, which extend to third order in $`q`$.) We do not actually compute the largest correction, which is the “center-of-energy” correction. This correction arises because potentials and kinetic energies should be included in the definition of the center-of-momentum frame, but have not been. The center-of-energy correction becomes important when the leading order LWA operator vanishes, as in $`d(\alpha ,\gamma )^6\mathrm{Li}`$, but it should amount to only +2.4% of the $`E1`$ cross section in $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ and +3.1% of the $`E1`$ cross section in $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$. We did not compute center-of-energy corrections (or include an estimate of them in the results presented below) because of the extra computation necessary to find energies during the capture calculation. Their omission is not serious because 1) their effect is to change the normalization, not the energy dependence of the $`E1`$ cross section, 2) our model calculation is only accurate to about 5–10% at best, and 3) the above estimates of the size of the effect should be quite accurate, being based only on the differences between using nuclear masses and using integer multiples of the mean nucleon mass in the factor $`[(Z_2m_1Z_1m_2)/(m_1+m_2)]^2`$ in the LWA cross section.

### A Matrix element integration

Actual computation of the matrix elements was performed with a modified version of the code described in Ref. , which is itself a modified version of a code developed to compute energies and other properties of light nuclei  for variational calculations. The method used to integrate over nucleon configurations is the Metropolis Monte Carlo algorithm, with a weight function proportional to the bound state wave function involved in the computed transition. As discussed below, this weight function was chosen to reflect in general detail the form of the matrix element integrands, and to obtain significant numbers of Monte Carlo samples over a broad range of cluster separations. The final calculation consisted of $`10^6`$ samples for each transition.

We have applied the approach of splitting the calculation into energy-dependent and energy-independent parts, as in Ref. , so that the reduced matrix elements of Eq. (21) are written

$`T_{\mathrm{}}^{LSJ_iJ_f}(q)`$ $`={\displaystyle \frac{\sqrt{2J_f+1}}{<J_iM_i,\mathrm{}\lambda |J_fm_f>}}`$ (26)
$`\times {\displaystyle _0^{\mathrm{}}}dxx^2\varphi _{\alpha \tau }^{J_iL}(x)\psi _7^{J_fm_f}|T_\mathrm{}\lambda (q)𝒜\{\delta (xr_{\alpha \tau })Y_L^{M_L}(\widehat{𝐫}_{\alpha \tau })){\displaystyle \underset{ij}{}}G_{ij}|\psi _\alpha \psi _\tau ^{m_S}\}_{LSJ_iM_i}`$
and computed using photon polarization $`\lambda =+1`$ for the multipole operators $`T_\mathrm{}\lambda (q)`$. The delta function is applied by accumulating the Monte Carlo integral in radial bins of thickness $`0.25`$ fm. The final integration over $`x`$ is performed by inserting the appropriate dependences on photon energy in each term (since the LWA expansion is in powers of energy, this dependence may be taken out of the integral), and computing $`\varphi _{\alpha \tau }^{JL}`$, at each energy. This allows the time-consuming Monte Carlo integration to be performed only once for each partial wave and operator, so that computation of RMEs for many energies is relatively inexpensive in computer time. After initially setting up the code and checking that selection rules were satisfied, we did not explicitly compute RMEs for parity-forbidden operators.

At cluster separations beyond about 10 fm, the RME densities were subject to considerable noise in the Monte Carlo sampling. This is because while Monte Carlo weighting schemes based on the ground state give good sampling along directions other than the cluster separation in the $`3A`$-dimensional configuration space, they provide small numbers of samples at large $`r_{\alpha \tau }`$ (due to exponential decay of the wave function at large distances). VMC work usually uses weighting based on the inner product in spin-isospin space of a simplified bound state with itself. Such weighting is good for minimizing Monte Carlo variance of integrands that resemble the bound state probability density closely, such as energies, but it does not provide enough samples in the asymptotic region of the wave function to do low-energy direct capture calculations. Our previous paper  used the square root of the usual weight function, extending the Monte Carlo sampling out to large cluster separations, but at the expense of greater sampling noise for a fixed number of samples at fixed cluster separation. In the mass-six problem, it was straightforward to run the Monte Carlo integration of the RMEs until the densities had “converged” to the expected asymptotic forms at large cluster separation. This required about $`2\times 10^6`$ total configurations for a given RME. At $`A=7`$, the spin-isospin space is larger, so the code is slower by a factor of about 10, and it was only practical to obtain $`10^6`$ samples for each RME with available computing resources. However, since the configuration space gains three dimensions with the addition of a particle, $`10^6`$ samples do not provide as thorough of a sample of the configuration space at mass 7 as at mass 6. The result is at best a few-percent measure of the asymptotic normalization at 7–10 fm. Many samples are obtained at larger cluster separation with the new weighting scheme, but the samples beyond 15 fm appear to have correlated noise. This is presumably because the samples in question correspond to only a few excursions of the sampling Markov chain into the region of large cluster separation, and are therefore not very independent from each other. They tend to be either mostly high or mostly low relative to the asymptotic forms explicitly built into the wave functions. See Fig. 4 for illustration of these difficulties. A much larger sample (by a factor of 10) would be expected to exhibit much less of this sort of correlated noise, but is prohibited by the large amount of computer time that would require.

The matrix elements were therefore integrated out to 7 fm using the Monte Carlo results for the integrand of Eq. (26); integration beyond 7 fm was carried out using the known asymptotic forms of the matrix element densities, normalized to the Monte Carlo output by least-squares fitting at 7–10 fm. This range was arrived at by comparing results for the ANC of the 7Li ground state arising from two different weighting schemes and varying numbers of samples. Reasonably consistent agreement was found by fitting at 7–10 fm in all cases. The accuracy in the cross section then depends on how accurately the ANCs can be determined in this region.

## IV Cross sections

### A $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$

The computed $`S`$-factor for $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ is shown broken down into contributions from various terms of Eq. (21) in Fig. 5, and in comparison with laboratory data in Fig. 6. The dominant processes are obviously $`E1`$ captures, with large contributions from captures into both the ground and excited states. Contributions from other partial waves and multipole operators are not present above the 1% level. Contributions from higher-order LWA corrections to the $`E1`$ operator are less than 0.02%. Our calculations are therefore limited in accuracy only by 1) the accuracy of the (bound- and scattering-state) wave functions, and 2) Monte Carlo statistics.

The present calculation is the result of sampling $`10^6`$ points in the seven-particle configuration space. Formal error estimation on the resulting cross sections is difficult because all 35 partitions of the seven nucleons into alpha and triton clusters were computed for each configuration. This saves computation time, and enforces antisymmetry of the initial state exactly; it also introduces correlations between values of the operator densities (integrands of Eq. (26)) at different cluster separations. The uncertainties in different transitions and partial waves are also correlated, because they are based on the same random walk of particle configurations. There is also an uncertainty from the wave function normalizations, because they are also derived from Monte Carlo integrations. This last uncertainty amounts to about 3% in the cross section.

We take the best indication of the Monte Carlo uncertainty to be the formal uncertainties on the asymptotic normalizations of the Monte Carlo matrix element densities, which we used to compute matrix element contributions beyond cluster separations of 7 fm. The asymptotic forms of the matrix element densities are the Whittaker-function asymptotic forms discussed above, multiplied by the radial dependences of the electromagnetic multipole operators. The normalizations were fitted to the matrix element densities at cluster separations of 7–10 fm. Using fitted asymptotic forms amounts to treating a large part of the matrix element as an external direct capture, and it fixes two problems. First, it removes the need for large numbers of Monte Carlo samples in the remote tails of the wave function. Second, the asymptotic normalizations are found by a weighted least-squares procedure, and formal error estimates on these normalizations are possible. In practice, the correlations between matrix element densities produce reduced $`\chi ^2`$ significantly less than unity. A common approach when confronted with such a problem in experimental data is to assume that uncertainties have been overestimated, and to reduce formal uncertainties accordingly. We have not done this, and we arrive at uncertainties of approximately 10% in the cross section, based on the formal error estimates (corresponding closely to the sizes of ANC errors in Table II). More detailed analysis is problematic, and it is not called for because the 10% estimate is already larger than any other contribution to the error budget.

The results themselves are best characterized in three ways:

1) Normalization. Our results are lower than the data of Brune et al.  by 0–20%, depending on the cluster potential. Although other data exist , those of Brune et al. are much more precise, and permit the best test of our results. Those data share a common 6% normalization uncertainty, not shown in Fig 6. The systematic discrepancy between some of our results and the data is probably small enough to ascribe to the combined uncertainties of the Monte Carlo integration and of the data. Taking the branching ratio to have its experimental value of $`R=0.453`$, the computed $`S`$-factors for transitions to the ground state match the measured low-energy $`S`$-factors, while those for the excited state do not. The normalization is affected significantly by the choice of potential used to generate the inter-cluster correlations $`\varphi _{\alpha t}^{JL}`$ for the scattering states. By applying five different potentials from the literature as described above, we obtain a variation of $`\pm 5\%`$ in the $`S(0)`$ (total range for the five potentials) about a mean of 0.90 keV$``$b – a full range equal in size to the Monte Carlo uncertainty. Summarizing our results in one number, we obtain using the inter-cluster potential A of Kim et al. (the best fit to the $`S`$-wave phase shifts), $`S(0)=0.095`$ keV b, similar to other estimates found in the literature.

2) Branching ratio. The branching ratio $`R`$, defined as the ratio of the cross section for capture into the excited state to that for the ground state, is shown in comparison with the Brune data in Fig. 7. A weighted least-squares normalization of the calculation to the laboratory data shows that our calculation of $`R`$ is lower than the data by 15%, which can be combined with the results above to infer that the excited-state ANC is low by 8%, within the range of the Monte Carlo sampling error. The energy dependence of the branching ratio matches the data with a $`\chi ^2`$ of 26.0 for 15 degrees of freedom, about as well as the straight line fits of Brune et al.

3) Energy dependence. The energy dependence of the $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ $`S`$-factor at low energy is almost completely independent of cluster-cluster potential, for the five potentials examined. After normalizing the computed cross section (with cluster correlations computed from potential A of Kim et al. ) to match the Brune data, we obtain a $`\chi ^2`$ of 38.7 for 16 degrees of freedom. A significantly better fit results if the highest 3 points in energy (where $`D`$-wave capture becomes important) are excluded. The residuals of $`1.8\%`$ to $`+8.5\%`$ compare favorably with other theoretical calculations, but the energy dependence of the calculation is systematically shallower than that of the data at the highest-energy points. Our calculation of the $`S`$ factor gives a logarithmic derivative at $`E=0`$ of $`0.972\mathrm{MeV}^1`$ – about equal to that found by Mohr et al. , but half of that found in other theoretical studies , and about 2/3 of that suggested in a recent compilation of astrophysical reaction rates .

### B $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$

The low-energy $`S`$-factor computed from the Kim A cluster-cluster potential for $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$ is shown in Fig. 8, along with the contributions of individual terms of Eq. (21). Total $`S`$-factors from five cluster-cluster potentials are shown along with the laboratory data in Fig. 9. The computed branching ratios are shown along with the corresponding data in Fig. 10. The discussions of small contributions to the $`S`$-factor and of the precision of the results for $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ above also apply here, again with a Monte Carlo error estimate of about 10% in $`S`$-factor normalization. We again break down the results for $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$ into normalization, branching ratio, and energy dependence:

1) Normalization. Our $`S`$-factors are more than 10% lower than the lowest data set (after applying the renormalization of the Kräwinkel data set  recommended by Hilgemeier et al. ), and nearly a factor of two lower than the highest data sets (note that this includes the data of Volk et al. , which are not shown because their results were published only as extrapolated $`S(0)`$ values). By applying five different phenomenological cluster-cluster potentials that are not in disagreement with the low-energy elastic-scattering data, we obtain a small variation in $`S(0)`$ about a mean of $`S(0)=0.40`$ keV b. Using the potential which best matches the low-energy $`S`$-wave scattering (potential A of Kim et al. ), we obtain $`S(0)=0.40`$ keV b. Although this is closer to matching the lower numbers found in capture photon experiments than the delayed activity experiments, it is not a close match in normalization to any of the experimental results. Our results are therefore not useful for addressing possible systematic problems in the data.

2) Branching ratio. It is seen from the branching ratios in Fig. 10 that our calculation is in reasonable agreement with the laboratory data with regard to relative strengths of transitions to the ground and first excited states of 7Be. This suggests that the ANCs of the VMC wave functions for both states are too small, and by about the same factor in each case. Using the results for $`S`$-factor normalization above, we conclude that this factor is in the range 5–25%.

3) Energy dependence. We renormalized our results to best fit each of the larger data sets separately, and computed chi-squared statistics in each case to determine goodness of fit. The results are shown in Table III, and indicate general agreement. For the logarithmic derivative of $`S(E)`$ at $`E=0`$, we obtain $`0.57\mathrm{MeV}^1`$, in reasonable agreement with other estimates in the literature , where published models fall in the range $`0.50`$ to $`0.72`$. It is about equal to the value presently used in solar neutrino work .

## V Implications and Applications

We have carried out the first calculation of the low-energy $`S`$-factors for the processes $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ and $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$ based on realistic nucleon-nucleon potentials. Seven-body wave functions for these potentials, constructed by the VMC method and constrained to have the correct asymptotic forms, produce $`S`$-factor energy dependences for the processes $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ and $`{}_{}{}^{3}\mathrm{He}(\alpha ,\gamma )^7\mathrm{Be}`$ which agree reasonably well with experiment.

This work indicates no serious problems in the reaction rates presently used in astrophysical models. In fact, the most important implication of these results for astrophysics is probably that the previous understanding of these reactions remains essentially unchallenged. For example, using the present calculation to extrapolate the Robertson et al. activity measurement  of the $`S`$-factor from 0.9 MeV to 0 MeV, we obtain essentially the same result as with the energy dependence currently used in the standard solar model . For big-bang nucleosynthesis, there is no low-energy extrapolation problem. The most useful result which a theoretical study could provide for cosmology would therefore be tighter constraints on cross-section normalizations than the current body of experimental data provides. The present calculations have not achieved that goal, but future first-principles calculations based on realistic nucleon-nucleon interactions might.

The only serious problem with the results presented here is the low normalization of the $`S`$-factors, and its principal cause is probably easy to identify. Because most of the cross section arises at large ($`>10`$ fm) cluster separations, the low normalizations most likely arise from the form of the seven-body bound states at large separations of the p-shell nucleons from the s-shell core. Part of the discrepancy may arise from the correlation $`f_{sp}`$. The correlation between $`\alpha `$ and $`\tau `$ clusters in the bound states is proportional to the twelfth power of $`f_{sp}`$, so that the long-range correlation between clusters is very sensitive to the choice of $`f_{sp}`$. However, it is hard to see how the $`f_{sp}`$ used in this study could affect the wave function beyond about 5 fm cluster separation. In the case of $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$, a more likely explanation is that the Monte Carlo uncertainty has been underestimated, and a (prohibitively) long integration would “converge” on values in better agreement with experiment. Other aspects of the present calculation which are not done as exactly as one may wish, and therefore may bear on this problem, include the fact that asymptotic forms have been imposed on the wave functions which are inconsistent with their energies, as well as the use of cluster-cluster interactions which were phenomenologically constructed and could conceivably be inconsistent with the rest of the calculation in subtle ways.

Despite these lingering difficulties, we have demonstrated the applicability of a new, almost ab initio approach for computing low-energy radiative captures in cases for which precise measurements exist, achieving the same accuracy as previous theoretical approaches. We also have presented additional evidence for uncertainty in the logarithmic derivative of the $`{}_{}{}^{3}\mathrm{H}(\alpha ,\gamma )^7\mathrm{Li}`$ $`S`$-factor at zero energy.

Regarding the present calculation as a first application of realistic potentials to radiative capture at mass seven, it has been successful, and it has taught important lessons about the role of the $`f_{sp}`$ correlations that were not apparent in previous VMC studies. This work clears the way for more refined models of these radiative captures based on realistic potentials. Specific improvements which will be possible in the near future include the use of improved three-body potentials now in development  and the use of essentially exact wave functions derived by the Green’s function Monte Carlo technique for the seven-body bound states. We note that the application of GFMC will require new wave functions, starting from VMC wave functions of the type used here to get the asymptotic forms correct. In the more distant future, it should be possible to perform the whole calculation self-consistently, constructing scattering wave functions from the potentials in a way similar to that used for the bound-state wave functions. VMC-based work would also profit by a more systematic effort to produce a Monte Carlo weighting scheme well-suited to computing the sorts of matrix elements encounted in direct capture calculations. Prospects for significant improvement on this initial investigation are very good.

###### Acknowledgements.

The author gratefully acknowledges R. Schiavilla, S. C. Pieper, V. R. Pandharipande, and M. S. Turner for many helpful discussions, and especially R. B. Wiringa for providing the bound state wave functions and the variational Monte Carlo code used to compute them, as well as extensive support in their use. Computations were performed on the IBM SP of the Mathematics and Computer Science Division, Argonne National Laboratory. This work was supported by the U. S. Department of Energy, Nuclear Physics Division, under contract No. W-31-109-ENG-38. This work was also performed at Argonne National Laboratory while a Laboratory-Graduate participant, a program administered by the Argonne Division of Educational Programs with funding from the U. S. Department of Energy.


text:
# Vortex bending and tightly packed vortex lattices in Bose-Einstein condensates

## I Introduction

A stirred coffee climbing up the cup walls, a whirl of water spiraling down a drain, or a tornado walking over the earth, are natural phenomena involving the rotation of a fluid, let it be coffee, water or dusty air. In all of them the fluid is subject to a balance between an external force and the centrifugal force, and this balance is responsible for a depression of the density along the axis of rotation. These structures have no trivial symmetries as is evident in the case of tornadoes, which have continuously changing bent shapes. None of these structures have long lives because viscosity, imperfections and other dissipative mechanisms play such an important role that rotation cannot be self sustained.

This behavior contrasts with that of the so called superfluids, which represent a state of matter of negligible viscosity or dry fluid . The lack of viscosity allows a superfluid to host rotation for long times with little or no external intervention. Besides, when a superfluid rotates it follows a special type of flow which is irrotational, $`\times 𝐯=0`$, except for a finite number of extended singularities, which are called vortex lines. These vortex lines, first predicted for 4He condensates, are the superfluid equivalent of robust whirls and tornadoes. However, due to centrifugal forces, the density of the superfluid  becomes zero along the vortex line something which is difficult to achieve in ordinary liquids.

Until recently only two examples of true superfluids were known: the superfluid phases of 4He and 3He. In both cases, the strength of the many–body interactions obscure the signatures of superfluidity and make their study difficult. The search for superfluid–type, weakly interacting systems has led to the experimental realizations of Bose–Einstein condensates (BEC) using dilute gases of alkali atoms . These condensates are ruled at low temperature by a mesoscopic quantum wavefunction, $`\psi (𝐫)=\sqrt{\rho (𝐫)}e^{i\theta (𝐫)}`$, whose phase determines a velocity field, $`𝐯=\mathrm{}\theta (𝐫)/m`$, that may host vortices along the lines where $`\theta `$ is not well defined. The existence of vortices in BEC has been directly confirmed in two sets of experiments. The first one is based on the manipulation of the internal degrees of freedom of an almost spherical condensate . In the second one, vortices are produced by “mechanical means” in rotating elongated traps with transverse asymmetries , a geometry that leads to intriguing effects on the vortex nucleation process.

In Ref.  we find several commonly accepted predictions about vortices in rotating traps. Namely, (i) a rotating trap leads to the production of one or more vortices of unit topological charge $`m`$, being $`m=(2\pi )^1\theta `$. (ii) Vortices appear in finite numbers, with straight shapes, forming lattices with p–fold symmetries. (iii) Certain critical rotation speeds $`\mathrm{\Omega }_p<\mathrm{\Omega }_{p+1}`$ should be surpassed before each new vortex is nucleated. (iv) Asymmetric states without rotational symmetry of any kind are found to be energetically unstable, which means that there are perturbations with no energetical cost that can destroy such states. However, Ref.  is based on a variational ansatz for the weak coupling limit, which lacks longitudinal degrees of freedom and implicitly induces the $`p`$–fold symmetry of the vortex lattice. Recent works abandoned these constraints but either focused on infinitely long condensates or studied almost spherically symmetric traps .

In this paper we study an elongated three–dimensional Bose–Einstein condensate subject to rotation. Our main result is that stationary vortex lines are bent even in completely symmetric setups and that the bending depends on the type of trap on the number of atoms of the experiment. The fact that vortex lines may bend was mentioned in Ref.  where longitudinal excitation modes of a vortex line induce vortex bending. However the physics in that paper differs drastically from our work, in which we consider the ground state and the generation of stable and stationary structures. Second, in a complementary paper  we have presented some evidences in favor of that bending and how it could be related to the experimental results of the ENS . In this work we extend that work to the analysis of arrays of bent vortices and its properties as well as the dependence of the bending of the vortices on the geometry of the trap.

Our work consists on two parts. In Section II we provide the mathematical foundations for the search of the rotating condensate most favorable configurations. To do so we first introduce the mean-field model used to describe the condensate together with some basic definitions and a new energy functional, whose minima are the possible ground states of the condensate. In this section we also develop an efficient minimization method for this functional.

All these tools are then applied in Sect. III to different configurations. We start with setups from current experiments  and show that in elongated traps with many bosons the vortex lines are bent and form Abrikosov lattices which are regular only in the core of the gas cloud and deform close to the boundaries. Next we study the dependence of the vortex bending with respect to the free parameters in our theory. Finally in Sect. IV we offer our conclusions together with some open questions.

## II Mathematical tools

### A The mean-field model

For current experiments it is an accurate approximation to use the zero temperature mean field theory of the condensate, in which the atomic cloud is described by a single wavefunction $`\psi (𝐫,t)`$ ruled by a Gross-Pitaevskii equation (GPE). In the case of rotating systems it is useful to consider the problem on the mobile reference frame that moves with the trap , in which the equation reads

$$i\frac{\psi }{t}=\left[\frac{1}{2}\mathrm{}+V_0(𝐫)+g\left|\psi \right|^2\mathrm{\Omega }L_z\right]\psi .$$
(1)
Here $`L_z=i\left(x_yy_x\right)`$ is the hermitian operator that represents the angular momentum along the z-axis and the effective trapping potential in given by

$$V_0(𝐫)=\frac{1}{2}\omega _{}^2[(1\epsilon )x^2+(1+\epsilon )y^2]+\frac{1}{2}\omega _z^2z^2].$$
(2)
In Eq. (1) we have applied a convenient adimensionalization which uses the harmonic oscillator length, $`a_{}=\sqrt{\mathrm{}/m_{Rb}\omega _{}}`$ and period, $`\tau =\omega _{}^1.`$ With these units the nonlinear parameter becomes $`g=4\pi a_\text{S}/a_{}`$.

### B Variational formulation

There are several conserved quantities associated to Eq. (1). The first one is the norm of the wavefunction, $`N[\psi ]=\left|\psi \right|^2𝑑𝐫`$, which is related to the number of bosons in the condensate. The second conserved equation is the energy of the gaseous condensate

$`E[\psi ]`$ $`=`$ $`{\displaystyle \overline{\psi }\left[\frac{1}{2}\mathrm{}+V_0\left(𝐫\right)+\frac{g}{2}\left|\psi \right|^2\mathrm{\Omega }L_z\right]\psi 𝑑𝐫}`$ (3)
$`=`$ $`E_0(\psi )\mathrm{\Omega }L_z(\psi ).`$ (4)
In current experiments with stirred condensates , the gaseous cloud reaches certain long lived configurations with one ore more vortices. In this work we are not interested on the precise dynamics of the nucleation process, but look for the final stationary configurations. Due to thermodynamical considerations we expect the condensate to achieve the least energy for given experimental parameters, $`\{N,g,\mathrm{\Omega }\}`$. Such configuration is called a *ground state*.

Within the framework of our mean-field model, a stationary state is represented by a wavefunction with the form

$$\psi _\mu (𝐫,t)=e^{i\mu t}\varphi (𝐫).$$
(5)
The first way to find such solutions is to introduce (5) in our mean–field model (1) and directly solve the resulting partial differential equation. This method has several disadvantages: (i) It is difficult due to the many degrees of freedom that it involves, (ii) this new equation is satisfied not only by the ground state but also by excited states, and (iii) there is no mathematical guarantee yet for the existence of solutions of this problem.

The second way to characterize a stationary is solution is by studying the energy functional and use the fact that each stationary configuration is a critical point of the energy (4) for given $`\{N,g,\mathrm{\Omega }\}`$,

$$\frac{E}{\psi }|_N[\psi _\mu ]=0.$$
(6)
Since we are indeed looking for the ground state, one should use a minimization procedure able to find the minima of the energy functional for given parameters.

### C Reshaping the energy functional

As we stated above, our objective is to find the solutions $`\psi _\mu `$ which are the minima of the energy subject to the restriction $`|\psi |^2N`$, for a fixed angular speed of the trap and for a given interaction. The existence of at least one minima or ground state for this variational problem has been proved elsewhere , they are guaranteed to be stable and represent the energetically most favorable configurations for given $`N`$ and $`\mathrm{\Omega }`$.

Such a problem is called a constrained optimization and we will refer to the solution, $`\psi _\mu `$, as a constrained minimum. The existence of constraints in an optimization problem poses serious difficulties to traditional descent methods, since it is difficult to design a efficient minimization algorithm which takes care of the constraints at each step . To take into account this restriction one should use Lagrange’s multipliers, i.e. just adding a fraction of the constraint $`\omega (\mu ,N)`$ (also called a “penalizer”) to the original functional

$$F[\psi ]=E[\psi ]+\omega (\mu ,N).$$
(7)
It is not difficult to show that $`F[\psi ]`$ and $`E[\psi ]`$ have the same stationary states, and that any absolute or relative minimum of $`F[\psi ]`$ is also a constrained minimum of $`E[\psi ]`$.

One expects that a wise choice of the penalizer will establish a one-to-one correspondence between the values of the chemical potential, $`\mu `$, the absolute minima of $`F[\psi ]`$ and the ground states of our condensate. The advantage of $`F[\psi ]`$ over $`E[\psi ]`$ is that in the new functional there is no need to care for constraints: the value of $`N`$ for the ground state is determined by the chemical potential, $`\mu `$, which is fixed throughout the minimization process.

The traditional choice for a Lagrange multiplier leads to the definition of the free energy

$$[\psi ]=E[\psi ]\mu N[\psi ].$$
(8)
However, this functional is not bounded below and thus $`\psi _\mu `$ is at most a local minimum of $``$. Consequently, this functional cannot be used to characterize the state of the condensate. This is an interesting result, since one could be tempted to think that the traditional definitions of the thermodynamical potentials are suitable to characterize all physical systems. We must remember, however, that our condensate is being described by a mean-field model, a simplification of a more complex model, and as such what works for the full problem needs not work for the simpler one.

In our search for suitable functionals we have found a simple one which we call the nonlinear free energy

$$F[\psi ]=E[\psi ]+\frac{1}{2}\left(N[\psi ]\lambda \right)^2.$$
(9)
First and most important, it can be proven that $`F[\psi ]`$ has at least one finite norm absolute minimum for each value of $`\lambda `$ , and that each of those minima corresponds to a constrained minima of the energy, $`E[\psi ]`$. And second, there exists a simple and invertible relation between the usual thermodynamic variables —the chemical potential and the number of particles, $`(\mu ,N)`$—, and our new variables —Lagrange’s constant and the number of particles, $`(\lambda ,N)`$—,

$$\mu =N[\psi _\mu ]\lambda .$$
(10)
An important feature of our new functional (9) is that it highly nonlinear with respect to $`\psi `$. This poses no new difficulty, since our original equations (1) were already nonlinear. Indeed Eq. (9) has proven to be the most natural choice for many other problems, such as the propagation of incoherently coupled laser beams through saturable media .

### D Optimal methods for minimization

To minimize the nonlinear free energy (9) we follow Ref. . First we choose the right function space, which in our case is the Sobolev space $`H^1(R^3)\{\psi /\psi ,\psi L^2\}`$ of functions which admit at least one spatial derivative. This space is equipped with a scalar product

$$\psi ,\varphi \left[\overline{\psi }(𝐫)\varphi (𝐫)+\overline{\psi }(𝐫)\varphi (𝐫)\right]d^nr.$$
(11)
and a norm $`\psi =\psi ,\psi `$. To obtain an explicit expression for the gradient of the functional in $`H^1`$ we perform a first order expansion of $`F[\psi ]`$ around a trial state $`\psi `$

$`F\left[\psi +ϵ\delta \right]`$ $`=`$ $`F[\psi ]+2ϵ\mathrm{Re}{\displaystyle \left(\delta ^{}\frac{E}{\overline{\psi }}+\delta ^{}\frac{E}{\psi }\right)}`$ (12)
$`+`$ $`O(ϵ^2).`$ (13)
We have to turn this expression into something which looks like Frechet’s definition of a derivative. This means that we have to find some $`\varphi `$ such that

$$\mathrm{Re}\left[\overline{\delta }\frac{E}{\overline{\psi }}+\overline{\delta }\frac{E}{(\overline{\psi })}\right]=\mathrm{Re}\left[\overline{\delta }\varphi +\overline{\delta }\varphi \right].$$
(14)
If we integrate by parts and impose that this equality be satisfied for all perturbations, $`\delta `$, the problem has a formal solution which is given by a Lagrange equation

$$\left(1\mathrm{}\right)\varphi =\frac{E}{\overline{\psi }}\frac{E}{(\overline{\psi })}.$$
(15)
In consequence, our formal expression for the Sobolev gradient of $`F[\psi ]`$ finally reads $`_SF\left(1\mathrm{}\right)^1F,`$ where $`_SF`$ stands for the Sobolev gradient, $`F`$ is the ordinary one, and $`(1\mathrm{})^1`$ represents the inverse of a linear and strictly positive definite operator.

We may now use the Sobolev gradient of our functional as the direction of descent for a minimization procedure

$`{\displaystyle \frac{\nu }{\tau }}(𝐫,\tau )=\left(1\mathrm{}\right)^1\left[{\displaystyle \frac{1}{2}}\mathrm{}+V+g|\nu |^2\mathrm{\Omega }L_z\right]\nu .`$ (16)
The preceding equation converges to some stationary state $`\varphi (𝐫)=lim_\tau \mathrm{}\nu (𝐫,\tau )`$. To grant convergence to the true ground state one needs to improve the descent method, using for instance the nonlinear conjugate gradient method instead of (16), and introducing some type of relaxation which helps avoid saddle node points. In any case, the minimization process must be performed on suitable space of functions that we built using a discrete Fourier basis with $`64\times 64\times 128`$ modes.

## III Results

### A Bent vortex lattices in current experiments

We have applied the numerical methods outlined above to different setups. The family of numerical experiments in this subsection resembles the experiments of Madison et al. . For these experiments with 87Rb condensates, the bosonic interaction is ruled by the scattering length $`a_\text{S}5.5\text{nm}`$. For several pictures in this paper, the number of atoms was chosen to match, $`Ng=9000`$, which corresponds to a few times $`10^5`$ Rb bosons but our results remain qualitatively valid for an ample range of $`gN`$ values as will be shown below.

The geometry of the trap is a very elongated one, $`\omega _{}>18.7\omega _z`$, and for the small transverse deformation of the trap we have tried $`\epsilon =0.0`$ as well as $`\epsilon =0.03`$, the latter being the closest one to the actual experiment . This means that we have studied both *axially symmetric setups and completely asymmetric setups*, and in both cases our study has provided essentially the same results. This is due to the fact that an intense nonlinear interaction $`(Ng10^4)`$ effectively prolongs the existence of all solutions from the symmetric setup to the asymmetric setup .

The first result is that an increase of the angular speed causes the minimum of the energy to move from a nodeless ground state into states with one, two, three and more vortices. Consequently, as it was already predicted , there exists a cascade of increasing angular frequencies, $`\{\mathrm{\Omega }_1<\mathrm{\Omega }_2<\mathrm{}\}`$ for the nucleation of one, two and more vortices —the larger the rotation frequency, the more vortices, and for our setup, $`\mathrm{\Omega }_10.4(0)`$, $`\mathrm{\Omega }_20.5(3)`$, $`\mathrm{\Omega }_30.5(8)`$, and $`\mathrm{\Omega }_40.6(4)`$. Figure 1 shows the structure of a ground state hosting two vortices at $`\mathrm{\Omega }=0.55`$.

The second and most important result, which is already evident in Fig. 1 is that vortices are nucleated with a stable and stationary bent shape even in axially symmetric traps. In Fig. 2 we show three–dimensional pictures of a condensate with up to four bent vortex lines.

These contorted shapes are absolutely stable configurations which lack rotational symmetry of any kind, even discrete, something which was thought to be forbidden according to Ref. . In that case asymmetric states were found only for discrete values of the angular speed, right on the transition between different number of vortices, and always exhibiting energetic instabilities. Furthermore, due to the bending, while the trunk of the condensate rotates, the caps remain almost still. This feature allows the gas to accommodate a fractional value of the angular momentum which is between 0 and 1, and which never reaches 1, as it was already conjectured in . As the speed of the trap increases, the angular momentum grows continuously by means of pulling the vortex line straighter and not only with discontinuous jumps.

### B Regularity of the vortex lattice

From old studies with liquid helium and more recent works with condensates of simpler geometry, it has been expected that vortices should form regular triangular lattices. On the other hand, our finding above seems to suggest that the vortex aggregate rather adopts an irregular shape. As we will show below, the condensate actually develops a medium range order, which dilutes in the vicinity of the cloud boundary.

In Fig. 3 we use the solutions of our previous section to simulate the pictures that should be seen in experiments. The first row of two–dimensional plots \[Fig. 3(a-d)\] are upper views of the condensate, and represent the column density of the bosonic cloud along the transverse directions. These pictures bear an extreme resemblance with the experimental photographs of Refs. , showing blurred clouds where the vortex cores seem partially filled.

In Fig. 3(i-l) we show the density of the condensate as seen at half the height of the condensate, i.e. $`|\psi (x,y,z=0)|^2`$. At the moment of revising this paper an experiment has been made which observes this type of transverse cuts of a stirred sodium condensate . In that experiment a thin slice of bosons is promoted to a different internal state and then imaged separately from the rest of the cloud. The resulting pictures are like those in Fig. 3(i-l), where the regularity of the Abrikosov vortex lattice is made evident. The blurring has disappeared and the holes arrange along the expected triangular lattice. In the above mentioned paper the bending of the vortex lines is also reported.

As further confirmation we suggest that using the setup of Ref.  and watching the condensate not from above, but from one side would allow a direct observation of bent vortex lines lead to pictures as those in Fig. (e-h).

### C What is the ultimate cause of bending?

The fact that the bending of the vortex line takes place not only on the asymmetric trap of current experiments, but also on an radially symmetric trap, $`\epsilon =0`$ represents a surprising novel type of symmetry breaking in which the superfluid not only chooses the sense of rotation, but also a plane for its bending.

It is important to emphasize the counterintuitive nature of the vortex bending in the case of an isolated vortex line. First, opposite to the case of 4He, the BEC is not constrained by any recipient and there are no asymmetric boundary conditions which could easily explain the deformation of the fundamental solutions. Second, although qualitative reasons for such bending may be found a–posteriori, they have never been reported before. And even though physical arguments may justify the bending of the vortex line, they can hardly support the fact that the bent vortex line is a stable stationary configuration. Both the stability and the stationarity can only be demonstrated with a direct study of the energy functional.

Nevertheless, it is legitimate to ask for the ultimate reasons of this rupture of symmetry. Since our search of states was based on energetical considerations, we have studied the dependence of the bending with respect to each of the free parameters in our problem: $`N`$, $`\mathrm{\Omega }`$, $`g`$, $`\epsilon `$, and the elongation $`\gamma =\omega _{}/\omega _z`$.

First, the influence of $`\epsilon `$ is discarded. As we mentioned above, bending exists in either symmetric or asymmetric traps. Rather, the asymmetry seems to take part only in the dynamics of the nucleation process, by changing the value of the critical frequencies and inducing a type of hysteresis .

A different role is played by the elongation of the trap, $`\gamma `$, and the effective interaction $`U=Ng`$. Starting with a configuration and lowering either the elongation or the number of atoms we see that the bending becomes smaller and eventually disappears \[Fig. 4(a-c)\].

We have tried to measure the bending so as to determine the minimal enlogation and interaction which are required to induce this phenomenon. The problem is that due to the small changes of $`L_z`$, the full minimization procedure does not allow us to bound this minimal elongation accurately. An alternative procedure is to study the relation between the angular speed which is required to stabilize a straight vortex, $`\overline{\mathrm{\Omega }}_1`$ and the angular speed at which the ground state acquires some angular momentum, $`\mathrm{\Omega }_1`$. The first value arises from the study of normal modes around a straight vortex, while the second value is the energy difference between a straight vortex and a vortexless state . When $`\overline{\mathrm{\Omega }}_1\mathrm{\Omega }_1`$ there are values of the angular speed where bending can be favorable. This phenomenon has been referred to in the literature as “anomalous modes” , and it is a signature of bending.

In Fig. 5 we plot the difference $`\overline{\mathrm{\Omega }}_1\mathrm{\Omega }_1`$, as a function of the elongation, $`\gamma `$, and of the effective interaction, $`U=Ng`$ for an axially symmetric trap. According to this study, a large elongation of the trap is required to make a bent vortex energetically favorable over a straight one.

## IV Conclusions

In this work we have studied vortex lines in a very elongated Bose–Einstein condensate in a rotating trap. By developing both a customized energy functional (Sect. II C) and an optimized descent method (Sect. II D) we have been able to find the ground state of a condensate in a family of different experimental configurations.

The first part of our study (Sect. III A) focused on realistic values of the experimental parameters, taken from the work by Madison et al. . The main conclusion of this work is that by increasing the rotation speed, the condensate achieves states with one, two, three and more bent vortices. Such vortices form a regular Abrikosov lattice at half the height of the condensate and deform not very far from the condensate core.

Appart from the observation reported in  and while current experiments at the ENS  are carried on a regime in which vortex lines should bend, this has not been yet observed there. The reason is that the imaging of a condensate which is a few $`\mu `$m large requires a previous expansion in which the cloud becomes disk shaped and bending cannot be appreciated. Furthermore, the study of the elongated cloud without vortices  shows that the nucleation of vortices is subject to hysteresis and one must actually exceed some rotation frequency $`\mathrm{\Omega }_m`$, which is rather large and makes it difficult to selectively produce one, two or more vortices.

On the other hand the bending of vortex lines has signatures which are also observed in current experiments. First, the bending by itself can explain the apparent filling of vortex lines when seen from above , and the regular pictures which appear in most recent observations with sliced sodium condensates . Second the bending also accounts for the continuous growth of the angular momentum with respect to the rotation speed, and the fractional values of the angular momentum, $`0<l<1`$, which arise in the indirect measures of  and . And finally, the nucleation of many bent vortices and their subsequent interaction during the expansion phase may lead to turbulent structures that should explain the lack of regular pictures above a certain angular speed .

The second part of this paper studied the dependence of the bending with respect to the parameters of the mean field model (1). Here we conclude that it is both the elongation and interactions which induce the bending of the vortex line, while the transverse asymmetry plays no important role.

Our results imply that past studies devoted to the quasilinear limit ($`U1`$) would become of no applicability for elongated traps due to the lack of bending in the simplified models. Second, along this line it would be nice to find analytically the geometry of the vortex line as a function of the interaction $`U`$ and the elongation $`\gamma `$.

Our study also reveals that current experiments are being developed in a regime which is qualitatively different from the setups that have been studied up to date . Values of the relevant parameters $`(\gamma ,U)`$ in current experiments are so far from most studies that new methods must be developed to accurately describe these amazing systems.

As an example, let us pose one of the technically difficult questions that arise in this work. From our numerical work it seems that the extremes of some of the vortex lines actually reconnect not too far away from the core of the cloud, forming what is called a vortex ring. Although it is not possible with current numerical methods to fully support this conjecture, vortex rings have already been observed in BEC . The limits of zero and infinite radius of a vortex ring lead to a dimensionless zero and an isolated vortex line, respectively. Therefore vortex rings would be a nice tool for explaining the nucleation of vortices and they would allow to interpret the states in this work as the result of an incomplete nucleation. The proof or refutation of this conjecture remains a mathematically challenging problem.

This work has been supported by grant BFM2000-0521.