text:
# Reconstruction and Analysis on Demand: A Success Story

## I INTRODUCTION

To deal with the vast amounts of data collected, HEP experiments usually write their own data access software. This software has to be used for a multitude of tasks: calibration, reconstruction, monte carlo generation, and analysis. In this paper we compare the ’standard’ data access system design to the CLEO III design which uses an ’on-demand’ mechanism known as lazy evaluation. In addition, we describe our experience of having used the CLEO III system for the past four years and in particular the response of the physicists who use this system.

## II STANDARD SYSTEM

### II.1 Description

The standard system design basf ; d0  that most HEP experiments use for their data access software is primarily designed to deal with the task of reconstruction. For this task, all data objects need to be created for each event that is being processed, since it is the ultimate goal of this task to write all those data objects to persistent storage.

In the standard system design, each processing step is partitioned into its own software Module. For example, the track finding algorithm and the track fitting algorithm are placed into their own software Modules. In addition, there are Input Modules responsible for reading objects from persistent storage and Output Modules responsible for writing objects to persistent storage.

To form an actual processing job, the software Modules are run in a user-specified sequence. When it is a Module’s turn to run, that Module executes its algorithm and then places its data into the event. In addition, if a Module for some reason rejects that event, then the Module can tell the system to stop processing that event (thereby skipping the processing of all Modules appearing later in the sequence) and to restart the sequence with the next event.

A simple example of a processing job run using the standard system design is shown in Figure 1. In this example four different Modules are used to process an event. The first Module is the Input Module which reads objects from a persistent object store (in this case it would be objects corresponding to hits that have been calibrated) and then inserts those objects into the event. The second Module in the sequence pulls the calibrated hits from the event and runs the track finding algorithm which then inserts the found tracks into the event. The third Module pulls the calibrated hits and the found tracks out of the event, refits the tracks and then inserts the refitted tracks into the event. The final Module is the Output Module which pulls objects out of the event and writes them to a persistent object store.

### II.2 Critique

The standard system design has worked well for many experiments for a number of reasons. One of the main reasons is it has a very simple conceptual model for how a job actually runs. The job simply runs a series of algorithms that are in a user defined sequence. This simple conceptual model makes the physicists feel confident that they know how the system works which in turn aids in acceptance of the system. Another reason experiments use the standard design is that jobs are fairly easy to debug since it is easy to determine which module had a problem.

Unfortunately, the simplicity of the system leads to a number of problems, particularly when the system is used for analysis. One of the main problems is the physicists using the system must know the inter-module dependencies in order to place them in the correct sequence. To avoid this problem, physicists often run jobs with many Modules they do not need in order to avoid missing a Module they might need. Another problem is that optimization of the Module sequence must be done by hand. For example, reconstruction jobs often have filter Modules which stop the processing of events that contain no interesting physics (e.g. accelerator related backgrounds). To get the best performance for the reconstruction job, only the Modules creating the data needed by the filter Module should be ahead of the filter Module in the Module sequence. Since the Module sequence must be explicitly set by the physicist it is up to the physicist to find the optimal Module sequence. A final problem is reading back from storage is almost always inefficient. The reason is that all the objects from storage must be created at the beginning of the event processing, since that is where the Input Module is in the sequence, even if the job does not use all objects.

## III ON-DEMAND SYSTEM

### III.1 Description

The on-demand system used by the CLEO III experiment cleo3  was primarily designed to be used for analysis batch processing. The main difference between reconstruction and analysis processing is that in analysis processing not all data objects need to be created every event. Therefore the on-demand system only creates data objects as they are requested.

Similar to the standard system, the on-demand system breaks the processing task into separate Modules. But in the on-demand system the types of Modules are further refined into two categories , each with two Module types.

return data when requested.

reads data from a persistent store.

creates data on demand by running an algorithm.

sequentially run for each new event.

analyzes and filters events.

writes data to a persistent store.

Since data objects are only created when requested, data providers register with the system what data they can provide. This registration is how the system knows which provider to request the data from when another Module requests that type of data. Since the method used to create a data object (either by reading from a persistent store or by running an algorithm) is activated the first time that data object is requested, the processing sequence is now set implicitly by the order of data requests and not explicitly by the physicist. Essentially the processing sequence is automatically determined by the system without the need for manual intervention.

While the data providers set what data is available in a job, what is done with that data is determined by the data requestors. Processors analyze events (e.g. fill histograms or ntuples) and/or filter events of interest by stopping further processing of certain events. In this system only Processors can stop the processing of an event. The order of the Processors must be explicitly set by the physicist since different Processor orders can be meaningful and only the physicist knows which order is proper for the job. For example, a physicist could place the event display either before or after an event filter to either see all the events going into the filter or to just see the events that pass the filter. In contrast, Sinks are always at the end of the processing order and therefore only write out data for events that pass all the filters.

Figure 2 shows how an event is processed in the on-demand system. A Source determines the next event to be processed (but does not fill that event with data). This event is passed to the first Processor in the sequence. If that Processor does not reject that event, the event is passed to the following Processor in the chain. This continues until one of the Processors rejects the event, at which time the system goes back to the Source to get a new event, or all Processors have been run. If no Processor rejects the event, the event is given to the Sinks for storage.

### III.2 Data Model

A data access system provides not only a mechanism for accessing data but also a data model that describes how the data is organized. The CLEO III data model is unusual because it provides a unified model for all data so, for example, event data and calibration data are treated in the same manner c3data-model .

Figure 3 gives a graphical view of the CLEO III Frame/Stream data model. In this data model, all data items are held in Records. Data items are grouped into Records based on the ’life-time’ of the data. For example, since both fitted tracks and electro-magnetic showers are only relevant for the instant an event is recorded, both those data items are placed in the event Record. In contrast, the energy of the beams in the accelerator’s storage ring is constant through-out a data run so the beam energy is held in the begin run Record rather than the event Record. Streams are a time ordered sequence of Records of the same type. For example, the event Stream holds all the event Records in the proper time order. The last concept in the data model is the Frame. A Frame is a collection of Records that describe the state of the detector at an instant in time. For example, to understand the data in event Record ’F’ requires the data at begin run Record ’B’, geometry Record ’B’ and calibration Record ’A’.

Because all data are accessed in the exact same way, the physicist’s learning curve is reduced. For example, accessing fitted tracks or detector alignment data are done in exactly the same manner. First the physicist must get the appropriate Record from the Frame. Second, he must ask the Record to return the data in question. Also, since the event Record is no different from any other Record, physicists can study any type of data in the exact same way they study event data. For example, a physicist could analyze the way the alignment of the detector changes with time by creating a Processor that analyzes detector alignment Records rather than event Records.

Since the Frame contains all data relevant for any processing job, the system is designed so that the different component Modules that make up a job can only communicate with each other via the data they have registered in the Frame. In a sense, this makes the Frame the communication bus for the system. Using the Frame as the communication bus allows us to avoid explicit dependencies between Modules. Modules are only dependent on the data objects they use or produce, and not on how those data objects are created.

### III.3 On-Demand Mechanism

Since the CLEO III system uses on-demand processing, we need a mechanism to communicate between the data requestors and the data producers. The mechanism we use employs Proxies. A data provider registers a Proxy for each data type the provider can create. These Proxies are placed in the appropriate Record and are indexed by a key. The key is composed of three different tags:

the class type of the object returned by the Proxy

an optional compile-time string describing the use of the object

an optional run-time settable string.

The two string tags used in the key allow a Record to contain many objects of the same class type. The Production tag allows us to compare the results of two providers which create data with the same Type and Usage tags.

Physicists access the data via a type-safe templated function call. E.g.,

```
List<FitPion> pions;
extract(iFrame.record(kEvent), pions);

```

In this example, the first line defines a variable named ’pions’ to hold a list of FitPions. In the second line, the templated ’extract’ function is given the event Record (which is obtained from the Frame via the member function ’record’) and the ’pions’ variable. The Type tag is determined at compile time based on the type of the ’pions’ variable. The Usage and Production tags are set to their default values. Once the key is built the extract call asks the Record for the Proxy. After the Proxy is obtained, the extract call tells the Proxy to run the algorithm to deliver the data. If the Proxy’s algorithm runs successfully, the data is cached in case another request for the data occurs later in the job. If an error occurs while trying to obtain the data, an appropriate exception is thrown.

Two examples of the on-demand system in action are shown in Figure 4. In the left hand example, the Processor SelectBtoKPi gets FitPions by running the full tracking reconstruction algorithm. In the example, The FitPionsProxy extracts Tracks and CalibratedHits. The TracksProxy also extracts the CalibratedHits. The CalibratedHitsProxy then extracts the RawData (which is obtained from a file) and the pedestal and alignment calibrations which are stored in the Calibration database. In the right hand example, the exact same SelectBtoKPi, run in a different job, now gets its FitPions directly from the event database. Through the use of dynamic loading, the exact same SelectBtoKPi binary object can be run in the two different jobs without the need to relink the main executable.

### III.4 Critique

The on-demand design has many positive features. First, this design supports data access for all HEP jobs: online software trigger, online data quality monitoring, online event display, calibration, reconstruction, MC generation, offline event display and physics analysis. Second, in contrast to the standard design, the physicist does not need to explicitly set the proper procedural order of the job. The proper order is determined automatically simply by the order in which data is requested. Third, the design optimizes access from storage. At the beginning of processing a Record, a source only needs to say when a new Record is available. Decoding or retrieval of the data associated with that Record can be deferred until the first request for the data.

However, the added complexity of not explicitly setting the procedural order does make debugging and profiling more challenging. Good encapsulation of Modules through the Frame mechanism ensures that problems are isolated within a Module and are not from Module to Module interactions. With such isolation, you only need to know which Module you are in when a problem arises in order to know where to start debugging. For extremely serious problems that cause a signal to be thrown by the operating system (e.g., a memory access violation), looking at the last few entries in a function call trace-back in a debugger is usually sufficient to isolate the error. For standard run-time problems we have found that the use of exceptions is critical in understanding the system. The first implementation of our system used null pointers to data to signify that a problem had occurred during an attempted data access. This lead to a number of problems:

* it was impossible to know why a problem occurred. E.g., if FitPions could not be accessed was it because track fitting failed or because no calibrated hits were available?
* developers had to propagate any error encountered while accessing data needed for their algorithm.
* the return value of a data access had to be checked to be certain that data was obtained.

Using exceptions solved all of these problems. Now when an exception occurs, the type of and message in the exception explain the problem. Additionally, only in blocks of code that catch the exception and continue would it be necessary to check to see if the returned data was valid. To further aid in understanding what happened when an exception is thrown, we created a stack that holds the present data request chain. When an exception is caught by the system (because no Module caught the exception), we print a message such as

