text:
# GOTPM: A Parallel Hybrid Particle-Mesh Treecode

## 1 Introduction

As we enter the era of precision cosmology, the need for higher resolution cosmological simulations has never been clearer. Recent determinations of the cosmological parameters to accuracies of a few percent by the WMAP mission (Bennett et al 2003) will entirely transform the industry of cosmological N-body simulations. Since the initial conditions are now well-defined, the focus will be on seeking detailed quantitative consistency with the prevailing cosmological paradigm with observations of the large-scale structure through clustering and lensing studies and detailed analysis of the structure and dynamics of galaxy clusters and galaxies themselves. During the last decade, there has been a huge growth in computational power brought on by steady increase in processor speed and the integration of processors into massively parallel architectures. Numerical cosmologists and dynamicists have developed parallel, N-body codes to take advantage these architectures and there are now a variety of codes that are being used to simulate the formation of large-scale structure and galaxies. In this paper, we describe a new parallel code that is adaptable and scalable to the next generation of parallel supercomputers.

Two classes of cosmological N-body codes are now used widely and are based on either the Particle-Mesh (PM) algorithm (Hockney and Eastwood 1981) or the Barnes-Hut Tree algorithm (1986). PM codes are the fastest available taking advantage of the fast Fourier transform algorithm (FFT) to solve Poisson’s equation. While useful for some cosmological applications, pure PM codes fail to resolve structures on scales smaller than the mesh size. Different approaches have been taken increase the force resolution on the sub-mesh scale. The simplest approach is the P3M algorithm where the mesh softened force is corrected at short range using direct summation of pairwise forces between particles (e.g. Hockney and Eastwood 1981, Efstathiou et al 1985). As systems become clustered, the P3M algorithm scales as O(N$`{}_{}{}^{2}{}_{h}{}^{}`$) where Nh is the size of the largest halos limiting its application to small N or very large cosmological volumes where halos contain modest numbers of particles. Another approach is to lay down refined sub-meshes on regions of high density to compute force refinements. This can be done hierarchically and leads to accurate force determination and code speed up. Couchman’s (1991) AP3M algorithm used in the Hydra code has been highly successful. A slight disadvantage of this approach is that a density threshold criterion is necessary to lay down sub-meshes and there may be discrete changes in force resolution across arbitrarily set boundaries that have unpredictable effects. An alternative mesh-refinement scheme is provided by Kravtsov et al. (1997) in the ART algorithm that lays down refined meshes as particles begin to crowd into mesh cells. A completely different approach uses the BH tree algorithm modified to include periodic boundary conditions using Ewald’s method (Bouchet, Hernquist and Suto 1991). This method has the advantage that forces are accurately determined at all stages of the calculation on all scales. The penalty is that tree codes generally have much larger memory requirements and require significantly more computation to achieve a similar force accuracy than PM based algorithms. Despite these problems, tree codes based on Ewald’s method have been easier to parallelize then AP3M are currently very competitive in many applications (e.g. Davé et al 199?, Springel et al. 2001).

A natural evolution of cosmological codes is the developments of hybrids that can incorporate the high efficiency of PM codes for long range force determination and the advantage of tree codes for short range force corrections. Xu (1995) developed a hybrid Tree/PM (TPM) algorithm that is in the spirit of AP3M and the algorithm has since been enhanced by Bode et al. (2000). BH trees are used in place of refining sub-meshes and constructed around dense clustered regions to determine short range force corrections. This method differs from the Ewald-method treecodes in that short range forces drop to zero at a few mesh lengths so trees only need to be built locally to calculate the force corrections. The algorithm has also been parallelized and is scalable to large processor numbers (e.g. Bode and Ostriker 2003). Bagla has also introduced a hybrid TREEPM algorithm that breaks the forces up into short and a long range parts using Ewald’s original formalism. The main difference from TPM is that short range forces are corrected using the BH tree method throughout the simulation volume with no need to specify a density threshold for laying down a tree.

In this paper, we describe a new Tree/PM hybrid that is in the same spirit as Bagla’s method but with parallelization. We call this new algorithm GOTPM standing for Grid-of-Oct-Trees-Particle-Mesh. In §2, we describe the details of the algorithm including a new parallel PM method based on the FFTW fast Fourier transform package (Frigo et al. 1998) and a method for computing the force refinements using localized BH trees. Since tree refinements are confined to a few mesh lengths, the process of building trees within the simulation volume can be done sequentially leading to large savings in memory compared with Ewald tree codes. We also discuss how to load balance the short range force correction calculations. In §3, we present measures of the force accuracy and timing tests on a variety of machines to measure the code’s performance and parallel scalability. In §4, we describe some preliminary results of various simulations with the largest containing N=$`1024^3`$ particles. In §5, we discuss the future incorporation of multi-timestepping and smoothed particle hydrodynamics and algorithmic improvements that will permit simulations containing $`N=864G`$ particles on next generation parallel machines.

## 2 Parallel TREEPM

### 2.1 Parallel PM

The PM algorithm is well documented so we highlight here the features of our implementation and method of parallelization. We use the PM code developed by Park (1990, 1997) as our starting point and it follows the usual steps:

1. We determine a density field on a mesh from the particle distribution using the triangular-shaped cloud (TSC) mass assignment scheme.
2. A Fast Fourier transform (FFT) of the density field is calculated and then multiplied by the appropriate Green’s function for the TSC scheme to determine the transform of the potential field.
3. An inverse FFT determines the potential field in real space.
4. We interpolate the potential field to determine the gravitational force at each particle position using a 4-point finite difference approximation (FDA).

We have parallelized this algorithm using a simple, slab domain decomposition of the simulation cube using the Message-Passing Interface (MPI) library. Particles are first distributed in slabs of equal width with an integral number of mesh divisions. We determine the density field using the TSC method in each slab in each processor. Since the TSC method weights the particle mass to each of the surrounding 27 grid points from the nearest grid point, we need to communicate the contribution to the density field from particles near the boundaries of the slabs. We achieve this by binning the mass on ghost slices of mesh that are communicated to their neighbours and summed on to the corresponding local slices to find the correct density field quantities.

Once the density fields on slabs in each processor are determined we then use the FFTW library (Frigo et al. 1998) with MPI extensions to calculate the FFTs used to obtain the gravitational potential. The FFTW package requires that the complete mesh be divided in slabs of equal width so we can determine the gravitational potential immediately with a parallel FFT. With the potential, we then calculate the particle forces in each slab using FDA interpolation. Again, we need to import images of mesh slices at the boundaries to determine forces on particles. We communicate the 3 or 4 mesh slices from the front and back sides of the slab and then calculate the force on each particle in each processor in parallel. Our implementation assumes all neighbouring mesh slices are in adjacent slabs so the slab width can be no smaller than 4 mesh divisions. This sets a limits on the number of processors that can be used in a parallel run e.g. only 128 processors can be used with a $`512^3`$ mesh. Although the code will work with thin slabs, we find in practice below that problems can arise with memory imbalance in the tree algorithm if slices are too thin so slabs with more than 8 mesh divisions are preferred.

### 2.2 Force Shaping

In the complete algorithm, we treat the PM force as the long range contribution to the total force and determine the short range contribution using the tree algorithm. The PM method softens the Newtonian potential within a few mesh divisions and also introduces some anisotropy to the force on the mesh size scale. We follow the lead of other codes and smooth the PM potential using a Gaussian filter with an e-folding radius of 0.9 times the mesh size. To characterize the resultant PM force, we calculate the force many times at a point from a single particle placed at random positions within the mesh in a otherwise homogeneous density field. Figure 1 shows usual behaviour with an asymptotic Newtonian force with significant scatter due to anisotropy at the maximum around $`r=1.4\mathrm{\Delta }x`$ where $`\mathrm{\Delta }x`$ is the mesh cell size. Our strategy is to fit these points with a smooth function of radius which asymptotically behaves as $`1/r^2`$ with a best fit shown in Figure 1. After some experimentation, we settled on the rather complex fitting function,

