text:
# Structure and energetics of Ni clusters with up to 150 atoms

## Abstract

We present a method (the Aufbau/Abbau method) for optimizing the structure of a whole series of clusters without making any assumptions on the structure. Subsequently, the method is combined with the embedded-atom method in determining the structure of the two energetically lowest isomers of NiN clusters with $`N`$ up to 150. Finally, various analytical descriptors are introduced that are used in studying the overall shape of the clusters, their structure and stability, and possible growth and dissociation processes.

About two decades ago the experimental observation and theoretical description of so-called magic numbers in clusters of simple metals (see \[deheer, ; brack, \]) provided an excellent example of how experiment and theory could supplement each other in the development of the understanding of materials properties. The theoretical studies for any but the smallest clusters were based on the jellium model, i.e., the internal structure of the cluster was completely neglected and the core electrons and the nuclei were smeared out to a homogeneous jellium in whose field the valence electrons moved. Moreover, in the first studies the jellium was supposed to be spherical and to have sharp boundaries, although these assumptions in later studies have been relaxed partially (see \[brack, \]).

For theoretical studies of any other metal clusters MN than those of simple metals, one encounters the serious problem that when the valence electrons are anything but completely delocalized, the (spherical) jellium model may not provide a realistic approximation. This is, e.g., the case for transition-metal clusters where the nearest-neighbor interactions between the atoms are important, so that the precise structure of the cluster as well as the inhomogeneities of the electron density are to be included in any realistic description. Then, furthermore, the fact that the number of structural degrees of freedom ($`3N6`$ or $`3N5`$) together with the fact that the number of metastable structures grows very fast with $`N`$ makes detailed systematic studies of how the properties depend on $`N`$ for any but the absolutely smallest clusters extremely demanding. On the other hand, the further development of this field requires a precise understanding of how the properties of interest depend on the structure, size, and composition of the clusters, which, unfortunately, is very difficult to achieve theoretically or experimentally.

Here we shall report results of a theoretical study of the structural and energetic properties of NiN clusters for all $`N`$ up to 150. Our calculations have made no assumptions about the structure of the clusters whereby they differ from most other studies of these systems (see, e.g., \[doye, ; luo, ; mc, ; boy, ; wetzel, \] and references therein) of which the previously most detailed and unbiased studydoye  considered only up to $`N=80`$. Since we have had to perform often several 1000 structure-relaxations for a given $`N`$ without making any assumption on structure or symmetry, parameter-free methods cannot be applied, and, therefore, we have employed the embedded-atom method (EAM) in the calculation of the total energy for a given structure. Moreover, we have developed an Aufbau/Abbau method for the structure optimization. This method is partly inspired by experimental conditions where clusters often are produced by adding atoms individually to a seed, although the experiments often take place at conditions far from equilibrium and from $`T=0`$ K. A preliminary description of our approach has been presented elsewhere,pccp ; psik  but compared with that work the present one is much more accurate (through a much more detailed structure optimization that also includes more isomers), and, in addition, we shall here present various theoretical tools with which we can extract the fundamental structural and energetic properties of these clusters.

In our Aufbau/Abbau method we consider clusters with $`L`$ and $`L+P`$ atoms with $`P510`$. For each of those we study a set of randomly generated structures. Using a quasi-Newton method relaxed local-total-energy minima are identified and the structures of the lowest total energy selected. One by one atoms are added to the structure with $`L`$ atoms (many hundred times for each size), and the structures are relaxed. Parallel, one by one atoms are removed from the structure with $`L+P`$ atoms. From this two series of structures for $`LNL+P`$ those structures of the lowest energies are chosen and these are used as seeds for a new set of calculations. First, when no lower total energies are found, it is assumed that the structures of the global-total-energy minima have been identified, and we proceed to larger clusters. In some few cases, extra high-symmetry structures were added as possible structures, which in one case ($`N=75`$) resulted in an otherwise unidentified lowest-energy structure. According to our experience, for a given $`N`$ we have to consider about 500 randomly generated structures for the smallest values of $`N`$ and about 1000 for the largest values. Moreover, about 3 Aufbau/Abbau cycles are needed for the smaller $`N`$ and up to 10 for the larger $`N`$.

In comparison with our earlier study, the method has been extended so that by keeping track of more structures of the lowest total energies, we obtain the energetically lowest isomers for each size, and a much more detailed search has led to the identification of new structurals for several $`N`$. Since our approach (as any other unbiased method for structure optimization of larger, low-symmetric systems) requires a large number of total-energy calculations for different, larger structures, it cannot be combined with accurate parameter-free electronic-structure methods, but has to be based on approximate schemes for calculating the total energy. Here, we combine the Aufbau/Abbau method with the EAM method for the latter purpose.

The basic idea of the EAMeam1 ; eam2 ; eam3 ; eam4  is to consider every atom as an impurity embedded in a host provided by all other atoms. The tails of the electron densities of the host atoms inside the guest atom lead to a modification in the energy of this atom, $`F_i(\rho _i^h)`$, where $`i`$ labels the atom, and $`\rho _i^h`$ is the electron density from the host atoms. Usually this is evaluated at the site of the nucleus of the $`i`$th atom. The remaining part of the binding energy of the system of interest is approximated through pair potentials, $`\frac{1}{2}_{ji}\varphi _{ij}(R_{ij})`$ with $`R_{ij}`$ being the distance between atoms $`i`$ and $`j`$. The embedding energies $`F_i`$ and the pair potentials $`\varphi _{ij}`$ have been obtained by fitting to theoretical and experimental information on different infinite systems.eam3 ; eam4  The EAM is conceptually simple, but it grasps the major parts of the interatomic interactions and gives accurate structural and energetic information on various extended and truncated systems.eam3 ; eam4 ; mishin

Before discussing our results in detail we shall address the question of the accuracy of the EAM method. As such, it is constructed first of all for describing cohesion in extended metallic systems like crystals without or with defects and surfaces. Due to the complexity of the present study, it is not possible to apply a parameter-free method, and simpler methods like the EAM have to be used, but nevertheless it may be asked whether the EAM can be used for the, after all, relatively small systems of our present study. We shall therefore compare our calculated structural properties for nickel clusters with $`N=`$ 2, 3, 4, 5, 6, 7, 8, and 13 with those of accurate electronic-structure calculations using different density-functional methods.calleja ; michelini ; desmaris ; krueger ; reuse  Table 1 list some of the key quantities in the structural characterization of the different clusters as obtained with the different methods, and in Fig. 1 our optimized structures are shown. It is seen in Table 1 that our results agree well with those of the other studies, both concerning the structure itself and concerning the interatomic distances. Most of the discrepancies are related to smaller distortions that may be due to electronic effects (e.g., Jahn-Teller distortions) that are not included in the EAM method. Thus, we conclude that our approach is accurate.

In Fig. 2 we show our calculated binding energy per atom for nickel clusters with $`N`$ ranging from 1 to 150. As is generally found, the binding energy for clusters is first of all a monotonously increasing function of $`N`$ with, however, some extra structure. In order to identify the latter more clearly, we show in the figure also the stability function $`S(N)=E_{\mathrm{tot}}(N+1)+E_{\mathrm{tot}}(N1)2E_{\mathrm{tot}}(N)`$ which has maxima for particularly stable structures. These are also marked in the figure. As we shall see, some of the particularly stable structures correspond to particularly symmetric ones. Particularly stable structures are found for $`N`$ = 13 (the most pronounced one), 19, 23, 39, 55, and 147.