```
ERROR: caught an exception:
"Starting from SelectBtoKPi extracted
[1] type "List<FitPion>"
     usage ""
     production ""
[2] type "List<Track>"
      usage ""
      production ""
[3] type "List<CalibratedHit>"
     usage ""
     production "" <== exception occurred

No data "List<CalibratedHit>" "" ""
   in Record event.
 Please add a Source or Producer to your
  job which can deliver this data."

```

and then either terminate the job or skip this event and begin processing the next Record.

## IV WHAT WE HAVE LEARNED

We have been using a version of this system since September 1998. During that time we have gained a great deal of experience with an on-demand system. One of the most important findings is that the on-demand mechanism can be made fast. We have found that the Proxy lookup takes less than 1 part in $`10^7`$ of the CPU time on a simple job that processed 2000 events/s on a moderately powerful computer. Another finding is that cyclical dependencies (a Proxy that winds up extracting its own data) are easy to find and fix. We only have had one cyclical dependency and it showed up on the first test of the program by causing a stack overflow. In addition, we have found that we do not need to modify data once it has been created, our Records hold immutable data. In the cases where we need a new version of the data, we put the new data into the Record with a different key. One surprising finding is the on-demand system automatically optimizes the performance of reconstruction jobs. For example, it was trivial to add a filter to reconstruction which removed junk events by using found tracks rather than fitted tracks. The on-demand system made sure only the algorithms needed for the junk filter are run so that the minimum time is spent processing the event before the filter decides whether to keep it. We have also found that the best way to speed up analyses is to store many small objects. That way the job only needs to retrieve and decode the data needed for the current job.

Because the on-demand design is so different from the standard design, we have been very concerned with physicists’ reaction to the system, particularly since CLEO’s previous data access system was a standard design implemented in FORTRAN. In general, their response has been very positive to the new system. The programmers of the reconstruction code like the system. We made it easy for them to get started by having code skeleton generators to make a Proxy, Producer or Processor. By using the generators, coders only have to write the code specific to their algorithm and not write the code required to work within the system. In addition, the reconstruction coders found it very easy to test their code since they could easily swap Modules to compare the results of the jobs. This is extremely easy in our system since we use dynamic loading to load in the different Modules a job uses. Therefore when developers are making changes they only need to recompile their own Module, not relink the entire system.

To ease the transition for our physicists from the old standard system to the new on-demand system, we made sure that they could still program their analysis the ’old way’, i.e. all analysis code in the ’event’ routine. However, some of our analysis coders are now pushing the bounds of the system. These physicists are placing selectors (e.g., cuts for tracks) in Producers. This allows them to reuse the same binary for many different analyses. We are even seeing physicists share these selector binaries. So once selections are done in Producers, the physicists only use Processors to fill histograms or ntuples and to filter events. Then if a physicist stores the results of her selections into a persistent store, any subsequent pass through the data only requires reading back the selections from storage (and not rerunning the algorithm) and rerunning the Processor.

## V CONCLUSION

Our five years of experience with using an on-demand system has shown us it is possible to build such a system so that it is efficient, debuggable, capable of dealing with all the different types of data (not just data in an event), easy to write components, good for reconstruction and acceptable to physicists.

We believe that several reasons contributed to our success. The first is our use of skeleton code generators, so that physicists only have to write new code, not infrastructure ’glue’. Second, users do not need to register what data they may request, only what data they provide. This makes life easier for physicists since data reads occur more frequently than writes. And third, we use a simple rule for when algorithms are run: if a Producer is added it takes precedence over a Source that can deliver the same data.

###### Acknowledgements.

This work was supported by the National Science Foundation.


format: text

text:
# Correlated 𝑁-boson systems for arbitrary scattering length

## I Introduction

Condensation of a macroscopic number of bosons in the same quantum state was predicted many years ago bos24 . Much later this was experimentally achieved in the laboratory for dilute systems of alkali gases and95 ; bra95 ; dav95 . The average properties of these gases are accounted for by the Gross-Pitaevskii equation edw95 ; bay96 . Exhaustive reviews of the theoretical developments after the experimental breakthrough can be found in dal99 ; pet01 .

Degrees of freedom beyond the mean-field are crucial for the stability of the condensates, e.g., recombination into bound dimer and trimer cluster states nie99 ; esr99 ; bed00 . The importance of such correlations is revealed in recent experiments rob01 ; don01 . One example is the collapse of a Bose gas with large scattering length don01  where the lack of atoms in the condensate challenged the mean-field description in terms of the time-dependent Gross-Pitaevskii equation adh02b ; adh02d .

The mean-field description is valid for $`n|a_s|^31`$, where $`n`$ is the density and $`a_s`$ is the two-body $`s`$-wave scattering length, when the particles on average are outside the interaction volume of the order of scattering length, $`a_s`$, to the third power pet01 . The mean-field method neglects all correlations and thus breaks down at larger densities where correlations become important. Going beyond the mean-field is often very complicated as exemplified by Jastrow theory jas55 ; cow01 ; mou01 , which leads to high-dimensional equations. A later formulation is contained in the Faddeev-Yakubovskiĭ equations fad61 ; yak67  where the wave function is expressed in terms of components describing the asymptotic behaviour of all kinds of clusters.

Comparison of different models is not always straightforward, since different degrees of freedom are treated and the two-body interactions must be renormalized accordingly. For example, using realistic potentials in selfconsistent mean-field calculations leads to disastrous results because the Hilbert space does not include correlations esr99b  as needed to describe both the short- and long-range asymptotic behaviour. In Gross-Pitaevskii calculations the $`\delta `$-function interaction is renormalized to give the correct scattering length in the Born approximation pet01 . However, this substitution is only valid in the low-density limit. When correlations are included a different, and more realistic, interaction must also be used. Furthermore, only average properties can be described in the mean-field approximation. Thus comparisons of correlation dependent quantities are meaningless.

Five years ago an interesting alternative study of a condensate was formulated in terms of hyperspherical coordinates without any two-body correlations boh98 . Using the same coordinate system a theoretical frame for describing correlations was given soon after bar99a . Detailed three-body calculations with zero total angular momentum were recently performed in the same framework blu02c . Here the scattering length is varied and excited three-body states and any number of bound two-body molecular states are allowed. The claim in blu02c  is that higher-lying Bose-Einstein condensed states in a trap of length $`b_\text{t}`$ do not collapse when $`N|a_s|/b_\text{t}0.5`$ as otherwise indicated by experiments don01 . The recombination takes place at distances several times the scattering length. They conjecture that the properties for $`N>3`$ are quantitatively similar to these three-body results. Another study in the same framework investigates the model-dependence of the three-body energy and finds that only the large scattering length enter jon02 . They also indicate that the energy is insensitive to possible higher-order correlations for systems with many particles in a trap.

In a further development using the adiabatic hyperspherical expansion we formulated a method to describe two-body correlations in many-boson systems sor02 ; sor02b . This method is a novel attempt to describe correlated systems of low density. The formulation heavily relies on an additive set of components of the wave function as in the Faddeev decomposition but in contrast to the Jastrow multiplicative formulation. The numerical studies were limited to Bose-Einstein condensation for 20 particles.

The purpose of the present paper is to extend the applications to arbitrary scattering lengths and large particle numbers. We want to extract the general properties of the solutions especially for large scattering lengths where mean-field computations are invalid. We obtain naturally self-bound many-body systems, even when the two- and three-body subsystems are unbound. The paper begins with a brief description of the hyperspherical adiabatic expansion method in section II. Then in section III the details of the properties of the angular eigenvalue is discussed as the decisive ingredient in the radial potential prociding the information about the two-body interaction. In section IV we discuss the radial potential and the properties of the corresponding solutions. Finally in section V we discuss stability criteria for condensates expressed in terms of various time scales and decay rates. Section VI contains the conclusions.

## II Hyperspherical adiabatic method

We use the hyperspherical adiabatic expansion method with finite-range two-body interactions and simplifying assumptions about the wave function. We shall briefly describe the method and the assumptions. Details are given in sor02b .

The system of $`N`$ identical particles of mass $`m`$ is in the center of mass frame described by hyperspherical coordinates, i.e., one length, the hyperradius $`\rho `$, given by

$`\rho ^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i<j}{\overset{N}{}}}r_{ij}^2={\displaystyle \underset{i=1}{\overset{N}{}}}r_i^2NR^2`$ (1)
and $`3N4`$ hyperangles $`\mathrm{\Omega }`$ bar99a ; sor01 . The $`i`$th single-particle coordinate is $`\stackrel{}{r}_i`$, $`\stackrel{}{R}`$ is the center of mass coordinate, and

$`r_{ij}=|\stackrel{}{r}_j\stackrel{}{r}_i|\sqrt{2}\rho \mathrm{sin}\alpha _{ij},`$ (2)
with $`\alpha _{ij}[0,\pi /2]`$. The atoms are trapped in an external field approximated by a spherically symmetric harmonic oscillator potential of angular frequency $`\omega `$:

$`V_{\text{ext}}={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{1}{2}}m\omega ^2r_i^2={\displaystyle \frac{1}{2}}m\omega ^2(\rho ^2+NR^2).`$ (3)
Without any two-body interaction between the particles the ground-state wave function is a Hartree product of Gaussian amplitudes:

$`\mathrm{\Psi }_{\text{total}}={\displaystyle \underset{i=1}{\overset{N}{}}}e^{r_i^2/(2b_\text{t}^2)}=e^{\rho ^2/(2b_\text{t}^2)}e^{NR^2/(2b_\text{t}^2)},`$ (4)
where the trap length is $`b_\text{t}=\sqrt{\mathrm{}/(m\omega )}`$. The second radial moments are $`r_i^2=3b_\text{t}^2/2`$ and $`R^2=3b_\text{t}^2/(2N)`$. For large $`N`$ the average hyperradius therefore approaches the average mean-field radial coordinate times $`\sqrt{N}`$, see eq. (1). The hyperangles $`\mathrm{\Omega }`$ determine the relative orientations of the particles.

With a two-body interaction term $`V(r_{ij})`$ the total Hamiltonian becomes

$`\widehat{H}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{\widehat{p}_i^2}{2m}}+{\displaystyle \frac{1}{2}}m\omega ^2r_i^2\right)+{\displaystyle \underset{i<j}{\overset{N}{}}}V(r_{ij}).`$ (5)
It separates into a center of mass part ($`\widehat{H}_{\text{c.m.}}`$), a radial part ($`\widehat{H}_\rho `$), and an angular part ($`\widehat{h}_\mathrm{\Omega }`$) depending respectively on $`\stackrel{}{R}`$, $`\rho `$, and $`\mathrm{\Omega }`$ sor01 :

