text:
# On the scaling of computational particle physics codes on cluster computers

## 1 Introduction

A compute-intensive application might have to be distributed over many compute nodes in order for it to run at a viable speed. One needs to know how well a platform can support the required level of granularity. In this paper we examine this in the case of cluster computers.

The use of clusters of computers in high performance computing has achieved some success with the price/performance advantage of using low-cost and yet high-speed commercial off-the-shelf (COTS) processors to build a fast machine with the budget available. The significant presence of clusters in the TOP500 list  attests to the power of this approach.

Cluster computers are also becoming more widely accepted with the growing ease with which they can be configured and maintained; there are now numerous open source and commercial “Cluster in a box” packages available. It can also be important that portable application code can be developed and run on these machines, using standard operating systems, libraries, tools and parallel APIs.

Essential for the use of clusters in many areas of high performance computing is the presence of high performance system area network technology. The question is whether this is good enough to make parallel computation on clusters viable in the face of very fast communications built into MPP or SMP/NUMA machines.

In this paper we intend to make some quantitative statements on the subject by examining the scaling of some codes from our parallel application, viz. Lattice QCD. In section 2 we give a brief description of Lattice QCD, emphasising the computational task involved. Then in section 3 we describe how our code implements this with particular reference to the parallelisation. We proceed to describe in section 4 the scaling tests we apply and present their results. These are summarised and discussed in section 5, after which we present our conclusions.

## 2 Lattice QCD

Lattice QCD is a numerical evaluation of Quantum Chromodynamics, the theory of the strong interaction, which takes its name from the use of a four-dimensional hypercuboid lattice to represent a discretised spacetime (see e.g.  for a recent overview of the subject). In the lattice theory, each lattice site (vertex of the lattice) represents a spacetime point.

These computations are demanding enough that they routinely make much use of supercomputers, and the grid-based formulation of the problem makes it most suitable for data decomposition on parallel machines.

One of the central defining features of a Lattice QCD calculation is the discretised representation of the Dirac fermion operator which describes the interaction of quarks in the theory. On the lattice, the discretised operator takes the form of a fermion matrix. The algorithms used for Lattice QCD calculations are such that the bulk of the computer time is spent in the inversion of this matrix. There is a variety of representations of the fermion matrix which have different properties on the lattice but they should all become identical to the continuum Dirac operator in the limit of zero lattice spacing. The particular formalism used in this work is that originally due to Wilson . The Wilson fermion matrix is

$$M_{xy}=\delta _{xy}\kappa D_{xy}$$
(1)
where the indices $`x`$ and $`y`$ refer to lattice sites and $`\kappa \text{}`$ is the hopping parameter, related to the bare quark mass.

The Wilson hopping matrix is

$$D_{xy}=\underset{\mu =0}{\overset{3}{}}(1\gamma _\mu )U_{\mu x}\delta _{yx+\widehat{\mu }}+(1+\gamma _\mu )U_{\mu x\widehat{\mu }}^{}\delta _{yx\widehat{\mu }}$$
(2)
where the lattice gauge fields $`U_{\mu x}SU(3)\text{}^{3\times 3}`$ carry colour indices and the Dirac $`\gamma `$-matrices $`\gamma _\mu \text{}^{4\times 4}`$: $`\{\gamma _\mu ,\gamma _\nu \}=2\delta _{\mu \nu }`$ carry spin indices.

Therefore $`M`$ is a matrix in lattice site $``$ spin $``$ colour space. If the number of sites in the lattice is $`V`$ then the fermion matrix is a complex matrix of size $`(12V)^2`$ acting on vectors consisting of $`12V`$ complex numbers.

The lattice is of finite extent, and we impose periodic or antiperiodic boundary conditions.

Note that the lattice site structure of this matrix is either local or connects nearest-neighbour sites. Therefore it is a sparse matrix and is implemented as an operator rather than stored explicitly

Other formulations of the fermion matrix include the Sheikholeslami-Wohlert discretisation which adds to the Wilson matrix a $`34\times 34`$ complex matrix diagonal in lattice sites  , variations of the Kogut-Susskind matrix  which have no spin components but can have a more than nearest-neighbour connectivity  and the Fixed Point formulations which have spin and colour indices and more than nearest-neighbour connectivity . The choice depends on the desired properties, and has consequences for the balance of computation to communication required; for example, the Sheikoleslami-Wohlert matrix has more local computation than the Wilson matrix with no extra communication and will therefore scale better.

The Wilson matrix has the favourable property of only requiring nearest-neighbour communications but since its computational requirements are comparatively modest it does put some pressure on the scaling abilities of a parallel machine. Although it is one of the older formulations it is still interesting not least because it can be used as the basis for the new breed of Ginsparg-Wilson lattice fermion formulations , which have very favourable field-theoretic properties but are still very expensive to compute with; a single multiplication by the GW fermion matrix involves many ($`𝒪(100)`$) multiplications by the Wilson matrix. For these reasons we believe that our test results with the Wilson matrix are relevant. However we also consider the effects of a fermion matrix requiring more than nearest-neighbour communications.

## 3 Description of Lattice QCD software

We give here some details of the implementation of our lattice QCD applications, describing the computational task itself, how the code is written and how the code is parallelised.

Much of the code used in these tests was actually developed for a specific Lattice QCD project called GRAL 111GRAL stands for “Going Realistic And Light”. . The GRAL project has as its aim the finite size scaling analysis of QCD physics close to the chiral regime of light dynamical quarks using the Wilson formulation. This requires a large computational effort, and although a cluster is just one of the machines that is used for the project, good performance of the Lattice QCD code on it is important to the project’s success.

The bulk of the code is written in C compiled with the Compaq C Compiler ccc with a set of optimisations enabled with the -fast flag. We use assembly language kernels  for certain low-level numerically intensive routines where performance is critical. Routines from the optimised BLAS library are also used where appropriate. Both 32-bit and 64-bit code is implemented, but since 64-bit accuracy is preferred for the GRAL project, we present here only 64-bit results.

The communications are entirely handled with message passing via MPI. The fast MPI-library on ALiCE is provided by the ParaStation™ cluster middleware .

### 3.1 Fermion matrix inversion

As stated in section 2 most computer time will be spent in those portions of this code which invert the fermion matrix. This involves solving the linear system

$$M\psi =b$$
(3)
where $`M`$ is the fermion matrix defined in section 2, $`b`$ is the source vector and $`\psi `$ is the solution to be found. A Krylov subspace algorithm is best suited to the solution of such a sparse linear system. We have implemented several such algorithms but here we show results using the BiCGStab  method which is usually the fastest. The communication requirements (described in section 3.3) of these different algorithms are comparable.

An additional refinement is the use of a parallel (locally lexicographic) SSOR preconditioning  to reduce the number of iterations of the BiCGStab solver. This requires a slightly different pattern of communication and so we present these results separately.

### 3.2 Grid decomposition

The four-dimensional lattice is decomposed over a grid of processing nodes. This grid can be from zero dimensions (i.e., a single node) to four dimensions (since the lattice is four-dimensional). The number of nodes along any dimension is restricted in that the local lattice on each node must have identical size and shape. This constraint is not strictly necessary on a MIMD machine but it makes the implementation simpler and ensures perfect load balancing.

### 3.3 Nearest-neighbour communication requirements

Every iteration of the solver algorithm requires three global summations of a complex number and performs two multiplications by the fermion matrix.

Each fermion matrix multiplication requires nearest-neighbour communication of the multiplicand vector. On a lattice boundary in a non-parallelised direction this is taken care of by the boundary condition. In a parallelised direction the required data is on another node. With the message passing communication paradigm, this means that the data on those lattice sites on the boundary of a local lattice is sent to the neighbouring node. Each node receives data from neighbouring nodes and stores it in buffers. This communication pattern is illustrated in figure 1.

The amount of data to be sent depends on the size of the lattice and the grid decomposition. If these are such that on each node of a grid of dimensionality $`d>0`$ the local lattice has $`V`$ sites, then the number of sites on the boundaries from which data must be sent is

$$2V\underset{\mu =0}{\overset{d1}{}}\frac{1}{N_\mu },$$
(4)
where the sum is over parallelised directions $`\mu `$ and $`N_\mu `$ is the length of the local lattice in direction $`\mu `$. Recall that a vector consists of 12 complex numbers on each lattice site, which is 192 bytes at 64-bit precision.

For the normal BiCGStab algorithm all the data on a boundary can be sent at once as a single message. We send each boundary before the multiplication by the fermion matrix begins. The order in which they are sent is essentially unimportant.

The SSOR preconditioning requires most of the data to be sent incrementally as a succession of smaller messages at various points within the matrix-vector multiplication; the order is important. This arrangement is potentially less favourable because of a possibly inefficient use of available bandwidth and the effect of multiple latency penalties.

#### 3.3.1 Data layout

The parallelisation efficiency is affected by the way the data is stored along with the design of the communication buffers.

The vectors involved in the calculation are stored as site-major arrays. Each array element contains the 12 complex numbers and the $`n`$th element corresponds to the local lattice site with cartesian coordinates $`(x_0,x_1,x_2,x_3)`$ via $`n=x_3+N_3x_2+N_3N_2x_1+N_3N_2N_1x_0`$ where $`N_\mu `$ is the size of the local lattice in the $`\mu `$th direction. The gauge fields are stored in a similar way.

This storage order is illustrated in figure 1. Note that because of the layout of the data in lattice, the data on the direction 1 boundaries is strided, which takes more time to copy than contiguous data.

The communication buffers are, for convenience, attached to this array, which increases its length by an amount the size of the local lattice boundaries. Two arrangements of the buffers were tried:

The “halo” layout

The “halo” layout places the buffers on either side of the array body, e.g. for a local lattice of size $`N_0\times N_1\times N_2\times N_3`$, the first two buffers would be

$`x_1=1`$$`x_0=1`$$`n=0\mathrm{}N_0N_1N_2N_31`$$`x_0=N_0`$$`x_1=N_1`$

The “tail” layout

An alternative is the “tail” layout:

$`n=0\mathrm{}N_0N_1N_2N_31`$$`x_0=1`$$`x_0=N_0`$$`x_1=1`$$`x_1=N_1`$

where the buffers are attached to the end of the array. In section 4.4.1 we test the differences between these two layouts.

### 3.4 Next-to-nearest-neighbour communication
requirements