There exists some previous experimental and theoretical studies of the properties of NiN clusters, of which many focus on the inner structural aspects (including symmetry, see, e.g., \[stave, ; boy, ; wetzel, ; doye, ; parks1, ; parks2, \]). The most detailed study of the energetics of NiN clusters is the theoretical one of Doye and Walesdoye  who used the Sutton-Chen potential in calculating structural and energetic properties of NiN clusters with $`N`$ up to 80. They find particular stability for $`N=`$ 13, 38, 50, 55, 64, 71, 75, and 79, which is in partial agreement with our results. As discussed by Parks $`et`$ $`al.`$parks1 ; parks2  the structures given by Doye and Wales seem to deviate somewhat from experimental findings, suggesting that the Sutton-Chen potential is not sufficiently accurate for such studies. In a tight-binding molecular-dynamics study, Luo studied clusters with up to 55 atoms.luo  He found most stable structures for $`N=`$ 21, 31, 35, 38, 40, 47, 50, and 55, only partly in agreement with the results of Doye and Wales and of us. On the other hand, Montejano-Carrizales $`et`$ $`al.`$mc  used the EAM in analysing the so-called umbrella growth models for icosahedral Ni clusters. Although they thereby make certain assumptions on the structure their $`N`$ values for particularly stable structures are in general in good agreement with ours and with those of Doye and Wales. This leads support to the umbrella growth model. This conclusion is further supported below.

The special stability of the cluster sizes mentioned above is supported by comparing the energies of the two energetically lowest isomers. First for $`N=6`$ and onwards we identify two different isomers, and their total-energy difference is shown in the lowest panel in Fig. 2. Here, the most significant features correlate closely with those of the stability function.

Besides the stability a number of other issues have attracted significant interest in cluster physics. These include the structure of the individual clusters, whether growth modes can be identified, and how the clusters will dissociate. We shall now address these issues.

For a given cluster NiN we calculate the eigenvalues $`I_{\alpha \alpha }`$ of the matrix containing the moments of inertia $`_{i=1}^Ns_it_i`$ with $`s`$ and $`t`$ being the $`x`$, $`y`$, and $`z`$ in a coordinate system with the origin at the center of mass of the cluster. For a spherical jellium $`I_{\alpha \alpha }N^{5/3}`$, and therefore we show (Fig. 3) $`I_{\alpha \alpha }/N^{5/3}`$. If the three eigenvalues are identical the overall shape of the cluster is approximately spherical. Such clusters are marked with points in the lowest row in the uppermost part of the figure and are found for $`N`$ = 1, 4, 6, 13, 26, 28, 38, and 55, which partially correspond to particularly stable clusters; cf. Fig. 2. When two of the eigenvalues are small and only one is large, the cluster is essentially prolate (cigar-shaped), whereas it with two large eigenvalues and one small one is oblate (lens-shaped). In the figure we have marked which of the clusters have these two types of overall shapes. Except for the absolutely smallest clusters, we have that, if the cluster structures can be considered as obtained by adding atom by atom to a largely fixed core, the clusters will either be close to spherical (so that several changes between prolate and oblate structures are possible), or one or the other of the two structure types will dominate in large $`N`$-intervals. However, the results of Fig. 3 show that except for $`N`$ above about 137, the eigenvalues $`I_{\alpha \alpha }`$ are very different, so that the clusters have pronounced prolate or oblate shapes. Moreover, several changes between these two structure types occur in the range of $`N`$ we have considered.

A more detailed way of analysing whether the cluster with $`N`$ atoms can be considered as formed by adding a single atom to the cluster with $`N1`$ atoms is provided by defining similarity indices as follows. We calculate the $`(N1)(N2)/2`$ interatomic distances $`d_{ij}^0`$ (in a.u.) for the $`(N1)`$-atom cluster and sort them. Subsequently, we consider $`all`$ the possible $`(N1)`$-atomic parts of the $`N`$-atom cluster, and for those we also calculate and sort the interatomic distances $`d_{ij}`$. The smallest value of $`q_1`$ with $`\frac{(N1)(N2)}{2}q_1^2=_{i>j=1}^{N1}(d_{ij}^0d_{ij})^2`$ defines then the similarity index, $`s_1=(1+q_1)^1`$, which is 1 for a cluster that is obtained by adding one atom to the previous one, and approaches 0 for very different structures. Another similarity index $`s_2`$ is obtained by instead of using the interatomic distances we use the distances of the individual atoms to the center of mass of the cluster. Denoting these so-called radial distances $`r_i^0`$ (in a.u.) for the $`(N1)`$-atom cluster and $`r_i`$ (in a.u.) for the $`(N1)`$-atom fragment of the $`N`$-atom cluster, we consider the smallest value of $`q_2`$ with $`(N1)q_2^2=_{i=1}^{N1}(r_i^0r_i)^2`$, and define $`s_2=(1+q_2)^1`$.

These two functions are shown in Figs. 4(a) and (b). Up to $`N`$ around 40 both functions possess a number of peaks indicating that no simple growth mechanism can be identified. This is in accord with the experimental observations of Pellarin $`et`$ $`al.`$pellarin  who moreover suggest that above that the umbrella growth model agrees with their observations, in agreement with the results of Montejano-Carrizales $`et`$ $`al.`$mc  However, both of our similarity indices find for larger $`N`$ peaks at $`N=`$ 39, 51, 63, 73, 80, 87, 105, and 138, implying that also here size-dependent structural differences occur. Thus, in total we find an overall support for the umbrella growth model although for many values of $`N`$ this model does not give an exact description of the structure.

In order to obtain a quantitative estimate of the similarity of the structure for a given $`N`$ with that of a fragment of the crystal we proceed as follows. We consider a large, but finite fragment of the fcc crystal structure and place a center (arbitrarily) at (0,0,0), ($`\frac{1}{2}a`$,0,0), and ($`\frac{1}{4}a`$,$`\frac{1}{4}a`$,0) (with $`a`$ being the length of the unit cell). The set of distances from the atoms to any of these centers, $`\{r_i^{\mathrm{fcc}}\}`$, is sorted and compared with the radial distances for a given $`N`$-atom cluster, $`Nq_{\mathrm{fcc}}^2=_{i=1}^N(r_i^{\mathrm{fcc}}r_i)^2`$. This defines the fcc similarity index, $`s_{\mathrm{fcc}}=(1+q_{\mathrm{fcc}})^1`$. A related icosahedral similarity index $`s_{\mathrm{ico}}`$ is defined by using the radial distances for the (relaxed) $`N=309`$ icosahedral cluster in comparing with the different clusters. The two functions are shown in Figs. 4(c) and (d). It is seen that in the range 73–79 the structure is related both to that of a fragment of the fcc crystal with, in addition, the center at (0,0,0) and to the icosahedral structure. This indicates that in exactly this size range more different structural motifs are present. Two further peaks can be identified in Fig. 4(c): at $`N=28`$ which corresponds to the center at ($`\frac{1}{2}a`$,0,0) and at $`N=38`$ which corresponds to the center at ($`\frac{1}{4}a`$,$`\frac{1}{4}a`$,0). The lower panel shows that the icosahedral clusters are important around $`N=13`$, between $`N=51`$ and $`N=55`$, and above $`N=138`$. We finally notice that the structural characteristics often are not isolated to singular values of $`N`$ but can be identified also around the peaks.

In order to get insight into possible fragmentation channels we analyse the calculated energies $`E_{\mathrm{tot}}`$ as a function of $`N`$. Neglecting kinetic effects (this is a crude approximation and may have severe consequences for our conclusions), we seek for a given $`N`$ that value of $`K`$ for which the dissociation energy $`E_{\mathrm{tot}}(N)[E_{\mathrm{tot}}(K)+E_{\mathrm{tot}}(NK)]`$ is lowest. It turned out that, without exception, this value of $`K`$ was always 0, indicating that any cluster is stabler than its dissociation products. This result is a consequence of the fact that $`E_{\mathrm{tot}}(N)`$ is upward curved, as indirectly can be seen in the upper part of Fig. 2.

Neglecting the $`K=0`$ solution we found the results shown in Fig. 5. Except for one case, $`K`$ was found to be either 1 or 2, indicating that the clusters upon fragmentation will dissociate in a highly asymmetric way by splitting off a single or a pair of atoms. The exception is that of $`N=16`$ for which $`K=3`$ was found as a consequence of the particular stability of the $`N=13`$ cluster. The dissociation energies, also shown in Fig. 5, are seen to lie around 4 eV for $`K`$ = 1 and 2 and around 5 eV for $`K=3`$ for the larger clusters.