$`\begin{array}{cc}F(x)=\hfill & a\mathrm{tanh}(bx)/x^2ab/(x\mathrm{cosh}^2(bx))\hfill \\ & +cx(1+dx^2)\mathrm{exp}(dx^2)\hfill \\ & +e(1+fx^2+gx^4)\mathrm{exp}(hx^2),\hfill \end{array}`$ (4)
and find it successfully fits the data. The short range force is simply determined as the difference between this smooth empirical function and the Newtonian gravitational force. We also determined the scatter in the expected relative force error by laying down a single particle randomly within the grid and measuring the acceleration due to the mesh plus correction at random points. Figure 2 show that the scatter peaks at a separation of about 2 mesh cell sizes revealing the anisotropy introduced by working on a regular cubic mesh. These errors are intrinsic to PM methods and we see below in more general situations they introduce a tail in the acceleration error distribution in the 1-2% range.

### 2.3 Correcting Forces with BH Trees

The parallel PM force is very efficient at determining long range gravitational forces but because of the finite mesh size the forces from nearby particles are highly softened. The standard way to overcome this deficiency is to define a new correction force law for local particles that is the difference between the desired Newtonian force and the force law determined by the PM method, $`F_c=F_NF_{PM}`$. The correction force has the property that it is nearly Newtonian on the sub-mesh scales but then falls to zero beyond four mesh cell sizes where the PM forces become nearly Newtonian. The gravitational potential is also softened using a Plummer law on the subgrid scale. We typically use a softening radius of $`5\times 10^5L`$, where L is the boxsize i.e. so typical simulations are equivalent in resolution to a PM method with a $`20000^3`$ mesh.

Three approaches have been used to determine $`F_c`$ in PM methods. Direct force evaluation using the particle-particle (P2) method (Hockney and Eastwood 1981), adaptive particle sub-meshes with P2(Couchman 1991), and tree methods (Xu 1995, Bode et al. 2000, and Bagla 2000). Direct force evaluation is the most general method and is easily parallelizable. Unfortunately, when clustering is strong the method grinds to a halt because of the $`O(N^2)`$ nature of the method. Couchman’s (AP3M) method relies on laying a local refined mesh on regions where strong clustering is occuring and calculating correction forces from this sub-mesh. It is highly efficient in serial codes but has proved difficult to parallelize. Finally, Xu’s hybrid TPM algorithm works in a similar way to AP3M but a Barnes-Hut tree is built locally on clustered regions. This method has been parallelized and performs well (e.g. Bode et al. 2003). A disadvantage of both these methods is that one requires an arbitrary criterion to establish whether a mesh volume contains clustered material that should be refined. Furthermore, boundaries between finely resolved and more coarsely resolved volumes are sharp and it is not clear if there are physical consequences to assigning these boundaries.

Bagla (2000) provides a way to characterize both the short and long range forces in his TREEPM algorithim and is completely analogous to P3M in that all local forces everywhere within the volume are corrected. We describe here a parallel version of this method that builds on the parallel PM algorithm above and incorporates most of the efficient tree-building and tree walking software of Dubinski’s (1996) parallel tree code PARTREE. We call our new parallel algorithm GOTPM for Grid of-Oct-Trees-Particle-Mesh. This new code has the generality of a $`P^3M`$ code in that all forces throughout the simulation volume are calculated accurately on the submesh scale with no density criteria imposed for the addition of refinement. The equivalent of mesh refinement is done automatically by building trees everywhere within the simulation volume at every step. We show below that this algorithm also is easy to load balance using algorithms taken from pure parallel tree codes.

### 2.4 GOTPM: A New Tree/PM Hybrid Algorithm

It is easiest to describe this algorithm by comparing it directly to the P3M method. The $`P^2`$ correction force on a particle is determined by summing directly over particles within a sphere centred on the particle which extends out to 2-3 mesh cell sizes. Beyond this boundary corrections on the mesh force are small, generally less than 1%, and are so neglected. The problem with $`P^2`$ correction, however, is that the algorithm scales inherently as $`O(N_h^2)`$ and a code is therefore slowed down considerably when clustering occurs. In constrast, the tree algorithm is $`O(N\mathrm{log}N)`$ and generally requires about 2 to 3 times the amount of work when going from a homogeneous to a clustered state.

The strategy we employ is to replace $`P^2`$ force correction with a correction determined by the tree method. The details of the basic BH tree method as implemented here are presented in Dubinski (1996) with the only difference being that the force law is given as the difference between the Newtonian and PM force law given in equation (1). In brief, the BH tree algorithm works by arranging particles in a hierarchy of cubes or cells in an oct-tree structure. Forces are determined by descending the tree and summing contributions from interactions with distributions of particles in “distant” cells and nearby particles. Cell interactions are determined from the quadrupole order expansion of the potential of the particles within the cube rather than through direct summation so there are great computational savings. The trick is to introduce a cell opening criterion based on the ratio of the size of a cube to the distance of its centre of mass, i.e. an opening angle. If the cell is big and nearby, then the opening criterion is satisfied such that the cell is broken down into its 8 sub-cells each of which is examined in turn. If the cell is small and sufficiently distant, the opening criterion fails and the quadrupole order expansion of the particles in the cell are used to determine the force and the tree descent ends at this node. In practice the tree method requires only $`O(N\mathrm{log}N)`$ computations per step to determine forces with an error of 1% and therefore has a great advantage over direct summation methods. There are different choices available for the cell opening criterion as well as strategies for grouping particles to minimize the overhead in tree walks. The cell opening criterion and grouping strategy we use here is the same as that described in Dubinski (1996) with some slight modifications that we discuss below.

The application of the tree method to determine force corrections in a PM code requires a few general modifications. The force correction effectively drops to zero at a few mesh lengths. It is therefore not necessary to build a global tree for the entire particle distribution in the simulation cube. Since the forces are localized, we can build a grid of trees, with a dimension of a few mesh lengths. If we define $`R_{max}`$ as the radius where the correction forces drops to small fraction of a percent then we only need to calculate the contributions from particles within this radius. In P3M, one generally sets $`R_{max}=23`$ but we find in practice for the tree algorithm that $`R_{max}=4`$ is necessary to assure accuracy. This scale is set to the length of the cubes enclosing the trees. Each particle is incorporated into one tree so to determine the force correction on a particle we need only descend the particle’s home tree and the surrounding 26 trees (Figure 3).

In principle, P2 force corrections are necessary from all particles within the volume although particles beyond $`R_{max}`$ contribute little to the correction. In practice, particles beyond $`R_{max}`$ are simply ignored saving a lot of time. This is necessary since the number of particles surrounding the point where the potential is calculated grows rapidly with radius i.e. $`Nr^3`$ for a homogeneous system and $`Nr`$ when clusters with densities of $`\rho r^2`$ develop. In a tree code, however, the number of interactions does not grow rapidly with radius. Distant particles are grouped into larger and larger cells so the number of cell interactions only grows roughly as $`N\mathrm{log}r`$ under most circumstances. It is therefore not necessary to impose a hard spherical boundary of $`R_{max}`$ as is done in the $`P^3M`$ method. As long as the trees are $`R_{max}`$ in size then all effective contributions to the force correction will be accounted for and the extra forces from the irregular jagged region beyond the spherical boundary will be small and generally average out. The force corrections themselves will have an inherent error which depends on the value of the critical opening angle in the tree algorithm and the value of $`R_{max}`$. In practice, we show below that as clustering develops the number of particle cell interactions only grows slowly from the initial homogeneous state to the strongly clustered regime. Computationally, the CPU time per step only grows by a factor of three from the beginning to the end of a simulation.

Since the trees are localized and compact there is one further refinement that can be added to improve accuracy and efficiency. If the number of particles within a cube is small, the cost of building and walking the tree for force determinations can be greater than simply doing a direct forces summation. We introduce an additional parameter, $`N_{tree}`$, that defines the minimum number of particles required for a tree build. If the number of particles in a bin is less than $`N_{tree}`$ then we do a direct P2 summation for all particles within the radius $`R_{max}`$. Otherwise, we use the tree machinery for the force calculation. In this mode, a typical run behaves as a P3M code at early times while at late times trees are used to calculate forces wherever clusters develop.