Although the Wilson fermion matrix requires only nearest-neighbour communication of the gauge field $`U`$, for algorithmic reasons we also need next-to-nearest-neighbour access222 In a typical GRAL application the gauge field is communicated $`𝒪(100)`$ times less often than the spin-colour vector, therefore the efficiency of the gauge field communications is not so important for our considerations..

Furthermore, as mentioned in section 1, some fermion matrices connect more than nearest-neighbour lattice sites. We investigate the effects of this by considering next-to-nearest-neighbour communication. This is already possible on a 1-dimensional grid with the communication patterns already presented, but on a 2-dimsional grid the communication in direction 1 has to be changed as shown in figure 2. Additional data from the corners of the lattices has to be sent.

For a lattice with $`V`$ sites on each node of a grid of dimensionality $`d>1`$, there are

$$4V\underset{\mu =0}{\overset{d1}{}}\underset{\nu =0}{\overset{\mu 1}{}}\frac{1}{N_\mu N_\nu }$$
(5)
lattice sites on the corners from which data must be sent, where the sum is over parallelised directions $`\mu `$ and $`\nu `$ and $`N_\mu `$ and $`N_\nu `$ are the lengths of the local lattice in those directions. This is in addition to the amount of boundary data described in section 3.3. We wish to test the effect of this additional communication.

#### 3.4.1 Data layout

The next-to-nearest neighbour communication pattern in more than one dimension requires extra buffers on each array for the corner data. The possibilities for different data layouts are therefore increased. We illustrate three possibilities for the layout on a two-dimensional grid (which can be generalised for a greater number of dimensions):

Layout 0

The corner buffers are added as halos to the $`x_1=1`$ and $`x_1=N_1`$ buffers;

$`x_1=1,x_0=1`$$`x_1=1`$$`x_1=1,x_0=N_1`$

$`x_0=1`$$`n=0\mathrm{}N_0N_1N_2N_31`$$`x_0=N_0`$

$`x_1=N_1,x_0=1`$$`x_1=N_1`$$`x_1=N_1,x_0=N_1`$

layout 1

The corner buffers are added to the $`x_1=1`$ and $`x_1=N_1`$ buffers as tails;

$`n=0\mathrm{}N_0N_1N_2N_31`$$`x_0=1`$$`x_0=N_0`$

$`x_1=1`$$`x_1=1,x_0=1`$$`x_1=1,x_0=N_1`$

$`x_1=N_1`$$`x_1=N_1,x_0=1`$$`x_1=N_1,x_0=N_1`$

layout 2

The corner buffers are added as a tail at the end of the array;

$`n=0\mathrm{}N_0N_1N_2N_31`$$`x_0=1`$$`x_0=N_0`$$`x_1=1`$$`x_1=N_1`$

$`x_1=1,x_0=1`$$`x_1=1,x_0=N_1`$$`x_1=N_1,x_0=N_1`$$`x_1=N_1,x_0=1`$

Layouts 0 and 1 have the advantage that the data in the $`x_0=N_0`$ and $`x_0=1`$ buffers with $`x_1=0(N_11)`$ can be sent in the negative (positive) direction 1 along with the data with $`0x_0<N_0\text{}x_1=0(N_11)`$ as a single message.

Layout 2 requires that the corners buffers be sent seperately. However, since it is merely the tail layout of section 3.3.1 with extra buffer space added to the end of it, a single layout can be used for fields both with and without next-to-nearest-neighbour interaction without wasting space.

## 4 The scaling tests

Our goal is to test how well the code introduced in section 3 scales on a cluster. We run tests to see how this is affected by various factors; the data layout, the locality of the fermion matrix and the solver algorithm.

In this section we first describe the essential details of all the scaling tests. We then describe the hardware and software on which these tests are run.

### 4.1 The test specifications

In order to test the scaling of our lattice QCD applications, we choose to test the fermion matrix inversion described in section 3.1 since it is the area in which many QCD programs and in particular the GRAL applications will spend most run-time and it has the considerable communication requirements discussed in section 3.3. All the tests are run in 64-bit precision.

We choose to time the inversion on a $`12^4`$ lattice. This rather arbitrary size is smaller than that used in most realistic lattice QCD calculations. As our machine size is 128 processors, this lattice size should give an idea of how the performance is affected by the size and surface/volume ratio of the sub-lattices on each node when extrapolated to larger computers with 512 or 1024 processors.

Each solve is performed on a random $`SU(3)`$ gauge field with a random gaussian source vector. The solve is repeated many times until a good estimate of the average run-time per solve can be obtained. The run-time on a single node is very well reproducible, but on multiple nodes there is some variation, perhaps due to effects of different physical grid topologies or the load on the network and switches, so the entire run is repeated at least five times to average over such variations.

Table 1 shows the grid topologies used in the tests. Results from grids of the same dimensionality and number of nodes are averaged. The ratio of the communicable surface of the local lattice to the local lattice volume is also shown.

Our main metric of the parallel performance is the speedup;

$$\text{speedup}=\frac{\text{speed on }\text{n}\text{ nodes}}{\text{speed on 1 node}}.$$
(6)
It can also be useful to know the parallel efficiency (scaled speedup);

$$\text{parallel efficiency}=\text{scaled speedup}=\frac{\text{speedup}}{\text{n}}.$$
(7)
### 4.2 The test platform

ALiCE, the Alpha-Linux Cluster Engine at Wuppertal University is a general purpose high-performance computer consisting of 128 Compaq DS10 servers, each with a 616MHz Alpha 21264 processor, 64 kB 2-way set-associative level 1 data and instruction caches and a 2 MB level 2 cache . The Alpha 21264 is 4-way superscalar, sustaining 2 floating-point operations per cycle. The operating system is SuSE Linux .

ALiCE is clustered using Myrinet with the ParaStation3 cluster middleware. This arrangement takes advantage of the open architecture of the Myrinet communication hardware to replace the standard firmware supplied by Myricom with that supplied by Par-Tec, together with kernel extensions and communications libraries. The ParaStation system endeavors to achieve high performance by implementing a reliable communications protocol in firmware and keeping communications in user space. The two-way 1.28 Gbit/s Myrinet network is configured as a multistage crossbar and the Myrinet cards have a 64 bit/33 MHz PCI buses. We use MPICH 1.2.3 as the interface between our application code and the ParaStation middleware. On ALiCE, the latency for ParaStation/MPI is $`16\mu `$s; the bi-directional bandwidth is about 170 MB/s .

### 4.3 Single node performance

In evaluating the following results it will be helpful to know the performance of the solvers on a single node. Figures 3 and 4 show how the speed of the BiCGStab solver and the SSOR preconditioned BiCGStab solver vary with the number of sites in the lattice. Memory cache and memory bandwidth effects mean that performance deteriorates as the lattice volume increases.

BiCGStab sustains a speed of 237 Mflops on our $`12^4`$ test lattice described in section 4.1, while SSOR preconditioned BiCGStab sustains a speed of 220 Mflops. The latter solver is preferred because it still finds the solution faster due to the smaller number of iterations.

### 4.4 The results of the scaling tests

We now examine what happens to the performance of our solvers as we move from a serial to a parallel implementation. There are numerous factors which might affect the parallel performance; we make a preliminary investigation into the effects of data layout and communication semantics before presenting our final results for the speedup and parallel efficiency.

#### 4.4.1 Investigation of data layout

We compare how the two data layouts described in section 3.3.1 affect the scaling of the data. It might be expected that the halo layout would have better spatial locality properties (and therefore use the cache better), but the performance, shown in figure 5, seems to be much the same. This could be because of the data prefetching written into our assembler kernels.

We choose to use the tail format. It is easily adapted to the case of next-to-nearest-neighbour communications (which we need for the gauge field) by simply adding extra buffer space at the end of the array.

#### 4.4.2 Investigation of communication semantics

We can hope to improve performance using non-blocking (asynchronous) communications (e.g MPI\_Isend) rather than blocking (synchronous) communications (e.g MPI\_Send).

Figure 6 shows that this improves BiCGStab performance on a 1-dimensional grid but degrades it on a 2-dimensional grid. We try a semi-blocking approach where non-blocking communications are used in direction 0, then a block allows this communication to complete before non-blocking communications are used in direction 1; this aproach does not work either. We have verified that this is also the case for grids of dimension greater than 2. Henceforth with BiCGStab we use non-blocking communications on one-dimensional grids and blocking communications on grids of more than one dimension.

Figure 7 shows that the situation is different with SSOR preconditioned BiCGStab: non-blocking communications are faster even on the two-dimensional grid. Henceforth we use non-blocking communications in connection with SSOR preconditioned BiCGStab. We suspect that the communication of BiCGStab, where long messages are sent, might negatively interfere with caching and compute processes. A better separation of communication directions with smaller message lengths, as is the case for SSOR, might lead less often to a situation in which data used for compute operations competes for cache lines with data used simultaneously for communication

### 4.5 Results for solver speedup and parallel efficiency

Having determined that the data layout is not very significant and when it is better to use blocking or non-blocking communications, we can present the results for the best scaling of our solvers on one- to four-dimensional grids.

What we see is that for the BiCGStab solver on this test lattice (figures 8 and 9) the one- and two-dimensional grids become quite inefficient on more than a few nodes. On the three- and four-dimensional grids the drop in efficiency on large numbers of nodes is not as dramatic; a reasonable speedup is achieved due to the more favourable surface/volume ratio. The speedup on a four-dimensional grid is not convincingly better than that on a three-dimensional grid, presumably because the data on the direction 3 boundaries is so finely strided as to negate the advantage of better surface/volume ratios.

The SSOR results (figures 10 and 11) show that the speedup is less affected by the dimensionality of the grid i.e. the surface/volume ratio. However the fact that the data on the direction 3 boundaries are so finely strided adversely affects the four-dimensional grid performance, since in this case it means that very many send/receives have to be done, which is evidently too inefficient due to latency effects.

A main result of these investigation is that a small $`12^4`$ lattice still performs with a scaled speedup of 0.45 on 64 processors. We will comment on this in section 5.1.

To shed further light on these results we plot in figures 12 and 13 the percentage of the total wall-clock execution time required for communication against the ratio of the communicable surface of the local lattice to the local lattice volume. The correlation does not usually depend very strongly on the dimensionality of the grid; a notable exception is the SSOR preconditioned solver on the four-dimensional grids. Note also how the dependence is weaker when the solver is SSOR preconditioned.

### 4.6 Results for next-to-nearest-neighbour communication