Concluding we have presented results of a theoretical study of structural and energetic properties of NiN clusters with all $`N`$ up to $`N=150`$. We stress that this study is the first where a systematic study of energetics, growth, dissociation, and structure of not only the absolutely smallest clusters has been performed without more or less severe assumptions of the structure. Our structure optimization was performed using our own Aufbau/Abbau method, and due to the large number of total-energy calculations that was necessary we used the approximate EAM for this purpose. As the main result we obtained the structure and the total energy as functions of $`N`$, and these results were then subsequently used in identifying particularly stable clusters, in classifying the overall shape of the clusters, in discussing whether a simple growth mechanism could be applied (this was found not to be the case), and in discussing possible fragmentation channels. In order to extract useful information from the obtained structural and energetic properties we introduced a number of quantitative measures and subsequently applied them on our results. For the first time it was shown that for $`73N79`$ the clusters possess more different types of structural characteristics, and, moreover, our theoretical study is the first such to show that for $`N=147`$ a magic number occurs for a multilayer icosahedron structure. This result is in perfect agreement with mass-spectroscopical experiments of Pellarin $`et`$ $`al.`$pellarin

The authors are grateful to Fonds der Chemischen Industrie for very generous support. This work was supported by the SFB 277 at the University of Saarland.


format: text

text:
# Twelve Ways to Build CMS Crossings from ROOT Files

## I Introduction

The CMS experiment CMS  is one of the two multipurpose experiments being under construction to operate at the future Large Hadron Collider (LHC) at CERN. A particularly important aspect of the CMS core software is the database system that will be used to handle the petabytes of data that the experiment will produce. The experiment has recently decided to move away from its Objectivity based system in favor of an hybrid solution based on ROOT I/O ROOT . This paper describes the work done in this context so to evaluate the performances of using the very specific ROOT classes (especially TTree and TClonesArray) for the data storage. In the meantime, CMS also started a more direct replacement of the existing Objectivity implementation Bill .

We did not choose to explore the many use-cases of CMS, but rather focused on a single representative one and studied the many ways to implement it with ROOT I/O.

The selected use-case is the last step of the simulation chain Production : starting from the events produced by the detector simulation, we simulate the raw data produced by the detector. Due to the high luminosity of the LHC machine, this involves firstly the superimposition to the signal event of a number of pileup events. The resulting crossing is then digitized, that is the effect of the front end electronics is simulated so to produced a digitised crossing or raw data. From the applications currently developped by CMS, this step is the most critical from the I/O point of view: for the simulation of each raw data event of $`2MBytes`$ size, one must load in memory about one hundred and fifty minimum bias events of $`300KBytes`$ size. This requires huge memory and is an intensive data reading process. Moreover, this is an unusual use of ROOT I/O on several aspects:

* The CMS code contains C++ templates, STL containers and it uses external packages whose classes cannot be instrumented for ROOT. The ROOT support for templates, standard containers and external classes is quite new and still not fully mature.
* The pileup events, stored in ROOT files, are not read sequentially. Ideally, such events should be taken from a large enough statistic, so to produce uncorrelated sequences of pileup to be added on each signal event. Since we have only a finite statistic of pileup events, we chose to select them in a random fashion so to limit the possible effect of correlation between different sequences. Nevertheless, the selection method is not fully random (real randomness would imply a change of file for each new pileup event and would have a strong impact on performances).
* The events are not processed one by one in memory, such as in standard analysis jobs as described in the ROOT documentation. As explained before, we rather have at the same time in memory all the hits from the signal and from the many minimum bias events that we want to pile up.

In this context, we have compared the performance of four kinds of containers (section II.2), combined with three different ways (section II.3) to store them in ROOT files. This leads to the twelve strategies advertized in the title. Also worth to be mentionned, we really focused on the huge event data to be transfered: we did not store neither meta-data nor the links or pointers between the elements of the events. Surely, the latter should be integrated in a later version of the testbed.

The section II gives a more detailed description of the testbed, what we have tested and on which platform. Then, the section III presents the most interesting results.

## II Testbed

### II.1 Main Use-Case

The testbed is basically able to build crossings, as described in figure 1. As input, we have a file of five hundred signal events, and one hundred files of five hundred minimum bias events (actually, for the testbed, one hundred times the same original file). A typical job is the production of five hundred crossings. Building a crossing and digitizing it consists of the following steps:

1. A persistency manager loads in memory the next signal event from the signal file
2. Another persistent manager loads 153 minimum bias events from the minimum bias files (this number corresponds to the luminosity of $`10^{34}cm^2s^1`$). These events are selected pseudo-randomly: we take X consecutive events, then we do a random jump between 0 and Y, etc. X is called ”burst” and Y is called ”jump”. Their typical values in CMS applications are 3 and 10. Also, when loading an event, we select randomly which rank it will have in the collection in memory.
3. The digitizer collects all the hits from the signal and minimum bias events in memory, then simulates the detector front end electronics and produces the correspondant digis.
4. A third persistent manager takes the digis and stores them as the next digitized crossing in the digis output file.

The expected bottleneck and area of interest is the loading of the minimum bias events. Below follows the mean number of objects and their size for such events:

* 351 instances of RtbGenParticle (whose raw size is 46 bytes).
* 584 instances of RtbSimVertex (whose raw size is 34 bytes).
* 169 instances of RtbSimTrack (whose raw size is 38 bytes).
* 3282 instances of RtbCaloHit (whose raw size is 20 bytes).
* 1871 instances of RtbTrackHit (whose raw size ist 56 bytes).

This leads to an event size of 208 Kbytes, if we just consider the pure data, without any adjunction for the support of persistency. As one can see, the size is mainly dominated by the detector simulated hits. We also computed that the mean event size would be 392 Kbytes if each numerical attribute would be of type double (this information will prove useful in section II.3.2).

In the results section, we will especially study the size of the minimum bias files, and the time necessary to load 153 pseudo-random events from them.

### II.2 Crossing Data Model

We chose to put all our persistent data in the shared hierarchy of folders proposed by ROOT (tree of instances of TFolder). This permits to decouple completely the persistency mechanism from the digitizing code.

The folder called //root/crossing/digis is the output ”event” of the current crossing: it contains a container of instances of RtbCaloDigis and a container of instances of RtbTrackDigis.

The 153 folders which we have called //root/crossing/minbias\*, plus the single one called //root/crossing/signal represent the input events composing the current crossing. Each of these folders contains a container for each kind of input event objects: RtbTrackHit, RtbCaloHit, RtbSimTrack, RtbSimVertex and RtbGenParticle.

Each time we run the testbed, we chose a given kind of container which is used for all the input and output data. All those containers inherit from a common abstract class:

```
  template <class T>
  class RtbVTArray<T>
   {
    public:
    // write interface
    virtual void clear() =0 ;
    virtual void add( const T & ) =0 ;
    // read interface
    virtual UInt_t size() const =0 ;
    virtual const T & operator[]
     ( UInt_t i ) const =0 ;
   } ;

```

This (simplified) class header shows that we rely on size() and operator\[\] to read the objects of a collection. On the other side, when creating a collection, we do not pre-allocate the size of the container to be filled: in the CMS code we have imported in the testbed, the size of the created collections is rarely known in advance. This write interface is not the most efficient, but it fits the user needs and what matters most for this testbed is the read efficiency. Four kinds of concrete containers have been implemented, which are described below.

#### II.2.1 Standard STL vector

The class RtbStlArray\<T\> wraps a standard std::vector\<T\>. The support of this kind of container is quite recent within ROOT. The main benefit of this container is that any kind of T can be collected (it is not needed that T inherits from TObject). The associated disadvantage is that ROOT does not apply any of its optimizations when storing/retrieving the objects, expecially its attribute-wise serialization.