Since force corrections are local, we see that this algorithm is easily parallelizable within the context of the PM code we have described. We outline the implementation of this algorithm below.

### 2.5 Slice domain decomposition

The parallel PM code we describe above distributes the particles into domain slabs of equal width and assigns them to independent processors on a parallel computer. We retain this slab domain decomposition scheme for force corrections with trees but relax the constraint of equal widths for the slices to achieve better load balance. As the system evolves, particles cross these boundaries and must be moved between processors by exchanging messages.

The inner loop of the algorithm is as follows:

1. Particles are redistributed into domain slices of equal width by exchanging messages for the next PM force evaluation.
2. The PM force is determined using the slice domain decomposition parallel method described above.
3. Particles are now redistributed into slices of variable width such that the amount of treecode computation per slice is equal. The boundaries of these domains are determined iteratively using the number of tree force interactions from the last step. It is a one-dimensional version of the orthogonal recursive bisection (ORB) domain decomposition method described in Dubinski (1996) and used in other parallel tree codes like GADGET (Springel et al. 2001). At the end of this step, each processor contains a list of particles which requires the same amount of computation to determine tree force corrections.
4. We are now ready to build a local grid of trees for the force correction. To do this, we need to import “ghost” particles from neighbouring domains extending a distance $`R_{max}`$ from the current slice boundaries. This is accomplished by exchanging particles with neighbouring processors. These ghost particles along with the local particles are then binned into cubes of dimension $`R_{max}`$ in the slice. To save memory, we only import the positions (and masses if variable masses are present) for the ghost particles.
5. For each particle bin, we then build a BH tree.
6. We now loop through the local particle list determining the force correction using tree descents from the surrounding 27 trees centred on each particle. We remember the number of interactions required for each particle for load balancing purposes in step 3 above.
7. Force corrections are now added to the PM force to determine the corrected Newtonian force in a periodic universe.
8. Particles are pushed forward according to the equations of motion and a leapfrog integration scheme.

Figure 4 illustrates the domain decomposition scheme needed for the FFT’s and achieve load balance for short range tree forces.

There are a few refinements to this algorithm that improve speed and minimize memory overhead. First, step 3 is generally expensive in terms of communication time since the domain boundaries must be determined iteratively in $`\mathrm{log}_2N_{proc}`$ message exchanges each possibly requiring the communication of large numbers of particles. To speed this step up, the full ORB decomposition is only done every 10 steps. At intermediate steps the boundaries are held fixed according to the last determination. Since the system does not change dramatically from step to step, only a small fraction of the particles actually cross these boundaries between steps thereby minimizing communication and preserving load balance. Second, we save on memory when building the trees in steps 5 and 6 by only building a small subset of the trees at a time and discarding them after we are done with them. Tree structures inherently require large amounts of data to describe the relational aspects of the tree through pointers and the cells themselves require lots of extra information including the quadrupole components as well as cell dimensional information. Since all force corrections for a particle are local to the surrounding 27 trees, we can build the trees in small groups, determine the correction forces for corresponding particles then discard the trees before building a new set. In practice, we build only $`8\times 8\times n_z`$ grid columns of trees and loop through the volume where $`n_z`$ is the number of trees in the vertical direction. The slice height $`n_z`$ varies from step to step because of the domain decomposition. For example, in a typical simulation on a 32 processor machine with a $`1024^3`$ mesh with $`R_{max}=4.0`$, there are on average $`256\times 256\times 8=512K`$ trees that need to be built in each processor. We only need to build 512 trees at a time with this method saving a factor of 1000 in memory for tree data structures. The only information one needs to save at each step is the basic particle data, i.e. positions, velocities, particle index and computational load. Extra memory is required to handle neighbour particles and memory imbalance induced by the load-balancing algorithm itself. For small cosmological volumes with a great deal of clustering, the memory imbalance leads to roughly a factor of 2 variation in the number of particles per processor. In the worst case, we require about 3 times the memory required for basic particle data (masses, positions and accelerations) to handle this problem.

### 2.6 Memory budget

The memory budget for a pure n-body simulation can be broken down into particle memory and mesh memory. A typical mesh requires $`4\times N_{mesh}^3`$ bytes for single precision and twice that for double precision. The particles themselves require memory for positions, velocities, accelerations, masses and softening lengths, a computational work measurement, a particle index to retain identity and a pointer for constructing linked lists when binning particles in cubes. The total budget is the equivalent of 56 bytes per particle in single precision or 112 bytes in double precision. Because of imbalance in particle numbers in the domain slabs the worst case that occurs with the smallest box sizes leads to domains with up to 3 times the average number of particles. Conservatively, approximately 150 bytes are required per particle in a typical single precision run. The ghost particles can potentially eat up significant amounts of memory when the domain slabs become thin. This deficiency can be removed in future versions by using fully 3-D domain decomposition in the tree correction phase though considerable bandwidth may be required to move particles back and forth from slabs in the PM phase and cuboids in the tree phase.

In summary, a typical run with $`N_{mesh}=2048^3`$ and $`N_p=1024^3`$ requires $``$200 GB in single precision in the worst case of a small box. This new algorithm is considerably more memory efficient than pure cosmological treecodes that require lots of memory to describe global trees.

## 3 Code Verification

### 3.1 Accuracy Tests

To validate the correctness and accuracy of the code, we calculated the accelerations from a 2M particle dark matter simulation in a $`50h^1`$ Mpc box at $`z=0`$ in a standard cosmology with $`\mathrm{\Omega }=0.3`$ and $`\mathrm{\Lambda }=0.7`$. This resulting particle distribution is strongly clustered and so represented a good test for the method’s accuracy. We used several methods to measure the accelerations. To start, we used a cosmological parallel treecode based on the Ewald method (Davé et al. 1997) to determine accelerations with small errors by setting the opening angle parameter $`\theta =0.1`$ and computing cell-particles interactions to quadrupole order. We used a Plummer softening length of $`5\times 10^5`$ times the box size. This is a costly calculation but the mean relative error in the accelerations $`|\delta a/a|`$ is estimated to be $`<0.01`$%. We consider these accelerations to be exact for the purpose of comparison to other methods. We then calculated acceleration errors using the Ewald method with a more standard value of $`\theta =0.8`$ and the new GOTPM algorithm with $`\theta =0.0`$ (making it a $`P^3M`$ method) and $`\theta =1.0`$ with $`N_{tree}=4`$ and 32. The size of the mesh is $`256^3`$. In the GOTPM algorithm, we expect errors to arise both from the treecode approximation as well as the PM method itself. The technique of fitting a smooth potential to the PM force function minimizes but does not eliminate errors induced by close range anisotropies in the mesh force.

Figure 5 shows the distribution of relative acceleration errors for various methods and parameter combinations. The first result is that the values of accelerations between the Ewald-tree method and GOTPM algorithm agree with each other on the percent level. This result reassures us that the computed accelerations are physically correct since the two numerical methods are very different in their mathematical formulation: the first method is based on a combination of the Ewald expansion of a periodic point mass potential and the second method based on the solution of Poisson’s equation using Fourier methods. We can clearly see the convergence to smaller errors when going from the P3M to GOTPM method with $`\theta =1.0`$ and smaller values of $`N_{tree}`$. The distribution of errors in the P3M represents the “best” one can do to reduce errors in the accelerations using a hybrid PM method of this form. The distribution of errors between the P3M method and GOTPM methods are almost indistinguishable when $`N_{tree}=32`$ while the time taken to compute these forces using P3M requires about 10 times as much CPU time. For larger N simulations, the discrepancy in CPU time becomes even greater. The distribution of errors for a Ewald-method treecode using a typical value of $`\theta =0.8`$. Although, the mean error is slightly less than the GOTPM method there is a broader tail in the distribution to larger errors precluding the use of a larger value of $`\theta `$. The Ewald-method treecode method requires about 3 times as much CPU time to compute accelerations of comparable accuracy to the GOTPM method in the clustered state. We note, however, the situation is very different when the systems are less clustered. The GOTPM method is several times faster and more accurate when the system is less clustered since it operates essentially like a pure PM method. In constrast, the Ewald tree method is only slightly faster in the less clustered case and acceleration errors tend to be much larger for a given value of $`\theta `$.