We examine the overhead introduced by the need to send the extra data from the boundary corners by implementing next-to-nearest-neighbour communications and timing exactly the same BiCGStab solve of the same fermion matrix as before on two-dimensional grids. Any difference is therefore entirely due to the extra comunication. Again it was found that for BiCGStab, blocking communications were fastest, so these tests are performed with blocking communications.

#### 4.6.1 Investigation of data layout

We compare in figure 14 the three layouts defined in section 3.4.1. Because of the seperate messages required, the third layout is the least scalable, but the differences only become significant in the region of poor scaling. Because of this and the fact that for our applications the next-to-nearest-neighbour communications are relatively unimportant, we adopt the third layout for the gauge field in the GRAL code, utilising the benefits of a unified data layout as described in section 3.4.

#### 4.6.2 Results for the speedup

In figure 15 we take the most scalable next-to-nearest-neighbour layout and compare its speedup with the nearest-neighbour communications on a two-dimensional grid, demonstrating the effect of sending the extra data from the corners of the lattice.

This is perhaps a favourable comparison since the amount of extra data that must be sent for next-to-nearest-neighbour communications increases with the dimensionality $`d`$ of the node grid as described in section 3.4

## 5 Summary and conclusions

We have presented the results of tests of the scaling of the Wilson matrix solver for a $`12^4`$ test lattice on ALiCE . Here we summarize those results and discuss their implications.

### 5.1 Summary

On small numbers of nodes the BiCGStab solver scales without much dependence on the grid decomposition. If a minimum effiency of 75% percent is demanded then figure 9 indicates that no more than 4 nodes should be used. However it still maintains around 50% efficiency on about 16–36 nodes on a three- or four-dimensional grid, corresponding to local lattices of $`6^4`$$`4^2\times 6^2`$ or surface/volume ratios of 1.33–1.67. Thereafter the efficiency degrades, but not too sharply.

The SSOR preconditioned solver has similar small grid characteristics; figure 11 suggests that the 75% efficiency bound is reached at around 4–9 nodes. The four-dimensional grids rapidly become very slow but an efficiency of around 50% is maintained on two and three-dimensional grids of around 36 nodes, i.e. a $`4\times 4\times 3\times 12`$ local lattice with surface/volume ratio 1.67 or a $`2^2\times 12^2`$ lattice with surface/volume ratio 2.00. The performance degrades slowly thereafter; we still have 45% efficiency on 64 nodes with a $`3^3\times 12`$ local lattice.

That the scaling performance need not be seriously affected by the data layout is a heartening result for code developers; we are free to organise data in whichever way is convenient.

A curious feature of our system is the performance of the non-blocking communications. It appears that if the communications in different directions are well separated then non-blocking semantics are an improvement. But they are a hindrance if those communications are done at the same time, even when separated by an explicit barrier as in the semi-blocking case. We have argued above that cache effects might be responsible for these differences.

One gratifying result is how well the SSOR preconditioned solver scales. Its algorithmic performance is excellent but its more complicated pattern of communication can in practice lead to scaling problems due to the communication latency. The vice of needing many smaller messages rather than a few large ones appears to become a virtue, since the separation of sends allows the faster non-blocking communications to be used. This demonstrates that a more complicated algorithm need not be deleterious to efficiency and performance.

### 5.2 Discussion

We see that in general the scaling ability is limited by the additional time that has to be used for the communications and this is essentially determined by the surface/volume ratio. Therefore the results presented here can serve as a guide for planning future runs on realistic lattices of larger global size, bearing in mind that increasing the global size introduces additional overheads from the collective communications.

For example, the discussion in the previous section concluded that our $`12^4`$ lattice runs with 75% efficiency on 4 nodes with BiCGstab. According to table 1 this corresponds to a surface/volume ratio of $`\frac{2}{3}`$. Therefore if we require 75% efficiency on any lattice we should aim for grid decomposition which gives a local surface/volume ratio of about $`\frac{2}{3}`$. From equation (4) we obtain the constraint on the dimensions of the local lattice in parallel directions

$$\underset{\mu =0}{\overset{d1}{}}\frac{1}{N_\mu }\frac{1}{3}.$$
(8)
This is achieved for a $`12^3\times 24`$ lattice with a grid of $`1\times 2\times 2\times 2`$ (8 processors) or possibly $`1\times 1\times 2\times 6`$ (12 processors). A $`16^3\times 32`$ lattice needs a $`1\times 2\times 2\times 2`$ (8 processors) or $`1\times 2\times 2\times 4`$ (16 processors) grid. A $`24^3\times 48`$ lattice needs a $`2\times 2\times 2\times 4`$ (32 processors) grid.

With SSOR preconditioning we might be able to use 1.0 as a target surface/volume ratio, and use grids of $`1\times 2\times 2\times 4`$ (16), $`1\times 2\times 4\times 4`$ (32) and $`1\times 4\times 4\times 4`$ (64) respectively. A similar analysis suggests that a lattice of size $`24^4`$ can be expected to run with a scaled efficiency of about 0.45 on 512 processors.

While there must be an optimum grid decomposition for a given problem and problem size, the complicated nature of the relationship between latency, bandwidth and message length makes finding this “sweet spot” a difficult optimisation problem.

In general there is a trade-off between speed and efficiency which must be determined according to the needs of each project and the resources available; one would also take into account the single-node speed when deciding the local lattice size.

The GRAL project typically requires a number of concurrent simulations at different physical parameters which will run for some considerable length of time. The former requirement will probably mean using a modest number of nodes on ALiCE for each run, it being a shared multi-purpose facility, and the latter requirement will favour running at a reasonable efficiency. Fortunately these considerations are entirely compatible and we have shown here that the required performance is available. We also now know that should a run need to proceed at greater speed, albeit less efficiently, then this is also feasible.

The target surface/volume ratio will in general change when a different sort of fermion matrix is used. Our measurements indicate that a fermion matrix with next-to-nearest-neighbour communications might introduce an additional communications overhead of about 10%.

### 5.3 Outlook

It is reasonable to speculate about what these results mean, or how they would change, for future cluster computer platforms. A faster CPU will increase the proportion of time spent in communications relative to computation, but it should be remembered that a faster CPU will speed up the communications too, so the overall effect might not be too different.

One should be aware that network and communications technology is also developing: we can look forward to high-bandwidth PCI buses and faster networks . In fact much of this technology exists already, and as cluster computing becomes more widespread and new network standards (e.g. InfiniBand) become available, current higher-performance technologies (e.g. Gigabit System Network) will become affordable.

A greater network speed helps scalability, but also of importance for the developer, in that it could change the way application code is designed, are implementations of communications APIs that make good use of DMA technology to provide support for simultaneous sends/receives of messages to multiple nodes, or for overlapping communications with computations.

Of great importance to performance of scientific codes are factors like data cache sizes and memory bandwidths. If these are greater then it is more difficult to achieve efficient scaling, but the benefits to the performance on each node might make decomposition still worthwhile. Conceivably, there could then be a situation where a finer grained decomposition actually improves speedup by moving into a régime where all local data fit into cache.

## Acknowledgements

ZS acknowledges the financial support provided through the European Community’s Human Potential Programme under contract HPRN-CT-2000-00145, Hadrons/Lattice QCD. We thank Guido Arnold for support with the ALiCE cluster. The project has been supported by the Deutsche Forschungsgemeinschaft under contract Li701/3-1, Optimierung von Leistung und Datenorganisation von Clusterrechnern der 150 Gflops-Klasse.


format: text

text:
# Ground-state energy fluctuations in the Sherrington-Kirkpatrick model

(August 26, 2003)

## Abstract

The probability distribution function (PDF) of the ground-state energy in the Sherrington-Kirkpatrick spin-glass model is numerically determined by collecting a large statistical sample of ground states, computed using a genetic algorithm. It is shown that the standard deviation of the ground-state energy per spin scales with the number of spins, $`N`$, as $`N^\rho `$ with $`\rho 0.765`$, but the value $`\rho =3/4`$ is also compatible with the data, while the previously proposed value $`\rho =5/6`$ is ruled out. The PDF satisfies finite-size scaling with a non-Gaussian asymptotic PDF, which can be fitted remarkably well by the Gumbel distribution for the $`m`$-th smallest element in a set of random variables, with $`m6`$.

Models with quenched disorder describe a wide class of systems in physics and other fields books ; YOUNG . In thermal systems, the effects of disorder are especially important at low temperatures, where the physics is dominated by the low-energy states. These states are also of direct interest in many applications, for instance to combinatorial search problems books . Understanding the properties of the low-energy states is therefore a central task.

In this context, an important aspect concerns the statistical fluctuations of the energy of the low-lying states with respect to the disorder. Consider a system described by an energy function $`H_N(𝐒,𝐉)`$, where $`𝐒`$ represents $`N`$ degrees of freedom and $`𝐉`$ a set of quenched random variables signifying the disorder. The ground-state energy per degree of freedom, $`e_NN^1\mathrm{min}_𝐒H_N(𝐒,𝐉)`$, when considered as a function of $`𝐉`$, is also random variable, and typically one would like to know its mean $`e_N`$, standard deviation $`\sigma _N=(e_N^2e_N^2)^{1/2}`$, and possibly its entire PDF, $`p_N(e_N)`$. In the thermodynamic limit, $`N\mathrm{}`$, usually the properties of extensivity ($`e_Ne_0`$, with $`e_0`$ a quantity of order one) and self-averaging ($`\sigma _N0`$) hold, so one simply has $`p_N(x)\delta (xe_0)`$. Thus, the interesting question is to determine the finite-$`N`$ scaling behavior. More precisely, one may ask for example:

(i) How does $`e_N`$ scale with $`N`$?

(ii) How does $`\sigma _N`$ scale with $`N`$?

(iii) Given the natural scaling variable $`y_N=(e_Ne_N)/\sigma _N`$, does its PDF $`\stackrel{~}{p}_N(y_N)`$ converge for large $`N`$ to a nontrivial asymptotic PDF, $`\stackrel{~}{p}_{\mathrm{}}(y_N)`$?

To leading order, we expect $`e_N=e_0+bN^\omega `$ and $`\sigma _NN^\rho `$, so (i) and (ii) amount to determine the exponents $`\omega `$ and $`\rho `$. Question (iii) amounts to determine whether the scaling relation