#### II.2.2 Dynamic C Array

The class RtbCArray\<T\> contains a simple C array dynamicaly allocated. Each time the array is full, the size is multiplied by two and the array reallocated.

The objects which are collected in the C array must be instrumented with ClassDef, otherwise ROOT I/O will not be able to save them. Since we do not want to impose that T is instrumented for ROOT, we wrote a template class RtbClassDef\<T\>, which inherits from T and is instrumented with ClassDef. Each time we want to add an instance of T into the C array, we first change its type from T to RtbClassDef\<T\>. So to evaluate the eventual cost of this process, we also kept the possibility to compile the testbed with event data classes directly instrumented with ClassDef (this is done by unsetting the macro RTB\_FOREIGN).

The class RtbCArray\<T\> should not be seen as a container to be used in a real application. It is rather a toy container, just written for comparison and to evaluate the performance of a good old simplistic C array.

#### II.2.3 TObjArray

The class RtbObjArray\<T\> wraps a TObjArray.

As for any ROOT collection, we cannot add an object to the TObjArray when its class does not inherit from TObject. Again, we do not want to impose that T is instrumented for ROOT, so we wrote a template class RtbObj\<T\> which inherits both from TObject and T. Each time we want to add an object into the collection, we first transform it into an instance of RtbObj\<T\>. So to evaluate the eventual cost of this process, we also kept the possibility to compile the testbed with top event data classes directly inheriting from TObject (this is done by setting the macro RTB\_TOBJECTS).

As for RtbCArray\<T\>, RtbObjArray\<T\> should not be seen as a container usable in a real application. It is here for comparison, and should have the worst performances of all the containers, because it is handling the objects by pointers (and not by value).

#### II.2.4 TClonesArray

The class RtbClonesArray\<T\> wraps a TClonesArray. As for RtbObjArray\<T\>, the instances of T must be transformed into instances of RtbObj\<T\> before they are added to the collection.

The class RtbClonesArray\<T\> is the real alternative to RtbStlArray\<T\>. Since it is based on TClonesArray, which is the official optimized container of ROOT, it should show the best perfomance. The counterpart is that it is harder to use and the collected objects must be kind of TObject.

### II.3 Persistency Managers

The task of a persistency manager is to transfer an event from memory (a TFolder and its contents) to disk (an entry in a TFile) and vice-versa. Three flavors have been implemented, which are described below.

#### II.3.1 RtbPomKeys

This persistency manager implements the simplest approach, which only uses the base ROOT I/O level and the class TKeys: one directly writes the TFolder to the TFile, each time with a different meaningful name: the name of the original folder plus an incrementing rank.

For example, in our main use-case, the output of the first digitized crossing will be saved as digis0, the second one as digis1, etc.

#### II.3.2 RtbPomTreeMatrix

The central idea of the second manager is to avoid the use of ROOT instrumentation and dictionaries: this is achieved by transfering the data of each container into an instance of TMatrixD and storing the matrix instead of the container. More precisely, each row of the matrix is the copy of one object from the original container, and each column corresponds to a given attribute of the class of the objects.

In this manager, we also chose to use a TTree. Each entry of the tree is an event. Each branch is dedicated to one kind of event objects, and attached to the corresponding matrix.

In this approach, we do not take profit of the ROOT persistency features, such as the generated streamers, and we must write by ourselves the code which transfers the data between the containers and the matrices (for each persistent class). Also, since any number is transformed into a Double\_t, we expect files to be twice bigger than their normal size.

On the other hand, we do not suffer from the ROOT parser limitations and bugs. We are quite confident here in the data retrieved from the files, and this manager is primarily used in the testbed so to check that the others managers are also correctly retrieving the data.

#### II.3.3 RtbPomTreeDirect

This third kind of manager is the ROOT recommended approach. Each file contains a single TTree whose entries are the events. Similarly to the previous manager, there is one top level branch for each kind of event object, but this branch is directly attached to the corresponding container in memory.

The level of split is a parameter of the testbed, but we generally use the recommended default of 99. Concerning the size of buffers, we rather tend to reduce them to 8000 KBytes (empirical best value).

#### II.3.4 Main Use-Case

Once a persistent manager is built, and before writing or reading an event from a file, one must connect the manager to a given folder in memory and to a given file.

As one can see in the main use-case (see figure 1), there is a manager for the signal events, and another one for the digis, because they are always connected to the same files and folders.

With regards to the minimum bias events, we were not able to build a manager for each of the 153 events in memory (each manager has some internal buffers and this would require a large amount of memory space). Thus, we use a single manager, which must be reconnected to a new memory folder and eventually to a new file after reading each event. This connexion time is also something we have closely looked at during the analysis of the results.

### II.4 Implementation issues

Provided one uses the option -p of rootcint, ROOT has very greatly improved its support of foreign classes, templates and std containers. It is now also possible to enforce the respect of ANSI C++ when compiling. However, there are still some issues with the use of ROOT that we discuss below.

#### II.4.1 Documentation

We really lack a central place where would be documented which subset of C++ is supported in the interpreter, which subset can be made persistent, which one can be used within a TTree, and which one can be used with a TClonesArray within a TTree. Since the ROOT team consider as a bug whatever is not supported, they try to fix any such case rather than report it. As a result, each user must rediscover by himself the unsupported cases, when they do not do invisible damage.

#### II.4.2 LinkDef

It has proved painful to write the configuration files for the generation of dictionaries. One must explicit all the classes which must be parsed, and in the right order. Even with only seven top classes to be made persistent, we felt the need to write a perl script for the generation of the LinkDef file. A key cause is that when one parses a given class, one must have parsed before all the classes of all the attributes, including each instanciated template. We wonder if one could not find a way to automatize this within rootcint.

#### II.4.3 Tuning of TChain branches

Since we handle a very large number of minimum bias files with the same internal tree, it was rather logical to use an instance of TChain. Actually, the fact that TChain inherits from TTree is misleading. In particular, if you get the branches and customize them, all your changes will disappear when the chain move internally to a new file. Thus, you must detect yourself any change of file and do again the branch customizations. In our use-case of random events and detailed branch tuning, TChain has finally not proved helpful.

#### II.4.4 TBranch attachment

When one attaches a branch to a given variable, it does not give the address of a variable, but the adress of a pointer to a variable. This will let any C++ programmer think that he can change later the pointer value, so to fill another variable. This is not true! We cannot imagine any technical reason for this, but if such an obstacle exists, the signature of the attachment method should be changed for the address of a variable.

#### II.4.5 TClonesArray

This class has really turned out to be hard to understand , with a least six size-like methods. We can understand that this comes from backward compatibility constraints, but still it is a problem as this class is the central piece of the persistency service. Our proposal to have a new ROOT collection class for the persistency (without gaps !) has not been supported by the ROOT team. Actually, much work is currently done for the efficient support of std::vector\<T\>, and we guess that this class could become what we would like to see.

### II.5 Parameters of the Testbed

A few options can be set before compiling, thanks to macros in the central header file. They have been kept as compilation option because it was hard to make them runtime options, and not necessarily useful. Here they are:

* RTB\_FOREIGN: if set (the default), the persistent classes are not instrumented with the macro ClassDef, and they will be considered as foreign classes by ROOT.
* RTB\_TOBJECTS: if set (not the default), the top persistent classes inherit from TObject and are instrumented with ClassDef (whatever the value of RTB\_FOREIGN). If not set, it implies the use of the RtbObj\<T\> when appending the objects to ROOT containers and the use of RtbClassDef\<T\> when putting them into a dynamic C array.
* RTB\_RESET: if set (not the default), the empty constructors of the persistent classes set all their attributes to 0 (we were expecting an eventual impact on the compression performance).

At runtime, one must choose within four kinds of containers and three kinds of persistent managers. This leads to twelves base strategies. All the testbed results will be displayed as an array of twelve cells corresponding to these strategies (see Table 1).