$`\widehat{H}`$ $`=`$ $`\widehat{H}_{\text{c.m.}}+\widehat{H}_\rho +{\displaystyle \frac{\mathrm{}^2\widehat{h}_\mathrm{\Omega }}{2m\rho ^2}},`$ (6)
$`\widehat{H}_{\text{c.m.}}`$ $`=`$ $`{\displaystyle \frac{\widehat{p}_R^2}{2Nm}}+{\displaystyle \frac{1}{2}}Nm\omega ^2R^2,`$ (7)
$`\widehat{H}_\rho `$ $`=`$ $`\widehat{T}_\rho +{\displaystyle \frac{1}{2}}m\omega ^2\rho ^2,`$ (8)
$`{\displaystyle \frac{\mathrm{}^2\widehat{h}_\mathrm{\Omega }}{2m\rho ^2}}`$ $`=`$ $`\widehat{T}_\mathrm{\Omega }+{\displaystyle \underset{i<j}{}}V_{ij},`$ (9)
where $`\widehat{T}_\rho `$ and $`\widehat{T}_\mathrm{\Omega }`$ are radial and angular kinetic energy operators. Then the center of mass motion always separates from the relative motion since the $`V_{ij}`$-terms are independent of $`R`$.

We remove the center of mass motion and study the Schrödinger equation for relative coordinates

$`(\widehat{H}\widehat{H}_{\text{c.m.}})\mathrm{\Psi }=E\mathrm{\Psi }.`$ (10)
The adiabatic hyperspherical expansion of the wave function is

$`\mathrm{\Psi }(\rho ,\mathrm{\Omega })`$ $`=`$ $`\rho ^{(3N4)/2}{\displaystyle \underset{\nu =0}{\overset{\mathrm{}}{}}}f_\nu (\rho )\mathrm{\Phi }_\nu (\rho ,\mathrm{\Omega }),`$ (11)
where $`\mathrm{\Phi }_\nu `$ is an eigenfunction of the angular part of the Hamiltonian with an eigenvalue $`\mathrm{}^2\lambda _\nu (\rho )/(2m\rho ^2)`$:

$`\widehat{h}_\mathrm{\Omega }\mathrm{\Phi }_\nu (\rho ,\mathrm{\Omega })`$ $`=`$ $`\lambda _\nu (\rho )\mathrm{\Phi }_\nu (\rho ,\mathrm{\Omega }).`$ (12)
Then eq. (10) leads to a set of coupled radial equations. Neglecting couplings between the different $`\nu `$-channels yields the radial eigenvalue equation:

$`\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{d^2}{d\rho ^2}}+U_\nu (\rho )E_\nu \right)f_\nu (\rho )=0,`$ (13)
$`{\displaystyle \frac{2mU_\nu (\rho )}{\mathrm{}^2}}={\displaystyle \frac{\lambda _\nu }{\rho ^2}}+{\displaystyle \frac{(3N4)(3N6)}{4\rho ^2}}+{\displaystyle \frac{\rho ^2}{b_\text{t}^4}},`$ (14)
where $`E_\nu `$ is the energy and the adiabatic potential $`U_\nu `$ acts as an effective mean-field potential as a function of the hyperradius. This potential consists of three terms, i.e., the external field, the generalized centrifugal barrier, and the angular average of the interactions and kinetic energies. The neglected non-diagonal terms are typically about 1% of the diagonal terms for attractive Gaussian potentials.

We have so far no restriction on the many-body wave function, but include in principle any structure of the system. To choose a convenient form we follow the philosophy in the Faddeev-Yakubovskiĭ formulations fad61 ; yak67 , i.e., the additive decomposition of the wave function reflects explicitly the possible asymptotic large-distance behaviour of cluster subsystems. We expect that two-body correlations are most important and we select the corresponding terms in the decomposition. Higher-order correlations are then essentially neglected. This procedure assumes a very different starting point compared to the Jastrow factorization into products of two-body wave functions bij40 ; din49 ; jas55 . The traditional Jastrow form is expected to be more efficient for large densities while our method is well suited for the low densities encountered for Bose-Einstein condensates.

Emphasizing two-body correlations we therefore decompose the angular wave function $`\mathrm{\Phi }`$ in the symmetric Faddeev components $`\varphi `$

$`\mathrm{\Phi }(\rho ,\mathrm{\Omega })={\displaystyle \underset{i<j}{\overset{N}{}}}\varphi _{ij}(\rho ,\mathrm{\Omega }){\displaystyle \underset{i<j}{\overset{N}{}}}\varphi (\rho ,r_{ij}),`$ (15)
where the last approximation assumes that only relative $`s`$-waves between each pair of particles contribute. Then the coordinate dependence reduces to the distance $`r_{ij}=\sqrt{2}\rho \mathrm{sin}\alpha _{ij}`$. Neglecting higher-order partial waves is justified when the large-distance properties are decisive. The capability of this assumption for large scattering length has been demonstrated for $`N=3`$ by describing the intricate Efimov effect fed93 ; jen97 .

The angular eigenvalue equation (12) can by a variational technique be rewritten as a second order integro-differential equation in the variable $`\alpha _{12}`$ sor01 . For atomic condensates the interaction range is very short compared to the spatial extension of the $`N`$-body system. Then this equation simplifies even further to contain at most one-dimensional integrals. The validity of our approximations only relies on the small *range* of the potential, whereas the scattering length can be as large as desired.

We shall use the finite-range Gaussian potential $`V(r)=V_0\mathrm{exp}(r^2/b^2)`$. Thus we have either overall attractive or overall repulsive potentials depending on the sign of the strength $`V_0`$. It is convenient to measure the strength of the interaction in units of the Born-approximation of the scattering length

$`a_\text{B}{\displaystyle \frac{m}{4\pi \mathrm{}^2}}{\displaystyle d^3\stackrel{}{r}_{kl}V(\stackrel{}{r}_{kl})}={\displaystyle \frac{\sqrt{\pi }mb^3V_0}{4\mathrm{}^2}},`$ (16)
where the last expression is for the Gaussian potential. We use the sign convention that the scattering length $`a_s>0`$ for a purely repulsive potential, such that $`a_sa_\text{B}`$ for $`|a_\text{B}|/b1`$. Thus $`a_s>0`$ for purely repulsive potentials while purely attractive potential can lead to any, positive or negative, value of $`a_s`$ depending on $`V_0`$ and $`b`$. In appendix A is collected the connections between the Gaussian strength measured in $`a_\text{B}/b`$ and the scattering length $`a_s/b`$ for the cases applied in this work. In most of the numerical work we have $`|a_\text{B}|/b`$ close to unity.

## III Angular potentials

The key quantity in the radial equation (13) is the angular eigenvalue $`\lambda `$ obtained from eq. (12). This eigenvalue depends on the number of particles, on the size of the system through the hyperradius, and on the two-body potential through the scattering length. The behaviour of $`\lambda `$ is decisive for the effective potential in eq. (14) which in turn determines the properties of the solutions to eq. (13). We shall therefore first study the dependence of $`\lambda `$ on the parameters in the model. We use the method described in sor02b . The two-body interaction is a simple Gaussian either purely attractive or purely repulsive. This finite-range interaction never produces the collapse at short distance arising from an attractive $`\delta `$-force fed01 . Thus we can as well use attractive potentials with one or more bound states.

### III.1 General eigenvalue behaviour

The angular eigenvalue spectrum coincides with the free spectrum (without interaction) at both small and large hyperradii; for $`\rho =0`$ because all interactions are multiplied by $`\rho ^2`$, see eqs. (6) and (9), and at $`\rho =\mathrm{}`$ because the short-range interaction has no effect at large distances. Thus, perturbation theory for small $`\rho `$ for a Gaussian potential shows that the eigenvalues change from their hyperspherical values $`\lambda _\nu (0)=2\nu (2\nu +3N5)`$, $`\nu =0,1,\mathrm{}`$, as

$`\lambda _\nu (\rho )\lambda _\nu (0)={\displaystyle \frac{mV_0}{\mathrm{}^2}}N(N1)\rho ^2.`$ (17)
If the two-body potential is attractive, but too weak to support any bound state, the eigenvalues reach a minimum as function of $`\rho `$ and then return to one of the finite hyperspherical values. For more attractive potentials there is a one-to-one correspondence between one given two-body bound state of energy $`E^{(2)}<0`$ and one eigenvalue $`\lambda `$ diverging with $`\rho `$ as $`\lambda =2mE^{(2)}\rho ^2/\mathrm{}^2`$. The corresponding structure describes, appropriately symmetrized, one pair of particles in that bound state and all others far apart from the pair and from each other. In addition to this finite number of such negative eigenvalues the hyperspherical spectrum emerges at large distances.

To illustrate we show in fig. 1 a number of possible angular eigenvalues $`\lambda `$ as functions of hyperradius for different potentials. The entirely positive (solid) curve corresponds to a repulsive Gaussian. The diverging (dotted and thick dot-dashed) curves correspond to potentials with one bound two-body state. For our purpose the curves approaching zero for large $`\rho `$ (dashed and thin dot-dashed curves) are the most interesting, since they are crucial for the later description of the condensate. This is true even when the potential has lower-lying bound states corresponding to diverging $`\lambda `$ (thick dot-dashed curve).

The convergence of $`\lambda `$ as $`\rho 0`$ is due to the finite range of the potential and the behaviour depends on the interaction range $`b`$. The deep minima in fig. 1 at small to intermediate distances depend strongly on both the number of particles and the strength of the attraction. They are substantially deeper than reported in sor02 ; sor02b  where one term inadvertently was used with the wrong sign in the numerical examples. This error is only significant for small and intermediate $`\rho `$. Increasing the strength of the attraction always leads to larger negative values of $`\lambda `$. However, at some point one more bound two-body state reveals its presence by changing convergence to zero into a parabolic divergence with $`\rho `$.

For distances much larger than the range of the potential the eigenvalues could as well be computed from a zero-range interaction, i.e., $`4\pi \mathrm{}^2a_s\delta (\stackrel{}{r})/m`$. The hyperharmonic angular wave function should then be appropriate for $`\mathrm{\Phi }`$ and the eigenvalue $`\lambda `$ obtained as the corresponding expectation value. The lowest hyperharmonic is a constant independent of angles and the result is boh98