$$p_N(e_N)=\frac{1}{\sigma _N}\stackrel{~}{p}_{\mathrm{}}\left(\frac{e_Ne_N}{\sigma _N}\right)$$
(1)
holds for large $`N`$, and the shape of $`\stackrel{~}{p}_{\mathrm{}}(x)`$. The answer to (ii) and (iii) is simple for systems with short-range interactions, where the volume can be subdivided so that $`e_N`$ is the sum of $`O(N)`$ nearly independent contributions. For the central limit theorem, then, $`\sigma _NN^{1/2}`$ (as rigorously proven for short-range spin glasses aw ) and $`\stackrel{~}{p}_{\mathrm{}}(x)`$ is a Gaussian for moderate values of $`x`$ (we will not be concerned with large deviations here). For sufficiently long interaction range, in general one has $`\rho 1/2`$ and a non-Gaussian $`\stackrel{~}{p}_{\mathrm{}}(x)`$ foot\_diluted .

In this paper, the questions above are addressed numerically for the infinite-range Sherrington-Kirkpatrick (SK) spin-glass model sk . Much is known about this model, following Parisi’s replica-symmetry-breaking (RSB) solution rsb ; books , but many issues remain open, finite-size scaling being a prominent one. The exponents $`\omega `$ and $`\rho `$ remain unknown analytically, while previous numerical work indicates $`\omega 2/3`$ thesis ; krzakala ; boettcher  an $`\rho 3/4`$ cabasino ; krzakala . Determining the full ground-state energy PDF is a much more challenging task, which so far has received limited attention krzakala . Here, we will address these issues with a much higher statistics than in previous studies.

The interest of this problem lies also in its relationship bouchaud  with extreme value statistics gumbel . The goal of extreme value statistics is to determine whether the PDF of the $`m`$-th smallest value $`X`$ in a set of random variables $`X_1,\mathrm{},X_M`$ (the energy levels in our case) approaches a scaling form $`p_M(X)=1/b_M\stackrel{~}{p}_{\mathrm{}}[(Xa_M)/b_M]`$ for large $`M`$, with suitable $`a_M,b_M`$. The case of independent, identically distributed (i.i.d.) $`X_i`$’s is well understood gumbel : if scaling holds, $`\stackrel{~}{p}_{\mathrm{}}(x)`$ is one of three universal distributions, depending on the PDF of the $`X_i`$’s. If the latter is unbounded and decays faster than a power law for large $`|X_i|`$, then $`\stackrel{~}{p}_{\mathrm{}}(x)`$ is the Gumbel PDF (defined below). An example of a disordered system with i.i.d. energy levels is the Random Energy Model rem , but in general the levels are correlated, and this case is much less understood (an exception is when $`e_N`$ is the sum of spatially uncorrelated variables, as mentioned above). The Gumbel form is preserved for short-range correlations and violated for sufficiently long-range correlations gumbel ; bouchaud . It is therefore interesting to ask if it applies, at least approximately, to the SK model, where full RSB leads to hierarchically correlated energy levels (see, however, Ref. dean,  for a model with hierarchical correlations giving a non-universal PDF).

Numerical details – We study the standard SK model with energy function

$$H_N=\underset{i<j}{}J_{ij}S_iS_j$$
(2)
where the $`S_i`$’s are $`N`$ Ising spins ($`S_i=\pm 1`$) and the $`J_{ij}`$’s are i.i.d. random variables drawn from a Gaussian PDF with zero mean and variance $`1/N`$. The sizes studied are $`N=19,39,59,99`$, and $`199`$ with, respectively, 250000, 250000, 180000, 120000, and 64000 independent realizations of the disorder (samples). The ground state of each sample was determined using a genetic algorithm pal ; py . The $`l`$-th order cumulants of $`p_N(e_N)`$, denoted by $`k_{l,N}`$, and their statistical errors were then estimated with a bootstrap method foot\_jack  up to $`l=5`$. With large statistics, it is especially important to control the systematic errors due to occasionally missing the ground state. This was done by performing much longer runs for a fraction of the samples, checking that the change in the cumulants was smaller than the desired statistical error.

Scaling of the mean and standard deviation – Previous numerical studies indicate a value $`\omega 2/3`$ for both the SK model thesis ; krzakala  and diluted spin glasses thesis ; krzakala ; boettcher ; liers . The analytical result $`\omega =2/3`$ was obtained for the finite-size scaling of the internal energy near $`T_c`$ and the Almeida-Thouless line slanina , but no result is available deep in the ordered phase. Our data for $`e_Ne_0N^{2/3}`$, using the exact Parisi result rsb  $`e_0=0.7633\mathrm{}`$, are shown in the inset of Fig. 1. The small residual dependence on $`N`$ foot\_error  either indicates that $`\omega `$ is greater than 2/3 or, if $`\omega =2/3`$, that the data are accurate enough to resolve subleading corrections. A fit with $`e_N=0.7633+bN^\omega `$ gives $`\omega =0.673\pm 0.002`$ ($`\chi ^2=11.2`$, $`\text{ndf}=3`$, ndf being the number of degrees of freedom), where the error gives the range for which the goodness-of-fit $`Q`$ is larger than $`10^3`$. A fit with $`\omega =2/3`$ (leaving $`b`$ as the only free parameter) gives $`Q<10^{16}`$, which is clearly not acceptable. Nevertheless, since corrections of just a few percent could accommodate the value $`\omega =2/3`$ (note the vertical scale in the inset of Fig. 1), this value cannot be ruled out. Fits with corrections (for example $`e_Ne_0N^\omega =b+cN^{\omega _2}`$ with fixed $`e_0`$ and $`\omega `$) are not conclusive since the corrections are small and thus the parameter range is large.

Turning now to the fluctuations, in the replica formalism one can write the cumulant generating function of the (quenched) free-energy as the annealed free energy of $`n`$ replicas, $`f(n)`$. Crisanti et al. crisanti  argued that if the first nonlinear term in the small-$`n`$ expansion of $`f(n)`$ is $`n^i`$ (in the thermodynamic limit) then $`\rho =11/i`$. Kondor kondor  found $`f(n)=nfcn^6`$ from the “truncated model”, which would imply $`\rho =5/6`$. Recently, Aspelmeier and Moore aspel  and De Dominicis and Di Francesco dedominicis  found a higher free-energy (hence better) solution $`f(n)=nf`$, from which the conclusion $`\rho =5/6`$ no longer follows, and Refs. aspel2, ; krzakala,  gave qualitative arguments indicating $`\rho =3/4`$. On the numerical side, Cabasino et al. cabasino  found a result compatible with $`\rho =3/4`$, but with large statistical errors, and Bouchaud et al. krzakala  found $`\rho 0.76`$, also compatible with $`3/4`$, but they could not rule out $`\rho =5/6`$.

Our numerical data for $`\sigma _N`$ are displayed in Fig. 1, multiplied by $`N^{1/2}`$ to stress that $`\sigma _N`$ decays faster with $`N`$ than in short-range models. An exponent $`\rho =3/4`$ follows reasonably well the data, but deviations can be noticed. These are seen more clearly in the plot of $`N^{3/4}\sigma _N`$ in Fig. 2(a). Here the curvature in the data demonstrates the presence of subleading corrections. The dashed line represents a power law fit $`\sigma _N=aN^\rho `$ for $`N59`$, which gives $`\rho =0.765\pm 0.01`$ ($`\chi ^2=3.93,\text{ndf}=1`$), in agreement with $`\rho =3/4`$ (the error corresponds again to $`Q>10^3`$). However, the decreasing trend for large $`N`$ foot\_error  and the negative curvature suggest that $`\rho `$ may in fact be greater than $`3/4`$. As for the mean, fits with correction terms (for fixed $`\rho `$) are not very illuminating. Fig. 1 also shows that $`\rho =5/6`$ clearly does not fit the data, as seen also in Fig. 2(b) where $`N^{5/6}\sigma _N`$ is far from constant. Since the sizes considered are already well in the scaling regime (at least as far as Eq.(1) is concerned, see below), crossover to a constant for larger $`N`$ is rather unlikely, hence we conclude that $`\rho =5/6`$ is ruled out. Finally, if $`\sigma _NN^\rho `$ and Eq.(1) holds, then $`k_{l,N}N^{l\rho }`$. Figs. 2(c) and 2(d) show that also for $`l=3,4`$ the data agree with $`\rho 3/4`$ but not with $`\rho =5/6`$ (note that the statistical error is much larger than that of $`\sigma _N`$).

Scaling of the PDF – Scaling plots of $`p_N(e_N)`$ according to Eq.(1) are shown in Fig. 3 and, on a logarithmic scale, Fig. 4. The data collapse is remarkably good, indicating that we are already well in the scaling regime. Comparison with a Gaussian PDF, represented by the dashed lines in Figs. 3 and 4, clearly indicates a non-Gaussian behavior. Table I displays the scaled cumulants $`\stackrel{~}{k}_{l,N}k_{l,N}/\sigma _N`$ for $`l=3,4,5`$ (since $`\stackrel{~}{k}_{2,N}=1`$ by definition, $`\stackrel{~}{k}_{3,N}`$ is the skewness and $`\stackrel{~}{k}_{4,N}`$ the kurtosis of $`\stackrel{~}{p}_N`$): the cumulants are independent of $`N`$ within the errors, again strongly supporting scaling (small deviations exist for $`N=19`$ and $`l=3,4`$, presumably due to corrections to scaling). Note that scaling of the PDF is independent of the values of $`\omega `$ and $`\rho `$.

Next, we compare the scaled PDF to known theoretical distributions. The Gumbel PDF for the $`m`$-th smallest in a set of i.i.d. random variables has the form gumbel ; foot\_gumbel

$$g_m(x)=w\mathrm{exp}\left[m\frac{xu}{v}m\mathrm{exp}\frac{xu}{v}\right]$$
(3)
where $`u`$ and $`v`$ are rescaling parameters and $`w`$ a normalization constant. To compare this to our data, we impose zero mean and unit variance which, together with normalization, fixes $`u`$, $`v`$, and $`w`$ as a function of $`m`$. Expressions for the cumulants of $`g_m(x)`$ are easily obtained. These are decreasing functions of $`m`$, in absolute value, and vanish for large $`m`$, thus $`g_m`$ provides an interpolation between the standard Gumbel ($`m=1`$) and a Gaussian PDF ($`m=\mathrm{}`$). As shown in Table I, the cumulants are fitted well for $`m=6\pm 0.3`$ ($`\stackrel{~}{k}_{5,N}`$ shows small deviations from the $`N39`$ data, possibly an effect of corrections to scaling). Table I also shows that $`m=\pi /2`$ (see below) is far from agreeing with the data, hence even more so is $`m=1`$, as also noted in Ref. krzakala, . This shows that the correlations in the energy levels give rise to a more Gaussian-like PDF with respect to the uncorrelated ($`m=1`$) case. The PDF $`g_6(x)`$ is displayed with solid lines in Figs. 3 and 4, and the corresponding cumulative distribution in the inset of Fig. 3: the agreement with the data is excellent everywhere. Although it is not clear why a Gumbel distribution with $`m1`$ should work here, $`g_6(x)`$ can certainly be used to accurately represent the true PDF up to several standard deviations. It would be interesting to reach higher statistics to check whether the agreement breaks down. Incidentally, the Gumbel distribution with a non-integer value $`m=\pi /2`$ has been conjectured to describe bramwell , at least to a good approximation, the PDF of spatially averaged quantities in many correlated systems, even when no extreme value is apparently involved. A modified Gumbel form with non-integer $`m`$ was also used to describe the order-parameter PDF in three-dimensional spin glasses berg .