In the table, the third column is just there to see how bad is TObjArray, the second column is just there to see the beavior of a good old dynamic C array, and the second line is mainly a way to counter-check the validity of the retrieved data. So, what matters most are the corners of the array, especially the top left corner (Keys/Stl) which is the strategy used by CMS for the direct replacement of its Objectivity implementation, and the bottom right corner (Tree/Clones) which is the solution advertized by ROOT team and evaluated in this work.

On top of the twelve strategies mentionned above, we have additional parameters, whose value can affect the performances of these strategies differently:

* Compression level: it should slow down the writing of objects, also slightly the reading, and reduce the size of files. A value of 1 is expected to be the best compromise.
* Split level: TTree and ROOT containers use some sort of attribute-wise storage mechanism. The split level is the depth of the decomposition.
* Size of tree buffers: amount of data which is read in one bunch from the file, for a given branch.
* Randomness: as described in the main use-case (see section II.1), the burst and jump can be changed so to be close to a sequential use of minimum bias events, or on the contrary close to a really random access pattern.
* Size of containers: this parameter permits to reduce the size of events by a given factor, so to measure the effect of this size on the performance. The default is 1 (takes all the data). If the value is 10, when the input files are prepared, only 1 from 10 elements is taken.
* Number of crossings: the use-case specifies that each job must build 500 crossings. This number can be lowered to see if the mean performance is the same or to shorten the execution.

### II.6 Platform used for the tests

All the results which are given below have been obtained with a PC where the testbed was the only application running. Its characteristics and software environment are the following:

* Processor: Pentium 4, 1.8 GHz.
* Memory: about 512 Mbytes.
* Disk: IDE.
* System: RedHat Linux 7.3 .
* Compiler: gcc 3.2 .
* Root release: 3.05/03 .

## III Results

We can hardly give here all the results we have obtained when tuning the parameters of the testbed. We will rather start with the best performance we have obtained, then show and comment the effect of changing the value of some relevant parameters.

### III.1 Best Results

In Table 2, you will find the size of the minimum bias files (divided by the numbers of events, i.e. 500) and the mean time to read all the events of a crossing (153 minimum bias event plus a signal event). We made the 500 crossings and used all the input data. We also used the default compression level (1) which appeared to be always the best compromise.

As one can expect, TObjArray is always the worst choice for the read performance.

Also expected, all the strategies with matrices give files of the same size. The read time differs from one container to the other, because of the final read step where the data is taken from the matrices and transformed into new objects added to the containers. It can be seen as a measurement of the efficiency of the container add() method.

Let’s now compare std::vector versus TClonesArray. Within a TTree, TClonesArray is by far the fastest, and the files are incredibly smaller. The key reason here is the split mechanism: each attribute of each persistent class is given its own branch and buffers (some sort of ”attribute-wise” storage). This, combined with the compression of data, proves very efficient. We do not have such performance with std::vector\<\> because ROOT does not support yet the splitting in such a case (yet one can notice a small improvement when the instances of std::vector\<\> are stored in a TTree rather than directly in the file). It seems that a maximum split is always worth, even when one reads back all the branches from the tree. With regards to the size of buffers, it appeared very complex to predict the best value, depending on the split level, the type of object attributes and the value of burst: it is useless to read much in advance, if there is a random event jump coming. We proceeded empirically and finished with a size of 8000 KBytes, largely smaller than the ROOT default.

The storage of TClonesArray directly in TFile (top right cell, Keys/Clones) exhibits a rather poor performance. The reason is surely because we had to switch off the ByPassStreamer option, apparently buggy in such a context. As a result, std:vector\<\> is the quickest alternative when not using a tree.

### III.2 Remove compression

When switching off the default compression of the data (table 3), one can measure how much it is efficient for the size of file! Without compression, the read time is globally smaller, without changing the classification of the corners.

Surprisingly enough, the strategy Matrix/Stl is really fast, and one can wonder what it would be if its implementation was improved. Also a surprise, but a bad one, only the Tree/Clones and the Matrix/\* strategies have file sizes which match the predictions (see section II.1). Other strategies have files largely too big, and we did not fully investigate why.

### III.3 Then increase randomness and reduce data volume

Our next step has been to reduce the parameter burst to 1 and to increase jump to 1000, i.e. to take the minimum bias events almost fully randomly.

In table 4, one can see that increased randomness reduces more the performance of the \*/Clones strategies than the \*/Stl, and reduces more the performance of Tree/\* than Keys/\*. This makes Keys/Stl clearly the best choice.

On top of that, if we reduce by a factor of ten the number of elements we have in the events (table 5), the effect on the specialized ROOT classes is even worse: the strategy Tree/Clones becomes the worst strategy (apart from Tree/C and \*/ObjArray) !

The lesson is quite clear: in our specific CMS use-case, the use of TTree and TClonesArray is by far the most efficient strategy, but this cannot be generalized. It highly depends on the volume of data and the amount of randomness.

### III.4 Other Results

Resetting the attributes to 0 in the empty constructors of the event data does not appear to help compression. So after trying it, we went back to an implementation where the empty contructors let undefined values in the attributes.

Unsetting RTB\_FOREIGN or setting RTB\_TOBJECTS has not greatly improved the performance of the ROOT collections, so we turned back to the use of our templated wrappers RtbClassDef\<T\> and RtbObj\<T\>.

For what concerns the write cpu time of the different strategies, the twelve strategies compares almost similarly: only the Tree/\* strategies are found slightly slower.

At last, we must confess we did not systematically build the whole set of 500 crossings that is specified by the use-case. When we did, we always noticed that the overall performance was slightly better.

## IV Conclusion

We have succeeded to read pseudo-random entries from a TChain and to dispatch them to a few hundred TFolders (despite the fact thar the tuning of the TChain branches has not been straightforward). Support for foreign classes, templates and C++ standard library has greatly improved in the recent releases of ROOT.

The magic couple TTree/TClonesArray has proved very efficient for our use-case, yet it requires top level TObjects and the benefits can become losses with smaller data volume or random access pattern. One can simply use STL vectors and store them directly into root files. Their integration in a TTree is not yet as good as a TClonesArray, but this could change in a future release of ROOT.

If this testbed were to be improved, one major step would be to add real pointers or TRefs between the objects (instead of the current indexes) and measure the impact on performance.

You can obtain the testbed source code (for linux with gmake and gcc) by contacting the authors.

###### Acknowledgements.

We would like to thank the ROOT team who has always quickly answered to any of our question or bug report, and for the discussion we had concerning our results. We are also grateful to Pascal Paganini for his contribution to the implementation of the CMS use-case in the testbed.


text:
# Understanding of the Vortex Lattice Behavior Induced by Selective Removal of Atoms from a Bose-Einstein Condensate

## 1 Introduction

It is well known that quantized vortices play an important role in understanding many phenomena of superfluid, superconductivity and even astrophysics. In the latest several years a remarkable experimental development has been made in the rapidly rotating Bose-Einstein condensate(BEC). A large amount of angular momentum is imparted into the system and its rotation frequency becomes faster and faster. The ground state of the fast rotating BEC was observed to be a regular triangular vortex lattice at nearly zero temperature. This provides a new stage of studying both the equilibrium and the dynamic properties of the vortex matter.

Recently two novel experiments were reported by Cornell’s group in JILA. Atoms in the hyperfine spin $`|F=1`$ state, were rotated rapidly by the evaporative spin-up mechanism. About $`180`$ vortices were observed in the cloud and they were arranged into a regular lattice structure. And then, an atom removal laser, whose frequency was tuned to the transition energy between the $`|F=1`$ state and the $`|F=0`$ state, was applied to the center of the condensate, and the recoil from a spontaneously scatted photon blasts atoms out of the condensate. In this way, the atoms with lower angular momentum were selectively removed from the system.