$`\lambda _\delta (N,\rho )`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3N3}{2}\right)}{\mathrm{\Gamma }\left(\frac{3N6}{2}\right)}}N(N1){\displaystyle \frac{a_s}{\rho }}`$ (18)
$`\stackrel{N1}{}`$ $`{\displaystyle \frac{3}{2}}\sqrt{{\displaystyle \frac{3}{\pi }}}N^{7/2}{\displaystyle \frac{a_s}{\rho }}.`$
This zero-range result is inversely proportional to $`\rho `$ for all hyperradii and consequently with a non-physical divergence when $`\rho 0`$. The only length scale arises from the strength of the $`\delta `$-function. In mean-field calculations this strength is chosen to reproduce the correct scattering length $`a_\text{B}`$ in the Born approximation pet01 ; esr99b . To reach this limit with a Gaussian potential then requires that the $`\delta `$-function is approached while $`a_\text{B}`$ is maintained equal to the desired value of $`a_s`$.

This artificial construction is due to the lack of correlations in mean-field computations where the effective interaction is adjusted to the available Hilbert space. We use finite-range Gaussian potentials and include two-body correlations. Then we expect the large-distance asymptotic behaviour to be described by eq. (18) with the correct scattering length. This tests the efficiency of the simplified structure of the wave function in eq. (15). Mathematically this should result from the structure of the second order integro-differential angular eigenvalue equation sor01 ; sor02b .

Numerically we investigate the asymptotic behaviour of $`\lambda `$ in this context by comparing to the zero-range result $`\lambda _\delta `$ in fig. 2. The convergence to the limiting value is fastest for the smallest value of $`|a_s|`$ (dashed and solid curve) reflecting that the correlations arising for large scattering lengths (dotted line) cannot be accounted for by the zero-range result. This is well understood for three particles where the Efimov effect (very large $`a_s`$) extends correlations in hyperradius to distances around four times the average scattering length fed93 ; jen97 . These effects are not present in the zero-range expectation value contained in $`\lambda _\delta `$. When $`\rho `$ exceeds $`|a_s|`$ by a sufficiently large amount the Efimov effect disappears in $`\lambda `$ and $`\lambda _\delta `$ is approached.

A stronger attraction corresponding to one two-body bound state produces one diverging eigenvalue (figure 1) while the second eigenvalue converges towards $`\lambda _\delta `$ (dot-dashed curve). In fact $`\lambda _1(10,1)`$ almost coincides with the lowest eigenvalue $`\lambda _0(10,0)`$ for the same scattering length but for a potential without bound two-body states (dotted curve).

The numerical deviations from $`\lambda _\delta `$ at large distance is in all cases less than 10%. The asymptotic behaviour is very smooth but still originating in systematic numerical inaccuracies.

These results demonstrate that the scattering length entirely determines the asymptotic behaviour of the potentials. The radial shape of the two-body potential could be Gaussian, square-well, Woods-Saxon, or Yukawa, still the same $`a_s`$ would produce the same angular eigenvalue at sufficiently large distance.

### III.2 $`N`$-dependence

The angular eigenvalues increase rapidly with $`N`$ as seen already from the $`N^{7/2}`$-dependence in $`\lambda _\delta `$. The major variation in magnitude is then accounted for by using this large-distance zero-range result as the scaling unit. We show in fig. 3 a series of calculations for the same two-body interaction for different numbers of atoms. All curves are similar with a systematic increase in the characteristic hyperradius $`\rho _a`$ where they bend over and approach the zero-range result. We then numerically determine this characteristic length $`\rho _a`$ to be proportional to the scattering length and a particular power $`7/6`$ of $`N`$, i.e.,

$`\rho _a(N)|a_s|N^{7/6}.`$ (19)
The quality of this scaling is illustrated in fig. 4, where all curves essentially coincide for distances smaller than $`\rho _a`$. At larger hyperradii the zero-range result of $`+1`$ should be obtained. However, systematic deviations from a common curve is apparent. For each $`N`$ one smooth curve is followed at small and intermediate distances implying that the numerical inaccuracies here are systematic until random fluctuations set in at large $`\rho `$.

The smooth numerical curves can be rather well reproduced by the function

$`\lambda ^{()}(N,\rho )=\lambda _\delta (N,\rho )g^{()}(\rho /\rho _a),`$ (20)
$`g^{()}(x)=g_{\mathrm{}}\left(1e^{x/x_a}\right)\left(1+{\displaystyle \frac{x_b}{x}}\right),`$ (21)
where $`g_{\mathrm{}}=1`$ in accurate calculations. The exponential term is introduced to reproduce the rather fast approach to the asymptotic value as seen in fig. 4. The behaviour at smaller distance, depending on the range of the interaction, is simulated by the $`x_b`$-term. The extreme limit of $`\rho 0`$ is attempted reproduced on the function in eq. (21).

The two groups of computations in fig. 4 are reasonably well reproduced by the parameter sets $`x_a0.74`$, $`x_b2.310^3`$, and $`g_{\mathrm{}}0.8`$ or $`g_{\mathrm{}}0.4`$. These parameters may depend on the scattering length, and we therefore repeated the computation for various $`a_s`$. The best choice of parameters are shown in table 1. We notice that $`g_{\mathrm{}}`$ and $`x_a`$ both are of order unity, and that the fraction $`g_{\mathrm{}}/x_a`$ is almost constant, except for the smallest scattering length. The parameter $`x_b`$, introduced to account for the finite interaction range, is almost equal to $`b/|a_\text{s}|`$. At large hyperradii, where $`\rho \rho _a`$, $`\lambda ^{()}`$ approaches $`g_{\mathrm{}}\lambda _\delta `$. The rather accurate results for $`N=100`$ displayed in fig. 2 confirm that $`g_{\mathrm{}}1`$ by deviating less than $`10\%`$ from $`\lambda _\delta `$ at large hyperradii.

The angular eigenvalue is given by $`g^{()}(x)g_{\mathrm{}}x/x_a`$ for $`x_b\rho /\rho _ax_a`$. Numerical calculations in this intermediate region of hyperradii therefore rather accurately determines the fraction $`g_{\mathrm{}}/x_a1.08`$ as given in table 1. With $`g_{\mathrm{}}=1`$ this implies that $`x_a1/1.080.92`$. The parameters of $`g^{()}(x)`$ in eq. (21) are then given by

$`g_{\mathrm{}}=1,x_a0.92,x_b0.92{\displaystyle \frac{b}{|a_s|}}.`$ (22)
We can compare with the rigorous result for $`N=3`$ jen97  where the angular eigenvalue at large-distance coincide with $`\lambda _\delta `$ in eq. (18), i.e.,

$`\lambda _\delta (N=3,\rho )={\displaystyle \frac{48a_s}{\sqrt{2}\pi \rho }}.`$ (23)
Thus also for $`N=3`$ the universal function $`g^{()}`$ asymptotically approaches $`1`$ for all scattering lengths. More accurate results for $`N=3`$ have been calculated with the zero-range model from fed01b  and are shown in rescaled form in fig. 4. The behaviour is similar to the behaviour of the eigenvalue for the $`N`$-body systems, which confirms the schematic model.

The accuracy of the parametrization in eq. (21) is seen in figs. 5a-d, where the angular eigenvalues are shown in units of $`\lambda ^{()}`$ with the individual set of parameters from table 1. Good agreement is found for $`\rho /\rho _a>x_b`$, except at large hyperradii where the numerical inaccuracy increases with increasing scattering length. Fortunately, the large-distance behaviour is known from analytic considerations and we do not need to rely on numerical computations at these distances. The remaining deviations occur at small hyperradii for $`\rho /\rho _a<x_b`$ or equivalently for $`\rho <N^{7/6}b`$, where the result depends on the radial shape of the two-body interaction.

### III.3 $`N`$-dependence with bound two-body states

In the presence of a bound two-body state of energy $`E^{(2)}`$ one angular eigenvalue eventually diverges at large hyperradii as nie01

$`\lambda ^{(2)}(\rho )={\displaystyle \frac{2m\rho ^2}{\mathrm{}^2}}E^{(2)},E^{(2)}<0.`$ (24)
In the limit of weak binding, or for numerically large scattering lengths, the energy of the two-body bound or virtual state is given by

$`E^{(2)}={\displaystyle \frac{\mathrm{}^2}{ma_s^2}}c,`$ (25)
where $`c`$ approaches unity for large scattering lengths.

We now parametrize the angular eigenvalue by an expression similar to eqs. (20) and (21). The effect of the bound two-body state is only expected to show up at large distances where the behaviour corresponds to eq. (25). The small and intermediate distances resemble the behaviour when no bound state is present. Therefore we arrive at the parametrization

$`\lambda ^{(+)}(N,\rho )=\lambda _\delta (N,\rho )g^{(+)}(\rho /\rho _a),`$ (26)
$`g^{(+)}(x)=x\left(1+{\displaystyle \frac{x_b}{x}}\right)\left({\displaystyle \frac{g_{\mathrm{}}}{x_a}}+c{\displaystyle \frac{4}{3}}\sqrt{{\displaystyle \frac{\pi }{3}}}x^2\right),`$ (27)
with the notation and estimates from eq. (22).

We compare in fig. 6 the parametrization in eqs. (26) and (27) with the computed angular eigenvalues for a potential with one bound two-body state. For the large scattering length ($`a_s/b=100`$) in fig. 6a one smooth curve applies for all the particle numbers; numerical inaccuracies set in at larger hyperradii, which is most obvious for the largest particle numbers. This smooth curve is in a large interval of hyperradii at most deviating by 20% from the parametrized form, and even less than 10% at large hyperradii, before the numerical instability sets in.

For smaller $`a_s`$ ($`a_s/b=+10`$) the deviation at large hyperradii is less than $`1\%`$. The deviation at intermediate distances would decrease by inclusion of a linear term in eq. (27). The smooth curve at small hyperradii is outside the range of validity of the parametrization, i.e., within the range of the two-body potential and then depending on details of the interaction.

The lowest eigenvalue $`\lambda _0`$ diverges at large hyperradius as described by eq. (26). If the two-body potential only has one bound state the second eigenvalue $`\lambda _1`$ is expected to approach zero at large distances as $`\lambda _\delta `$. This pattern should be repeated with more than one bound two-body state, i.e., the first non-divergent angular eigenvalue would behave as $`\lambda _\delta `$ for large $`\rho `$.

We therefore in fig. 6b compare the computed first excited angular eigenvalue with $`\lambda _\delta `$ for different $`N`$. As in fig. 4 we obtain smooth and almost universal curves at large $`\rho `$, where the approach to unity sets in exponentially fast depending on $`N`$, but now much later when $`\rho 10^2\rho _a`$. Clearly a parametrization would also here be possible.

The large-distance asymptotic behaviour of $`\lambda _1`$ now corresponds to an effectively repulsive potential. However, at small and intermediate hyperradii the potential is still effectively attractive ($`\lambda _1<0`$). This attractive region may support a self-bound system located at distances far inside and independent of the confining external field.

This feature is absent in the mean-field description of Bose-Einstein condensation. For overall repulsive potentials corresponding to positive scattering lengths no attractive part is possible. For attractive potentials either a zero-range potential would produce a collapsed wave function and a finite-range potential would not give repulsion at large distance.

### III.4 The properties of $`\lambda ^{(\pm )}`$

The functions $`\lambda ^{(\pm )}`$ coincide when $`\rho \rho _a`$ and depends only on the absolute magnitude of the scattering length $`|a_s|`$. For $`\rho \rho _a`$ the functions differ qualitatively, i.e., $`\lambda ^{()}`$ approaches zero as $`\lambda _\delta `$ while $`\lambda ^{(+)}`$ diverges as $`\rho ^2`$.

At intermediate hyperradii, $`x_b\rho /\rho _ax_a`$, when

$`b{\displaystyle \frac{\rho }{N^{7/6}}}|a_s|,`$ (28)
the angular eigenvalue $`\lambda ^{(\pm )}`$ approaches a constant value $`\lambda _{\mathrm{}}`$, i.e.,

$`\lambda _{\mathrm{}}\lambda _\delta {\displaystyle \frac{\rho g_{\mathrm{}}}{\rho _ax_a}}={\displaystyle \frac{3g_{\mathrm{}}}{2x_a}}N^{7/3}\sqrt{{\displaystyle \frac{3}{\pi }}}1.59N^{7/3}.`$ (29)
This numerical result is in agreement with the following derivation.

The angular eigenvalue for large scattering length $`a_s`$ is independent of hyperradius $`\rho `$ when $`\rho `$ is large compared to the range $`b`$ of the potential but small compared to $`a_s`$. The plateau value $`\lambda _{\mathrm{}}`$ can be estimated as the intersection between two curves at the point $`\rho _a`$. The first curve is the parabollically decreasing $`\lambda (\rho )`$ corresponding to a bound two-body state, i.e., $`\lambda (\rho )=2m\rho ^2E^{(2)}/\mathrm{}^2=2\rho ^2/a_s^2`$, where $`E^{(2)}`$ is given by eq. (25) with $`c=1`$. The second curve is the increasing $`\lambda _\delta (\rho )`$ for an attractive potential ($`a_s<0`$), see eq. (18). Thus $`\lambda _\delta (\rho _a)=\lambda (\rho _a)`$ gives

$`\rho _a`$ $``$ $`\sqrt[3]{{\displaystyle \frac{3}{4}}}N^{7/6}|a_s|,`$ (30)
$`\lambda _{\mathrm{}}(N)`$ $``$ $`\sqrt[3]{{\displaystyle \frac{9}{2}}}N^{7/3}1.65N^{7/3},`$ (31)
which is very close to the numerical results in eqs. (19) and (29).

The symbol $`\lambda _{\mathrm{}}`$ is chosen for this constant, since the $`\rho `$-region where $`\lambda =\lambda _{\mathrm{}}`$ increases proportional to $`|a_s|`$, see eq. (28), and thus extends to infinity for $`|a_s|=\mathrm{}`$. With no bound two-body states ($`a_\text{s}<0`$) the lowest angular eigenvalue approaches zero at larger hyperradii, whereas it diverges towards $`\mathrm{}`$ as $`\rho ^2`$ when a bound two-body state is present ($`a_s>0`$). On the threshold for a two-body bound state $`a_s=\pm \mathrm{}`$ and the angular eigenvalue therefore remains constant.

In sor02  $`\lambda _{\mathrm{}}5N^2`$ was estimated by courageous extrapolation of calculations for $`N=10,20,30`$ and the analytic result for $`N=3`$. The much better estimate in eq. (29) of the large-$`N`$ asymptotics of $`\lambda _{\mathrm{}}(N)`$ increases with a slightly higher power of $`N`$ but with a smaller proportionality factor.

## IV Radial potentials and solutions

The radial equation is the next step in the process of obtaining knowledge about the physical properties of the many-boson system. The angular potentials found in the previous section now enter the effective radial potential in eq. (14) and infer information about the interactions to quantities like energy, size, and structure of the system.

### IV.1 Properties of the radial potential

The radial potential in eq. (14) consists of three terms where the repulsive centrifugal barrier and the confining external field both are positive. The interaction term $`\lambda _\nu `$ can be either repulsive or attractive depending on hyperradius and which eigenvalue we consider. The combination has structure depending on the interaction. For a purely vanishing or repulsive two-body potential we arrive at a simple behaviour qualitatively similar to the non-interacting (dashed) curve shown in fig. 7 for $`N=100`$. All solutions are confined to the region between the infinitely large potential walls at small and large hyperradii.

For a moderately attractive two-body potential a different structure already appears for the lowest angular potential (solid curve in fig. 7). The large-distance behaviour is determined by the trap and is roughly as without interaction, but the barrier at intermediate distance is now finite both in height and width. The barrier height is small compared to the potential at both small and large hyperradii. At smaller hyperradii a rather deep and relatively narrow minimum is present outside a hard core repulsion. The minimum occurs for $`N=100`$ at about 150 times the range of the interaction which corresponds to a mean distance $`(2\rho ^2/N)^{1/2}`$ between each pair of particles of about 15 times the interaction range $`b`$.

With this potential we solve the diagonal radial equation. The solutions can be divided into groups related to either the first or the second minimum. The lowest-lying of the first group of solutions have negative energies. In the model they are truly bound states as they cannot decay into continuum states at large hyperradii sor02 . Their properties are independent of the external trap which only has an influence at much larger distances. These self-bound $`N`$-body states can decay into lower-lying states consisting of various bound cluster states, e.g., diatomic or triatomic clusters. The possibility of self-bound many-body systems, even though the two- and three-body sub-systems are unbound, is also discussed by Bulgac bul02 , who, however, considers the three-body interaction strength as a determining parameter for the properties of the self-bound many-boson system.

The group of states in the higher-lying minimum at larger distance all have positive energies. They are only stable due to the confining effect of the external trap potential. The lowest of these is interpreted as the state of the condensate and indicated by a horizontal line in fig. 7b. This second minimum almost coincides with the minimum of the radial potential arising without any two-body interaction. Thus the structure of the condensate is similar for both positive and negative scattering lengths arising from either attractive or repulsive interactions. However, an attraction produces in addition a series of lower-lying states at smaller hyperradii.

Increasing $`N`$ leaves semi-quantitatively the same features for pure repulsion, whereas an unchanged attraction leads to decreasing barriers at intermediate hyperradius and at some point this barrier vanishes altogether. At the same time the attractive minimum at smaller hyperradius becomes deeper. This in turn leads to an increasing number of bound states in this minimum as function of $`N`$.

As the scattering length increases, the barrier disappears and the effective potential inside the trap has the $`\rho ^2`$-behaviour characteristic for Efimov states, see eq. (14) with $`\lambda =\lambda _{\mathrm{}}`$ of eq. (29). The lowest-lying states are influenced by the details of the two-body interaction and without Efimov features. However, the higher-lying states, located for $`\rho `$-values obeying eq. (28), exhibit the Efimov scaling. They easily become very large and located far outside the minimum responsible for the binding. Only a finite number of bound states is possible due to the confining external field.

These states are many-body Efimov states arising when the two-body scattering length is large. This is also precisely the condition for the three-body Efimov states efi70 ; fed93 . Therefore the many-body Efimov states are embedded in the continua of dimer, trimer and higher-order cluster states. They could be artifacts of the model where only special degrees of freedom are treated. However, these states may also be distinguishable resonance structures which are relatively stable because the particles are very far from each other and the couplings to the continuum states therefore are very weak. So far this remains an open question.

### IV.2 Interaction energy

The total energy of a state in the first minimum are independent of the external field as these states are located at small distances. These states have no analogue in mean-field calculations. In contrast, total energies of the states in the second minimum are dominated by the contribution from the confining field and therefore are rather insensitive to anything else than this field and the corresponding harmonic oscillator quantum numbers. It is then much more informative to compare the interaction energies where the large external field contribution is removed.

In fig. 8 is shown the interaction energy per particle as a function of the particle number for a relatively weak attraction corresponding to a small scattering length. The Gross-Pitaevskii solution exists and the related interaction energy is negative due to the attraction between the particles. A nearly linear behaviour is observed at small particle numbers, since each particle interacts with the $`N1`$ other particles. As $`N`$ increases the mean-field attraction increases and the Gross-Pitaevskii solution becomes unstable when $`N|a_s|/b_\text{t}>0.58`$. This corresponds to $`N=1000`$ with the present parameters of $`|a_s|/b_\text{t}=5810^5`$.

These mean-field interaction energies are in fig. 8 compared to the results obtained with the correlated wave functions from the present formulation. Only a few of the large number of bound states for each $`N`$-value resemble the Gross-Pitaevskii solutions with a radius corresponding to the second minimum. For $`N=20`$ the lowest six states are located in the first minimum. Their interaction energies are large and negative outside the scale of fig. 8. The seventh state is located in the second minimum with an interaction energy very close to the Gross-Pitaevskii result, while the eighth has a positive interaction energy. This feature is repeated for increasing $`N`$, i.e., the lowest state located in the second minimum is similar to the mean-field result and the higher-lying states in this second minimum are less bound. When the mean-field solutions collapse, the correlated solutions remain stable due to the use of a finite-range potential.

The correlated and mean-field interaction energies are remarkably similar when both exist. It may at first appear odd that the mean-field interaction energy is marginally lower than by use of the better suited form of the correlated wave function. The reason is that we compare the mean-field result for an effective interaction which has the correct scattering length in the Born approximation while the correlated solution is obtained for an interaction with the correct scattering length. The mean-field interaction is more attractive to compensate for the limited mean-field Hilbert space. The more revealing comparison is to use the same interaction in both calculations.

We can then compare results for the same $`a_\text{B}/b=0.5`$, i.e., a Gross-Pitaevskii calculation with $`a_s/b=0.5`$ and a Gaussian of $`a_\text{B}/b=0.5`$ corresponding to $`a_s/b=0.84`$. As seen in fig. 8 (dashed curve) now the mean-field energies are much smaller. However, it is remarkable that the correlated solution essentially reproduces the energy of the mean-field calculation where the interaction is renormalized to reproduce the correct energy, but with the wrong wave function. The implication is that the correlated wave function is sufficient to describe the correct structure with the correct interaction. The large-distance average properties are at best obtained in mean-field computations, but all features of correlations are absent by definition.

### IV.3 Definition and size of a condensate

The lowest-lying positive energy solutions located in the second minimum have properties similar to the condensates obtained in Gross-Pitaevskii calculations. The present formulation also provides lower-lying negative-energy states. It is therefore necessary to discuss how to distinguish a condensate state from other (perhaps very unstable) $`N`$-body states.

In mean-field treatments, with repulsive two-body potentials and confining trap potentials, the condensate is uniquely defined as a statistical mixture of single-particle states with the ground state dominating pet01 ; pou02 . A condensate has on average many particles in the lowest single-particle state. This many-body state is only stable against total fragmentation due to the trap and as such first of all determined by the properties of the trap. Even with the trap the many-body state is still at best only approximately stationary due to the neglected degrees of freedom which allow energetically favored di-, tri-, and multi-atomic cluster states. This instability is also an experimental fact seen by permanent loss of trapped atoms, e.g., in recombination processes don01 .

Without any two-body interaction the properties of the many-body system is determined by the thick, dashed potential curve in fig. 7. Then we can easily identify the condensate as a state in this potential where the dominating component for finite temperature is the ground state. Including attractive two-body interactions (full curve) the deep minimum at small hyperradius is produced. Then the corresponding ground state, located in this minimum, has nothing to do with a condensate. The density is so high that couplings to other degrees of freedom would develop higher-order correlations and processes like three-body recombinations would quickly destroy the single-atom nature of the gas. This $`N`$-body ground state does not show the signature of a Bose-Einstein condensate, where many particles occupy one single-particle level.

The formulation in the present work does not use the concept of single-particle levels. Therefore we cannot talk about a statistical distribution of particles with the majority in the lowest state. However, we can talk about a many-particle system described as a superposition of many-body eigenstates, where the lowest states are favored in thermal equilibrium. To clarify we can think of a quantum state $`\mathrm{\Psi }`$ as a superposition of different eigenstates $`\mathrm{\Psi }_n(\rho ,\mathrm{\Omega })`$ in eq. (11) given by

$`\mathrm{\Psi }(\rho ,\mathrm{\Omega })={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}c_n\mathrm{\Psi }_n(\rho ,\mathrm{\Omega })`$
$`=\rho ^{(3N4)/2}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}c_n{\displaystyle \underset{\nu =0}{\overset{\mathrm{}}{}}}f_{\nu ,n}(\rho )\mathrm{\Phi }_\nu (\rho ,\mathrm{\Omega }),`$ (32)
with the normalization $`_n|c_n|^2=1`$. A condensate must be sufficiently large to exceed a certain minimum interparticle distance, $`d_\text{c}`$, below which the atoms are too close and recombine very fast. This distance depends on the scattering length and on the number of particles. Therefore, in our formulation the stationary states cannot be characterized as a condensate if $`r_{ij}^2d_\text{c}^2`$ for this wave function.

We define one of the stationary states in this model as the “ideal condensate” state, i.e., the state of lowest energy with

$`r_{ij}^2_{n_\text{c}}d_\text{c},`$ (33)
characterized by $`n=n_\text{c}`$. This state is dominated by the component in the lowest adabatic potential although not necessarily the states of lowest energy, which might have an average particle distance less than $`d_\text{c}`$. The appropriate of these excited states depends on the number of particles and on the scattering length. The ideal condensate is then characterized by one dominating component, with $`|c_{n_\text{c}}|1`$ and $`|c_{nn_\text{c}}|1`$.

If $`d_\text{c}`$ is significantly smaller than the trap length $`b_\text{t}`$, then the state of lowest energy located in the second minimum can be identified as the condensate. This state is characterized by a radial wave function $`f(\rho )`$ with the root mean square radius $`\rho ^2`$ approximately equal to the hyperradius at the second minimum of the lowest adiabatic potential $`U_0(\rho )`$.

To be specific we show in fig. 9 the root mean square interparticle distance given by $`r_{ij}^2_n=2\rho ^2_n/(N1)`$ for the lowest excited states (labeled by $`n`$) in the potential of fig. 7. All states with $`n58`$ have negative energy and $`20b(r_{ij}^2_n)^{1/2}100b`$, which implies that the particles are separated more than their interaction range. Whether these average distances allow qualification as condensates depends on the decay rate of these states.

When $`n59`$ the energies are positive and the average particle distance now suddenly exceeds $`2000b`$. In fact we now find $`r_{ij}^23b_\text{t}^2`$ which approximately is obtained in the limit of a non-interacting gas. These states probably qualify as condensates. Their interaction energies are in fig. 8 compared to the Gross-Pitaevskii values. The discontinuity at $`n=58,59`$ is due to the intermediate barrier. Decreasing and eventual removal of the barrier would smear out this abrupt change of size. Some of the negative-energy states could then extend very far out and in fact have sizes comparable to the trap length. This investigation could be repeated for the higher adiabatic potentials, still neglecting the couplings. The same pattern is obtained with fewer states of small interparticle distance.

## V Decay rates

The condensate is unstable due to the neglected couplings into other degrees of freedom. The condensate therefore has to be located at relatively large distances. The decisive radial potentials are sensitively depending on the scattering length. In fig. 10 we illustrate the different behaviour by using the angular eigenvalues parametrized through eqs. (14), (20), (21), (26), and (27). In fig. 10a the scattering length is relatively small and a large barrier separates the outer minimum from the inner region. By increasing the scattering length the barrier decreases first into a relatively flat region as in fig. 10b and then disappears completely as in fig. 10c when the trap length is exceeded. With these potentials we can now discuss various decay processes, i.e., three-body recombination into dimers, macroscopic tunneling through the barrier and macroscopic collapse after sudden removal of the barrier.

### V.1 Three-body recombination

Bound state dimers can be formed by a three-body process where the third particle ensures conservation of energy and momentum. The number of these three-body recombination (rec) events per unit volume and time can be estimated by the upper limit given in nie99 ; bed00 :

$`\nu _{\text{rec}}=67.9{\displaystyle \frac{\mathrm{}|a_s|^4n^3}{m}},`$ (34)
where $`n`$ is the density of the gas. This expression can be converted into an estimate of the recombination rate for a given hyperradius $`\rho `$. With the volume $`𝒱=N/n`$, the relation between density and mean distance $`1/n=4\pi r_{ij}^2^{3/2}/3`$, and $`r_{ij}^2=2\rho ^2/(N1)`$ obtained from eq. (1), the total recombination rate becomes

$`{\displaystyle \frac{\mathrm{\Gamma }_{\text{rec}}}{\mathrm{}}}=\nu _{\text{rec}}𝒱0.5{\displaystyle \frac{\mathrm{}|a_s|^4N^4}{m\overline{\rho }^6}},`$ (35)
where the mean square average $`\overline{\rho }`$ is defined as $`\overline{\rho }^2\rho ^2`$. In the spirit of the adiabatic hyperspherical expansion method we use $`\overline{\rho }`$ as a classical parameter. The recombination rate increases rapidly with decreasing $`\overline{\rho }`$, as indicated by the vertical arrows in fig. 10.

The recombination time $`T_{\text{rec}}`$ is defined by $`N(t)=N(0)\mathrm{exp}(t/T_{\text{rec}})`$, where $`N(t)`$ is the number of remaining atoms. We then obtain $`\mathrm{\Gamma }_{\text{rec}}/\mathrm{}dN/dt=N/T_{\text{rec}}`$ and

$`T_{\text{rec}}={\displaystyle \frac{N\mathrm{}}{\mathrm{\Gamma }_{\text{rec}}}}={\displaystyle \frac{2m\overline{\rho }^6}{\mathrm{}|a_s|^4N^3}}={\displaystyle \frac{m\overline{r}_{ij}^6}{4\mathrm{}|a_s|^4}},`$ (36)
where we used $`\overline{r}_{ij}=\sqrt{2/N}\overline{\rho }`$. The final expression for $`T_{\text{rec}}`$ is independent of $`N`$. Since the condensate has to form in the external trap it is reasonable to define stability against recombination by $`T_{\text{rec}}T_\text{c}2\pi /\omega `$, where $`T_\text{c}`$, is the oscillator time. With $`1/\omega =mb_\text{t}^2/\mathrm{}`$ and eq. (36) we get stability when $`\overline{r}_{ij}\sqrt[6]{8\pi }|a_s|^{2/3}b_\text{t}^{1/3}=d_\text{c}`$, which provides $`d_\text{c}`$ introduced in section IV.3. In units of $`b_\text{t}`$ we obtain

$`{\displaystyle \frac{d_\text{c}}{b_\text{t}}}=\sqrt[6]{8\pi }\left({\displaystyle \frac{|a_s|}{b_\text{t}}}\right)^{2/3}.`$ (37)
Thus for $`|a_s|/b_\text{t}1`$ also $`d_\text{c}/b_\text{t}1`$. The average particle distance $`\overline{r}_{ij}`$ for a state located in the second minimum is of the order $`b_\text{t}`$ and therefore $`\overline{r}_{ij}d_\text{c}`$, i.e., for these states $`\overline{r}_{ij}`$ is larger than the critical stability length $`d_\text{c}`$. These states then qualify as condensates. For <sup>87</sup>Rb atoms ($`a_s100`$ a.u.) trapped in a field of $`\nu 100`$ Hz we obtain $`T_{\text{rec}}7`$ days.

### V.2 Macroscopic tunneling

The second decay process is related to macroscopic tunneling through the barrier, as indicated in fig. 10b. The present model provides stationary eigenstates (within the allowed Hilbert space) which by definition are time independent. Thus, strictly the states do not tunnel through the barrier, but an exponentially small tail extends to small hyperradii. All particles in this tail would immediately recombine into molecular clusters, because the density is very large in the inner region (both $`\overline{\rho }`$ and $`r_{ij}`$ are small). The rate of this two-step decay, i.e., tunneling through the barrier and subsequent recombination, can be computed as the knocking rate multiplied by the transmission coefficient, which is a measure of the ratio of the probabilities at the turning points inside and outside the barrier. The rate of recombination due to macroscopic tunneling can then be estimated semi-classically as in boh98  by

$`{\displaystyle \frac{\mathrm{\Gamma }_{\text{tun}}}{\mathrm{}}}`$ $``$ $`{\displaystyle \frac{N\nu }{1+e^{2\sigma }}},`$ (38)
$`\nu `$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}\sqrt{{\displaystyle \frac{1}{m}}{\displaystyle \frac{d^2U(\rho )}{d\rho ^2}}|_{\rho _{\text{min}}}},`$ (39)
$`\sigma `$ $`=`$ $`{\displaystyle _{\rho _{\text{in}}}^{\rho _{\text{out}}}}𝑑\rho \sqrt{{\displaystyle \frac{2m}{\mathrm{}^2}}\left[U(\rho )E\right]},`$ (40)
where the factor $`N`$ is needed to give the total number of recombined particles. Here $`\rho _{\text{min}}`$ is the position of the second minimum of $`U`$ and $`\rho _{\text{in}}`$ and $`\rho _{\text{out}}`$ are the classical turning points of the barrier.