Another example of extreme value of correlated random variables is the smallest eigenvalue, $`\lambda _M`$, of an $`M\times M`$ random matrix. For three well-known ensembles, the scaling PDF of $`\lambda _M`$ is the Tracy-Widom (TW) distribution tw , $`f_\beta (x)`$, with $`\beta `$ a parameter dependent on the ensemble. We tested it against our data finding a poor agreement, although better than for $`g_1(x)`$. Table I shows that the cumulants for $`\beta =1`$ (and thus for $`\beta =2`$), computed from the tabulated $`f_\beta (x)`$ spohn , are well outside the statistical errors. The PDF $`f_1(x)`$ is displayed with dotted lines in Figs. 3 and 4: one can see small deviations from the data for small $`x`$ and large deviations in the tails.

In conclusion, we presented accurate numerical data consistent with $`\omega =2/3`$ and $`\rho =3/4`$ when considering subleading corrections, although both exponents may be in fact slightly larger than these values. The PDF of the ground-state energy satisfies finite-size scaling with an asymptotic PDF empirically well described by a Gumbel distribution.

After this work was virtually finished, we learned that $`p_N(e_N)`$ was also studied numerically in Ref. andreanov, .

Acknowledgements – The author thanks A.P. Young and M.A. Moore for very useful exchanges, and the authors of Ref. andreanov,  for communicating some of their findings prior to publication. The author also thanks M. Mézard for a useful discussion and M. Prähofer and H. Spohn for making their computation of the TW distribution publicly available spohn . Support from STIPCO, EC grant number HPRN-CT-2002-00319, is acknowledged.


text:
# Measurement of phosphorus segregation in silicon at the atomic-scale using STM

## I Introduction

If electronic devices continue to shrink at the same rate as predicted by Moore’s law Moore65  they will within the next 10 - 15 years reach the atomic-scale. The fabrication of silicon-based devices at the atomic-scale requires the ability to control the position of dopant atoms, such as phosphorus with atomic precision. Dopant atoms are usually introduced into the silicon crystal either by ion implantation or by diffusion from highly doped surface layers Sze85 . However, these techniques do not allow for atomically precise control of the position of dopant atoms. An alternative way is to dope a 2D silicon surface and subsequently overgrow the sheet of dopant atoms with epitaxial silicon Oberbeck\_APL02 ; O'Brien\_PRB01 ; Gossmann85 . One concern with this method is that dopant atoms are known to diffuse at elevated temperature Fahey89  and to segregate during epitaxial silicon growth Nuetzel96 . It is therefore important to monitor the precise position of dopant atoms after high temperature processes such as crystal growth or annealing to ensure segregation and diffusion are kept to a minimum.

Secondary ion mass spectrometry (SIMS) is a commonly used measurement technique to investigate dopant segregation in Si Nuetzel96 ; Friess92 ; Kasper00 . However, it only has the depth resolution of a few nanometers depending on the sample and the sputtering conditions Wilson89 . In the case of phosphorus in Si, additional problems arise with SIMS measurements due to the high mass resolution required to distinguish between 31P (mass 30.9738 webelements ) and 30SiH (mass 30.9816 webelements ). Finally, ion beam mixing during the sputtering process can broaden the dopant profiles thereby influencing any quantitative analysis of dopant segregation.

In this paper, we compare three different methods to quantify the segregation of P atoms in silicon. We use SIMS to study segregation of P atoms in Si at the nanometer scale, Auger electron spectroscopy (AES) for studies at the sub-nanometer scale, and scanning tunnelling microscopy (STM) to measure the segregation and diffusion of phosphorus atoms to the Si surface at the atomic-scale. In each case a Si(001) surface is $`\delta `$-doped using PH3 gas as the dopant source followed by an anneal at 600C to incorporate the P atoms into the surface. In order to determine the amount of segregation that can occur during Si encapsulation the sheet of dopants is then overgrown with several monolayers (ML) of epitaxial silicon either at room temperature (RT) or 250C.

We demonstrate that reducing the silicon growth temperature to room temperature can minimize the segregation length of P in Si to sub-nanometer values, the lowest segregation length values for P in Si shown so far. Subsequent sample annealing, necessary to enhance the structural quality of the low temperature grown epitaxial Si layer, causes the encapsulated P atoms to diffuse to the surface. However, the density of P atoms that reach the surface can be reduced to a few percent of the original P atom density (even after annealing at 500C) if silicon encapsulation is performed at RT.

## II Experimental

Experiments were performed using an Omicron variable temperature STM with Auger electron spectroscopy and a Si sublimation cell. The base pressure of the two-chamber analysis/growth ultra-high vacuum (UHV) system was $`<`$ 5 $`\times `$ 10-11 mbar. Phosphorus doped 1 $`\mathrm{\Omega }`$cm Si(001) wafers were cleaved into 2 $`\times `$ 10 mm2 sized samples, mounted in sample holders, and then transferred into the UHV system. Sample preparation was performed in vacuum without prior ex-situ treatment by outgassing overnight at $``$550C using a resistive heater element followed by direct current heating at the same temperature for several hours. The sample was then flashed for 30 to 60 s to 1150C by passing a direct current through the sample. After flashing, the sample was cooled slowly ($``$3C/s) from 900C to room temperature. The clean sample surface at room temperature was then directly exposed for 15 min to a beam of phosphine (PH3) molecules at a chamber pressure of 10-9 mbar resulting in a surface saturation dosing with phosphine molecules of 0.37 ML Lin99 . An annealing step at 600C for 1 min incorporated P atoms into the Si surface and desorbed the H atoms resulting from dissociation of the PH3 molecules giving a total density of 0.25 ML of P atoms at the surface Lin99 .

For the fabrication of $`\delta `$-doped samples for SIMS measurements, the sheet of P atoms was then overgrown by an $``$24 nm thick epitaxial intrinsic Si ($`\rho `$ $`>`$ 1 k$`\mathrm{\Omega }`$cm) layer deposited from a Si sublimation cell at a substrate temperature of either 250C or RT, and a growth rate of 2 ML/min. Note that RT growth means without intentional sample heating. However, during deposition the sample was radiatively heated by the glowing Si arc used as a sublimation source for Si deposition. The sample temperature is estimated to increase to a few 10C above RT during Si deposition. For growth at 250C the Si sample was heated by direct current flow through the sample and the sample temperature was controlled by an infrared pyrometer. SIMS measurements were performed on the $`\delta `$-doped samples to determine the level of surface segregation of the P atoms during epitaxial silicon overgrowth. Using an ATOMIKA 6500 system with Cs+ primary ion beam the P depth profile was recorded at primary ion energies of 2, 5.5, and 15 keV. As the mass resolution of the ATOMIKA system was not sufficient to distinguish between the masses of 31P and 30SiH, additional SIMS data was acquired at 14.5 keV using a CAMECA SIMS system with a higher mass resolution. The SIMS results were then used to quantify the segregation of phosphorus atoms at the nanometer scale.

In addition, Auger electron spectroscopy and STM were used to study the segregation and diffusion of P atoms at the sub-nanometer scale and atomic-scale, respectively, during epitaxial Si growth and subsequent annealing steps. As before, a Si(001) surface was phosphorus $`\delta `$-doped and then overgrown by several monolayers of intrinsic Si at a growth rate of 2 ML/min and at substrate temperatures of either $``$250C or room temperature. After growth the sample was transferred to the analysis chamber where the Auger electron spectroscopy and STM measurements of the density of P atoms segregated to the Si surface were carried out. For Auger electron spectroscopy measurements the sample was then annealed in successive steps at temperatures between $``$300 and 900C in 150C steps for 1 min at each step. Auger electron spectroscopy measurements were then performed after cooling down of the sample to RT. For separate STM experiments sample annealing was carried out between 350 and 950C in $``$50C steps for 5 s at each step. After each process step the sample surface was investigated at room temperature and the number of P atoms at the surface counted.

## III Results and Discussion

### III.1 SIMS measurements to study segregation at the nanometer scale

Figure 1 (a) shows SIMS measurements of phosphorus $`\delta `$-doped Si(001) samples overgrown at 250C and RT, respectively. Using the ATOMIKA system, two distinct mass-31 peaks in the two curves at a depth of $``$10 and $``$24 nm were measured. Since it is difficult to distinguish between the 30SiH and 31P peaks additional measurements were carried out using a high mass resolution CAMECA system which displayed only one 31P peak at $``$ 20 nm (see inset of Fig. 1 (a)). This additional SIMS data confirms that the broad peak at $``$10 nm is not due to 31P, but is 30SiH arising from H adsorption at the Si surface. The peak at 24 nm can however, be ascribed to the P $`\delta `$-doped layer at the interface between the epitaxial layer and the substrate. The 31P peaks of the samples grown at 250C and RT, respectively, display a full width at half maximum (FWHM) of 6 and 4 nm, respectively, for a primary ion energy of 5.5 keV demonstrating the smaller amount of P segregation occurring during RT Si overgrowth. The segregation length Nuetzel96  $`\mathrm{}`$ for P segregation in Si was obtained by fitting the exponential decrease of the P concentration towards the sample surface displayed in Fig. 1 (a), giving values of $`\mathrm{}`$ = 2.3 nm and 1.5 nm for the samples grown at 250C and RT, respectively.