Many exotic features were observed when this technique is applied. After the atom removal laser was applied continuously, the vortex number contained in the cloud increased form $`180`$ to about $`230`$, and the vortices began to aggregate to the center of the condensate, eventually they formed a giant vortex. If the laser was applied for a short limited time, it was reported in Ref. that a periodic formation and disappearance of the giant vortex core was observed after applying a strong atom removal pulse, and in Ref. that the Tkachenko oscillation of the vortex lattice would be found after a weak pulse. Both kinds oscillation were found damped.

A number of questions are raised from these experiments, such as how to understand the vortex number increase and the vortices aggregating to the center of condensate, why the strong and weak laser strength will cause two different kinds of oscillations, how to understand the relation between the oscillation frequency and the initial conditions of the condensate, and which mechanism leads to the damping. We notice that recently there are lots of theoretical efforts around these experiments, some of them try to understand these observations numerically, and some of them focus on the Tkachenko wave of the vortex lattice state. However, in this paper, we will present an explanation to these questions from another point of view, namely using a dual $`2+1`$-dimensional Maxwell electrodynamic (MED) description, and we will show that the main features observed can be understood naturally in this framework, although qualitatively in some aspects.

This dual electromagnetic description of the interacting bosons model was obtained by the duality transformation early in the study of Helium superfluidity, and then in the past decades it has been discussed in the contents of superconductivity and superfluid film by various authors. Compared to the model discussed before, we would like to emphasize some important differences of the system studied here. First, a large amount of angular momentum have been imparted into the atoms and the rotational invariance is restored before the formation of condensate, therefore the total angular momentum of the condensate should be conserved to a non-zero value. The second, the atom removal laser violates the current conservation condition. And the third, the system is dilute and the interaction between atoms is weak, the sound velocity is consequently relatively small. We will show that these differences lead to these exotic phenomena in these experiments, in other words, these observed phenomena reflect some intrinsic properties of the rotating BEC.

## 2 The Dual Picture

In this section we will firstly make a brief description of the dual picture for the interacting bosons system, for the details we refer to the Ref.. Beginning with the coherent state path integral formulism, we write down the quantum mechanical propagator

$$Z=D[\psi ,\overline{\psi }]\mathrm{exp}\{i𝑑td^2\stackrel{}{r}\},$$
(1)
where the Lagrangian

$$=\mathrm{}\overline{\psi }\frac{_t}{i}\psi \frac{\mathrm{}^2}{2m}\left|\frac{\stackrel{}{}}{i}\psi (r)\right|^2\frac{g}{2}\left|\psi (r)\right|^4.$$
(2)
$`g`$ is the interaction strength and equals to $`4\pi \mathrm{}^2a_{sc}/m`$ for ultracold atomic system, where $`a_{sc}`$ is the $`s`$-wave scattering length between atoms. Here $`\psi `$ is a complex field describing the bosonic many-body system, and $`\overline{\psi }`$ is its complex conjecture. We can factor $`\psi `$ as $`\sqrt{\rho }e^{i\theta }e^{i\theta _v}`$. The field $`\rho `$ describes the density of bosons. The phase of $`\psi `$ contains two parts, $`\theta `$ is the smooth part free from singularity and $`\theta _v`$ is that caused by the presence of vortices.

First of all, we can introduce an auxiliary field $`\stackrel{}{J}`$ by performing a Hubbard-Stratanovich (H-S) transformation, the physical meaning of which is the superfluid current. By defining the three vectors $`\stackrel{~}{J}(\mathrm{}\rho ,\mathrm{}\stackrel{}{J})`$ and $`\stackrel{~}{}=(_t,\stackrel{}{})`$, in the low frequency and long wavelength limits one can neglect the second order gradient term $`|\rho |^2/\rho `$, and obtain that

$`Z={\displaystyle D[J]D[\rho ]D[\theta ]D[\theta _v]\mathrm{exp}\left\{i𝑑td^2\stackrel{}{r}\left(\stackrel{~}{J}\stackrel{~}{}\theta \stackrel{~}{J}\stackrel{~}{}\theta _v+\frac{m}{2\rho }|\stackrel{}{J}|^2\frac{g}{2}\rho ^2\right)\right\}}.`$
(3)
Integrating over $`\theta `$, one can naturally obtain a $`\delta `$-function constraint $`\stackrel{~}{}\stackrel{~}{J}=0`$, which is just the conventional current conservation condition

$$\frac{\rho }{t}+\stackrel{}{}\stackrel{}{J}=0.$$
(4)
Owing to this constraint, we can introduce a three-component gauge field $`\stackrel{~}{b}`$ defined as $`(b_0,\stackrel{}{b})`$, which satisfies

$$\stackrel{~}{J}=\stackrel{~}{}\times \stackrel{~}{b},$$
(5)
the density and the current of the bosons are thereby determined by the gauge field. Hence, $`\rho =\stackrel{}{}\times \stackrel{}{b}`$ is compared to a magnetic field perpendicular to the 2-dimensional plane. We can expand $`\rho `$ as $`\overline{\rho }+\delta \rho (\stackrel{}{r},t)`$, $`\overline{\rho }`$ is the average density which is compared to an external magnetic field, while $`\delta \rho (\stackrel{}{r},t)`$ is the density fluctuation which is compared to the electromagnetic wave. The electric field $`\stackrel{}{E}`$ is defined as $`\stackrel{}{}b_0_0\stackrel{}{b}`$, and it is perpendicular to the boson current $`\stackrel{}{J}`$. Therefore the Eq.(4) can be rewritten as

$$\frac{B}{t}+\stackrel{}{}\times \stackrel{}{E}=0,$$
(6)
which is just the Maxwell equation in a monopole free space-time.

With replacing the $`\rho `$ by $`\overline{\rho }`$ in the denominator, the last two terms in Eq.(3) can also be reexpressed in terms of the gauge field $`\stackrel{~}{b}`$, which is just the dynamic term $`F^{\mu \nu }F_{\mu \nu }`$ of the gauge field as in the conventional MED. It determines the dynamic properties of boson density fluctuation, i.e. phonon, and recovers the well-known result for the sound velocity $`s=\sqrt{g\overline{\rho }/m}`$, which plays the same role as the light velocity in $`2+1`$-dimensional MED.

The coupling between the vortices and bosons can be illustrated as

$`\stackrel{~}{J}\stackrel{~}{}\theta _v=\stackrel{~}{b}\left(\stackrel{~}{}\times \stackrel{~}{}\theta _v\right)=\stackrel{~}{b}\stackrel{~}{J}_v.`$ (7)
Here $`\stackrel{~}{J}_v`$ is defined as $`\stackrel{~}{}\times \stackrel{~}{}\theta _v`$. It can be verified that the zero-component of $`\stackrel{~}{J}_v`$, which is denoted by $`\rho _v`$, represents the vortex density, and the other two spatial components represent the vortex current. The Eq.(7) indicates that the vortices can be viewed as charged particles, which are coupled to an electromagnetic field via gauge coupling. The coupling constant, that is the charge $`q`$ of the charged particles, is unit. In the Coulomb gauge, the scalar potential $`b_0`$ is instantaneous, one can integrate it out and obtain the mass of vortex $`m_v`$, which is equal to $`\overline{\rho }\mathrm{}^2\mathrm{ln}\frac{R}{\mathrm{\Lambda }}/(4\pi s^2m)`$ for a uniform condensate, and the mutual logarithmic interaction between vortices described by

$$_{int}=\frac{\overline{\rho }\mathrm{}^2}{4\pi m}d^2\stackrel{}{r}d^2\stackrel{}{r^{}}\rho _v(\stackrel{}{r})\mathrm{ln}\frac{|\stackrel{}{r}\stackrel{}{r^{}}|}{\mathrm{\Lambda }}\rho _v(\stackrel{}{r^{}}).$$
(8)
Here $`R`$ is the radius of the condensate and $`\mathrm{\Lambda }`$ is the radius of vortex core.