The barrier depends strongly on the combination $`N|a_s|/b_\text{t}`$ boh98 ; sor02 . When $`N|a_s|/b_\text{t}1`$ the barrier is large and the very small rate can be estimated through eqs. (38), (39), and (40). The WKB action integral is

$`\sigma {\displaystyle \frac{3}{2}}N\mathrm{ln}\left({\displaystyle \frac{b_\text{t}}{N|a_s|}}\right).`$ (41)
The barrier is absent when $`N|a_s|/b_\text{t}0.53`$. Close to, but before reaching, this threshold of stability the WKB-exponent can be approximated by

$`\sigma 1.7N\delta _s,\delta _s1{\displaystyle \frac{N|a_s|/b_\text{t}}{0.53}},`$ (42)
which is valid when $`\delta _s`$ is close to zero.

The barrier is observed to vanish when $`N|a_s|/b_\text{t}0.53`$ rob01 ; don01 , which due to the factor of $`N`$ implies that $`|a_s|/b_\text{t}1`$. Therefore close to this threshold we have for a condensate in the second minimum that $`\overline{r}_{ij}b_\text{t}d_\text{c}`$, i.e., the three-body recombination does not limit the stability. In the limit $`\sigma 1`$ we get explicitly

$`{\displaystyle \frac{\mathrm{\Gamma }_{\text{rec}}}{\mathrm{\Gamma }_{\text{tun}}}}{\displaystyle \frac{1}{7.0N^4}}1`$ (43)
implying that the macroscopic tunneling process dominates. With $`\sigma 1`$ we obtain that $`\mathrm{\Gamma }_{\text{tun}}/\mathrm{}=N/T_{\text{tun}}0.5N\nu `$, which for $`\nu 100`$ Hz corresponds to a macroscopic tunneling time of $`10`$ ms. This is much faster than the three-body recombination time when the barrier is small ($`\sigma 1`$), i.e., $`T_{\text{rec}}T_{\text{tun}}`$, see eq. (43).