Figure 1 (b) highlights the limitations of the SIMS technique. Here we present SIMS measurements of the phosphorus $`\delta `$-doped Si(001) sample overgrown at 250C for various Cs+ primary ion energies of 2, 5.5, and 15 keV. Note that the apparent reduced P peak height of the 2 keV data results from a quantification of the data using the 5 keV signal as a reference. The P peak is observed to decay into the substrate with the decay length increasing from $``$1 decade/6.5 nm to $``$1 decade/20 nm as the primary ion energy is increased from 2 to 15 keV. This variation of the decay length with ion energy is a well-known measurement artefact arising from ion beam mixing, i.e. primary Cs+ ions pushing P atoms further into the substrate during the sputtering process. The SIMS data of the $`\delta `$-doped sample shows a FWHM of the 31P peak of 5, 6, and 14 nm for Cs+ primary ion energies of 2, 5.5, and 15 keV, respectively. While it can be seen that measurement artefacts due to ion beam mixing decrease with decreasing primary ion energy, it is also obvious that a small primary ion energy of 2 keV results in low count rates and therefore a stronger noise in the measurement signal and a higher detection limit. Most importantly the variation of the shape of the P peak with primary ion energy demonstrates that the resolution of the SIMS technique is only of the order of a few nm. As a result the segregation of P atoms during growth can only be detected at this level using SIMS highlighting that the technique is not sensitive enough to measure the segregation at the atomic-scale.

### III.2 AES measurements to study segregation at the sub-nanometer scale

As a consequence we also carried out studies of phosphorus segregation and diffusion at the sub-nanometer scale using Auger electron spectroscopy. Fig. 2 (a) shows Si and P signals at 107 eV and 120 eV Auger electron energies respectively, in the differentiated AES spectrum of the Si sample surface obtained from (i) the clean Si surface, the surface after (ii) PH3 saturation dosing and P incorporation, (iii) 5 ML Si overgrowth at RT, and (iv) various 1 min annealing steps at the temperatures displayed. Note: the intensity of the Auger signal decreases with sample depth and is therefore limited to approx. the first three ML in Si Seah79 .

In order to quantify the amount of P segregation/diffusion we have analyzed the AES data of the integrated spectrum. The number of Auger electrons calculated from the integrated spectrum is directly proportional to the P surface density Jacobson98 . Fig. 2 (b) shows the P density in the first three ML of the sample relative to the P density of the saturation dosed surface from two different growth experiments at 250C and RT. If we consider the figure in more detail we see that the clean Si surface after flashing does not show any discernable P concentration, as expected. After P saturation dosing and incorporation of the P atoms, the maximum P density normalized to 1 is visible on the surface giving the largest Auger electron peak at 120 eV. The P surface density then decreases to $``$0.74 and $``$0.12 after encapsulation with 5 ML of Si grown at 250C and RT, respectively, and increases again with sample annealing at various temperatures up to 600C due to P diffusion to the surface. After annealing at 750C the P density decreases markedly due to P desorption from the surface as P2 Jacobson98 . Finally after annealing to 900C the P peak disappears completely. Phosphorus is still detectable after the 750C annealing step as P diffusion from subsurface layers to the surface is still occurring during the anneal.

To calculate the segregation length $`\mathrm{}`$ of P in Si at 250C and RT from the AES data we use the relation Nuetzel96

$$p_{inc}=\text{a}_0/4\mathrm{}$$
(1)
with the incorporation probability $`p_{inc}`$ and $`\text{a}_0/4`$ = 0.1358 nm the distance between two subsequent Si(001) monolayers. The relation between the incorporation probability and the segregation length is only valid for first-order kinetic processes and assumes negligible desorption and bulk diffusion. To analyze our data we know that in Fig. 2 (b) we find $``$74$`\%`$ of P atoms at the surface after 5 ML Si growth at 250C. From this value we get an incorporation probability $`p_{inc}`$ = 0.06. According to equation 1 this value corresponds to a segregation length of $`\mathrm{}`$ = 2.3 nm. For RT growth we find a value of $`\mathrm{}`$ = 0.39 nm. In section D we discuss these values of segregation length obtained using AES measurements in relation to both the values in the literature and to segregation lengths obtained from SIMS and STM measurements. We note however that AES is sensitive to subsurface, non-segregated P atoms and therefore we expect these values to be slightly on the high side of the real $`\mathrm{}`$. In addition, in the case of RT Si growth due to the high surface roughness of the as-grown surface parts of the surface are covered by less than 5 ML of Si which makes it possible that AES detects non-segregated P atoms in their original lattice positions.

### III.3 STM measurements to study segregation at the atomic-scale

Finally in order to gain a better quantitative estimate of the P segregation/diffusion at the atomic-scale we used the STM directly. We imaged the clean Si(001) surface, the surface after PH3 dosing/P incorporation and after Si growth and subsequent annealing steps to determine the density of P atoms at the surface. Annealing of the sample after growth at 255C or RT is necessary to flatten the surface in order to identify the characteristic asymmetric appearance of the Si-P heterodimers which form after P incorporation from adsorbed PH3 molecules Oberbeck\_APL02 .

Figures 3 (a) - (f) show filled state STM images taken at room temperature of two separate growth experiments. For each experiment STM images of the clean Si(001) surface are shown after initial sample flashing (a); P incorporation (b); Si overgrowth at different temperatures (c); and successive annealing steps for 5 s at $``$350C, $``$500C, and $``$600C (d) - (f). Note that for space reasons only selected STM images of the annealing sequence are shown here. The complete annealing process is however summarized in Fig. 5 (a) and (b). The two experiments only differ in the Si overgrowth step (Fig. 3 (c)): either 5 ML of Si were grown at 255C or at RT.

Both experiments started with the preparation of a low defect density clean Si(001) surface, see Fig. 3 (a). The surface was then phosphine dosed and annealed at 600C for 1 min to incorporate the phosphorus atoms into the surface and desorb the associated hydrogen Pietsch95 . The resulting image in Fig. 3 (b) shows dimer rows on the Si surface that are strongly buckled due to the high phosphorus density at the surface. After 5 ML growth at 255C the surface shows small 2D islands and short Si dimer chains (c). Subsequent annealing of the surface at temperatures of 345 (d), 498 (e), and 600C (f) causes the Si surface to flatten due to island coarsening and diffusion of Si atoms to step edges. Simultaneously, the density of bright asymmetric features that result from segregated P atoms forming Si-P heterodimers at the surface increases in accordance with Auger electron spectroscopy data in Fig. 2 (b).

The difference between the two experiments is in the Si growth step (Fig. 3 (c)), where at RT 3D Si islands are formed due to the small mobility of Si atoms on the surface Herman96 . Successive annealing steps flatten the surface as with the 5 ML growth at 255C. Interestingly, the Si surface morphology in the different experiments looks very similar for the same annealing temperature, even for the first annealing step at $``$350C for 5 s (d). However, Figs. 3 (d) - (f) clearly show that the density of Si-P heterodimers at the Si surface is much lower for a given anneal temperature if the Si encapsulation occurs at RT compared with 255C Si growth.

Figure 3 (c) highlights an important limitation in the use of STM to investigate P segregation/diffusion: the high surface roughness that results from low temperature growth, especially room temperature growth, makes identification of Si-P heterodimers at the surface difficult. This occurs for two reasons: firstly, STM images of rough Si surfaces usually have a lower quality compared to images of an atomically flat surface. Secondly, to identify a Si-P heterodimer a sufficient brightness contrast between the heterodimer and the surrounding Si surface is necessary, which can only be obtained for atomically flat surfaces. Therefore, Si-P heterodimers only become clearly visible in our experiments after a short anneal at $``$350C for 5 s. As a consequence the P density at the Si surface after the first anneal is a result of not only P segregation that occurs during growth but also arises from diffusion of P atoms that occurs during the subsequent anneal. Nonetheless it is reasonable to assume that P diffusion, with an energy barrier of 3.66 eV Fahey89 , is negligible for short anneals at such low temperatures. Thus, the P density observed at the surface after the first anneal is most likely a result of P segregation that occurs during Si growth.

To further minimize the surface P density we here also encapsulated the sheet of P atoms with 10 ML, instead of only 5 ML, of Si grown at RT. Figures 4 (a) - (d) show a comparison between the two experiments. For both experiments the surface after Si growth consists of small 3D islands and flattens during subsequent annealing steps as before. As expected, the STM images show that the density of Si-P heterodimers at the Si surface is higher for 5 ML encapsulation compared to 10 ML Si growth.

To quantify the density of P atoms at the Si(001) surface observed in the STM images of Figs. 3 and 4 we have counted the number of Si-P heterodimers after each anneal. Figure 5 (a) shows the increase in the density of P atoms at the surface following subsequent annealing steps for the three different experiments. The relative density is obtained by comparison of the phosphorus density after Si growth and sample annealing, with the initial density of incorporated phosphorus atoms after phosphine dosing of the clean Si surface.

If the encapsulation of the sheet of P atoms is undertaken at 255C, already $``$25% of the initial number of phosphorus atoms have segregated during Si growth to occupy surface lattice sites after the first 5 s annealing step at the lowest applied annealing temperature ($``$350C), see Fig. 5 (a). After subsequent annealing at 400, 450, 500, 550, to 600C nearly 60% of the phosphorus atoms are present at the surface. These results demonstrate that even with a short, low temperature anneal encapsulation of P atoms in epitaxial Si grown at $``$250C results in significant P segregation.

In contrast, if P atoms are overgrown with 5 ML of Si deposited at RT and annealed at 350C for 5 s (Fig. 5 (a)), only $``$5% of the initial number of P atoms are present at the surface (a factor of 5 less). During subsequent annealing at 400, 450, 500, 550, and 600C the P density only increases slightly to $``$10% due to diffusion of phosphorus atoms from subsurface layers to the surface. The reduced density of P atoms at the Si surface compared to the 255C Si growth experiment is a direct consequence of the strongly suppressed segregation of P atoms during Si overgrowth at RT. Figure 5 (a) also shows that the surface P density can be further reduced by growth of a thicker Si layer of 10 ML at RT, as seen in the STM images of Figs. 4 (b) - (d). Note that the relative P concentration after various annealing steps obtained from Auger electron spectroscopy measurements shown in Fig. 2 (b) is higher than the values displayed in Fig. 5 (a). This is a direct result of the longer annealing time of 1 min for each step in the AES experiments and the fact that AES also detects subsurface P atoms, highlighting the limitations of this technique.

Figure 5 (b) shows the decrease of the surface phosphorus density for the case of 5 ML Si growth at 255C after successive anneals at temperatures between 600 and 950C due to phosphorus desorption from the surface. This figure also demonstrates that P atoms are still present at the surface even after annealing for 5 s at 900C due to the continual diffusion of P atoms from subsurface layers to the surface.