Thus we have obtained a dual electrodynamic description of quantized vortices, which is briefly summarized in the following table. In the dual picture the vortex lattice state corresponds to the Wigner crystal state of 2-dimensional electrons in the external magnetic field.

| Original Picture | Dual Picture |
| --- | --- |
| vortex | charged particle |
| average | external magnetic field |
| superfluid density $`\overline{\rho }`$ |  |
| superfluid current $`J`$ | electric field $`E_\mu =ϵ_{\mu \nu }J_\nu `$ |
| density fluctuation | electromagnetic wave |
| sound velocity | light velocity |
| the coupling between vortex | gauge coupling |
| and sound wave |  |

Such a description can be directly applied to the confined BEC when it is in the Thomas-Fermi regime. The average density here is therefore determined by the trapping potential, as well as the centrifugal repulsive force caused by rotating, and we will take it as inverted parabola in this regime. Here we would like to point out that the main assumption leading to the electrodynamic description is that the vortices are treated as point-like particles, the physics inside the vortex core is not considered, as a consequence we take the mass of vortex as its electromagnetic mass throughout this paper. Besides, in the following treatment we view the system as quasi-two dimensional and neglect the density inhomogeneity when estimating the vortex mass. In the quantum Hall mean field regime, it is shown that the physics of vortex core will become important. However, before entering this regime, and when the Thomas-Fermi theory is still valid, these neglected effects will only slightly renormalize the vortex mass as well as the interaction between vortices.

## 3 Explanation of the Experiments

### 3.1 Vortex Number Increasing

In JILA experiments the rotational invariance is restored before condensation, therefore the condensate should obey the total angular momentum conservation constraint. Such a system can be described by adding a Lagrange multiplier term to the Lagrangian.

$$Z=𝑑\psi 𝑑\mathrm{\Omega }\mathrm{exp}\{id^2\stackrel{}{r}𝑑t(\psi ,\overline{\psi })\mathrm{\Omega }(\overline{\psi }\widehat{L}\psi L_0)\}.$$
(9)
It is obvious that one will obtain a $`\delta `$-function which means that the total angular momentum must always equal to $`L_0`$ after integrating $`\mathrm{\Omega }`$ out. The physical meaning of $`\mathrm{\Omega }`$ is in fact the effective rotation frequency. Following the procedure discussed in above section, it is found that only three additional terms, namely $`2mb_0\mathrm{\Omega }+\frac{1}{2}m\mathrm{\Omega }^2\rho r^2+\mathrm{\Omega }L_0`$, will be added to the Lagrangian. One can see that only the zeroth component of the gauge field responds to the rotation. Integrating $`\mathrm{\Omega }`$ out and neglecting the density fluctuation, we obtain

$$=\frac{(2mb_0+L_0)^2}{2m\rho r^2}+\text{other terms}.$$
(10)
Following the same step which leads to Eq.(8), integrating out $`b_0`$, we can obtain a local chemical potential term in the following Lagrangian

$`_{int}={\displaystyle d^2\stackrel{}{r}d^2\stackrel{}{r^{}}\rho _v(\stackrel{}{r})G(|\stackrel{}{r}\stackrel{}{r^{}}|)\rho _v(\stackrel{}{r^{}})}+{\displaystyle d^2\stackrel{}{r}\frac{2L_0}{\rho r^2}\left(d^2r^{}G(|\stackrel{}{r}\stackrel{}{r^{}}|)\right)\rho _v(r)}`$
$`=_{int}+{\displaystyle d^2\stackrel{}{r}\mu (\stackrel{}{r})\rho _v(\stackrel{}{r})}.`$ (11)
Here $`G(|\stackrel{}{r}\stackrel{}{r}^{}|)`$ is the Green’s function for the $`b_0`$ field acquiring a very small mass $`2m/\rho r^2`$, which will vanish in the large atom number limit. The first term in the Eq. (11) indicates that the vortex system can be viewed as an interacting Coulomb gas system, and the second term defines a local chemical potential $`\mu (r)`$ for vortices emerging from the angular momentum constraint, which is proportional to $`L_0/\rho r^2`$. The atom removal laser, which removes atoms with lower angular momentum, decreases the atom density $`\rho `$ with little change to the total angular momentum $`L_0`$, leads to the increase of the vortex chemical potential. Notice that the interaction energy $`_{int}`$ depends quadratically dependence on vortex density $`\rho _v`$, and the interaction between vortices are repulsive, so

$$\frac{\mu }{\rho _v}=\frac{\delta ^2_{int}}{\delta \rho _v^2}>0.$$
(12)
This shows the increase of vortex density.

### 3.2 Aggregate

For simplicity we assume that the atom removal laser acts only at the center point. Extension to the Gaussian shape laser is straightforward and the qualitative results will not be changed. When the laser is applied, a centripetal atom current will be induced and the current conservation condition will be violated at the center. In the dual picture, a circular electric field will be generated and the Maxwell equation should be modified in the following way

$$\frac{B}{t}+\stackrel{}{}\times \stackrel{}{E}=I\delta (\stackrel{}{r}).$$
(13)
Here the value of $`I`$ equals to the number of atoms removed from the condensate per unit time. Therefore it is modelled as a monopole current running through the plane. The magnitude of the electric field is approximately $`I/(2\pi |\stackrel{}{r}|)`$ nearby the center. For the case of continuous laser the vortices as charges will be accelerated azimuthally by the electric field, and the Lorentz force will drive the vortices to aggregate toward the center, and eventually to form a giant vortex.

### 3.3 Two Kinds Oscillation

After a pulse laser which only acts for a short time $`\mathrm{\Delta }T`$, the vortices will acquire a velocity $`\stackrel{}{v}`$ and leave the balance positions of the lattice. For the vortices closest to the center, the magnitudes of their velocities can be approximated by

$$v=Eq\mathrm{\Delta }T/m_vI\mathrm{\Delta }T\mathrm{}/(2\pi m_vd),$$
(14)
$`d`$ is the vortex lattice spacing. Therefore, there are two forces which can cause the oscillation motion, one is the Lorentz force $`f_1`$, and the other is the mutual interaction force between vortices $`f_2`$. It is easy to obtain the ratio of the two forces:

$$\frac{f_1}{f_2}=I\mathrm{\Delta }T\frac{m}{m_v}.$$
(15)
For these experiments we can naturally assume that the number of atoms removed from the condensate $`I`$ is an increase function of the power of removal laser. Thus Eq.(15) indicates that the Lorentz force is dominant for a strong laser pulse, such as the pulse used in Ref.. In this limit, the vortices will execute circumnutation. When the vortices move very close to each other, they will merge into a giant vortex, and separate into several single vortices again when they oscillate away from each other. The mutual interaction force will be dominant for the weak pulse case. In this limit, it will cause a collective oscillation of vortex lattice known as Tkachenko oscillation.

The oscillation frequency $`\omega `$ depends on the initial conditions of the condensate, such as the number of atoms and the rotation frequency $`\mathrm{\Omega }`$. In the strong pulse limit, for the motion is mostly determined by the Lorentz force, $`\omega =\mathrm{}\rho /m_v`$. Because the density $`\rho `$ is proportional to the number of atoms $`N`$ and is a decreasing function of the rotation frequency $`\mathrm{\Omega }`$, the oscillation frequency is linearly dependent on $`N`$ and increases with $`\mathrm{\Omega }`$ decreasing.

In the weak pulse limit, the oscillation frequency $`\omega `$ is proportional to $`\sqrt{\rho }/d`$. It is not obvious to tell whether $`\omega `$ will decrease or not when the condensate rotates faster, because both $`\rho `$ and $`d`$ will decrease with the rotation frequency $`\mathrm{\Omega }`$ increasing. We make a rough estimate as following: For a fast rotating BEC confined in a harmonic trap with trapping frequency $`\omega _{}`$, the effective confining potential is $`(\omega _{}^2\mathrm{\Omega }^2)r^2/2`$. Using the Thomas-Fermi approximation we can obtain that the condensate density nearby the center $`\overline{\rho }`$ is proportional to $`\sqrt{1\frac{\mathrm{\Omega }^2}{\omega _{}^2}}`$, and then using the formula $`d^21/\overline{\rho }_v=\mathrm{}/(2m\mathrm{\Omega })`$, we can obtain that