The three-body recombination rate is in fig. 11 shown as a function of hyperradius (solid curve) and compared with the macroscopic tunneling rate (dashed curve) where all particles in the condensate simultaneously disappear. At small hyperradii the three-body recombination rate is clearly much larger than the macroscopic tunneling rate, whereas the opposite holds for large hyperradii. For the parameters in fig. 11 we find that the two time-scales are roughly equal around the second minimum where the condensate is located.

However, the tunneling rate depends strongly on the barrier through the combination $`N|a_s|/b_\text{t}`$. Varying either of the three quantities would then move the tunneling rate up or down in fig. 11. For a larger barrier the condensate would only decay by direct recombination. For a smaller barrier macroscopic tunneling would dominate and the condensate would decay by “collective” recombination of all particles in a very short time interval.

When a few particles recombine into dimers and leave the condensate, the remaining system is no longer in an eigenstate of the corresponding new Hamiltonian. An adiabatic adjustment of Hamiltonian and wave function could then take place. Since fewer particles and unchanged $`a_s`$ and $`b_\text{t}`$ means a larger barrier, the stability against macroscopic tunneling of the new system is therefore increased.

This stabilization by particle “emission” could also be the result of the recombination in the macroscopic tunneling process if the time-scale for recombination at the relevant small distances is longer than the adiabatic adjustment time. In a possible development first a number of particles are emitted, the adjustments follow, and a larger barrier appears which traps and stabilizes the part of the initial wave function in the second minimum. However, now the condensate contains fewer particles.