Using equation 1 we calculate the segregation length $`\mathrm{}`$ of P in Si from the STM data. We know that in Fig. 5 (a) we find about 25$`\%`$ of P atoms at the surface after 5 ML Si growth at 255C and the first annealing step at $``$350C. This means that 75$`\%`$ of the P atoms were incorporated. From this value (assuming that only segregation during growth occurred and that the diffusion of P atoms during the first annealing step at $``$350C is negligible) we get an incorporation probability $`p_{inc}`$ = 0.24. According to equation 1 this value corresponds to a segregation length of $`\mathrm{}`$ = 0.58 nm. For the 5 and 10 ML Si growth at room temperature, the same analysis gives an average segregation length of $`\mathrm{}`$ = 0.34 nm.

### III.4 Comparison of segregation length data

In this section we compare the segregation length values obtained from our SIMS, AES and STM data with values from the literature. Segregation length values for P in Si have been published by the Abstreiter group Friess92 ; Nuetzel96  and by Hobart *et al.* Hobart96  In these studies electrochemical capacitance-voltage profiling (eCV), SIMS Friess92 ; Nuetzel96 , and spreading resistance profiling Hobart96  are used to determine the segregation length at growth temperatures in the range of 320 to $``$900C and growth rates in the range of 0.06 Å/s to 2 Å/s. To compare the literature data with our experimental results we first re-scale the literature data to the growth rate of 0.045 Å/s = 2 ML/min used in our experiments using the relation Nuetzel96

$`\mathrm{}\sqrt{{\displaystyle \frac{1}{R}}},`$ (2)
where $`R`$ is the Si growth rate. The growth rate dependence of the segregation length arises from the fact that the growth rate determines the time available for dopant atoms to exchange sites from subsurface to surface lattice sites before they are encapsulated in Si. Since the literature data was obtained at higher Si growth temperatures we have extrapolated their data to RT, as shown in Fig. 6.

From this figure we can see that *three* observations can be made. Firstly, the experimental values for P segregation in Si determined by the Abstreiter group (Friess and Nützel) and Hobart *et al.* differ significantly. This difference can result for a number of different reasons, including a difference in the temperature calibration of the growth systems, different measurement techniques (eCV and SIMS vs spreading resistance profiling), and different P doping profiles and doping techniques (MBE growth of P doped layers vs P diffusion) used by the two groups. As a consequence the extrapolation of experimental data from Friess, Nützel, and Hobart *et al.* to lower temperatures also results in a large difference in the low temperature segregation lengths.

Secondly, there is a noticeable difference in the segregation length values obtained from our SIMS, AES, and STM experiments at the same temperature, see Table 1. The STM values are consistently lower than the respective values determined from SIMS and AES measurements. This can be attributed to limitations of these other techniques, most notably: (i) The depth resolution of the SIMS technique is limited to the nanometer range Wilson89  as a result of the sputtering process, where P atoms in the sample are measured from several monolayers simultaneously. This broadens the P signal leading to an increase in the segregation length obtained. (ii) AES not only detects segregated P atoms at the surface but also detects subsurface, non-segregated P atoms. The analysis of AES data therefore displays a segregation length which is higher than the real $`\mathrm{}`$. It is anticipated therefore that the segregation lengths obtained from STM measurements are more accurate than both SIMS/AES analysis.

Thirdly, a final observation that can be made is that our experimental data at these much lower temperatures lies at the limits of the data range expected from extrapolation of the literature values (see the grey area in Fig. 6). There could be two main reasons for this. (i) The sample structures measured are somewhat different: while we used P $`\delta `$-doped layers, Nützel Nuetzel96 , Friess Friess92  and Hobart Hobart96  used bulk doped layers. (ii) The strong temperature dependence of the segregation length can mean that a small difference in the temperature calibration of the different growth systems can give significant differences in the segregation length.

Despite the small differences in the segregation length between our data and extrapolated data from the literature, our study shows that low temperature Si encapsulation of P atoms with a few ML of Si, can suppress the segregation of P atoms very effectively. Our results also demonstrate that STM, in contrast to SIMS and AES, allows for a better quantitative determination of the segregation length at the atomic-scale.

Future experiments will focus on STM imaging of the buried phosphorus dopant atoms once they are encapsulated by epitaxial Si overgrowth. Buried dopant imaging has already been demonstrated for boron and arsenic dopants in Si using hydrogen terminated Si(001) surfaces Liu01 . Such a technique would be very important to detect buried P atoms in the Si matrix, not only to allow atomically precise determination of the lateral position, but also the vertical position of the P atoms.

## IV Conclusions

We have used SIMS and, for the first time, AES and STM to analyze P segregation/diffusion in Si at very low temperatures (250C and room temperature). We show that SIMS is unable to detect segregation lengths of phosphorus atoms in silicon below the nanometer scale. In contrast, AES and STM can be used to monitor the segregation and diffusion of phosphorus atoms in silicon at the sub-nanometer and atomic-scale, respectively. We highlight that Auger electron spectroscopy is limited as a quantitative technique since it detects subsurface, non-segregated P atoms. As a result we demonstrate that STM measurements can be used directly to quantify the surface density of segregated P atoms at the atomic-scale. Using STM we obtained segregation lengths below 1 nm for phosphorus encapsulated in silicon grown at 250C and room temperature respectively, which are the lowest values recorded in the literature. Differences between the segregation length obtained and extrapolated experimental data from the literature are accounted for by differences in sample preparation and experimental techniques. These results highlight the effectiveness of room temperature silicon encapsulation to minimize dopant segregation in silicon and how it is possible to quantify this directly using STM.

###### Acknowledgements.

The authors thank P. Priestley for analysis of STM data. LO acknowledges a Hewlett-Packard Fellowship. MYS acknowledges a Federation Fellowship. This work was supported by the Australian Research Council and the Army Research Office (ARO) under contract number DAAD19-01-1-0653.


text:
# The quantum dynamics of atomic magnets, co-tunneling and dipolar-biased tunneling

## Abstract

Multi-spins tunneling cross-relaxations in an ensemble of weakly-coupled Ho3+ ions, mediated by weak anisotropic dipolar interactions, can be evidenced by ac-susceptibility measurements in a high temperature regime. Based on a four-body representation, including the rare-earth nuclear spin, two-ions tunneling mechanisms can be attributed to both dipolar-biased tunneling and co-tunneling processes. The co-reversal involving entangled pairs of magnetic moments is discussed with a particular emphasis, giving new evidences to elucidate the many-body quantum dynamics.

Since the discovery of quantum tunneling of the magnetization in single molecule magnets , the effects of environmental degrees of freedom on tunneling of mesoscopic spins, especially due to other spins , have been investigated in detail within the single spin framework . Indeed, most features of the staircase-like hysteresis loops, observed in an *ensemble* of weakly-interacting mesoscopic spins, were rather well understood within this simple picture of a single spin submitted to external and internal magnetic fields. In the presence of an axial anisotropy, hysteresis occurs at low temperatures due to the energy barrier $`E_B`$ hindering the magnetization reversal, when thermal activation over the barrier requires a very long relaxation time $`\tau =\tau _0e^{E_B/k_BT}`$. Nevertheless, non-stationnary quantum mechanics allows for a faster time-evolution at avoided level crossings within the single-spin Zeeman diagram (tunneling), leading to magnetization steps for well-defined values of the magnetic field applied along the easy axis (resonance conditions). Each mesoscopic spin experiences a time-dependent distribution of internal fields which is determinant for the analysis of the quantum relaxation of weakly-coupled spins . These fields being generated by environmental spins, coupled by dipolar and eventually exchange interactions, one should ask the question to what extent the measured tunneling rate is related to quasi-isolated spins or to larger entities such as two (or more) spin centres.

Recent magnetization measurements of mesoscopic Ho3+ ions, highly diluted in a LiYF4 single crystal, evidenced the coherent rotation of both the electronic angular momentum and the nuclear spin of a nearly-isolated Ho3+ ion . However, some additional small magnetization steps were also observed at fast field-sweep rates, and still remain ambiguous because these steps have no explanation in term of level crossings within the single-ion model, including the rare-earth nuclear spin. In this letter, we demonstrate that new tunneling mechanisms originate from an *unused* feature of this simple picture: equally-spaced energy levels allow for non-dissipative spin–spin cross-relaxation tunneling processes. Ac-susceptibility measurements in a high temperature regime give us a unique mean to assess the *various* multi-spins tunneling mechanisms which may affect the quantum dynamics of an assembly of mesoscopic spins. In particular, we focus on the tunneling of pairs of Ho3+ ions, weakly coupled by anisotropic dipolar interactions, leading to new resonance conditions that are an intrinsic property of the two-ions picture (no matter is the coupling strength). We call this phenomenon a *co-tunneling* process, and show that pairs of Ho3+ ions may co-tunnel with interactions as weak as a few mK. This suggests that, under these specific field conditions, any environmental Ho3+ ion may co-tunnel with a ‘central’ spin, making the quantum spin-dynamics much more collective than expected before.

At low temperatures, a nearly isolated Ho3+ ion in LiYF4 behaves as an *atomic magnet* due to a large axial crystal-field anisotropy along the tetragonal easy $`c`$-axis , thus leading to very long spin–phonon relaxation times. Therefore, the magnetic moment is oriented along the $`c`$-axis and constitutes a two-level system (electronic doublet of Ising-type), with the moment being frozen either up or down . However, the latter may still evolve due to quantum tunneling at avoided electro-nuclear level crossings of a single Ho3+ ion, the rare-earth nuclear spin being $`I=7/2`$. The single-ion tunneling process thus results from the coherent rotation of both its electronic momentum and nuclear spin, the crossing fields being only related to the hyperfine constant $`A_J`$ . In this *strong* hyperfine coupling case, decoherence only occurs due to nuclear spins of diamagnetic ions (that is, mainly the nuclear spins of eight fluorine ions, with $`I=1/2`$). Furthermore, the effect of a many-body quantum dynamics was also observed in holmium-based spin glasses , but the way how quantum fluctuations arise in this more concentrated regime still remains unclear.