$$\omega =C\sqrt{\frac{\mathrm{\Omega }}{\omega _{}}\sqrt{1\frac{\mathrm{\Omega }^2}{\omega _{}^2}}}.$$
(16)
Here $`C`$ is a complicated constant, which is related to the atom mass, vortex mass, trap frequency and the density distribution along the third spatial dimension. The Eq.(16) shows that $`\omega `$ is a decreasing function of $`\mathrm{\Omega }`$ when $`\mathrm{\Omega }/\omega _{}>1/\sqrt{2}`$, which is shown in the left side of Fig.(2).

From above discussion we can see the atomic current toward the center is of essential importance for causing oscillations. Hence if there is any other ways to generate a similar atomic current, such as using optical dipole force to draw atoms into the middle of the condensate in the Ref., the similar oscillation of vortex lattice was observed.

### 3.4 Damping Effect

One feature of these oscillation motions shown in these experiments is that the oscillation amplitudes are heavily damped, especially for the Tkachenko oscillation reported in Ref.. It is quite natural in the view of the dual picture. Because the coupling between the vortex current and the phonon field is a gauge coupling, an accelerated vortex will radiate sound wave in the same way as an accelerated electron radiating electromagnetic wave. In the view of hydrodynamic picture, when a vortex is accelerating, the circular current $`\stackrel{}{J}`$, as well as the singularities of field $`\stackrel{}{J}`$, can not be stationary in any inertial frame. Due to the current conservation condition Eq.(4), a time dependent density fluctuation $`\delta \rho (\stackrel{}{r},t)`$ will be generated. Hence the oscillating vortices will induce sound wave, lose their energy and result in the damping of amplitude. In the dual picture this process can be calculated through a standard MED calculation.

The retarded Green’s function for $`2+1D`$ MED is

$$G(\stackrel{~}{r},\stackrel{~}{r}^{})=\frac{\mathrm{\Theta }(tt^{})}{2\pi }\frac{\mathrm{\Theta }(s^2(tt^{})^2(\stackrel{}{r}\stackrel{}{r}^{})^2)}{\sqrt{s^2(tt^{})^2(\stackrel{}{r}\stackrel{}{r}^{})^2}}.$$
(17)
And the gauge potential $`\stackrel{~}{b}_{rad}`$ for the radiation field is

$$\stackrel{~}{b}_{rad}=\underset{\mathrm{}}{\overset{t|\stackrel{}{r}\stackrel{}{r}(t^{})|/s}{}}𝑑t^{}\frac{\stackrel{~}{V}(t^{})}{\sqrt{s^2(tt^{})^2(\stackrel{}{r}\stackrel{}{r}(t^{}))^2}}.$$
(18)
Here $`\stackrel{~}{V}(t^{})`$ is the vortex $`3`$-current. Unlike the $`3+1`$D case, the Green’s function contains the step function $`\mathrm{\Theta }`$ instead of Dirac’s $`\delta `$ function, so one can not easily find a general expression for the radiation energy of a charge with arbitrary acceleration. Here we only consider two special cases, the linear oscillation $`\stackrel{}{V}(t^{})=(v\mathrm{cos}(\omega t^{}),0)`$, and circular motion $`\stackrel{}{V}(t^{})=(v\mathrm{sin}(\omega t^{}),v\mathrm{cos}(\omega t^{}))`$. Following the approximation made in Ref., we replace $`|\stackrel{}{r}\stackrel{}{r}(t^{})|`$ with $`R`$ and calculate the energy radiation along large circle $`|r|=R`$. Then we can perform an analytical calculation for the average energy radiated over one period,

$$P=\frac{1}{T}\underset{0}{\overset{T}{}}\underset{0}{\overset{2\pi }{}}𝑑t𝑑\theta r\stackrel{}{\sigma }\widehat{r},$$
(19)
where $`\stackrel{}{\sigma }`$ is the Poynting vector. For these two special cases, the result is

$$P=\frac{\rho \mathrm{}^2\omega v^2}{2ms^2}=\frac{\mathrm{}^2\omega v^2}{2g}.$$
(20)
The radiation rate is inverse proportional to the interaction strength $`g`$, so the damping effect is remarkable for the weak interacting boson gases. Recalling the approximate expression for the vortex mass $`m_v=\overline{\rho }\mathrm{}^2\mathrm{ln}(R/\mathrm{\Lambda })/(4\pi s^2m)`$, Eq.(20) can be reexpressed as

$`P={\displaystyle \frac{1}{2}}{\displaystyle \frac{\rho \mathrm{}^2}{s^2m}}\omega v^2={\displaystyle \frac{1}{2}}{\displaystyle \frac{4\pi \omega }{\mathrm{ln}\frac{R}{\mathrm{\Lambda }}}}m_vv^2={\displaystyle \frac{4\pi \omega }{\mathrm{ln}\frac{R}{\mathrm{\Lambda }}}}E.`$
Hence

$$\frac{dE}{dt}=\frac{4\pi \omega }{\mathrm{ln}\frac{R}{\mathrm{\Lambda }}}E,$$
(21)
this indicates that the oscillation energy exponentially decays as $`\mathrm{exp}\{t/\tau \}`$, and the ratio of the exponential damping time $`\tau `$ to the oscillation period $`T`$ is equal to $`\mathrm{ln}\frac{R}{\mathrm{\Lambda }}/(8\pi ^2)`$. A typical oscillation curve in the current experiment conditions is schematically shown in the right side of Fig.3, one can see that the oscillation amplitude will decay rapidly within one period, and this is qualitatively coincident with the experimental observation. However, there are also a little quantitatively disagreement between the theoretical damping curve and the experiments, in Ref. the amplitude vanishes within about three periods and in Ref. it is about one and a half periods. This disagreement occurs because in the real case the vortex mass should be modified due to the spatial density inhomogeneity, the third spatial dimension effect and the edge effect. The mass of vortex, which is proportional to the self-energy, is therefore changed with the position of vortex. It is expected that the damping curve will agree better with the experiments when these subtle effects are carefully considered.

## 4 Summary

According to the effective Maxwell electrodynamic model the vortices in the BEC can be viewed as an ensemble of classical charges moving in an external magnetic field. So far we have applied this description to understand the observed phenomena of vortices induced by selective removal atoms from a fast rotating BEC. The atom removal laser applied to the center of the condensate will cause two leading effects in the view of dual picture. One is the increase of the vortex chemical potential and then subsequently the increase of vortex density. The other is that it will induce a circular electric field, and each vortex will acquire an azimuthal velocity. We clarified that there are two forces that can cause oscillation motion, and discussed the the relation between the oscillation frequencies and the number of atoms together with the effective rotation frequency. We showed that the dispersion mechanism for the vortices oscillation comes from the emission of phonon, which is the same as the radiation of the electromagnetic wave of an accelerated charge. And the damping effect is more remarkable in these experiments because the interaction strength is weak.

In this paper the dual description of vortices has been used to explain current experiments. Although it is not yet a detailed quantitative description, it provides a natural physics picture of most key experimental observations, which convinces us to believe that this description is valid to the extent that the BEC is fast rotated but still not entering the quantum Hall mean field regime. We hope this will start a new route to study the issues in the rapidly rotating BEC system.

Acknowledgements: HZ would like to thank Professor C. N. Yang for encouragement. And the authors would like to acknowledge Z.Y. Weng, T.L. Ho, L. Chang, R. L$`\ddot{u}`$ and X.L. Qi for helpful discussions. This work is supported by National Natural Science Foundation of China ( Grant No. 10247002 )

## References