### V.3 Macroscopic collapse

These decay scenarios are open for direct experimental investigations since the interaction can be changed in an experiment by using the Zeeman splitting to tune to a Feshbach resonance ino98 ; rob01 ; don01 . An initial value of the scattering length (corresponding to a stable condensate in the second minimum) can almost instantaneously be changed to a value where the barrier is removed. The initial wave function for the condensate is now no longer a stationary state in the new potential.

If we assume that the only excitations are the degrees of freedom contained in the lowest new hyperspherical potential with $`s`$-waves, we can use the sudden approximation and expand on the corresponding eigenfunctions. The most important levels are then the lowest-lying positive energy states with energies comparable to the initial condensate. The time-scale for the time evolution of the initial state in the new potential is then determined by the energy differences between such levels. These states of positive energy and large spatial extension confined by the trap are roughly separated by the oscillator quantum of energy $`\mathrm{}\omega `$. The corresponding rate for populating smaller distances with the consequence of immediate recombination is then crudely estimated to be

$`{\displaystyle \frac{\mathrm{\Gamma }_{\text{col}}}{\mathrm{}}}{\displaystyle \frac{1}{T_{\text{col}}}}{\displaystyle \frac{\omega }{2\pi }}.`$ (44)
Experimentally don01  this macroscopic collapse time is verified to be of the order $`1/\omega `$, typically a few milliseconds, as given by the external trapping field.

This macroscopic collapse time is shorter than the macroscopic tunneling time for the parameters of the system in fig. 11. The motion in the potential is fast or slow compared to the recombination time for distances in the first or second minimum, respectively. The time evolution after the sudden removal of the barrier could then be a macroscopic collapse towards smaller hyperradii where dimers and trimers are “emitted” and the barrier begins to appear. The part of the wave function trapped at large distances in the second minimum can then stabilize into a condensate with fewer particles. The time-scale for these processes should then be between the macroscopic collapse time and the recombination time at the second minimum. Possibly other time scales due to the neglected degrees of freedom (angular momentum, clusterization, etc.) could be present in the full study of the dynamics of a many-boson system.

## VI Summary and conclusion

Correlations in a system of $`N`$ identical bosons of low density are described by use of a hyperspherical adiabatic expansion. The wave function is decomposed in additive Faddeev-Yakubovskiĭ components, where each term is related to one pair of particles and only $`s`$-waves are included. The adiabatic potentials are only weakly coupled and we investigate structures where only the lowest contributes. We use a finite-range purely attractive or purely repulsive Gaussian interaction and extract general properties of the lowest angular eigenvalue.

We establish universal scaling relations for the radial potential for arbitrary scattering length and particle number. These scaling rules are valid for large and intermediate distances where the particles on average are outside the range of the interaction. Only the short-distance behaviour is influenced by the choice of interaction potential.

We parametrize the model-independent part of the effective radial potential in a simple form with an interaction part, a centrifugal barrier term and a contribution from the external field. This potential diverges at small distances due to the centrifugal barrier and at large distances due to the confining external field. The two minima are generally separated by a barrier. The deepest minimum at small to intermediate distances supports self-bound $`N`$-body systems where the density is much larger than for a Bose-Einstein condensate. The second minimum at a much larger distance allows solutions with properties characteristic of a condensate. We distinguish by formulating a definition of a condensate in this context.

We compare properties of the correlated structures with those of the zero-range mean-field solutions. The large-distance asymptotic behaviour is found numerically to reproduce the mean-field result for a zero-range interaction renormalized to give the correct scattering length in the Born approximation. This is remarkable since the correct scattering length for the Gaussian potential is far from the Born approximation. Thus the different terms in the second-order integro-differential equation conspire to produce this large-distance result, which is rigorously established for three particles and on general grounds also expected for many particles. The choice of wave function is then a posteriori shown to be sufficient.

The stability of the condensate is limited by decay into lower-lying many-body cluster states reached by processes where three-body recombination resulting in bound dimers is very prominent. We compute various rates of decay and discuss the time-scales involved. The bare three-body recombination process is strongly scattering length and density dependent and therefore increases dramatically when the wave packet moves from the second minimum to smaller distances. An intermediate barrier would only allow quantum tunneling followed by a macroscopic collapse. When this barrier is very small by choice of parameters the macroscopic tunneling rate would dominate. When the interaction is changed during an experiment and the barrier is totally removed the already created condensate would collapse and a number of cluster configurations would appear. Stability may subsequently be automatically restored and a new condensate created with fewer particles.

In conclusion, we have discussed properties of condensates and extracted universal scaling relations. We have focused on the effects of correlations for large scattering lengths where the mean-field approximation breaks down. Finally we investigated time-scales for various decay mechanisms limiting the stability of the condensate. The parametrized potentials allow independent investigations without the full numerical machinery. More general $`N`$-body structures are studied than the simple condensates.

## Appendix A Numerical details

The angular equation can be scaled by using the potential range $`b`$ as the unit length sor02b . The only interaction parameter is then the Born approximation to the scattering length in this unit $`a_\text{B}/b`$. The only length coordinate is then $`\rho /b`$. All physical quantities are functions of such dimensionless ratios.

The $`s`$-wave two-body scattering length is the node of the zero-energy solution to the two-body Schrödinger-equation, i.e., $`u(r)(ra_s)`$. Table 2 shows the scattering length $`a_s`$ for different potential strengths $`a_\text{B}`$, see eq. (16). The Born-approximation equals the correct scattering length only in the limit of weak attraction, where the scattering length $`a_s`$ is much smaller than the range of the interaction $`b`$.

To exemplify, in experimental work <sup>87</sup>Rb atoms with a scattering length of $`a_s=100`$ a.u. are trapped in an external trap of frequency $`\nu =100`$ Hz boh98 . Assuming an interaction range around $`b=1`$ nm we obtain $`a_s/b=5.29`$, $`b_\text{t}/b=1442`$. This can be modelled by a Gaussian two-body interaction with $`a_\text{B}/b=1.5`$, where the lowest solution corresponds to a two-body bound state and the next accounts for the properties of the condensate.


text:
# On the center-vortex baryonic area law

## I Introduction

In this brief report we correct an error in an earlier work of the author co96 , claiming that in center-vortex QCD (that is, for gauge group $`SU(3)`$) the area law for the $`SU(3)`$ baryonic Wilson loop (Fig. 1) involved the sum of three extremal areas, each spanning a pair of quark world lines. We call this sort of area law a $`\mathrm{\Delta }`$ law, and the corresponding area $`A(\mathrm{\Delta })`$. We now claim that the correct law is a $`Y`$ law, with an extremal area $`A(Y)`$ composed of three areas spanning the quark world lines and a central Steiner line. (A Steiner line is a line where three surfaces meet at angles, generically $`2\pi /3`$, chosen to minimize the total area of the three surfaces.) Either the $`\mathrm{\Delta }`$ law or the $`Y`$ law accords with the principle, derived co03  from center-vortex theory, that area laws in general involve extremal areas. The $`Y`$ law leads to an absolute minimum of the area spanning three fixed quark world lines, but the $`Y`$ law always has co96  higher action than that of the $`\mathrm{\Delta }`$ law, because in $`SU(3)`$ the areas have different coefficients in the VEV $`W`$ of the Wilson loop :

$$W\mathrm{exp}[\frac{K_FA(\mathrm{\Delta })}{2}]\mathrm{or}\mathrm{exp}[K_FA(Y)].$$
(1)
It is also worth noting that, for group-theoretical reasons, lowest-order perturbative gluon exchange leads to (short-range) forces solely of the $`\mathrm{\Delta }`$ pattern co77 .

Recent lattice simulations of the heavy-quark baryonic area law, as far as the author knows, are in accord with the $`Y`$ law suga ; afj ; ibss . (Some of these workers had also found aft  evidence for a $`\mathrm{\Delta }`$ law, but this work is apparently superseded by that in afj .) The numerical distinction between the $`Y`$ law and the $`\mathrm{\Delta }`$ law is small, and so high accuracy is needed; it could be that the last word in lattice simulations has not yet been said. Still, the present lattice evidence seems to favor the $`Y`$ law.

Section II points out the error in co96  and corrects it for $`SU(3)`$, using the baryonic Stokes’ theorem of co96 . Section III extends this work to $`SU(N)`$. Section IV contains conclusions.

## II The baryonic area law for SU(3)