Using a larger mesh size also increases the accuracies of the accelerations and as we see below can also requires less computation so if enough memory exists it can be advantageous to use a finer mesh for the PM part of the calculation.

Also apparent is a tail of the distribution of errors that remains as $`\theta `$ decreases. The cause of this tail is due to the anisotropy of the PM force on the mesh scale which introduces a random error on the determined acceleration for particles. This error is inherent to the PM method itself and all techniques using meshes will have this problem. We also note that the acceleration can be near zero because of the periodic nature of the force so relative errors are amplified in this case and do not necessarily indicate low accuracy. A plot of the acceleration versus the relative error indeed shows that the largest acceleration errors occur for the smallest accelerations. For $`\theta =1`$, the acceleration errors are found to be $`0.43\%\pm 0.56\%`$, the range typically tolerated in treecode simulations. There does not seem to be a large gain in accuracy in going to smaller $`\theta `$ because of mesh-induced force errors. All methods base on PM forces plus correction will have these errors so we recommend using $`\theta =1`$ as a compromise between speed and accuracy since generally 0.5% errors are usually tolerated in collisionless n-body calculations. We note, however, that the relative errors are much smaller for particles feeling large accelerations in the centre of clusters in these simulations so the 0.5% error is a conservative estimate.

### 3.2 Code Performance