Ac-susceptibility of a 0.1$`\%`$ at. holmium-doped LiYF4 single crystal was studied as a function of a quasi-static field applied along the easy $`c`$-axis. The measurements were performed in a high temperature regime ($`T1.75`$ K) with a Quantum Design MPMS SQUID magnetometer, using a 4 Oe-amplitude excitation field. In Fig. 1, a low-frequency measurement, as obtained with $`f_m=163`$ Hz, reveals peaks (dips) in the in-phase (dissipative) response when the dc-field is changed, allowing us to investigate the role of the various quantum relaxation mechanisms. In every cases, the results can be qualitatively understood when considering a field-dependent relaxation time to describe the system dynamics. As shown in Fig. 2, if the constant frequency measurement $`f_m`$ is smaller than the characteristic frequency $`f_0`$ of spin fluctuations (case - *iii*-), an increase in the in-phase response, or a decrease in the dissipative response, is associated with a faster relaxation rate (that is, $`f_0`$ increases). The large deviations to the monotonous response, as observed in Fig. 1, are clear signatures of the *single-ion quantum dynamics* in a thermally-activated tunneling process through the energy barrier ($`E_B10`$ K, due to the first excited electronic singlet). Indeed, the faster relaxation of the magnetization observed for all the resonant fields $`H_n`$ (see solid lines in Fig. 1), with $`H_n=n\times 230`$ Oe ($`7n+7`$, $`n`$ being an integer), is in agreement with the crossing-fields of the electro-nuclear Zeeman diagram calculated within the single-ion picture (shown in Fig. 4a). Moreover, their width is temperature-dependent. Besides, similarly to non-expected magnetization steps observed in hysteresis loops at very low temperatures, additional peaks/dips are also evidenced for integer values of $`|n|`$ larger than 7, as well as for half-integer $`n`$ values with $`132n+13`$, which do not correspond to any level crossing in Fig. 4a). These strongly suggest that *two* new tunneling mechanisms occur. Because of the nature of the intra-atomic hyperfine coupling, all these fields involve a unique property within the single-ion representation: some energy levels are equally separated, so that simultaneous multi-spins transitions may conserve the total energy. Furthermore, part of these cross-spins relaxation (CSR) processes may *not* conserve the total magnetization of the system, some electro-nuclear states being associated with opposite directions of the magnetic moment, thus allowing for a *multi-ions quantum dynamics*. Among these transitions, we can clearly distinguish between the *co-tunneling* and the *biased tunneling* processes. In particular, we discuss the role of co-tunneling transitions for the first time in mesoscopic magnetism. As shown below, ac-susceptibily anomalies for half-integer $`n`$ are only related to co-tunneling events involving a two-ions process, which only occur within the dipolar window, their width thus being temperature-independent. In addition, this width is smaller than the one related to thermally-activated single-ion processes, because at $`T1.75`$ K the homogeneous broadening due to spin–phonon couplings is larger than the inhomogeneous dipolar broadening. Note that thermal activation over the energy barrier still requires longer times than the CSR *tunneling* events; the spin–phonon relaxation time can be mesured for an applied field located out of any resonance . In this strongly diluted limit, $`N`$-ions co-tunneling transitions with $`N>2`$, allowed at fields $`H_p=p\times H_{(n=1)}/N`$ (integer $`p`$, with $`|p|7N`$), are less probable. If not, $`N=3`$ transitions should be observed, the first resonance condition $`H_{(n=1)}`$/3 being out of the dipolar width of the zero-field resonance.

Finally, to clearly distinguish between every tunneling mechanisms, we used a higher measurement frequency $`f_m`$ to evidence the role of dipolar-biased tunneling processes. Indeed, these lead to additional anomalies in the dissipative response. At $`f_m=800`$ Hz and $`T=1.75`$ K, the relaxation rate associated with a thermally-activated single-ion tunneling process is smaller than $`f_m`$ (case - *i*\- in Fig. 2), whereas it is larger for the dipolar-biased tunneling process (case - *iii*-). Both rates increase being closer to resonance conditions related to integer $`n`$ values with $`|n|7`$. However, the dissipative response now increases close to a resonance in the former case, whereas it decreases in the latter one. As the phonon-induced broadening is larger than the dipolar-induced one in this temperature range, the response first *increases* when the applied field gets closer to the resonance, but still being out of the dipolar window so that CSR do not occur. When this field enters the dipolar window centred around the resonance, the relaxation is quickly dominated by the CSR rate and the dissipative response now *decreases*. This feature allows us to clearly evidence, for the first time, the dipolar-biased tunneling mechanism as the occurence of a dip dug in the dissipative response, as shown in Fig. 3b). Again, the width of this dip corresponds to the narrow dipolar window, but it is now temperature-independent only below $`T2.2`$ K. When the temperature is increased further, the thermaly-activated tunneling rate indeed gets faster than the measurement frequency. The large maximum envelope of the response then turns into a minimum and we cannot distinguish between these two processes anymore.

To underline the role of spin–spin interactions and discuss the physics of CSR in LiYF4:Ho3+ in detail, we go beyond the single-ion picture and introduce a two-ions representation energy scheme. As shown in Fig. 4b), within this four-body picture (two electronic angular momenta and two nuclear spins) the former situation of equally-spaced levels now corresponds to the existence of energy level crossings of pair states, associated either with a co-tunneling event (that is, two spins, initially in the same state, flip simultaneously to the opposite state), a biased tunneling process (that is, only one spin flips) or a flip-flop process (that is, two spins with an opposite initial state flip simultaneously). Only the first two mechanisms contribute to the relaxation of the magnetization, the latter being associated with the so-called spin-diffusion phenomenon. The two-ions Zeeman diagram was calculated on the basis of the single-ion electro-nuclear states shown in Fig. 4a). Details of the calculation will be published elsewhere. Without taking any interaction between Ho3+ ions into account, we underline in Fig. 4b) the intrinsic new features of this four-body picture: some additional crossing-fields are only related to the energy-level separations within the single-ion picture (see dotted lines). The effect of an additional *isotropic* coupling is shown in the inset of Fig. 4b), using a large amplitude for the sake of clarity. Part of the degeneracies is removed so that, in an ensemble of weakly-coupled mesoscopic spins, some crossings are shifted by a distribution of equivalent molecular fields $`H_{bias}`$, leading to biased tunneling. Co-tunneling crossing-fields are not affected. Note that only *anisotropic interactions* allow for CSR tunneling transitions, leading to finite tunneling gaps, such as the non-secular part of the dipole–dipole interaction. All the features observed in Fig. 1 now have an interpretation in terms of level crossings involving a pair of Ho3+ ions, that is, within the four-body picture. The same ideas can be extended to weakly-coupled single molecule magnets (SMMs), such as the well-known molecular nanomagnets Mn12-ac and Fe8. Indeed, co-tunneling could explain specific heat anomalies observed at intermediate temperatures in Fe8 by Gaudin and co-workers . Moreover, isotropic exchanged-biased tunneling was also recently observed in a molecular dimer \[Mn4\]2 but in a strong antiferromagnetic coupling regime , thus preventing co-tunneling events, as well as in weakly-interacting Mn4 SMMs , still with a larger ferromagnetic exchange coupling than dipolar interactions.

In LiYF4:Ho3+, two-ions energy levels are strongly degenerate due to the high degree of symmetry of the single-ion Zeeman diagram, the strong dipolar hyperfine coupling resulting in many equally-separated energy levels. This unique property makes diluted Ho3+ ions in LiYF4 a *model system* to evidence, and describe in detail, all the various CSR tunneling mechanisms that can occur in an ensemble of weakly-coupled mesoscopic spins. Generally speaking, resonance conditions associated with the multi-spins tunneling dynamics can be easily deduced from energy-conserved multi-transitions events within the single-spin picture. In particular, co-tunneling processes are a pure *intrinsic property of the many-body picture*, even for vanishing interactions, the crossing fields being independent of the coupling strength. More importantly, co-tunneling events are only allowed by anisotropic interactions, which are also expected to be a key ingredient to understand quantum fluctuations at larger scales in complexe systems, such as in the quantum spin glass problem.

LIST OF FIGURES CAPTIONS

* (a,b) Ac-susceptibility measured at $`T=1.75`$ K with a frequency $`f_m=163`$ Hz. The non-monotonous ac-response, as a function of $`H_z`$, gives evidence for a faster relaxation at well-defined values $`H_n=n\times 230`$ Oe. Large and broad peaks or dips are clearly observed for integer $`n`$ values, with $`7n+7`$, whereas smaller and narrower ones are measured for half-integer $`n`$ values, with $`132n+13`$, as well as for $`n=\pm 8`$ and $`n=\pm 9`$. The solid lines correspond to electro-nuclear level crossings calculated within the single-ion representation (see Fig. 4a), whereas the dotted lines are associated with level crossings calculated within the two-ions representation (see Fig. 4b).
* Variation of the in-phase (a) and dissipative (b) responses in the ac-susceptibility for three different relaxation rates $`f_0`$, described within the Debye model assuming a single (dominant) relaxation time. The constant measurement frequency $`f_m`$ is suggested by the vertical solid line, and $`\chi _0`$ is the static susceptibility. Depending on the value of $`f_m`$, the dissipative reponse can be either increased (case - *i*-) or decreased (case - *iii*-) if the dynamics gets faster. In any case, the in-phase response follows the variations of $`f_0`$.
* (a,b) Ac-susceptibility measured at a larger frequency $`f_m=800`$ Hz. Dipolar-biased tunneling is clearly evidenced, at low-enough temperatures, due to a decrease in the dissipative response for every integer $`n`$ values. The large arrows in (b) show the dips observed for $`n=`$ 0 and $`n=`$ 1. The width of these dips is constant below $`T2.2`$ K and corresponds to the average amplitude of dipolar interactions.
* (a) Electro-nuclear energy levels calculated within the single-ion representation. Equally-spaced levels in the Zeeman diagram allow for non-dissipative CSR processes in the tunneling regime. As suggested by red arrows, *three* various kinds of CSR may occur at well-defined values $`H_n=n\times 230`$ Oe. (b) Within the two-ions Zeeman diagram, each cross-relaxation process involving a pair of magnetic moments now corresponds to a level crossing. Note that the finite slopes are twice larger than in (a). Their crossing describes a co-tunneling process (integer or half-integer $`n`$, with $`142n+14`$). Other crossings only involve the flipping of one magnetic moment (integer $`n`$, with $`14n+14`$). The inset shows how an isotropic interaction, with a given coupling stength, removes part of the degeneracies for integer $`n`$ values, leading to biased tunneling, as well as an avoided level crossing due to both Ho3+ nuclear spin and first electronic singlet.