The correct parts of the arguments of co96  will be briefly summarized, and the error of this reference corrected. The baryonic Wilson loop $`W`$ is:

$$W=\frac{1}{6}ϵ_{abc}ϵ_{a^{}b^{}c^{}}U(x,y;1)_{aa^{}}U(x,y;2)_{bb^{}}U(x,y;3)_{cc^{}}$$
(2)
where each $`U`$ is an ordered path integral from the point $`x`$ where the three quark world lines originate to the point $`y`$ where they rejoin:

$$U(x,y;j)_{kk^{}}=[P\mathrm{exp}_{\mathrm{\Gamma }(j)}𝑑z^\alpha A_\alpha (z)]_{kk^{}}.$$
(3)
Here $`\mathrm{\Gamma }(j)`$ is the path from $`x`$ to $`x`$ along quark world line $`j`$, and we use the standard anti-Hermitean matrix representation of the gauge potential, which has the gauge coupling incorporated:

$$A_\alpha (z)=\frac{g\lambda _j}{2i}A_\alpha ^j(z).$$
(4)
This Wilson loop is to be averaged over all center-vortex gauge-potential configurations. We idealize the center vortices to random closed two-surfaces ($`d=4`$) or closed lines ($`d=3`$) of infinitesimal thickness, finite persistence length, and finite density per unit area. (In reality the vortices have finite thickness; this gives rise only to perimeter-law corrections.) For $`SU(3)`$ there is only one kind of vortex, plus its anti-vortex; these can be distinguished by assigning an orientation to the closed surface (line) representing the vortex in $`d=4(3)`$. This orientation can be considered as equivalent to the direction of Abelian magnetic field lines lying in the vortex surface (line), as indicated in the figures. Just as in co03  this leads to the argument that the VEV $`W`$ has a form as in Eq. (1) which is exponential in an extremal area. The extremal area for the $`\mathrm{\Delta }`$ law is the sum of three areas, each extremal for the three pairs of quark lines that it spans, and for the $`Y`$ law it is the absolute minimum area spanning the three lines jointly.

For an infinitesimally-thick center vortex, we can choose the corresponding gauge potential as a pure singular gauge:

$$A_\alpha (z)=V(z)_\alpha V^1(z).$$
(5)
In particular, up to a regular gauge transformation $`V`$ can be chosen to be of the form

$$V(z)=\mathrm{exp}[2\pi iQ\beta (z)]$$
(6)
where $`Q`$ is the diagonal matrix (1/3,1/3,-2/3) and $`\beta (z)`$ gives rise to an Abelian gauge potential:

$$\stackrel{~}{A}_\alpha (z)=_\alpha \beta (z)$$
(7)
that is singular on a closed surface and has unit magnetic flux. For the anti-vortex, $`QQ`$.

For this pure gauge potential the path integrals of Eq. (3) are:

$$U(x,y;j)=V(x)V^1(y)\mathrm{exp}[2\pi iQ\mathrm{\Lambda }(j)];\mathrm{\Lambda }(j)=_{\mathrm{\Gamma }(j)}𝑑z^\alpha \stackrel{~}{A}_\alpha (z).$$
(8)
This path integral may or may not depend on the path $`j`$; if it does not, then the Wilson loop of Eq. (2) is just $`detU=1`$. If it does, then the vortex is linked to the Wilson loop and further analysis is needed.

Possible path dependence can be understood from the non-Abelian Stokes’ theorem for baryonic Wilson loops co96 . As explained there, the appropriate surface to integrate over, for $`SU(3)`$, is three-bladed. Fig. 2 shows an intermediate step in the standard nast  conversion of a non-Abelian contour integral into a surface integral and Fig. 3 shows how this loop is changed into a surface by using the construction of Fig. 2 on all three blades of a baryonic Wilson loop. The central line in Fig. 3, marked 78, gives a factor which can be co96  reduced to $`detU(78)=1`$. Because Stokes’ theorem expresses a line integral as a surface integral, it is necessary that the surface form of the integral not depend on the detailed form of the spanning surface, as long as it spans the Wilson loop contour and has no other boundaries. This will be true in our center-vortex picture because of quantization of the magnetic flux carried by the center vortices.

As co96  shows, for a given center vortex of the form in Eqs. (6, 7) the baryonic Wilson loop is, from Eq. (2):

$$W=\frac{1}{3}\mathrm{exp}[\frac{2\pi i}{3}\{\mathrm{\Lambda }(1)+\mathrm{\Lambda }(2)2\mathrm{\Lambda }(3)\}]+\mathrm{cycl}.\mathrm{perm}.$$
(9)
We now deviate slightly from the path of co96 , where the development of Eqs. (28-33) of this reference are in error. Evidently the $`\mathrm{\Lambda }(i)`$ are not integers, although the difference of any two of them is, and it is useful to express the contents of the curly brackets in Eq. (9) as differences of integers. Consider Fig. 4(a), showing (in $`d=3`$, for clarity) a center vortex linked to quark line 1. This linkage is expressed as a (signed) intersection by the vortex of a surface which we will call $`\mathrm{\Gamma }(01)`$, which spans quark line 1 and the central line, labeled 0. The link number, denoted $`Lk(01)`$, is the sum of intersection numbers; as shown in the figure, the link number is 1. As mentioned above, the line integral of the Wilson loop must be independent of the position of the central line, and so precisely the same topology must be expressed in Fig. 4(b), where $`Lk(01)=0,Lk(02)=Lk(03)=1`$, leading to a total of $`Lk(02)+Lk(03)`$=-2. Although the link numbers for these configurations are different, they are equal mod(3), and this is all that counts for calculating the Wilson loop, as we will soon see. With these link numbers we can express differences of the $`\mathrm{\Lambda }(i)`$ as differences of link numbers. For example, $`\mathrm{\Lambda }(1)\mathrm{\Lambda }(2)=Lk(01)Lk(02)`$, since either expression is the link number of the closed loop $`1\overline{2}`$, in which line 2 is traced in the opposite direction. This link number is found by summing the signed intersection number of a vortex with this loop $`1\overline{2}`$. So then we express the Wilson loop phase in terms of integer differences:

$$\mathrm{\Lambda }(1)+\mathrm{\Lambda }(2)2\mathrm{\Lambda }(3)=Lk(01)+Lk(02)2Lk(03)[Lk(01)+Lk(02)+Lk(03)]\mathrm{mod}(3).$$
(10)
We are almost through. From Eqs. (9,10) the Wilson loop VEV is:

$$W=\mathrm{exp}[\frac{2\pi i}{3}\{Lk(01)+Lk(02)+Lk(03)\}].$$
(11)
We claim that each of the phase factors is statistically independent of the other two, as simple geometrical arguments reveal, so that:

$$W=\mathrm{exp}[\frac{2\pi i}{3}Lk(01)]\mathrm{exp}[\frac{2\pi i}{3}Lk(02)]\mathrm{exp}[\frac{2\pi i}{3}Lk(01)].$$
(12)
But each of the factors in brackets is itself the expectation value of a $`q\overline{q}`$ Wilson loop, where, for example, the 01 factor involves the Wilson loop formed from quark line 1 and the central line 0 in Stokes’ theorem (taken to be directed from $`y`$ to $`x`$). As such, by previous arguments co03 , each bracketed factor is given by:

$$\mathrm{exp}[\frac{2\pi i}{3}Lk(01)]=\mathrm{exp}[K_FA_{min}(01)]$$
(13)
where $`A_{min}(01)`$ is the minimum area associated with the 01 loop. So the baryonic Wilson loop is just:

$$W=\mathrm{exp}[K_F\{A_{min}(01)+A_{min}(02)+A_{min}(03)\}].$$
(14)
It only remains to point out that the central line 0, so far undetermined, must be a line which minimizes the sum of areas in Eq. (14), that is, a Steiner line. And so we find the $`Y`$ area law instead of the $`\mathrm{\Delta }`$ law.

## III Extension to SU(N)

Stokes’ theorem for these higher $`N`$ is the obvious generalization co96  of the $`SU(3)`$ theorem. We construct $`N`$ sheets, each bounded by a quark world line and a central line, and form plaquettes on them as we did for $`SU(3)`$. This construction must be reconciled with the requirement that the baryonic area must be extremal. It is possible to find a local minimal area consisting of these $`N`$ sheets meeting along a single line, by appropriate choice of the line and sheet areas; the minimum is with respect to all areas spanning the quark lines and the mid-line. If this is the appropriate extremal area, then the Wilson-loop VEV is of the form $`\mathrm{exp}[K_FA_{min}(0i)]`$, generalizing Eq. (14). However, the global minimum is constructed from surfaces joined by Steiner lines. An example is shown in Fig. 5(a) for $`N=5`$. Fig. 5(b) shows the deformation of the 5-bladed Stokes’ theorem surface to conform to the extremal area. Note that the surfaces whose projection is the line AB or the line BC are double, so their contribution to the phase factor in the Wilson loop for a given vortex is doubled. That is, if a vortex of magnetic flux $`2\pi J/5`$, for $`J=\pm 1,\pm 2`$ is linked to quark line 1, as expressed in piercing the surface 1A of Fig. 5(b), the phase factor is $`\mathrm{exp}[2\pi iJ/5]`$. The same topology is expressed through piercings (with opposite signs) of surfaces 5B, BC, and 2A. The piercing is counted twice for surface BC, so the alternative phase factor is $`\mathrm{exp}[4\pi iJ/5]`$, which is equal to the first phase factor. Unfortunately, in this case it is not possible to give an explicit formula for the Wilson-loop VEV, since in $`SU(N)`$ there are \[$`N/2`$\] different vortices, each with an anti-vortex, where \[$`n`$\] is the integer part of $`n`$ and in general, these vortices occur with different areal densities. The baryonic Wilson-loop formula then is connected to the problem of so-called k-string tensions (string tensions for antisymmetric representations formed from the product of $`k`$ fundamental representations) kstring ; we will discuss it further elsewhere.

## IV Conclusions

We have corrected an error in co96 , leading to the conclusion that the $`SU(3)`$ baryonic area law is the $`Y`$ law, not the $`\mathrm{\Delta }`$ law. The error lay in an incomplete specification of the link numbers $`Lk(0i)`$, and an improper specification of the statistical independence of these link numbers.

Our present result can be extended in principle to other $`SU(N)`$ gauge groups, at least formally, although the final result may depend on unknown areal densities of the several different types of center vortices that occur for $`N>3`$. This extension uses the baryonic Stokes’ theorem for larger $`N`$, where $`N`$ sheets span the quark world lines and a central line. There is a local minimum area for such a configuration, and a simple baryonic area law which generalizes the $`Y`$ law in $`SU(3)`$. It is also possible to deform this $`N`$-bladed surface to become the global minimal area spanning the $`N`$ quark world lines in a baryon, with the aid of Steiner lines; no explicit form for the baryonic area law is known at present.