The GOTPM code has been run on a variety platforms including a 32 processor IBM SP3 in Seoul, an HP/Compaq GS320 and a Beowulf cluster based on Intel Xeon 2.4 GHz chips with gigabit networking both at CITA. We present results of the code performance and parallel scaling up to 128 CPUS on the Beowulf cluster though runs with up to 256 CPUS have been done and perform well. Our reference simulation computes the evolution of $`N=256^3`$ particles using a $`512^3`$ mesh in a 200 h-1 Mpc box in a standard cosmology with $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Lambda }=0.7`$.

To study scaling of the code speed with processor number, we follow the computation for 5 equal time steps from $`z=1`$ when the system is moderately clustered and time the PM and tree phase of the code separately to determine a speed. These timing tests are done on a Beowulf cluster with gigabit networking and dual Xeon 2.4 GHz processors with code compiled using the Intel C and Fortran compilers. Figure 6 shows the code speed versus the number of processors $`N_p`$ from 4 to 128 processors. The speed scaling is not perfect since the amount of communication time grows as $`N_p`$ increases. As $`N_p`$ grows, the domain slabs become thinner and the number of particles crossing boundaries through both particle redistribution for load-balancing purposes grows significantly as well as the number of ghost particles imported for tree builds. Nevertheless, the code scales well to $`N_p=128`$ requiring only 11.6 second per step and clocking at about 1.4 million particles/s. The PM portion of the calculation only requires 1.8 seconds compared with 9.8 seconds for the tree correction force so the code is dominated by determination of the short range forces.

As part of J. Kim’s doctoral work, we simulated $`1024^3`$ particles in LCDM and SCDM cosmological models in a cosmological box of size $`L=512h^1`$ Mpc. The simulations were run on a IBM SP3 (Regatta) system at the Korea Institute of Science and Technology Information (KISTI). Over the course of two months, we used one node of 32 CPU’s with 208 Gbytes of memory for two GOTPM simulations. Memory limited us to a $`1024^3`$ mesh. Each model was simulated with 460 time steps in 16 wall-clock days and the peak memory usage was about 170 Giga bytes. The simulation details will be presented in Kim’s thesis and a forthcoming publication. The CPU timing of the LCDM simulation is displayed in Figure 7. The CPU time per step grew from 40 minutes to 55 minutes largely due to the increase of tree calculation overload caused by the clustering of particles. Periodic spikes in DD and the total CPU time plot are due to the particle backup at every 15 time steps and pre-halo finding calculations every 5 time steps. The increase of CPU time in the FDA part was partially due to the steady increasing rates of cache miss of particle data in L2 cache. As the system evolves, the memory locality of the particle data differs significantly from the spatial positions of the particles. In Figure 8, we show the inhomogeneity of particle distributions of PM domain slabs in the LCDM simulation. We expect the density fluctuation$`\mathrm{\Delta }(r)`$ at length scale $`r`$ is related to the power index $`n`$ as (Longair 1998)

$$\mathrm{\Delta }(r)=\left(\frac{r}{8h^1\mathrm{Mpc}}\right)^{\frac{n+3}{2}},$$
(5)
where $`r`$ is defined to be the radius of a sphere having equal volume of the local domain slab. In the case of $`n=1,2`$, $`\mathrm{\Delta }(M)`$ of the simulation domain slabs is $`0.16,0.03`$, respectively. Hence the static domain decomposition of the PM part does not matter seriously to memory budget overload for cosmological volume.

## 4 Mass Functions

Another useful tool for the verification of the code is the comparison of simulated halo mass functions with other analytic halo mass functions. There have been several fitting functions proposed by many authors for halo mass functions. Press and Schecter(1974,hereafter PS), Lee and Shandarin(1998,hereafter LS) Sheth and Tormen(1999,hereafter ST) and Jenkins et al.(2000) have presented analytic and fitting mass functions. On the basis of spherical collapses and Gaussian fluctuations in the initial density field, Press and Schecter derived the mass function as,

$`\begin{array}{cc}\mathrm{d}N(M)=\sqrt{\frac{2}{\pi }}\frac{\rho _0}{M}\frac{\delta _c}{\sigma }\mathrm{exp}\left(\frac{\delta _c^2}{2\sigma ^2}\right),\hfill & \end{array}`$ (7)
where $`\delta _c`$ is defined $`1.686`$ and $`\sigma (M)`$ could be calculated from the initial power spectrum $`P(k)`$. Other analytic halo findings are rather empirical functions to fit to their simulated halo mass functions. We display the halo mass function for the 1G particle LCDM simulation with other analytic halo mass functions in Figure 9. In a forthcoming paper, we will describe a halo-finding method used to extract virialized halos from the simulation data during run time. The simulated halo mass function of mass above $`2\times 10^{14}h^1\mathrm{M}_{}`$ is well fitted by the PS predictions but underestimates halo numbers by about 30% below this mass scale. Even with the shortcomings at the high-mass end, our halo mass functions show global consistency with ST and Jenkins predictions.

We also note that the PM version of the code has already been used to do a series of $`512^3`$ particle simulations to model gravitational lensing convergence and shear for comparison to lensing surveys (Pen et al. 2003).

## 5 Conclusions

GOTPM is a new parallel hybrid Tree/PM cosmological N-body method that is fast, accurate, memory efficient and scaleable to hundreds of processors on parallel machines. We have verified the correctness and accuracy of the forces through comparison with other methods and produced a mass function consistent with other results in the literature. This new code is well-suited for the next generation of large-scale cosmological simulations on massively parallel machines.

At present the code has been used to compute $`N=1G`$ particle simulations but to meet the challenges of doing even larger calculations on future supercomputers with thousands of processors some modifications will be necessary. The current version uses a single leapfrog timestep. The incorporation of hierarchical multiple timestepping should be straightforward. Since the rapidly changing force field in dense regions will be mainly calculated as a tree correction, the PM part of the potential can be held fixed or extrapolated over the largest system time step. Since all near range forces are localized to particles in the algorithm only trees in the densest regions need to be updated at smaller fine-grained timesteps while remaining trees can be repaired using methods described by Springel et al. (2001). Furthermore, the communication needs will be small at the fine-grained timesteps since forces are all determined locally and few particles will cross domain boundaries. Work is proceeding to incorporate multiple timestepping.

One major weakness of the algorithm from the perspective of parallelization and memory needs is the use of slab domain decomposition for the tree correction phase. As the number of CPU’s increase slabs become thinner especially when the system becomes more inhomogeneous. This leads to dramatic increase in both communication to move ghost particles as well as memory usage since the imported ghost particles can actually dominate local memory. For $`N=1G`$ simulations with up to 256 processors the algorithm works well but slab domain decomposition will break down on future machines with thousands of processors. One solution to this weakness, is the use of 3-D cuboid domain decomposition instead of slabs for the tree phase. This will solve the ghost particle memory and communication problem. The PM phase will still require slabs since this is inherent to FFTW’s MPI implementation. Either particle data will have to be moved back and forth between cuboids and slabs every step or the mesh points themselves will have to be communicated. These various strategies will be examined in coming versions and will prepare the way for $`N=864`$G particle simulations.

Another weakness of the algorithm is the scatter in the PM mesh force anisotropy. We can minimize this problem by introducing a spherically symmetric mass allocation kernel. The extra cost of this scheme will still be small compared to the tree phase of the algorithm. The advantage of this method is that we can potentially reduce the size of $`R_{max}`$ and save on computation during the tree force corrections.

The code is also well suited for incorporation of smoothed particle hydrodynamics (SPH). The grid of localized oct-trees can be exploited immediately to form neighbour lists for SPH estimation of density and pressure fields. With the multiple timestepping scheme described above, GOTPM with SPH is a natural fit.

The code is currently available by request but will be released publicly in the near future.

## Acknowledgments

JD thanks NSERC and the Canadian Foundation for Innovation and acknowledges the CITA pre-doc programme that allowed JK’s visit. JK and CP thank the Brain Korea 21 program of the Korean Government and the Basic Research Program of the Korea Science & Engineering Foundation (grant no. 1999-2-113-001-5). JK and CP also acknowledge the support from the KISTI (Korea Institute of Science and Technology Information) under ‘Grand Challenge Support Program’ with Dr. S. M. Lee as the technical supporter. The use of the computing system of the Supercomputing Center is also greatly appreciated.


format: text

text:
# Detection of abundant solid methanol toward young low mass stars Based on observations obtained at the European Southern Observatory, Paranal, Chile, within the observing program 69.C-0441.

## 1 Introduction

The presence and origin of complex organic molecules in protostellar regions and their possible incorporation in protoplanetary disks is an active topic of research, both observationally and through laboratory simulations. Large organic molecules such as CH3OCH3 and CH2CH3CN have been detected with enhanced abundances in so-called ‘hot cores’ around massive young stars (e.g. Hatchell et al. 1998; Gibb 2001). In the laboratory, UV photolysis of ice mixtures of species prepared with interstellar abundances followed by thermal warm-up to 300 K is known to lead to a wealth of complex species such as carboxylic acids and esters (Briggs et al. 1992), hexamethylene-tetramine (Bernstein et al. 1995) and prebiotic molecules including amino acids (Munõz Caro et al. 2002; Bernstein et al. 2002). In the ice experiments, methanol is thought to be a key molecule in the production of these complex molecules. Understanding the methanol content of interstellar ice mantles and its variation in different circumstellar environments is therefore essential to test the validy of the basic assumptions that enter the chemical models and experimental approach.

So far, solid methanol has only been detected along lines of sight associated with intermediate- to high-mass protostars, with abundances ranging from 5% up to 30% with respect to the dominant ice mantle component, H2O (Dartois et al. 1999). Extensive searches toward low-mass protostars have resulted in typical upper limits of a few % (Chiar et al. 1996; Brooke et al. 1999). In contrast, gas-phase methanol has been found toward both high- and low-mass YSOs with evidence for abundance jumps of factors of 100–1000 in the inner part of the envelope, consistent with ice evaporation (e.g. van der Tak et al. 2000; Schöier et al. 2002; Bachiller et al. 1995). Still, the observed gas-phase column densities of methanol are usually $`1010^4`$ times less than those observed in the solid phase (e.g. Schöier et al. 2002). In hot cores, a rapid high-temperature gas-phase chemistry is triggered by the evaporation of methanol-rich ices, which forms high abundances of large organic molecules for a period of $`10^410^5`$ yr (e.g. Charnley et al. 1992; Rodgers & Charnley 2001). Understanding the origin of solid $`\mathrm{CH}_3\mathrm{OH}`$ and its variation from source to source is a key ingredient for these models.

As part of a large VLT-ISAAC 3–5 $`\mu `$m spectral survey toward young low mass stars we have detected an absorption feature at $`3.53\mu \mathrm{m}`$ towards at least three icy sources out of $`40`$, which is attributed to solid methanol. Two of the sources are part of the dense cluster of young low-mass stars, SVS 4, located in the Serpens cloud core (Eiroa & Casali 1989). The third source is a more isolated low-mass class I source in the Chameleon I cloud known as Cha INa 2 (Persi et al. 1999). The spectra presented here represent the first detections of solid methanol toward low mass YSOs. The observations and analysis are reported in §2, whereas the implications for the chemical evolution of circumstellar matter are discussed in §3.

## 2 Observations

### 2.1 The $`L`$-band spectra

The observations were performed using the Infrared Spectrometer and Array Camera (ISAAC) mounted on the Very Large Telescope UT1-Antu of the European Southern Observatory at Cerro Paranal. $`L`$-band spectra were obtained using the low resolution mode and the 0$`\stackrel{}{.}`$3 slit resulting in resolving powers of $`\lambda /\mathrm{\Delta }\lambda =1200`$. SVS 4-5 and SVS 4-9 were observed simultaneously on May 5, 2002. Cha INa 2 was observed on May 6, 2002. All spectra were recorded during good weather conditions. The data were reduced using our own IDL routines, which are described in Pontoppidan et al. (2003). Correction for telluric absorption features was done using the early-type standard stars BS 4773 (B5V) and BS 7348 (B8V) for Cha INa 2 and the SVS 4 sources, respectively. The wavelength calibration was obtained using the telluric absorption features with an estimated accuracy of $`0.001\mu \mathrm{m}`$. The spectra were flux calibrated using the standard stars with an estimated uncertainty of 15%, mainly due to uncertainties in slit losses.

### 2.2 Derivation of methanol abundances

The water band profile was obtained by adopting a blackbody continuum fitted to a $`K`$-band photometric point from the literature and the observed $`L`$-band spectrum from 4.0 to $`4.15\mu \mathrm{m}`$ as well as the $`M`$-band spectra, which are discussed in detail in Pontoppidan et al. (2003). The $`2.84.8\mu \mathrm{m}`$ spectra are shown in Fig. 1. The accuracy of the $`K`$-band points are crucial for estimating the column density of the water ice bands. However, the $`K`$-band fluxes have considerable uncertainties since many young stellar objects are known to be highly variable and the sources may have changed their flux level since the photometric measurement. Indeed, Cha INa 2 is known to vary significantly in the $`K`$-band (Carpenter et al. 2002) and has apparently continuously brightened with 1.7 mag between March 1996 and April 2000 (Persi et al. 1999; Kenyon & Gómez 2001; Carpenter et al. 2002) to $`K=11.044`$. It is found that at least $`K=10.2`$ is necessary to fit the continuum defined by the ISAAC $`L`$ and $`M`$-band spectra, consistent with the assumption that the source has continued to brighten until May 2002 with the same rate. The continua used for deriving the optical depths of the water bands are shown in Fig. 1 along with continua fitted using extreme $`K`$-band fluxes. For the Serpens sources a conservative estimate of the $`K`$-band variability is 0.2 mag when comparing with the photometry of Eiroa & Casali (1992). For Cha INa 2 continua using $`K=10`$ and $`K=12.75`$ are also shown.

A local third order polynomial was used to estimate the continuum around the methanol feature (see Fig. 2). The polynomial was fitted to the points between $`3.1\mu \mathrm{m}`$ and $`3.2\mu \mathrm{m}`$, and the region between $`3.63\mu \mathrm{m}`$ and $`3.70\mu \mathrm{m}`$. The confidence in the local continuum for extraction of the methanol profiles is high, partly due to the low order of the polynomial and partly due to the low curvature required for the three spectra. The uncertainty in the derived abundances is thus dominated by the main water band continuum. The measured water ice abundances are presented in Table 1. The values were derived by fitting a laboratory spectrum of pure amorphous water ice at 50 K to the observed water bands and integrating the scaled laboratory spectrum from 2.5 to $`3.7\mu \mathrm{m}`$. None of the observed water bands are saturated. However, due to the low S/N ratio in SVS 4-5, the fit in this source relies on the wings in the water band. The adopted band strength is $`2.0\times 10^{16}\mathrm{cm}\mathrm{molec}^1`$ (Gerakines et al. 1995). Note that laboratory measurements of solid state band strengths may in general only be accurate within 30%.

The final absorption profiles in the 3.3–$`3.7\mu \mathrm{m}`$ region are shown in Fig. 3 on an optical depth scale. Common to the profiles of SVS 4-9 and Cha INa 2 is a broad feature centered at $`3.47\mu \mathrm{m}`$ along with a narrow feature at $`3.53\mu \mathrm{m}`$ seen in all three sources. The former is most likely the same feature seen in this wavelength region toward many other young stars, generally known as the ‘$`3.47\mu \mathrm{m}`$ feature’ and thought to be due to diamonds or an ammonia hydrate. The 3.53 $`\mu `$m feature is seen toward much fewer objects and is assigned to the $`\nu _3`$ CH-stretching vibration band of solid methanol.

To estimate the true depth of the methanol band, it is necessary to remove the contribution from the $`3.47\mu \mathrm{m}`$ feature. A first-order approximation is to use a high signal-to-noise profile observed toward another YSO. We chose here the feature observed toward the intermediate mass YSO Reipurth 50 described in Dartois et al. (2003). The feature has a center position of $`3.478\mu \mathrm{m}`$ and a $`FWHM`$ of $`0.12\mu \mathrm{m}`$. $`3.47\mu \mathrm{m}`$ features from other sources give similar results (Brooke et al. 1999). The $`3.47\mu \mathrm{m}`$ band from Reipurth 50 was fitted to the observed spectra simultaneously with a laboratory methanol band.

It is well-known that the $`3.53\mu \mathrm{m}`$ methanol band is sensitive to the structure and composition of the ice (Ehrenfreund et al. 1999). Laboratory spectra are compared to the observed methanol bands in Fig. 3. A laboratory spectrum of pure methanol does not give a perfect fit of the exact profile to any of the observed spectra, mainly due to the relative weakness of the secondary band of solid methanol at $`3.40\mu \mathrm{m}`$ in the observed spectra. It is interesting to note that the ratio of the $`3.40\mu \mathrm{m}`$ and the $`3.53\mu \mathrm{m}`$ methanol features in low-mass sources are similar to that of the high-mass sources of Dartois et al. (1999). A laboratory spectrum rich in water and formaldehyde ($`\mathrm{H}_2\mathrm{CO}`$) can reproduce the spectra. However, the H2CO identification was found to be inconsistent for the massive stars W 33A and RAFGL 7009S (Dartois et al. 1999) due to the absence of bands at 5.8 and $`6.69\mu \mathrm{m}`$. Other mechanisms of reducing the strength of the $`3.40\mu \mathrm{m}`$ band cannot be ruled out. For example, irradiated pure methanol spectra annealed to temperatures higher than 60 K also gives good matches. For the estimate of the $`\mathrm{CH}_3\mathrm{OH}`$ column densities it is assumed that $`\mathrm{H}_2\mathrm{CO}`$ does not contribute to the $`3.53\mu \mathrm{m}`$ band.

The adopted band strength of the CH-stretching mode at $`3.53\mu \mathrm{m}`$ is $`5.3\times 10^{18}\mathrm{cm}\mathrm{molec}^1`$. This value is not found to vary significantly depending on the ice mixture (Kerkhof et al. 1999). The derived abundances are summarized in Table 1.

Brooke et al. (1999) and Chiar et al. (1996) find in general upper limits on the methanol abundance in the solid phase of a few % with respect to water. We find consistent results for the rest of our observed sample, i.e. methanol with an abundance of more than 5% is found toward only 3 out of a sample of $`40`$ low mass YSOs.

## 3 Chemical implications

The methanol can be produced in space and in laboratory experiments through several different pathways. Gas-phase chemistry produces only low abundances of CH3OH, of order 10-9 with respect to H2, so that simple freeze-out of gas-phase CH3OH cannot reproduce the observed solid state abundances of $``$$`10^5`$. A possible route is therefore thought to be grain surface hydrogenation reactions of solid CO with atomic H, leading to H2CO and eventually CH3OH (Tielens & Hagen 1982, Charnley et al. 1997). Experimental results on this reaction scheme at low temperatures are still controversial. Watanabe & Kouchi (2002) find efficient formation of $`\mathrm{H}_2\mathrm{CO}`$ and $`\mathrm{CH}_3\mathrm{OH}`$ by warm ($``$ 300 K) hydrogen addition to CO in an $`\mathrm{H}_2\mathrm{O}`$-CO ice mixture, but the temperature of the hydrogen is not relevant for quiescent dark clouds. Other authors find yields which show that the formation of methanol directly via hydrogenation of CO is inefficient in dark clouds in the 10–25 K regime (Hiraoka et al. 2002). Further experiments are needed to solve this discrepancy.

Theoretical models of grain surface chemistry indicate that the fraction of CH3OH ice with respect to H2O ice, as well as the H2CO/CH3OH ratio, depend on the ratio of H accretion to CO accretion on the grains and thus on the density of the cloud (e.g., Charnley et al. 1997, Keane et al. 2002). In the framework of these models, solid CH3OH abundances of 15–30% with respect to H2O ice such as found here can only be produced at low densities, $`n10^4`$ cm-3. Moreover, the region needs to be cold, $`T_{\mathrm{dust}}<20`$ K, to retain sufficient CO on the grain surfaces. The models thus suggest that solid CH3OH formation occurs primarily in the cold pre-stellar phase. However, this fails to explain the large range of methanol abundances observed along lines of sight toward low mass YSOs in the same star forming cloud, unless the duration of the pre-stellar phase varies strongly from source to source. As discussed by van der Tak et al. (2000), comparison of solid CH3OH with solid CO2 can further constrain the H/O ratio of the accreting gas and thus the density and duration of the pre-stellar phase. For low-mass YSOs, data on solid CO2 are not yet available, but will be possible with the Space Infrared Telescope Facility (SIRTF) via the $`15.2\mu \mathrm{m}`$ $`\mathrm{CO}_2`$ bending mode.

An alternative formation route of solid CH3OH may be through ultraviolet processing (Schutte 1988) or proton irradiation of ices (Hudson & Moore 1999). However, as shown by van der Tak et al. (2000), the CH3OH production rates in these experiments are orders of magnitude lower than those through grain surface formation. The fact that some low-mass sources show similarly high solid CH3OH abundances as many massive sources argues further against ultraviolet radiation playing a significant role.

It is clear that any formation scheme of solid CH3OH must be able to explain local ice mantle abundances of at least 15–30% compared to H2O ice, regardless of the potential differences in density, temperature and ultraviolet radiation field between high- and low-mass YSO environments. At the same time, such models need to explain the large variations in solid CH3OH abundances from object to object, with abundances of less than 3–5% found for both low- and high-mass YSOs (Chiar et al. 1996; Brooke et al. 1999). In this context, it is interesting that abundant solid methanol is now detected toward two sources which are members of the same dense cluster, SVS 4-5 and 4-9, while other nearby sources in Serpens such as CK 1 show solid methanol abundances of less that a few % . This suggests that the ice mantle chemistry may differ significantly between different regions even in the same star-forming clouds or that the formation and subsequent presence of methanol on grain surfaces is a transient phenomenon related to strong heating of the ices. In the last scheme the methanol may not exist for long in the solid phase before being desorbed by the continuous heating of the ice. Finally shock processing may produce a different ice chemistry (Bergin et al. 1998). Making a survey of the SVS 4 cluster lines of sight may help to decipher the relation of the methanol ice to the environmental differences between ‘quiescent’ environments as opposed to the massive and complex star-forming environments probed so far.

###### Acknowledgements.

The authors wish to thank the VLT staff for assistance and advice, in particular Chris Lidman for many helpful comments on the data reduction. This research was supported by the Netherlands Organization for Scientific Research (NWO) grant 614.041.004, the Netherlands Research School for Astronomy (NOVA) and a NWO Spinoza grant.


text:
# Ground State of a Heisenberg Chain with Next-Nearest-Neighbor Bond Alternation

## Abstract

We investigate the ground-state properties of the spin-half $`J_1J_2`$ Heisenberg chain with a next-nearest-neighbor spin-Peierls dimerization using conformal field theory and Lanczos exact diagonalizations. In agreement with the results of a recent bosonization analysis by Sarkar and Sen, we find that for small frustration ($`J_2/J_1`$) the system is in a Luttinger spin-fluid phase, with gapless excitations, and a finite spin-wave velocity. In the regime of strong frustration the ground state is spontaneously dimerized and the bond alternation reduces the triplet gap, leading to a slight enhancement of the critical point separating the Luttinger phase from the gapped one. An accurate determination of the phase boundary is obtained numerically from the study of the excitation spectrum.

Asymmetric spin-half ladders have recently attracted much theoretical interest voit ; cinesini ; sen ; saw ; naka  due to possible experimental realizations in delafossite cuprates such as $`\mathrm{YCuO}_{2.5}`$ (Ref. dela, ) and the unusual physical effects that the asymmetry in the leg exchanges could introduce both in the ground-state correlations and in the properties of the excitation spectrum. The simplest and most transparent Hamiltonian representation of these systems can be given in terms of frustrated ($`J_1J_2`$) spin-half Heisenberg chains with a spin-Peierls dimerization in the next-nearest-neighbor (n.n.n.) interaction, namely

$$\widehat{}=\underset{i}{}\left\{J_1\widehat{𝐒}_i\widehat{𝐒}_{i+1}+J_2[1+(1)^i\delta ]\widehat{𝐒}_i\widehat{𝐒}_{i+2}\right\},$$
(1)
where $`J_1,J_20`$ are the nearest- and next-nearest neighbor exchange couplings, respectively and $`0\delta 1`$ is the bond alternation parameter. In the following, all the energies are expressed in units of $`J_1`$. For $`\delta =0`$ the properties of the model and the existence of a quantum phase transition at $`(J_2/J_1)_c0.241`$ are well established. affleck ; j1j21d  Below the critical point the system is in a Luttinger spin-fluid phase with a gapless spectrum, no broken symmetry, and a finite spin-wave velocity. Spin correlations are characterized by a power-law decay defined by the same exponent of the pure Heisenberg chain ($`J_2=0`$). In the regime of strong frustration, instead, the ground state is spontaneously dimerized, two-fold degenerate and adiabatically connected to the Majumdar-Gosh exact solution for $`J_2/J_1=0.5`$ mg ; ss  with exponentially decaying spin correlations.

The effect of the bond alternation on the n.n.n. interaction has been recently investigated using field-theoretical approaches voit ; cinesini  with remarkable disagreement in the conclusions. In fact, Chen and collaborators voit  using bosonization and a renormalization group analysis claimed that in the limit of small frustration the asymmetry in the n.n.n. integrals destabilizes the isotropic Heisenberg fixed point leading to a new phase with gapless excitations and vanishing spin-wave velocity. On the other hand, with a similar bosonization approach Sarkar and Sen cinesini  found that bond alternation represents an irrelevant perturbation in the regime of small frustration. In this paper, using conformal field theory we show that the n.n.n. spin-Peierls operator, i.e., the one associated with the alternation in the n.n.n. exchange,

$$\widehat{O}_{n.n.n.}=\underset{l}{}(1)^l\widehat{𝐒}_l\widehat{𝐒}_{l+2},$$
(2)
represents indeed an irrelevant perturbation for the Heisenberg chain in the regime of weak frustration. Furthermore, using Lanczos diagonalization technique we find that in the dimerized phase the effect of the bond alternation is to reduce the spin gap leading therefore to a slight enhancement with $`\delta `$ of the critical value $`(J_2/J_1)_c`$ for the fluid-dimer transition. Finally, we present an accurate numerical determination of the shift introduced by the bond alternation on the phase boundary between the Luttinger and the dimerized phase. A preliminary account of this work was already presented in Ref. beto, .

In general, the relevance of the operator (2) for the Luttinger fixed point can be determined by calculating its scaling dimension $`X`$. This can be defined in terms of the size-dependence of its generalized susceptibility

$$\chi _O=\frac{2}{L}\psi _0|\widehat{O}(E_0\widehat{})^1\widehat{O}|\psi _0L^{2(1X)},$$
(3)
where $`L`$ is the number of sites of the ring, so that $`X<1`$ characterizes a relevant operator. bibbia  Here $`|\psi _0`$ represents the unperturbed Luttinger ground state. The calculation of $`X`$ can be done by identifying the quantum numbers $`\{j\}`$ (referenced to the ground state with energy $`E_0`$) of the intermediate states $`|j`$ appearing in the Lehmann representation of $`\chi _O`$bibbia  In fact, in a conformally invariant field theory, the scaling dimension $`X`$ of a given operator is related to the finite-size corrections of the energy of the lowest intermediate eigenstate $`|j`$ by

$$\mathrm{\Delta }E(L)=E_j(L)E_0(L)=2\pi v_sX/L,$$
(4)
where $`v_s`$ is the spin-wave velocity. In the case of the n.n.n. operator (2), the quantum numbers of the intermediate states, referenced to the ground state, are the total spin $`\widehat{𝐒}=0`$, momentum $`k=\pi `$ and even parity under spatial reflection ($`lLl`$), $`R=1`$. In all cases, the only difference with the nearest-neighbor (n.n.) spin-Peierls operator,

$$\widehat{O}_{n.n.}=\underset{l}{}(1)^l\widehat{𝐒}_l\widehat{𝐒}_{l+1},$$
(5)
the well-known relevant perturbation, is the spatial-reflection quantum number, $`R=1`$. The finite-size corrections can be computed either by Bethe ansatz (for $`J_2/J_1=0`$) or by bosonization. According to the standard notation, at the Luttinger fixed point, the effective Hamiltonian describing spin fluctuations is given by a free boson theory: affleckrev

$$H_L=\frac{v_s}{2}𝑑x\left\{K\mathrm{\Pi }^2(x)+\frac{1}{K}[_x\mathrm{\Phi }(x)]^2\right\},$$
(6)
where $`v_s`$ is the spin velocity, $`\mathrm{\Phi }(x)`$ is the spin fluctuation field and $`\mathrm{\Pi }(x)`$ is the conjugated moment. The dimensionless coupling constant $`K`$ is the key parameter of the bosonization theory, whose value can be analytically determined at $`J_2=0`$ by comparing the low-lying excitation spectrum of the Hamiltonian (6) with the Bethe ansatz solution of the lattice model: $`K=(22\mathrm{cos}^1\lambda /\pi )^1`$. Here $`0\lambda 1`$ is the coupling constant defining the easy-plane anisotropy: at the spin-isotropic point $`\lambda =1`$, leading to $`K=1/2`$. The finite-size corrections of the lowest excitation energies of this Hamiltonian are given by: bibbianoi

$$E_{n,m}(L)E_0(L)=\frac{2\pi v_s}{L}\left(\frac{n^2}{4K}+Km^2\right).$$
(7)
For each pair of integers $`n`$ and $`m`$ ($`n=m=0`$ being the ground state), a low-energy excited state with energy $`E_{n,m}`$ is defined by choosing appropriate boundary conditions for the fields:

$`\mathrm{\Phi }(L)\mathrm{\Phi }(0)`$ $`=`$ $`n\sqrt{\pi },`$ (8)
$`{\displaystyle \underset{0}{\overset{L}{}}}𝑑x\mathrm{\Pi }(x)`$ $`=`$ $`2m\sqrt{\pi }.`$ (9)
Using the relation between the boson fields and the physical electron operator ones, it was then derived in Ref. affleckrev,  that $`n`$ is the $`z`$ component of the spin, while $`m+n`$ is the momentum of the low-energy excitation in units of $`\pi `$. Following this work, it can be shown that the reflection symmetry $``$ acts on the fields as follows:

$`^{}\mathrm{\Phi }(x)`$ $`=`$ $`\mathrm{\Phi }(Lx),`$ (10)
$`^{}\mathrm{\Pi }(x)`$ $`=`$ $`\mathrm{\Pi }(Lx),`$ (11)
implying that the quantum number $`m`$ changes sign under spatial reflection $``$, whereas $`n`$ does not, and, therefore, the ground state has $`R=1`$. For $`m0`$, each pair $`\pm m`$ leads to two degenerate eigenstates with opposite reflections, $`R=\pm 1`$.

Clearly, the excitations of Eq. (7) do not exhaust the full energy spectrum of the Luttinger model because, on top of them, the usual boson spectrum is also present:

$$E_k(L)=v_s\underset{k0}{}|k|a_k^{}a_k,$$
(12)
where $`a_k^{}`$ is a boson field defined in terms of $`\mathrm{\Phi }(x)`$ and $`\mathrm{\Pi }(x)`$, which satisfies $`^{}a_k^{}=a_k^{}`$.

For both the spin-Peierls operators previously defined, the appropriate choices of quantum numbers are $`n=0`$ (i.e., $`S^z=0`$) and $`m=\mathrm{odd}`$ (i.e., momentum $`\pi `$). The two lowest excitations with $`m=\pm 1`$ and $`n=0`$ can be combined to construct two states with defined reflection parity, $`R=\pm 1`$. In the spin-isotropic case, $`K=1/2`$, the excitations with ($`m=\pm 1`$, $`n=0`$) and ($`m=0`$, $`n=\pm 1`$) are degenerate, the latter ones corresponding to a $`R=1`$ triplet (the spatial reflection does not change $`n`$). Therefore, the state with $`n=0`$ and $`R=1`$ is the $`S^z=0`$ state of the triplet, and the remaining state with $`n=0`$ and $`R=1`$ has to be a singlet. Analogously, all finite-$`m`$ excitations with $`n=0`$ are forbidden to be singlets with $`R=1`$ by the spin isotropy. In order to obtain a low-energy singlet excitation with $`R=1`$, we have to combine the $`m=1`$ ($`R=1`$) excitation with at least three $`a_k^{}`$ bosons with zero total momentum. In fact one boson cannot have zero momentum and two bosons cannot change the spatial reflection symmetry. The minimum energy of such state is readily evaluated \[with $`k_1=k_2=2\pi /L`$, and $`k_3=(k_1+k_2)`$\] as

$$\mathrm{\Delta }E(L)=2\pi v_s(K+4)/L,$$
(13)
leading to the result $`X=K+4`$ for any operator with the same symmetries of the operator (2).

Although the previous equation has been extracted in the isotropic case ($`\lambda =1`$), where $`K=1/2`$, it is indeed valid in general for $`0\lambda 1`$, as can be verified, analytically, in the XY model (i.e., at $`\lambda =0`$) or, numerically, for $`\lambda <1`$ (not shown). The reason is that, in the entire gapless regime, the low-energy excitations are described by Eqs. (7) and (12). Notice that, for $`\lambda 1`$, the total spin is no longer a good quantum number, nevertheless the parity under spin-reflection \[$`(S^x,S^y,S^z)(S^x,S^y,S^z)`$\] can still be used to label the low-lying states.

The previous analysis shows that the finite-size spectrum of the Luttinger chain does not contain low-lying singlets with $`k=\pi `$ and even parity under spatial-reflections, yielding to a critical exponent $`X>1`$ for the operator (2), and therefore proving its irrelevance. This is also illustrated in Fig. 1, where we report exact diagonalization results in the weak frustration regime of the isotropic chain. notesw  However, we stress that our conclusions follow directly from the bosonization technique and represent an exact result for the isotropic spin-half Heisenberg chain.

The irrelevance of the n.n.n. spin-Peierls perturbation in the regime of weak frustration can be also confirmed numerically from the Lanczos calculation of the associated generalized susceptibility. As explained in Ref. plaquetto, , this can be done by adding to the Hamiltonian the considered perturbation $`\eta \widehat{O}`$ and then computing the ground-state energy per site $`e(\eta )`$ for few values of $`\eta `$ in order to estimate numerically $`\chi =d^2e(\eta )/d\eta ^2|_{\eta =0}`$ as the limit $`\chi =lim_{\eta 0}2(e_0e(\eta ))/\eta ^2`$. For instance, for the $`28`$-site chain at $`J_2/J_1=0.1`$, this procedure gives $`\chi _{nn}6.4`$ and $`\chi _{nnn}10^6`$, for the n.n. and n.n.n spin-Peierls operator, respectively. This fact indicates that the response of the system to a n.n.n. dimerization is more than six orders of magnitude smaller than the response to the standard n.n. dimerization.

Entering the regime of strong frustration, the effect of the n.n.n. alternation is to reduce the spin gap $`\mathrm{\Delta }`$ as it is shown in Fig. 2 at $`J_2/J_1=0.5`$. As it has been already pointed out, voit ; cinesini  at the Majumdar-Gosh point ($`J_2/J_1=0.5`$) the presence of the n.n.n. spin-Peierls dimerization does not alter the ground-state manifold which is spanned by the two well-known dimerized states, with energy $`E_0=3/8N`$, independent of $`\delta `$. Nevertheless, the excited states are not known exactly so that we have calculated the spin gap using the Lanczos technique. By extrapolating the finite-size data according to the expected law, scaloppo  $`\mathrm{\Delta }(L)=\mathrm{\Delta }_{\mathrm{}}+a/L^2+b/L^4`$, we have obtained the $`\delta `$ dependence of the thermodynamic gap, which is shown in the same figure. Our results are in good agreement with the variational calculation of Refs. voit, ; sen, ; ss, , and indicate a $`8\%`$ spin-gap reduction from the Majumdar-Gosh model for $`\delta =0`$ ($`\mathrm{\Delta }=0.235`$) to the so-called saw-tooth chain for $`\delta =1`$ ($`\mathrm{\Delta }=0.216`$). note2

The reduction of the spin gap in the dimerized phase induced by the n.n.n. bond alternation leads to an increase of the critical value $`(J_2/J_1)_c`$ for the fluid-dimer transition. As originally suggested by Castilla and collaborators, emery1d  the phase boundary between the two phases can be calculated very accurately by exploiting the change in the excitation spectrum occurring at the critical point. In fact, below the critical point the system displays quasi-long-range antiferromagnetic order: in this case, gapless excitations are obtained by creating low-energy spinons both in triplet (spatially even) and singlet (spatially odd) states, the former having the lowest energy on finite size. In the dimerized phase, instead, the triplet excitations acquire a gap so that the ground state remains two-fold degenerate in the singlet sector, leading to the dimerized ground state which breaks the translational invariance. As a result, the critical coupling $`(J_2/J_1)_c`$ can be determined by performing a size scaling of the value $`J_2/J_1`$ where the level crossing between the low-lying singlet and triplet occurs (see Fig. 3). With this technique, using exact diagonalizations up to $`30`$ sites, it is easy to obtain very accurate determinations of the critical value $`(J_2/J_1)_c`$. For instance, in absence of n.n.n. bond alternation, we obtain a critical value $`(J_2/J_1)_c`$ in agreement with previously given estimates. j1j21d ; emery1d  The resulting phase boundary for the fluid-dimer transition is shown in Fig. 4, indicating that for $`\delta =1`$ the critical value of the frustration ratio for the fluid-dimer transition is only increased by about $`1\%`$.

In conclusion, we have studied the effects of a next-nearest-neighbor spin-Peierls dimerization on the ground state of the spin-half $`J_1J_2`$ Heisenberg chain, using conformal field theory and exact diagonalization. We have found that, despite the conclusions of a recent study based on bosonization and renormalization group approaches, voit  the bond alternation represents an irrelevant perturbation in both the regime of strong and weak frustration. In particular, as predicted by Sarkar and Sen, cinesini  the perturbation does not change the nature of the excitation spectrum, leading only to a weak reduction of the triplet gap in the dimerized phase, and, therefore, to a slight enhancement of the critical value $`(J_2/J_1)_c`$ associated to the fluid-dimer transition.

We would like to thank D. Sen for having attracted our attention to this model and for useful correspondence. L.C. and A.P. acknowledge kind hospitality at SISSA. This work was partially supported by MIUR (COFIN01) and INFM (PAIS-MALODI). L.C. was supported by NSF grant DMR-9817242.