text:
# 1 Introduction

## 1 Introduction

The world’s ecosystems are increasingly being subjected to stresses that result in large-scale changes in species population densities. These stresses often directly or indirectly arise from human activities and include pollution, over-exploitation, species invasions and habitat destruction (Carlton and Geller 1993, Milner-Gulland and Bennett 2003). Understanding how ecosystems respond to such perturbations is therefore highly important.

Here we will attempt to contribute to this understanding by focusing on species deletion, the complete removal of a species from an ecosystem community, using the theoretical framework of dynamical modelling. Species deletion is a large-scale perturbation of particular relevance, as it is a commonly used empirical tool to measure interactions strengths within real communities, and can be considered a reasonable approximation to other large perturbations (Paine 1980, Pimm 1980).

The theory of small perturbations in dynamical models of ecosystems is well developed. It began with May’s seminal work showing that the probability of an ecosystem with random interactions being locally stable decreases with both the number of species, the frequency of interactions and the strength of those interactions (May 1972, 1973). This result was an important contribution to the complexity-stability debate and contradicted earlier ideas that complexity should naturally lead to stability (Odum 1953, MacArthur 1955, Elton 1958). A crucial conceptual element to May’s work is that only a local knowledge of the dynamics, encapsulated in the “community matrix”, is necessary to determine the stability of the population equilibrium to small perturbations. This is also true of perturbations that actually alter the position of the equilibrium, provided that they are small enough (Yodzis 1989). The results of small perturbations can be determined using only “local models”, where the species growth rates are approximated by linear functions of the population densities (Yodzis 2001).

In contrast the study of large-scale perturbations requires a “global model” — one defined over the whole of phase space. The use of such a model will inevitably involve a modelling choice, but it is important that it incorporates phenomena such as non-linear functional responses and adaptive foraging, which will be likely to operate over these large changes in population density (Abrams 1996).

A global model with these features has already been developed as part of a larger model of community coevolution by some of us (Drossel et al. 2001). This model assembles ecosystem communities through the repeated addition of new species that are modified versions of those already present. It is therefore conceptually similar to community assembly models (Drake 1990, Law and Morton 1996, Morton and Law 1996, Lockwood et al. 1997), the difference being that new species are generated in situ rather than being taken from a species pool. In the model species are constantly being subjected to large perturbations in population densities as new species add and existing species go extinct. Crucially, species are allowed to alter their foraging strategies in response to these changes. Thus the population dynamics of the model are particularly well suited to the study of species deletion and will be used in this study. Since previous studies of deletion have used the Lotka-Volterra equations or the equivalent discrete time Ricker dynamics, this will give a unique perspective on the problem (Pimm 1979, Pimm 1980, Borrvall et al. 2000, Lundberg et al. 2000).

We will not only use the population dynamics of the model, we will also use it to generate the food webs from which species will be deleted. It might be preferable to use real food web structures, but then interaction strengths that possess a stable equilibrium for that structure would have to be determined, which is likely to be difficult for a large food web. The only studies of deletion that have used large empirical food webs circumvent this problem by adopting a network approach, without using an explicit population dynamics (Solé and Montoya 2001, Dunne et al. 2002). Their procedure reveals interesting features of the network topology, but by using the simple rule that species go extinct when none of their prey remain, they ignore the complex realities of the dynamics. Our model allows us to generate a data set of food webs with reasonably realistic topologies and interaction strengths (Drossel et al. 2001, Quince et al. unpublished), suitable for investigating the complex dynamics of deletion in multi-species communities.

The outline of the rest of the paper is as follows. We begin with a review of the model in Section 2, followed by the details of the generation of the food web data set and the deletion experiments themselves in Section 3. The next sections concern the results of these experiments. In Section 4 we investigate how the trophic relationship between a pair of species influences the outcome of deleting one of them, in Section 5 we go on to consider which species properties correlate with large changes in community composition following that species removal and in Section 6 we consider the stability of the communities as a whole to deletion and the extent to which this is determined by their food web structure. We finish with a brief discussion of the major results in Section 7.

## 2 The model

We will now give a short description of the model we will use to evolve the food webs to be used in our study. Further details are given in Caldarelli et al. (1998), Drossel et al. (2001), Quince et al. (2002) and Quince et al. (unpublished).

The dynamics of the model has a different appearance depending on the time-scales under consideration. On time-scales of the order of the lifetimes of individuals, the number of species is fixed and the dynamics is deterministic: it is given by a variant of the standard equations of population dynamics. This includes the capacity of all species to change their predator-prey characteristics with time, competition between predators and a varying amount of effort that a predator puts into catching a particular prey. The system is allowed to change according to this dynamics for what may be relatively long intervals until equilibrium is reached, that is, until the populations of the different species present remains unchanged. If during the population dynamics the density of a species falls below a value $`N^{\mathrm{min}}`$, usually taken to be 1, it is assumed to have become extinct, and is removed from the system.

Once an equilibrium of the populations densities is reached a speciation event is initiated: a new species is generated by changing one of the features of one of the individuals of a randomly chosen species. This new species is then added to the system with a population density $`N^{\mathrm{child}}`$, also taken to be 1 in these simulations. The system is then again allowed to develop under the deterministic equations of the population dynamics. The small population of the new species may give rise to a viable population, or it may die out, but eventually when a new equilibrium is reached, a new speciation event will take place. By repeating this procedure tens of thousands of times an entire food web can be evolved, using a combination of conventional population dynamics and stochastic speciation events. On these very long evolutionary time-scales the discrete time steps where speciation occurs are the main aspect of the dynamics.

In order to define the speciation process, and more generally to be able to characterise a species, we need to introduce features which when taken together make up a species. In our model, features are specified by integers: $`\alpha =1,\mathrm{},K`$. Any subset of $`L`$ of these features constitutes a species. It is assumed that the effectiveness of predator-prey relationships among species is due to the effectiveness of individual features against each other. Therefore the score of one species, $`i`$, against another, $`j`$, denoted by $`S_{ij}`$, is defined in terms of the $`K\times K`$ matrix $`m_{\alpha \beta }`$ which gives the score of feature $`\alpha `$ against feature $`\beta `$:

$$S_{ij}=\mathrm{max}\{0,\frac{1}{L}\underset{\alpha i}{}\underset{\beta j}{}m_{\alpha \beta }\}.$$
(1)
The matrix $`m_{\alpha \beta }`$ is antisymmetric. Its independent elements are random Gaussian variables with zero mean and unit variance chosen at the beginning of a simulation run and not changed during that particular run. This allows the score of one species against another to be calculated from (1): if $`S_{ij}>0`$ then species $`i`$ is adapted for predation against species $`j`$, if $`S_{ij}=0`$ then it is not. We will also need to define the overlap $`q_{ij}`$, between two species $`i`$ and $`j`$, as the fraction of features of species $`i`$ that are also possessed by species $`j`$. The external environment is represented by a species indexed 0. This is assigned a random set of $`L`$ features at the beginning of a run, and is not changed throughout the course of the run.

Having described the structure of a species in the model, and used this to define the score (1) and overlap, we will now use these quantities in the construction of the population dynamics that governs the changes in population sizes between speciation events. The rate of change of $`N_i(t)`$, the population size of species $`i`$ at time $`t`$, is given by

$$\frac{dN_i}{dt}=N_i+\lambda \underset{j}{}N_ig_{ij}\underset{j}{}N_jg_{ji}.$$
(2)
The function $`g_{ij}`$ is the functional response: the rate at which one individual of species $`i`$ consumes individuals of species $`j`$. The choice of $`g_{ij}`$ essentially defines the nature of the population dynamics. We will give the explicit form chosen below, but it is here that the dependence on the score and the overlap functions will enter. It will also depend on the $`N_k(t)`$, and so will change with time.

The terms on the right-hand side of (2) are now simply interpreted. The last factor represents the loss in resources for species $`i`$ due to predation by all of the other species (in the model the measure of resources and species number are synonymous). The factor $`_jN_ig_{ij}`$ on the other hand represents the gain to species $`i`$ from predation on the set of species $`j`$ including the environment, species 0. The environment is assigned a fixed population, $`N_0=R/\lambda `$, thus $`R`$ is a parameter of the model that controls the rate of input of external resources. If it is assumed that a fraction $`\lambda `$ of the resources gained through predation are used to create new members of species $`i`$, the second term on the right-hand side of (2) is obtained. Finally, the first term simply represents the rate of death of individuals in the absence of interaction with other species.

In order to briefly motivate the form of the functional response we will use, let us first discuss the case of a single predator $`i`$ feeding on a single prey $`j`$. In this case

$$g_{ij}(t)=\frac{S_{ij}N_j(t)}{bN_j(t)+S_{ij}N_i(t)},$$
(3)
where $`b`$ is a constant. We can gain more understanding of the structure of $`g_{ij}`$ by noting that when the predators are far more numerous than the prey ($`N_iN_j`$), $`N_jg_{ij}N_i`$: the feeding rate of the predators is limited only by the number of prey. In the other limit, when the prey is very abundant compared with the predators ($`N_jN_i`$), $`g_{ij}S_{ij}/b`$: each predator feeds at a constant maximum rate. This latter result also gives an interpretation to the constant $`b`$. Having introduced the basic form (3), we can now state general form for the functional response used in the model:

$$g_{ij}(t)=\frac{S_{ij}f_{ij}(t)N_j(t)}{bN_j(t)+_k\alpha _{ki}S_{kj}f_{kj}(t)N_k(t)}.$$
(4)
There are two new aspects present in (4) and absent in (3):

* Interference competition between predators of prey species $`j`$ is modelled by the factor $`\alpha _{ki}`$. We take $`\alpha _{ii}=1`$ and $`\alpha _{ki}<1,ik`$ to reflect the fact that competition between members of the same species is typically stronger than competition between different species. In fact, we expect that the more species are alike, the greater will be the competition between them, and therefore take

$$\alpha _{ij}=c+(1c)q_{ij},$$
(5)
where $`c`$ is a constant lying between zero and one which is the residual degree of competition that exists even if two competing predators have no features in common.
* Adaptive foraging is modelled using the factors $`f_{ij}`$. The effort $`f_{ij}`$ can be viewed as the fraction of time an individual of species $`i`$ spends predating on species $`j`$ or the fraction of the population of species $`i`$ dedicated to consuming only $`j`$. These efforts must satisfy $`_jf_{ij}=1`$ for all $`i`$. To determine the $`f_{ij}`$ it seems reasonable to assume that the gain which an individual of species $`i`$ makes in consuming individuals of species $`j`$ (that is, $`g_{ij}`$), divided by the amount of effort $`i`$ puts into this task (that is, $`f_{ij}`$), should be the same for all prey species $`j`$. Using this condition, together with the normalisation of the $`f_{ij}`$, leads to

$$f_{ij}(t)=\frac{g_{ij}(t)}{_kg_{ik}(t)}.$$
(6)
This choice of efforts can be shown to be an evolutionary stable strategy (Drossel et al. 2000), or in the terminology of foraging theory, an ideal free distribution of predators across prey (Fretwell and Lucas 1970).

The calculation of the efforts, through (6), effectively introduces a new behavioural timescale into the problem. We assume that the efforts change on a much shorter timescale than the population densities $`N_j(t)`$, and therefore that they may be found by iterating (4) and (6) assuming constant population densities. When this process has been completed, we may then move on to updating the population densities. Thus we do not treat the efforts as dynamical variables, instead we assume that they are a function of the population densities, even if we have no explicit form for this function. This describes the essential features of the model. In the rest of the paper we will explore the consequences of performing deletion experiments on food webs created by the model.

## 3 Deletion experiments

The set of model food webs used in this study was obtained by performing four hundred independent simulations of the model, each simulation lasting for 120000 speciations. These simulations were independent, in that different pseudo-random number sequences were used in their generation. Thus the simulations differed in their random matrices, $`m_{\alpha \beta }`$, environment features and speciation events. The same parameters were used in all the simulations these being $`R=1\times 10^5`$, $`b=5\times 10^2`$, $`c=0.5`$, $`\lambda =0.1`$ and $`N^{\mathrm{min}}=N^{\mathrm{child}}=1.0`$. The effect of altering the model parameters is investigated in Quince et al. (unpublished). The four hundred final food webs from these simulations, which after 120000 speciations will have structures drawn from a stationary distribution, constituted the ecosystem data set.

A typical food web is shown in Fig. 1. Each species in this diagram is represented as a circle, the sizes of which are the same for all species. This differs to the convention adopted in Quince et al. (2002), where the radii of the circles were proportional to the logarithm of the population densities. The arrows represent predator-prey interactions, with the arrow pointing from the prey to the predator. The intensity of the arrow is proportional to the fraction of the predators diet that consists of that particular prey. The vertical arrows originating from the base of the diagram, rather than from another species, indicate that the species is feeding off the environment. The species are positioned vertically according to trophic height, defined as the average path length from the species to environment, the average being weighted by predator diet fractions. The dashed lines show the position of integer values of the trophic height.

The four hundred food webs in the model data set span a range of structures. We shall quantify this variation with five food web statistics:

1. The total number of species in the food web denoted by $`S`$.
2. The links per species denoted by $`L/S`$. This quantity is simply the number of predator-prey interactions divided by the number of species, where we will use the convention of counting a link if it constitutes greater than 1% of a predator’s diet.
3. The fraction of omnivorous species in the food web denoted by $`O`$. We define omnivorous species as those which feed at more than one trophic level and define the trophic level of a species to be the shortest path from that species to the environment. The justification for this choice is the observation (Yodzis 1984) that the shortest path between a species and the environment tends to be the most important energetically.
4. The degree of cycling denoted by $`𝒞`$. The method we adopt for measuring the amount of cycling of energy in the food webs is based on both the ideas presented in Ulanowicz (1983) and on source code kindly provided by the author. Essentially a backtracking algorithm was first used to identify each cycle in a food web. Having done this, the amount of energy flowing in a cycle was identified with the strength of the weakest link of the cycle, exactly as in Ulanowicz (1983). To define $`𝒞`$ we measure the proportion of a predator’s prey obtained through cyclic flow averaged over the species in the web.
5. The ecosystem redundancy denoted by $`𝓇`$. The redundancy of an ecosystem is the proportion of species which can be considered superfluous to the functioning of the ecosystem as a whole (Walker 1992). We define a species to be redundant in our food webs if at least one other species possesses the same pattern of trophic links i.e. the same predators and prey. As for the calculation of $`L/S`$, only links forming greater than 1% of the predator’s diet are used in this calculation. Then the redundancy, $`𝓇`$, is the fraction of redundant species in the food web.

In Table 1 the means and standard deviations of the five statistics are shown for the four hundred food webs in the data set.

This data set was then used to investigate the effect of species deletion. A single species was removed from a web and the population dynamics iterated until a stable equilibrium was reached. If, during this process, the population of a species fell below $`N^{\mathrm{min}}`$, then it was removed (“went extinct”) in accordance with the criterion applied when evolving the communities. For each web every species was deleted independently, that is, the webs were returned to their original state between deletions. In total, over all the food webs, 25531 species were deleted. In Fig. 2 we show the frequency distribution of number of further extinctions for these deletions. This distribution does not decay exactly exponentially with further extinction number, but it does have a characteristic size of just a few species. The number of further extinctions seems bounded, the largest number is seventeen, and on no occasion is the whole web of typically sixty species close to collapse. The effect of deleting a species is localised in the webs.

## 4 Trophic relationships and species deletion

In this section we investigate how the trophic relationship between a pair of species influences the impact that deleting one of the pair will have on the other. We began by simply taking the four hundred webs described in the previous section and placing each species into four non-exclusive categories according to its trophic relationship with the deleted species. These categories were predators, prey, competitors and indirect predators. The definitions of the first two categories are obvious. A “competitor” was defined as any species that shared a prey with the deleted species. The final category of “indirect predators” requires a little more explanation. Consider any two species in an ecosystem, $`i`$ and $`j`$. If there exists one or more paths from $`j`$ to $`i`$ travelling only up trophic links, from prey to predator, and the minimum length of those paths is greater than one, then we define species $`i`$ to be an indirect predator of species $`j`$. In other words, species $`i`$ is an indirect predator of $`j`$ if some of the resources consumed by $`i`$ came originally from $`j`$, but species $`i`$ does not prey on $`j`$ directly. In calculating the categories all links which formed less than one percent of the predators diet were ignored.

Having defined these categories we calculated two quantities. The first was the mean number of species in each category that went extinct following species deletion, with the average performed over every deleted species in every web. The second was the probability that a species with a particular trophic relationship to the deleted species went extinct. This quantity was estimated by simply averaging the proportion of species with a particular relationship that went extinct. These results are shown in Table 2. The first row shows the results for all species, in order to aid comparison with the other results. We see that, on average, following the deletion of a species from a web, a further 1.366 species went extinct and that any given species had a probability of 0.021 of going extinct. A clearer picture is usually obtained by using probabilities, rather than mean values, when making comparisons between tropic relationships. This is due to the weaker dependence that probabilities have on the number of species with that particular relationship to the deleted species. From these results we see that that a predator of a deleted species is more than ten times as likely to go extinct following the deletion of its prey than an average species in the web. This makes obvious intuitive sense. More interestingly it is also found that a prey of the deleted species is twice as likely to go extinct as an average species in the web when its predator is removed. The most likely explanation for this is that an effect known as predator-mediated coexistence or keystone predation is operating. This is a mechanism whereby a predator allows inferior competitors to coexist with a superior competitor by predation of the competitively dominant species. Thus the deletion of the predator can lead to extinctions amongst its prey. This has been observed in real communities (Paine 1974, Lubchenco 1978, Navarrete and Menge 1996) and in theoretical studies of simple dynamical systems (Fujii 1977, Shigesada and Kawasaki 1988).

The plausibility of the above mechanism is supported by the fact that competitors of the deleted species are five times less likely to go extinct than an average species, suggesting that competition plays an important part in structuring these communities. The effect on indirect predators is also significant: they are five times as likely to go extinct following the removal of their indirect prey compared to the average. It is interesting to note that the effect on an indirect predator is less than that on a direct predator. We will return to this later.

### 4.1 Predator-prey interaction strengths and extinction probabilities

The results in Table 2 show that removing the prey, or indirect prey, of a species has a negative effect on the predator. Clearly the importance of a prey species to its predator will vary between predator-prey pairs. A measure of the positive effect that a prey species has on the predator will be the effort $`f_{ij}`$, corresponding to the fraction of the prey in the predators diet (Ulanowicz and Puccia 1990).

We can derive a similar quantity for the positive effect that an indirect prey species has on its indirect predator. Consider the square of the $`f`$ matrix $`f_{ij}^2=_{k=0}^Nf_{ik}f_{kj}`$, where the sum is over all species in the web and $`k=0`$ corresponds to the environment. Since $`f_{kj}`$ is the fraction of species $`k`$’s diet that comes from species $`j`$, then $`f_{ik}f_{kj}`$ is the fraction of species $`i`$’s diet that comes from species $`j`$ via species $`k`$. If we now sum $`k`$ over all species in the web, we obtain the fraction of species $`i`$’s diet that comes from species $`j`$ through all paths of length 2. Therefore if we define the matrix $`F=_{n=1}^{\mathrm{}}f^n`$, then $`F_{ij}`$ is the fraction of species $`i`$’s diet that comes from species $`j`$ via all possible paths. Since all resources originally derive from the environment, $`F_{i0}=1`$ for any species $`i`$. For the direct matrix $`f_{ij}`$, normalisation ensures that $`_{j=0}^Nf_{ij}=1.0`$. However for the indirect matrix,

$$\underset{j=0}{\overset{N}{}}F_{ij}=\underset{j=0}{\overset{N}{}}\left\{f_{ij}+\underset{n=2}{\overset{\mathrm{}}{}}(f^n)_{ij}\right\}\underset{j=0}{\overset{N}{}}f_{ij}=1.0,$$
(7)
with equality only when species $`i`$ is a basal species.

If the infinite series of matrices converges, it can be calculated using $`F=(1f)^1f`$. All eighty $`f`$ matrices in this study were such that the sum was finite, however we have not been able to prove the convergence of the sum for a general web.

We are interested in separating the effects of indirect and direct predation. This can be achieved by defining a further matrix $`I`$ such that

$$I_{ij}=\{\begin{array}{cc}F_{ij},\hfill & \text{ if }f_{ij}=0.0\hfill \\ 0,\hfill & \text{ if }f_{ij}>0.0\text{ .}\hfill \end{array}$$
(8)
Thus if $`I_{ij}>0`$, then species $`i`$ will be an indirect predator of species $`j`$ according to the definition given above.

We can now examine how the probability that a predator or indirect predator goes extinct, following removal of its prey, varies with the proportion of the prey in the predators diet. This was carried out by dividing all the predator-prey (indirect predator-prey) pairs in the four hundred ecosystems into bins of size 0.05 according to their associated $`f_{ij}`$ ($`I_{ij}`$) values. The proportion of prey removals that resulted in the extinction of the predator (indirect predator) was then calculated for each bin and plotted in Fig. 3.

The interpretation of Fig. 3 is quite straightforward. Ignoring fluctuations attributable to sampling effects, both curves are monotonic. As the fraction of the predators, or indirect predators, diet that is obtained from the deleted species increases, so does the probability of predator extinction. However the probability of extinction does not become large until the prey constitutes a significant fraction of the predators diet. For direct predators it does not reach 10% until $`f_{ij}0.65`$ and for indirect predators this occurs when $`I_{ij}0.45`$. This probably arises from incorporating adaptive foraging into the population dynamics: predators can survive events that remove a large portion of their prey. In both cases the extinction probability rapidly approaches, but does not quite reach, 100% as $`f_{ij}`$ ($`I_{ij}`$) approaches 1.0. In fact, when a prey constitutes greater than 99% of a predators diet its removal leads to predator extinction only 97% of the time. This again illustrates the effect of adaptive foraging: 3% of the time the predator must be surviving by exploiting a prey species that previously formed less than 1% of its diet.

The two curves in Fig. 3 have a similar form, but at any given $`f_{ij}`$ or $`I_{ij}`$ value the probability of extinction is greater for indirect predators. This seems to contradict Table 2, which shows a greater probability of deletion for direct predators than indirect predators. However this can be explained by the distribution of the non-zero elements of the $`I`$ and $`f`$ matrices. If there are more small values of $`I_{ij}`$ than $`f_{ij}`$, then defining a predator (indirect predator) as having $`f_{ij}>0.01`$ ($`I_{ij}>0.01`$) will lead to a lower probability of deletion for indirect predators than direct predators, even though for a given value of $`f_{ij}`$ ($`I_{ij}`$) the effect of removing the indirect predator is more significant. That this is the case can be seen from Fig. 4, where the probability distributions of $`f_{ij}`$ and $`I_{ij}`$ for all $`f_{ij}>0.01`$ and $`I_{ij}>0.01`$ are shown. This does not of course explain why an indirect predator is more likely to go extinct when prey constituting a given fraction of its diet is removed. This may arise from the assumption that the net importance of multiple paths to the predator can be obtained by simply summing their individual weights implicit in Eq. (7).

We showed in Table 2 that a prey species has an increased probability of extinction following removal of its predator. We might also expect that this effect will depend on the fraction of the prey in the predators diet. This is investigated in Fig. 5, where the probability of prey extinction following predator removal has been estimated as a function of $`f_{ij}`$, by placing all the predator-prey interactions from the four hundred food webs into bins of size 0.1, and calculating the proportion of deletions that led to extinction of the prey. From this graph we see that prey extinction probability peaks at an intermediate value around $`f_{ij}=0.25`$, where slightly less than 15% of prey are going extinct. This is what we would expect if the mechanism for prey extinction is predator-mediated coexistence as proposed above, since if a prey forms only a fraction of a predators diet, then that predator is also likely to be exploiting its competitors.

### 4.2 Trophic Levels

In Table 2 we categorised species according to their trophic relationship to the deleted species. Another way to categorise species in a food web is by trophic level. As mentioned in Section 3 we use the minimum path length definition of trophic level (Yodzis 1984). Having assigned a trophic level to each species in our ecosystems, we calculated the mean number of species on each level that went extinct as a function of the trophic level of the deleted species. These results are shown in Table 3, where the same averaging procedure was used as for Table 2, that is, the average was performed over all the deleted species in all eighty ecosystems.

The analysis of this data is aided by statistics on how the number of species $`S`$, the average prey number $`S_{\mathrm{prey}}`$, and average predator number $`S_{\mathrm{predator}}`$, vary between trophic levels. The number of species $`S`$ in each trophic level was calculated by averaging over the eighty webs in this study. The other statistics were averaged over the total number of species in each trophic level. These statistics together with standard deviations in brackets are shown in Table 4.

Several patterns emerge from Table 3. For example, if we remove a species from a particular level then this leads to extinctions on the level above and this effect diminishes as we increase the level of the deleted species. Thus the mean number of level 2 species going extinct upon removal of a level 1 species is 3.651, which is roughly seven times the mean number of level 3 species going extinct following a deletion on level 2 (0.560), which in turn is about forty times larger than the mean number of extinctions in level 4 as a result of a deletion on level 3 (0.015). We can explain this with reference to Table 4, from which we see that while the mean number of prey increases rapidly with trophic level, the mean number of predators decreases. These two effects will be complementary, so that as we increase the trophic level of the deleted species, the number of predators of that species decreases, and those predators become less specialised on the deleted species. This leads to less extinctions in the level above when the species is removed. As a concrete example, consider deletion of a level 1 species. This will, on average, affect about 4 species on level 2 and these 4 species will be almost entirely dependent on the deleted species (as mentioned previously, deletion of the only prey of predator leads to extinction of that predator almost 100% of time). Thus we would expect about 4 predators in level 2 to go extinct as a result of a deletion on level 1. This compares well with the actual value of 3.651. We also see from Table 3 that the effect of deleting a species on the level above propagates up trophic levels, so that removing a level 1 species leads to deletions on level 2, which in turn leads to extinction on levels 3 and 4. This is as we would expect, given that we have already illustrated that consequences of species deletion can propagate up food webs.

As well as extinctions occurring in the levels above a deleted species, we find that deleting a species can lead to extinctions in the level below it. The magnitude of this effect seems to decrease with the trophic level of the deleted species: the mean number of level 3 species going extinct given deletion on level 4 is 0.709, compared to 0.359 for the number of level 2 species going extinct as the result of a deletion in level 3. Similarly, the deletion of a level 2 species leads to an average of just 0.013 extinctions on level 1. These observations can be explained if we assume that predator-mediated coexistence causes the extinctions, since the more prey a predator has, the stronger we might expect this effect to be. Then the pattern of increasing extinctions in the level below, with increasing level of the deleted species, merely reflects the increase in average prey number with trophic level observed in Fig. 4. The question of why the model food webs show these changes in mean predator and prey number with trophic level is not addressed here, but is discussed in Quince et al. (unpublished).

## 5 Keystone species

Thus far the focus of this study has been on pairs of species and how the trophic relationship between them influences the effect that deleting one species will have on the other. It was found that the effects of deletion can propagate both up and down food webs. We now consider the related question of what factors determine the impact that deleting a species will have on the whole web. Specifically, we ask whether there are any consistent differences between species whose removal causes little change in the food web, and those which play a major role in structuring the community. The latter are sometimes referred to as “keystone species”, although this term can be restricted to those species whose importance is large relative to their population size (Power et al. 1996). Here we will simply be interested in identifying factors that statistically influence the number of further extinctions that follow the removal of a species.

The trophic level of the deleted species has an effect on the expected number of further extinctions. This is shown in the first column of Table 3, where the mean number of further extinctions in the web as a whole is seen to be largest for species on trophic level 1. This probably reflects the greater importance of effects propagating up, rather than down, the food web. An alternative way to categorise species according to trophic role, which is better at separating these two processes, is into the following three classes: ‘top’ (species with no predators), ‘intermediate’ (species with both predators and prey) and ‘basal’ (species with no prey). When this is done, it is found that removal of a basal species causes an average of 6.58 further extinctions as opposed 0.65 and 0.61 for intermediate and top species respectively. Comparing the figures for basal and top species reveals that top-down processes are indeed statistically less important than bottom-up in these communities.

The above categories are quite broad. One property which we can use to more finely discriminate between species is the number of other species they interact with through predator-prey links. We might expect this to correlate with the impact of deleting the species on the food web. That this is indeed the case is shown in Fig. 6, where the average number of further extinctions is plotted as a function of the node degree $`k`$ in network terminology. In the case of food webs the latter corresponds to the number of predators plus the number of a prey of a species. The results are shown for basal, intermediate and top species separately. For all three classes the average number of further extinctions increases with $`k`$, so the most connected species are the ones whose removal has the greatest effect on the food web structure. The results for top species provide further support, albeit circumstantial, that the top-down effect in our webs is predator-mediated coexistence, since we would expect the importance of this effect to increase with the number of prey of the predator. They can be contrasted with an earlier study that failed to find such a relationship (Pimm 1980).

These results have relevance for the studies of real food web robustness to deletion mentioned in the introduction. These obviated the need for an explicit dynamics by considering multiple removals, and judging a further extinction to have occurred when all the prey of a species are absent (Solé and Montoya 2001, Dunne et al. 2002). These studies found that removing the most connected species resulted in more secondary extinctions for the same number of species removed. Our results suggest that if a population dynamics was included in these studies, the food web structures would be even more sensitive to the removal of highly connected species.

The correlation between $`k`$ and further extinction number is very strong, especially given that it quantifies the number of direct interactions of a species, and does not give information on the strength of indirect effects, which we have already shown to be important (Jordan and Scheuring 2002). It is difficult to devise measures of species importance that do incorporate indirect effects. Potential candidates are the bottom-up and top-down keystone indices of Jordan et al. (1999). These were originally devised for binary food webs, but is easy to extend them to the model food webs including diet compositions considered here. In fact the bottom-up keystone index for a species $`i`$ is simply

$$K_b^i=\underset{j=1}{\overset{S}{}}F_{ji}.$$
(9)
Here the matrix $`F`$ is defined by $`F=_{n=1}^{\mathrm{}}f^n`$ and, as discussed in the previous section, its elements give the fraction of species $`i`$’s diet that comes from species $`j`$ via all possible paths. This makes the meaning of $`K_b^i`$ clear: it is the total number of species that depend on $`i`$ directly or indirectly for resources. The complementary quantity $`K_t^i`$ measures the strength of top-down effects it can be defined in a similar way for non-binary food webs, which we will not discuss here. There is in fact a very good correlation between $`K_b^i`$ and further extinction number for basal species. This shown in Fig. 7. This suggests that this quantity is effective at predicting the strength of bottom-up effects. We also found a correlation between $`K_t^i`$ and the number of further extinctions following the removal of a top species (Fig. 8) although not as strong as for basal species. Furthermore the composite quantity $`K^i=K_t^i+K_b^i`$ did have some success at predicting the outcome of removing an intermediate species. This suggests that although top-down effects might be more difficult to quantify with these network measures, and they are not yet as effective as simply counting the number of predators and prey, they may well prove a practical way of predicting the importance of species in food webs.

## 6 Stability of food webs to deletion

In the previous section we focused on the effect that deleting one species has on another in the model food webs. Here we will consider the stability of the communities as a whole. In particular we will ask whether that stability is determined by the topological structure of the ecosystem. We begin by defining deletion stability as the fraction of species that have the property that, when they alone are deleted, further extinctions do not take place. This definition is similar to that given by Pimm (1979), and has the advantage that it should not have any in-built dependence on species number. Its calculation is illustrated in Fig. 9.

We will use the four hundred model food webs detailed in Section 3 as our data set and describe their topological structure with the five food web statistics also defined there. These five statistics quantify four food web properties that it has been suggested may impact ecosystem stability: complexity in terms of number of species and linkage density ($`S`$, $`L/S`$), the amount of omnivory in the web ($`O`$), the importance of cycles ($`𝒞`$) and the redundancy ($`𝓇`$). We remind the reader that the means and standard deviations of these statistics for the data set are given in Table 1.

If we possessed a range of food webs that vary for each property independently, whilst the other properties remain fixed, then this analysis would be quite simple. Indeed this was the approach adopted in previous studies using small pre-defined food web structures (Pimm 1979, Pimm 1980, Borrvall et al. 2000). However by generating the food webs by evolving a number of communities at a particular set of parameter values, we are faced with a less straightforward situation. The food webs vary for all properties simultaneously, and some of the properties are significantly correlated with one another. In fact, the situation is much like analysing real food web structures. The advantage is that our range of structures are much more complex and realistic, the disadvantage is that we are forced to adopt a statistical approach.

We will commence our analysis by examining the correlation matrix for all the variables, both the dependent variable $`S_d`$, and food web properties. This is shown in Table 5. The deletion stability is significantly correlated with the fraction of omnivores (negatively) and the redundancy (positively). However because of the many significant correlations between the variables themselves, we can not conclude that a smaller proportion of omnivores or greater redundancy will be associated with higher deletion stability, all other properties being unchanged.

One way to shed some light on this problem is to use multivariate linear regression (Jobson 1991). This statistical procedure essentially assumes that the dependent variable, in this case $`S_d`$, can be described as a linear function of the independent variables, the five food web properties plus a noise term. The results of such an analysis are shown in Table 6. Examining the regression coefficients, $`\mathrm{b}_\mathrm{j}`$, for the individual variables we see that two, the redundancy again and the links per species, are deemed to have a strong influence on the deletion stability, and that both of these effects are positive. More pertinently however, examining the fit as whole we see from the $`^2`$ value that only 8% of the variation in stability is being explained by these variables.

We therefore conclude that for this data set a fairly robust positive relationship exists between the proportion of redundant species in the web and the stability of the web to deletion. This result makes sense. We defined a redundant species as one that is functionally equivalent, in the sense of possessing the same predators and prey, to at least one other species. Thus the removal of a redundant species is unlikely to cause further species to go extinct, since its functional equivalents should be able to increase their population sizes and compensate for the loss. This supports the hypothesis that increased redundancy in ecosystems will result in increased functional reliability (Walker 1992, Naeem 1998).

We find no evidence that complexity in terms of increased species number or links per species is destabilising. In fact there was evidence that increasing the latter actually reduced the probability of further extinctions when species are deleted. This can be compared to early work on deletion stability in small food webs modelled with Lotka-Volterra dynamics, where it was found that increasing either the number of species or connectance, $`L/S^2`$, rapidly decreased stability (Pimm 1980). The difference can probably be attributed to the more realistic structures and global dynamics used here, in particular incorporating adaptive foraging into the population dynamics.

This result can also be compared to more recent work examining the effect of deleting species from simple three level food webs, constructed such that all species on a trophic level were functionally equivalent (Borrvall et al. 2000). These authors found that stability to deletion increased with the number of species on each trophic level. Their set-up corresponds to keeping $`𝓇`$ constant and equal to one, whilst increasing the total number of species, $`S`$, a variable for which no effect was observed here. This study also used Lotka-Volterra dynamics and this may explain the discrepancy. In particular, they observed that removing a predator had no effect on its prey, in contrast to the results detailed in Section 4.

We also found that omnivory and cycling were unimportant in determining deletion stability. In fact, all our structural properties taken together explained very little of the variation between the webs. It may be that the measure of deletion stability we used is intrinsically noisy — it does seem sensitive to the presence of one or two vulnerable species — or it may be that our statistics are not capturing the properties that are important in determining robustness to deletion.

## 7 Conclusion

In this paper we have shown that deletion experiments, which are very difficult and time-consuming to carry out in real communities, can be easily implemented on model webs, which previous studies have shown have many of the characteristics of real webs (Drossel et al. 2001). We found a number of interesting results on species removal from communities which differed from some previous studies. The food webs as a whole were shown to be quite robust to deletions, and individual species were able to survive the loss of prey constituting a major fraction of their diet. Deletions were shown to cause further extinctions amongst species both above and below the deleted species in the food web. These phenomena arose out of the complex population dynamics used in the model which incorporates adaptive foraging. They illustrate the importance of using a realistic global dynamics when considering community responses to large scale perturbations such as deletion. This contrasts with studies that either lack an explicit dynamics (Solé and Montoya 2001, Dunne et al. 2000) or use Lotka-Volterra equations (Pimm 1979, Pimm 1980, Borrvall et al. 2000, Lundberg et al. 2000).

In addition to a realistic dynamics, effective studies of species removal require a realistic set of structures. By using a set of large webs composed of species spanning a range of trophic roles, we were able to show that removing the most connected species resulted in the most further extinctions and that recently developed ‘keystone indices’ were fairly effective at predicting species importance. A study using small simple webs failed to find these relationships (Jordan et al. 2002). The range of food web structures studied allowed us to show the important role that redundant species play in increasing food web robustness to deletion, and find that there was no correlation between increased complexity and decreased stability to deletion. This adds another component to the stability-complexity debate.

This paper represents the first attempt to join complex food web structures with realistic population dynamics to study species loss from communities. It would be interesting to see if the results would change if real, rather than evolved, food web structures were used with such a dynamics or if the particular choice of dynamics were changed. In any case, we believe that we have shown that it is crucial for models to display a degree of realism, if reliable deductions concerning the consequences of species deletions are to be made.

Acknowledgements: We would like to thank R. Ulanowicz for letting us have the source code for his program on cycling in food webs. We also wish to thank R. Law and M. Rattray for useful discussions. CQ thanks the EPSRC (UK) for financial support during the initial stages of this work.

## References

Abrams, P. A. 1996. Dynamics and interactions in food webs with adaptive foragers. - In: Polis, G. A. and Winemiller, K. O. (eds), Food webs: integration of patterns and dynamics. Chapman & Hall, pp. 113-121.

Borrvall, C., Ebenman, B. and Jonsson, T. 2000. Biodiversity lessens the risk of cascading extinction in model food webs. - Ecol. Lett. 3: 131-136.

Caldarelli, G., Higgs, P. G. and McKane, A. J. 1998. Modelling coevolution in multispecies communities. - J. Theor. Biol. 193: 345-358.

Carlton, J. T. and Geller, J. B. 1993. Ecological roulette – the global transport of nonindigenous marine organisms. - Science 261: 78-82.

Drake, J. A. 1990. The mechanics of community assembly and succession. - J. Theor. Biol. 147: 213-233.

Drossel, B., Higgs, P. G. and McKane, A. J. 2001. The influence of predator-prey population dynamics on the long-term evolution of food web structure. - J. Theor. Biol. 208: 91-107.

Dunne, J. A., Williams, R. J. and Martinez, N. D. 2002. Network structure and biodiversity loss in food webs: robustness increases with connectance. - Ecol. Lett. 5: 558-567.

Elton, C. S. 1958. The ecology of invasions by animals and plants. \- Methuen, London.

Fretwell, D. S. and Lucas, H. L. 1970. On territorial behaviour and other factors influencing habitat distribution in birds. - Acta Biotheor. 19: 16-32.

Fujii, K. 1977. Complexity-stability relationship of two-prey-one-predator species system model: local and global stability. \- J. Theor. Biol. 69: 613-623.

Jobson, J. D. 1991. Applied multivariate data analysis. \- Springer-Verlag.

Jordan, F., Takacs-Santa, A. and Molnar, I. 1999. A reliability theoretical quest for keystones. - Oikos 86: 453-462.

Jordan, F. and Scheuring, I. 2002. Searching for keystones in ecological networks. - Oikos 99: 607-612.

Jordan, F., Scheuring, I. and Vida, G. 2002. Species positions and extinction dynamics in simple food webs. - J. Theor. Biol. 215: 441-448.

Law, R. and Morton, R. D. 1996. Permanence and the assembly of ecological communities. - Ecology 77: 762-775.

Lockwood, J. L., Powell, R. D., Nott, P. and Pimm, S. L. 1997. Assembling ecological communities in time and space. - Oikos 80: 549-553.

Lubchenco, J. 1978. Plant species-diversity in a marine inter-tidal community: importance of herbivore food preference and algal competitive abilities. - Am. Nat. 112: 23-39.

Lundberg, P., Ranta, E. and Kaitala, V. 2000. Species loss leads to community closure. - Ecol. Lett. 3: 465-468.

MacArthur, R. H. 1955. Fluctuations of animal populations and a measure of community stability. - Ecology 36: 533-536.

Mathsoft Inc. 2000. S-Plus software: version 6.0. - MathSoft Inc., Cambridge, MA.

May, R. M. 1972. Will a large complex system be stable? \- Nature 238: 413-414.

May, R. M. 1973. Stability and complexity in model ecosystems. \- Princeton Univ. Press.

Milner-Gulland, E. J and Bennett, E. L. 2003. Wild meat: the bigger picture. - Trends Ecol. Evol. 18: 351-357.

Morton, R. D. and Law, R. 1997. Regional species pools and the assembly of local ecological communities. - J. Theor. Biol. 187: 321-331.

Navarrete, S. A. and Menge, B. A. 1996. Keystone predation and interaction strength: interactive effects of predators on their main prey. \- Ecol. Monogr. 66: 409-429.

Naeem, S. 1998. Species redundancy and ecosystem reliability. \- Conserv. Biol. 12: 39-45.

Odum, E. P. 1953. Fundamentals of ecology. - Saunders.

Paine, R. T. 1974. Intertidal community structure: experimental studies on the relationship between a dominant competitor and its principal predator. - Oecologia 15: 93-120.

Paine, R. T. 1980. Food webs : linkage, interaction strength and community infrastructure. - J. Anim. Ecol. 49: 667-685.

Pimm, S. L. 1979. Complexity and stability – another look at MacArthur’s original hypothesis. - Oikos 33: 351-357.

Pimm, S. L. 1980. Food web design and the effect of species deletion. - Oikos 35: 139-149.

Power, M. E., Tilman, D., Estes, J. A., et al. 1996. Challenges in the quest for keystones. - Bioscience 46: 609-620.

Quince, C., Higgs, P. G. and McKane, A. J. 2002. Food web structure and the evolution of ecological communities. - In: Lässig, M. and Valleriani, A. (eds), Biological Evolution and Statistical Physics. Springer-Verlag, pp. 281-298.

Shigesada, N. and Kawasaki, K. 1988. Direct and indirect effects of invasions of predators on a multiple-species community. - Theor. Pop. Biol. 36: 311-328.

Solé, R. V. and Montoya, J. M. 2001. Complexity and fragility in ecological networks. - Proc. Roy. Soc. Lond. B268: 2039-2045.

Ulanowicz, R. E. 1983. Identifying the structure of cycling in ecosystems. - Math. Biosci. 65: 219-237.

Ulanowicz, R. E. and Puccia, C. J. 1990. Mixed trophic impacts in ecosystems. - Coenoses 5: 7-16.

Walker, B. H. 1992. Biodiversity and ecological redundancy. \- Conserv. Biol. 6: 18-23.

Yodzis, P. 1989. Introduction to Theoretical Ecology. \- Harper & Row.

Yodzis, P. 2001. Must top predators be culled for the sake of fisheres? - Trends Ecol. Evol. 16: 78-84.


format: text

text:
# A prime prime primer

## 1 Introduction

Interest in prime numbers continues apace, e.g. computing, secure communications, quantum computing, statistical physics, engineering and, of course, mathematics. However, we are still far from finding a useful function $`p(n)`$ which generates all the primes in order. We present such a function next and discuss its usefulness . All variables below, save $`x`$, are integers.

## 2 Analysis

Consider the product function

$$p(n)=\underset{n^{}=2}{\overset{\sqrt{n}}{}}\left(\frac{n}{n^{}}\frac{n}{n^{}}\right),$$
(1)
where the integer floor $`x`$ (ceiling $`x`$) functions give the greatest (smallest) integer not greater (not smaller) than $`x`$: e.g. $`7=7=7`$, $`3.7=3`$, $`3.7=4`$, $`7.3=7`$, $`7.3=8.`$

Then, for $`n=2,3,4\mathrm{}`$,

$$p(n)=\{\begin{array}{cc}1,& \text{if n is prime}\hfill \\ 0,& \text{if not}\hfill \end{array}$$
(2)
Also, we have directly from (2) that the number $`\pi (n)`$ of primes not exceeding $`n`$ is

$$\pi (n)=\underset{n^{}=1}{\overset{n}{}}p(n^{})$$
(3)
## 3 Discussion

* That function (1) satisfies (2) may be seen by inspection of $``$ \- $``$:

$$xx=\{\begin{array}{cc}0,& \text{if }x\text{ is an integer}\hfill \\ 1,& \text{if not}\hfill \end{array}$$
(4)
Thus $`p(n)=1`$ only if all the divisors $`n^{}`$ of $`n`$ give non-zero remaninders, which is the definition of primality.
* What constitutes a “useful function” $`p(n)`$ is partly subjective , so we describe our $`p(n)`$ as: simple, compact, transparent, and using only two basic notions (floor $``$) and (ceiling $``$) of integer mathematics.
* We have programmed and checked algorithm (1,2,3) on a PC.

## Reference

 Hardy, G. H. & Wright, E. M. (1978) ’An introduction to the theory of numbers’, ( 5th edition). Clarendon Press, Oxford.


text:
# Probing the atmosphere of the bulge G5III star OGLE-2002-BUL-069 by analysis of microlensed H𝛼 line Based on observations made at ESO, 69.D-0261(A), 269.D-5042(A), 169.C-0510(A)

## 1 Introduction

Near extended caustics produced by binary (or multiple) lenses, the source star undergoes a large total magnification in brightness. Furthermore, its surface is differentially magnified because of a strong gradient in magnification. The relative lens-source proper motion is typically slow enough to allow the light curve to be frequently sampled. This translates to a high spatial resolution on the source star’s surface and hence permits its radial brightness profile to be inferred from the observations. Over the past four years, coefficients characterizing linear or square-root limb-darkening profiles have been obtained with the microlensing technique for several Bulge giants and sub-giants (Fields et al. 2003 and references therein) and a main sequence star (Abe et al. 2003). A new generation of stellar atmosphere models (Orosz & Hauschildt 2000) have revealed limb-darkening laws that are significantly different from the traditional analytic ones (Bryce et al. 2002; Claret & Hauschildt 2003). The centre to limb variation of spectral lines can show markedly different behaviour from that of the continuum. Most moderately strong and weak lines weaken towards the limb, but resonant scattering lines can vary in a much more pronounced way (Loeb & Sasselov 1995). In cool giants, $`{}_{}{}^{}\mathrm{H}_{}^{}\beta `$, being formed lower in the atmosphere, is more limb-darkened than $`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$. $`{}_{}{}^{}\mathrm{TiO}_{}^{}`$ and $`Ca\text{ii}`$ show strong variations, and some lines and bands may even be limb-brightened. Intensive spectroscopic monitoring of a caustic exit at high resolution with high S/N should reveal temporal variations in the equivalent widths of promising spectral lines that can be compared with predictions from stellar atmosphere models (Heyrovský et al. 2000).

We present here the first photometric and spectroscopic monitoring campaign that has successfully been performed at high resolution with dense sampling. Previously, Castro et al. (2001) obtained two KECK HIRES spectra of EROS 2000-BLG-5, but missed the trailing limb of the caustic where the effects are stronger, while the data in the qualitative analysis of Albrow et al. (2001) involved dense coverage but at low resolution. Afonso et al. (2001) provided a model that reproduces the photometric data and found an excess of $`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$ at the limb. They attribute the excess to chromospheric emission, but it could be due to shortcomings of the synthetic spectra. Using a determination of the spectral type of the source star, we have computed the limb-darkening profiles in $`\mathrm{R}`$ and $`\mathrm{I}`$ from appropriate PHOENIX synthetic spectra of the source star and fitted a fold-caustic model to our photometric data obtained during the caustic exit. This model has been used here to compute the synthetic spectrum for wavelengths around the $`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$-line. A more detailed report on the determination of the stellar parameters and the analysis of other spectral lines will be presented elsewhere (Beaulieu et al. 2004), as will the details of a full binary lens model fit to the PLANET observations of the event (Kubas et al. 2004).

## 2 OGLE 2002-BLG-069 photometry and spectroscopy

The PLANET collaboration, comprising six different telescopes, namely, SAAO 1.0m (South Africa), Danish 1.54m (Chile), ESO 2.2m (Chile), Canopus 1.0m (Australia), Mt. Stromlo 50in (Australia) and Perth 0.6m (Australia), commenced photometric observations of OGLE 2002-BLG-069 alerted by the OGLE-III early Warning System (Udalski 2003) in early June, 2002. From online data reductions on 25 June it was realised that the event involved a binary lens and a bright giant source star. With the source taking $`1.4`$ days to cross the caustic on entry, it appeared to be an excellent candidate for time-resolved spectroscopy of the caustic exit. Using the predictions based on modeling our photometry, and thanks to excellent coordination with the staff at La Silla and Paranal, very good coverage was obtained during the caustic exit together with post-caustic reference spectra. These observations were performed using the UVES spectrograph mounted on Kueyen (VLT UT2) as part of Target-of-Opportunity and Director’s Discretionary time. Thirty-nine spectra with exposure time of 20 min (1 hour for the post caustic observations) were obtained alternately in the so-called standard settings, 580 and 860, covering the full $`480010600`$ Å range at a resolution of 30000 with S/N ranging from 50 to 130. From an analysis of the curve-of-growth of 100 Fe i and Fe ii lines as done by Minniti et al. (2002), and independently from the study of the Ca ii and $`{}_{}{}^{}\mathrm{Mg}_{}^{}`$ lines (Jørgensen et al. 1992), a good fit to all the data was obtained with a plane-parallel model atmosphere having $`T_{\mathrm{eff}}=5000`$K, $`\mathrm{log}(g)=2.5`$, $`v_{\mathrm{turb}}=1.5\text{km}\text{s}^1`$ and $`[{}_{}{}^{}\mathrm{Fe}_{}^{}/{}_{}{}^{}\mathrm{H}_{}^{}]=0.6`$, The calibration of the MK spectral type gives, for a G5III star, $`T_{\mathrm{eff}}=5050`$K, $`\mathrm{M}_\mathrm{V}=+0.9`$, $`(\mathrm{V}\mathrm{I})=0.95`$, $`\mathrm{M}/\mathrm{M}_{}=1.1`$ and $`R/R_{}=10`$. At HJD=2452455.2754, the source was amplified by 17.7. The measured color from SAAO data are $`I=13.01\pm 0.01`$, $`(VI)=2.06\pm 0.02`$, so $`\mathrm{E}(\mathrm{V}\mathrm{I})=1.11\pm 0.02`$. We adopt a conservative $`\mathrm{A}_\mathrm{V}/\mathrm{E}(\mathrm{V}\mathrm{I})=2.2\pm 0.3`$, and then derive a distance of $`9.4\pm 1.4\mathrm{kpc}`$.

## 3 Modeling of the caustic exit

With $`\rho `$ denoting the fractional radius of the source star, we approximate its wavelength-dependent brightness profile $`I^\lambda (\rho )=\overline{I}\xi ^\lambda (\rho )`$ by concentric rings of constant intensity. A model atmosphere for the stellar parameters given in the previous section was computed with the PHOENIX grid v2.6 (based on the code described by Allard et al. (2001) but using spherical geometry and spherically symmetric radiative transfer). A synthetic spectrum was calculated from this model at a spectral resolution of 0.05 Å at 128 steps in source radius. The intensity profile was interpolated with cubic splines at 1000 points equally spaced in $`\mathrm{cos}\vartheta =\sqrt{(1\rho ^2)}`$ where $`\vartheta `$ is the emergent angle, so that the width of the rings of constant intensity decreases toward the stellar limb. For a broadband filter $`(s)`$, the intensity profile $`\xi ^{(s)}(\rho )`$ is obtained by convolving the source brightness profile with the filter, CCD transmission and atmosphere response functions. The observed magnitude $`m^{(s)}(t)`$ reads $`m^{(s)}(t)=m_\mathrm{S}^{(s)}2.5\mathrm{lg}\left\{A^{(s)}(t)+g^{(s)}\right\}`$ where $`m_\mathrm{S}^{(s)}`$ is the intrinsic source magnitude, $`g^{(s)}=F_B^{(s)}/F_S^{(s)}`$ is the ratio between background flux $`F_B^{(s)}`$ and intrinsic source flux $`F_S^{(s)}`$, and $`A^{(s)}(t)`$ the source magnification at time $`t`$.

The photometric data presented in the lower panel of Fig. 1 show a peak with the characteristic shape resulting from a fold-caustic exit. For a source in the vicinity of a fold caustic, the magnification $`A^{(s)}(t)`$ can be decomposed (e.g., Albrow et al. 1999a) as $`A^{(s)}(t)=A_{\mathrm{crit}}^{(s)}(t)+A_{\mathrm{other}}(t)`$, where $`A_{\mathrm{crit}}^{(s)}(t)`$ denotes the magnification of the critical images associated with the caustic, and $`A_{\mathrm{other}}(t)`$ denotes the magnification of the remaining images of the source under the action of the binary lens.

Let us consider a uniformly moving source crossing the caustic within $`2\mathrm{\Delta }t`$ and exiting the caustic by its trailing limb at $`t_\mathrm{f}`$. If one neglects the curvature of the caustic and the variation of its strength over the size of the source, the magnification of the critical images reads $`A_{\mathrm{crit}}^{(s)}(t)=a_{\mathrm{crit}}G_\mathrm{f}(\frac{tt_\mathrm{f}}{\mathrm{\Delta }t};\xi ^{(s)})`$ where $`G_\mathrm{f}(\eta ;\xi ^{(s)})`$ is a characteristic function (Schneider & Wagoner 1987) which solely depends on the intensity profile $`\xi ^{(s)}`$. For $`|\omega (tt_\mathrm{f})|1`$, the magnification of the non-critical images can be approximated by $`A_{\mathrm{other}}(t)a_{\mathrm{other}}\left[1+\omega (tt_\mathrm{f})\right]`$ where $`a_{\mathrm{other}}`$ is the magnification at the caustic exit at time $`t_\mathrm{f}`$ and $`\omega `$ measures its variation.

As the period during which our generic fold-caustic model is believed to be a fair approximation, we choose the range $`2463.45\text{HJD}24500002467`$. Restricting our attention to those data sets with more than two points in this region leaves us with 29 points from SAAO and 15 points from UTas in I as well as 98 points from the ESO 2.2m in R, amounting to a total of 142 data points. From a binary lens model of the complete data sets (Kubas et al. 2004) for these observatories and filters, we have determined the source magnitudes $`m_\mathrm{S}^{(s)}`$ and the blend ratios $`g^{(s)}`$. With the adopted synthetic spectra, we compute the stellar intensity profiles for I and R, and use these for obtaining a fit of the generic fold-caustic model to the data in the caustic-crossing region by means of $`\chi ^2`$-minimization, which determined the 5 model parameters $`t_\mathrm{f}`$, $`\mathrm{\Delta }t`$, $`a_{\mathrm{crit}}`$, $`a_{\mathrm{other}}`$, and $`\omega `$ as shown in Tab. 1. If a classical linear limb darkening is included in the list of parameters to fit (as in Albrow et al. 1999b), the best fit is obtained with $`\mathrm{\Gamma }=0.5`$.

The modeled light curve and the structure of the residuals of the two fits are given in the lower panels of Fig. 1. The fit with adopted PHOENIX limb darkenings clearly show systematic trends at the level of 1-2 %. The same trends in the residuals can be seen from the figure, which suggests some specific features lying in the real stellar atmosphere. We checked whether this discrepancy can be caused by the straight fold caustic approximation, but from running the already mentioned global binary lens model neither the very small curvature of the caustic nor the source trajectory can explain this systematic effect.

## 4 Comparison between synthetic spectra and spectroscopic data

The fold-caustic model derived in Sec. 3 is used to compute the source flux during the caustic exit in order to obtain a synthetic spectrum at each epoch in a consistent way. In Fig. 2 we show the observed and synthetic spectra at two epochs: 9 July at UT 00h04 and 9 July at UT 22h59 (trailing limb crossing the caustic).

A first analysis reveals good agreement between UVES and synthetic spectra for the wings part of $`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$ ($`65586561.8`$ Å and $`6563.86567.6`$ Å), whereas a clear discrepancy is observed in its core ($`6561.86563.8`$ Å). We note that the chromosphere is not included in the PHOENIX calculations. Althought it should have a minor influence on the broadband limb darkenings, it will have an effect on the core of lines like $`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$. We therefore divided the analysis in two parts: the wings and the core. In order to compare the observed equivalent widths of the $`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$-line with the predictions from our synthetic spectrum, we apply an overall scale factor to the equivalent width ($`1.035`$ for the wings part and $`1.495`$ for the core part of the line), so that the measurements derived from the fiducial model match the post-caustic observations.

The observed and predicted temporal variations of the equivalent widths of both the wings and the core of $`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$ are plotted in the upper part of Fig. 1. Note at the beginning a small short-term decrease in the model, when the leading limb hits the caustic, followed by little variation while the source centre crosses, but then a marked change when the trailing limb exits the caustic. The spectral profile of the core is not well matched by the model spectrum as seen in Fig. 2; both the width and the depth disagree. Furthermore, the differential variation in flux between the centre and the limb is not well reproduced, even after the scale factor has been applied. A better fit is obtained to the wings. Both parts of the line show larger equivalent width variations than predicted by the model.

Last but not least, in the spectrum corresponding to 10 July, UT 02h58 shown in Fig. 3, we clearly note a strong deviation in the core of the line but a much smaller one in the wing. For the purpose of the following discussion, we define the effective photospheric radius of the star as that at which the R band intensity has fallen to 5 % of the central intensity in R. Using the fit done with limb darkening derived from PHOENIX model, we infer that the spectrum started when the caustic was at a fractional radius $`\rho =0.996`$ and ended at $`\rho =1.015`$ while with a linear limb-darkening law, the spectrum started at a fractional radius $`\rho =1.007`$ and ended at $`\rho =1.026`$.

In both cases, it is clear that the spectrum was taken when the caustic was outside the photosphere in R. There is a clear emission peak in the core of the $`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$ absorption line, and this can only exist provided there is a temperature inversion in the atmosphere (i.e. a chromosphere). The doppler shift of the emission core is consistent with material moving outward through the giant source star’s chromosphere ($`{}_{}{}^{}\mathrm{H}_{}^{}\alpha `$ emission line) with a radial velocity of $``$ 7 $`\mathrm{kms}^1`$, the signature being amplified at the very end of the caustic exit.

###### Acknowledgements.

For this particular VLT follow up exercise, we are expressing our highest gratitude to the ESO staff at Paranal and La Silla whose flexibility and efficiency was absolutely vital to the success of the observations. We are very grateful to OGLE-III for their Early Warning System (EWS), to the observatories that support our science (European Southern Observatory, Canopus, CTIO, Perth, SAAO) via the generous allocations of time that make this work possible. We thank Andy Gould and Dimitar Sasselov for fruitful discussions. The operation of Canopus Observatory is in part supported by the financial contribution from David Warren, and the Danish telescope at La Silla is operated by IJAF financed by SNF. JPB acknowledges financial support via award of the “Action Thématique Innovante” INSU/CNRS. MD acknowledges postdoctoral support on the PPARC rolling grant PPA/G/O/2001/00475.


text:
# Fusion and breakup in the reactions of 6,7Li and 9Be

## 1 Introduction

The fusion reactions of neutron-rich nuclei around the Coulomb barrier provide a good opportunity to study the interplay between quantum tunneling and the breakup process. An important question here is: does the breakup process influence fusion reactions in a similar way as the inelastic process, which leads to subbarrier enhancement  of fusion cross sections over predictions for a single barrier? In addressing this question, measurements with weakly bound stable nuclei, such as 9Be, 6Li, and 7Li, have been proved to be useful. These beams are currently much more intense than radioactive beams, allowing more precise and extensive experimental studies. Also, these nuclei predominantly break up into charge fragments (9Be $``$ 2$`\alpha +n`$ or $`\alpha +^5`$He, 7Li $``$ $`\alpha +t`$, and 6Li $``$ $`\alpha +d`$), which are more easily detected. This allows a clean separation of the products of complete fusion from those of incomplete fusion, where only a part of the projectile is captured.

To describe the effect of breakup on fusion reactions theoretically, a model has to take into account the following three different processes: (i) the projectile as a whole is captured by the target, (ii) only one of the breakup fragments is captured, (iii) all the breakup fragments are absorbed after the breakup takes place near the target nucleus. Process (ii) is referred to as an incomplete fusion reaction, while both of processes (i) and (iii) lead to complete fusion. Calculations based on the continuum-discretized coupled-channels method may account for the total fusion cross section (the sum of all the processes) and/or the separate contribution of process (i), but it is not easy to compute cross sections for process (iii), that is, the breakup followed by complete fusion. In order to model this process, we need to follow the trajectories of the breakup fragments to determine whether one or both fragments are captured by the target nucleus. In this contribution, we present such a model, where many classical trajectories are sampled by the Monte Carlo method to compute the total complete fusion cross section. We note that a very similar model, called the three-body classical trajectory Monte Carlo (CTMC) method, has been developed in the field of atomic and molecular physics, and has been successfully applied to ion-atom ionization and charge transfer collisions .

## 2 Three-body classical trajectory Monte Carlo (CTMC) model

Let us assume the following Hamiltonian for a reaction of a projectile nucleus which consists of two cluster fragments $`P_1`$ and $`P_2`$,

$$H=\frac{𝒑_1^2}{2m_1}+\frac{𝒑_2^2}{2m_2}+\frac{𝒑_T^2}{2m_T}+V_{12}(𝒓_1𝒓_2)+V_{1T}(𝒓_1𝒓_T)+V_{2T}(𝒓_2𝒓_T),$$
(1)
where $`𝒑_1`$, $`𝒑_2`$, and $`𝒑_T`$ are the momenta of fragment $`P_1`$, $`P_2`$, and of the target, respectively. $`𝒓_1`$, $`𝒓_2`$, and $`𝒓_T`$ are the corresponding coordinates. The idea of the three-body classical trajectory model is to solve this Hamiltonian classically in the two-dimensional ($`x,y`$) plane . We assume that at $`t`$=0 the target is at rest at the origin, while the center of mass of the projectile is at $`(x,y)=(R_{\mathrm{max}},b)`$ moving towards the $`x`$ direction. Here, $`R_{\mathrm{max}}`$ is some large number, while $`b`$ is the impact parameter. The relative motion between the projectile fragments is confined at $`t=0`$ inside the potential $`V_{12}`$ at the energy corresponding to the breakup $`Q`$-value, $`Q_{\mathrm{bu}}`$, with a random orientation $`\theta `$ with respect to the $`x`$ axis. The initial conditions thus read,

$`x_1(0)`$ $`=`$ $`R_{\mathrm{max}}{\displaystyle \frac{m_2}{m_1+m_2}}r_0\mathrm{cos}\theta ,y_1(0)=b+{\displaystyle \frac{m_2}{m_1+m_2}}r_0\mathrm{sin}\theta ,`$ (2)
$`x_2(0)`$ $`=`$ $`R_{\mathrm{max}}+{\displaystyle \frac{m_1}{m_1+m_2}}r_0\mathrm{cos}\theta ,y_2(0)=b{\displaystyle \frac{m_1}{m_1+m_2}}r_0\mathrm{sin}\theta ,`$ (3)
$`v_{x1}(0)`$ $`=`$ $`v_{\mathrm{lab}}{\displaystyle \frac{1}{2m_1}}p_{\mathrm{rel}}\mathrm{cos}\theta ,v_{y1}(0)={\displaystyle \frac{1}{2m_1}}p_{\mathrm{rel}}\mathrm{sin}\theta ,`$ (4)
$`v_{x2}(0)`$ $`=`$ $`v_{\mathrm{lab}}+{\displaystyle \frac{1}{2m_2}}p_{\mathrm{rel}}\mathrm{cos}\theta ,v_{y2}(0)={\displaystyle \frac{1}{2m_2}}p_{\mathrm{rel}}\mathrm{sin}\theta ,`$ (5)
$`x_T(0)`$ $`=`$ $`y_T(0)=v_{xT}(0)=v_{yT}(0)=0,`$ (6)
where

$$v_{\mathrm{lab}}=\sqrt{\frac{2}{m_1+m_2}[E_{\mathrm{lab}}V_{1T}(r_{1T}(0))V_{2T}(r_{2T}(0))]},$$
(7)
is the initial velocity for the center of mass of the projectile, and $`p_{\mathrm{rel}}=\sqrt{2\mu [Q_{\mathrm{bu}}V_{12}(r_{12}(0))]}`$ with $`\mu =m_1m_2/(m_1+m_2)`$. The value of $`r_0`$ is restricted to $`0r_0r_t`$, where $`r_t`$ satisfies $`V_{12}(r_t)=Q_{\mathrm{bu}}`$.

With these initial conditions, we solve the Newtonian equations to follow the time evolution of the trajectories for each particle, $`x_i(t),y_i(t),v_{xi}(t)`$, and $`v_{yi}(t),(i=1,2,T)`$. At each time, we monitor the relative distance between the particles. We assume that the projectile fragment $`P_i(i=1,2)`$ is absorbed by the target nucleus when the relative distance between the fragment and the target is smaller than $`r_{iT}<1.1(A_i^{1/3}+A_T^{1/3})`$ fm, and that the breakup occurs when the relative distance between $`P_1`$ and $`P_2`$ is larger than the barrier radius for the potential $`V_{12}`$. By sampling many trajectories for random values of $`r_0`$ and $`\theta `$ with the Monte Carlo method, we calculate partial cross sections for a given impact parameter $`b`$. Total cross sections are then obtained by integrating them over the value of $`b`$.

## 3 Application to 6Li + 209Bi reaction

We now apply the three-body classical trajectory Monte Carlo method to the 6Li + 209Bi system and discuss the effect of breakup of 6Li into $`d+\alpha `$ on the fusion reaction near the Coulomb barrier. We assume the Woods-Saxon form in calculating $`V_{12},V_{1T}`$, and $`V_{2T}`$, using the same values for the parameters as in Refs. . We take $`R_{\mathrm{max}}`$ = 200 fm. Fig. 1(a) – 1(c) shows sample trajectories at $`E_{\mathrm{cm}}`$ = 40 MeV. These are for $`b`$=2.5, 4.65, and 5.5 fm, respectively. For relatively small values of the impact parameter (Fig. 1(a)), the projectile nucleus 6Li is captured by the target before breaking up into $`\alpha `$ and $`d`$, leading to complete fusion. For $`b`$=5.5 fm (Fig. 1(c)), only the deuteron is captured while the alpha particle escapes from the target. This is an example of incomplete fusion. For intermediate values of $`b`$, breakup takes place near the target, but both of the breakup fragments get absorbed (Fig. 1(b)). This is breakup followed by complete fusion. Figure 2 shows the partial fusion cross sections for these processes as a function of the impact parameter, for $`E_{cm}`$=40 MeV.

Figure 3 shows fusion cross sections as a function of energy in the center of mass frame. The dotted line denotes the sum of complete and incomplete fusion cross sections. The dashed line indicates the absorption of the whole projectile from the bound states, while the solid line shows the total complete fusion cross sections. The difference between the solid and the dashed line thus represents the contribution from breakup followed by fusion of both fragments. This makes a significant contribution to complete fusion at energies above the barrier.

## 4 Summary

We present a three body classical trajectory Monte Carlo method for fusion of weakly bound nuclei. Applying this method to the 6Li + 209Bi reaction, we have demonstrated that there is a significant contribution of breakup followed by complete fusion to the total complete fusion cross sections. This process has not been considered in any coupled-channels calculations so far, and the present result suggests that a consistent definition of complete fusion is necessary when one compares experimental data with theoretical calculations.


text:
# The spin-dependent nd scattering length - a proposed high-accuracy measurementnucl-ex/0401029. To be published in Nucl. Inst. Methods A as part of the Proceedings of the 9th International Workshop on Polarized Solid Targets and Techniques in Bad Honnef (Germany), 27th – 29th October 2003.

## 1 Introduction

In the past few years, a new strategy has been developed to describe nuclear forces at low energy. Chiral perturbation theory ($`\chi `$PT) is an effective field theory, which describes interactions of pions and between pions and nucleons (N). It leads to a systematic expansion of the scattering amplitude in powers of ratios of small momenta and low-energy input parameters like the pion mass over the breakdown scale of the theory. For the first time the accuracy of calculations can be estimated in a model-independent theory of nuclear interactions, providing reliable predictions of many important low-energy quantities. These are ground state properties of bound systems and processes involving external and exchange currents, as e.g. cross sections relevant for big-bang nucleosynthesis and stellar fusion . Also the determination of fundamental properties of the neutron from experiments on few-nucleon systems mandates a model-independent subtraction of nucleon binding and meson exchange effects.

Weinberg pointed out that chiral three-nucleon (3N) forces appear naturally in $`\chi `$PT . The most relevant processes are: a two-pion exchange, a 2N contact interaction with pion exchange, and a 3N contact interaction . The contact interactions parameterise the short-range physics. As in Fermi’s theory of weak interaction, they are characterised by effective couplings, called low-energy constants (LECs). They have to be fixed by measured data of two independent low-energy 3N observables. Very recently, a first complete analysis of nd scattering at next-to-next-to leading order has been performed with impressive results . All observables are expanded in powers of momenta and the pion mass over the $`\chi `$PT-breakdown scale of about $`800`$ MeV.

An even simpler approach to the nuclear few-body problem is an effective field theory without pionic degrees of freedom . This theory starts out from point-like interactions between nucleons, which only have to respect the symmetries of QCD. Like $`\chi `$PT, it describes phenomena in a systematic way, but is applicable only at energies well below a breakdown scale set by the pion mass. Again, only two LECs characterising 3N forces are required to predict observables with an accuracy of less than $`1`$ % in processes involving three nucleons.

However, this accuracy can only be achieved if the experimental inputs are known correspondingly well. The binding energy of the triton and the doublet nd scattering length $`a_2`$ are particularly well suited to determine the LECs for $`\chi `$PT and for the pion-free theory. First, there are no Coulomb effects to be considered. Second, $`a_2`$ is very sensitive to 3N forces: in the quartet channel, the incident neutron has its spins parallel to the one bound in the deuteron, so that the Pauli principle prohibits 3N forces to play any sizeable role at low momentum transfer. In contradistinction, the $`s`$-wave in the doublet channel allows for a momentum-independent 3N interaction. This turns out even necessary to achieve results which are not sensitive to physics at high energy scales, or respectively, at short distances beyond the range of applicability of the theory. While the triton binding energy is known with an accuracy of $`5\times 10^7`$, the experimental knowledge of the nd doublet scattering length is only $`6`$ %.

## 2 Present situation and accuracy goal

The scattering length of a neutron with spin $`𝐬`$ and a nucleus with spin $`𝐈`$ is given by

$$a=\frac{I+1}{2I+1}a_++\frac{I}{2I+1}a_{}+\frac{2\left(a_+a_{}\right)}{2I+1}𝐬𝐈$$

$$=a_\mathrm{c}+\frac{2a_\mathrm{i}}{\sqrt{I\left(I+1\right)}}𝐬𝐈,$$
(1)
where $`a_+`$ and $`a_{}`$ denote the scattering lengths in the state with total spin $`I+\frac{1}{2}`$, respectively, $`I\frac{1}{2}`$. Since one cannot prepare the latter state, it is not possible to measure $`a_{}`$ directly. Experimentally accessible are the spin-independent, coherent scattering length $`a_\mathrm{c}`$, and the factor $`a_\mathrm{i}`$ which parameterises the spin-dependence (sometimes called ”incoherent” scattering length). For the neutron-deuteron system, the doublet and quartet scattering lengths, $`a_{}=a_2`$, respectively $`a_+=a_4`$, are given as the linear combinations

$$a_2=a_{\mathrm{c},\mathrm{d}}\sqrt{2}a_{\mathrm{i},\mathrm{d}},$$

$$a_4=a_{\mathrm{c},\mathrm{d}}+\frac{1}{\sqrt{2}}a_{\mathrm{i},\mathrm{d}}.$$
(2)
The best experimental value of $`a_2`$ was obtained 30 years ago , using a combination of a measurement of the scattering cross section of the free deuteron

$$\sigma _{\mathrm{s},\mathrm{d}}=4\pi \left(\left|a_{\mathrm{c},\mathrm{d}}\right|^2+\left|a_{\mathrm{i},\mathrm{d}}\right|^2\right),$$
(3)
with a gravity refractometric measurement of the bound coherent nd scattering length $`b_{\mathrm{c},\mathrm{d}}`$ (through the relation of bound and free scattering lengths $`b`$, respectively $`a`$ by

$$a=\frac{M}{M+m}b,$$
(4)
with corresponding indices. $`M`$ is the mass of the nucleus and $`m`$ the neutron mass). The experimental values were $`\sigma _{\mathrm{s},\mathrm{d}}=3.390\pm 0.012\mathrm{barn}`$ and $`b_{\mathrm{c},\mathrm{d}}=6.672\pm 0.007\mathrm{fm}`$, leading to

$$a_2=0.65\pm 0.04\mathrm{fm}.$$
(5)
Recently, a new measurement of $`b_{\mathrm{c},\mathrm{d}}`$ was performed at the NIST interferometer at Washington DC . Including the result, $`b_{\mathrm{c},\mathrm{d}}=6.6649\pm 0.0040\mathrm{fm}`$, the present world average is

$$b_{\mathrm{c},\mathrm{d}}=6.6683\pm 0.0030\mathrm{fm}.$$
(6)
However, this improvement does not significantly reduce the experimental uncertainty of $`a_2`$, since this is dominated by the unsufficient knowledge of $`a_{\mathrm{i},\mathrm{d}}`$. On the other hand, the authors of ref. argue that, because $`a_4`$ should have a very small dependence on 3N forces, $`a_{\mathrm{i},\mathrm{d}}`$ may be derived from the experimental value stated in eq.(6) and a theoretical value of $`a_4`$, using eq.(2) and eq.(4). That way, they obtain a semi-experimental value $`a_2=0.645\pm 0.003\left(\text{expt}\right)\pm 0.007\left(\text{theory}\right)\mathrm{fm}`$.

The goal of the present experiment is a direct measurement of $`a_{\mathrm{i},\mathrm{d}}`$, which does not rely on any nuclear few-body theoretical input. As a minimum aim we hope to achieve an accuracy of $`10^3`$. Using eq.(2) together with eq.(6), this shall provide a new value of $`a_2`$ with an uncertainty of $`0.004\mathrm{fm}`$.

## 3 Method

The spin-dependent scattering length induces a spin-dependence of the neutron refractive index. As a result, $`b_\mathrm{i}`$ can be determined directly with a polarised neutron beam passing through a polarised target, via detection of pseudomagnetic neutron precession around the axis of nuclear polarisation . The pseudomagnetic precession angle is given by

$$\phi ^{}=2\lambda d\underset{k}{}\sqrt{\frac{I_k}{I_k+1}}P_kN_kb_{\mathrm{i},k}.$$
(7)
$`N_k`$ is the number density, $`I_k`$ the nuclear spin and $`P_k`$ the nuclear polarisation, with the sum index $`k`$ extending over the nuclear species with spin. $`\lambda `$ is the de-Broglie wavelength of the neutrons, and $`d`$ is the thickness of the sample. The angle $`\phi ^{}`$ can be measured with an accuracy of at least one degree, using the method described in  based on Ramsey’s well-known resonance technique with two separated oscillatory fields. The target is situated in the homogeneous magnetic field between the two high-frequency $`\frac{\pi }{2}`$ coils.

The uncertainty of the deuteron polarisation $`P_\mathrm{d}`$ would impose a severe limitation in accuracy if the deuterons were the only nuclei with spin in the target. This difficulty can be considerably relaxed in a measurement of $`b_{\mathrm{i},\mathrm{d}}`$ relative to $`b_{\mathrm{i},\mathrm{p}}`$ of the proton, which is known with the high accuracy of $`2.4\times 10^4`$. Using a single target which contains both deuterons and protons at hydrogen sites, absolute polarisation measurements can be avoided . Many materials are suitable to apply the method of dynamic nuclear polarisation (DNP), by which both isotopes can be polarised simultaneously under still rather convenient conditions .

The method combines several measurements described in the following. First, the sample is polarised via DNP. After freezing the nuclear polarisation one determines the pseudomagnetic precession angle. According to eq.(7) and including an additional, instrumental phase $`\phi _0`$, it is given by

$$\varphi _1=\phi _\mathrm{d}^{}+\phi _\mathrm{p}^{}+\phi _0$$
(8)
with

$$\phi _\mathrm{d}^{}=\sqrt{2}\lambda dP_\mathrm{d}N_\mathrm{d}b_{\mathrm{i},\mathrm{d}},$$

$$\phi _\mathrm{p}^{}=\frac{2}{\sqrt{3}}\lambda dP_\mathrm{p}N_\mathrm{p}b_{\mathrm{i},\mathrm{p}}.$$
(9)
Further, using hf-saturation, one can selectively depolarise either the protons or the deuterons, without significantly affecting the polarisation of the other spin species. The subsequent cross relaxation between the two spin systems can be held sufficiently slow by keeping the temperature sufficiently low. Saturating, e.g. first the protons, one can then measure

$$\varphi _2=\phi _\mathrm{d}^{}+\phi _0.$$
(10)
After a subsequent depolarisation of the deuterons one measures

$$\varphi _3=\phi _0.$$
(11)
Combining the measured values and using eq.(9), we obtain

$$b_{\mathrm{i},\mathrm{d}}=\sqrt{\frac{2}{3}}\frac{\varphi _2\varphi _3}{\varphi _1\varphi _2}\frac{P_\mathrm{p}N_\mathrm{p}}{P_\mathrm{d}N_\mathrm{d}}b_{\mathrm{i},\mathrm{p}}.$$
(12)
The method is completed with measurements of deuteron and proton NMR signals $`_k`$, taken as integral of the corresponding rf-absorption line before and after each measurement of a pseudomagnetic precession angle:

$$_k=C_kP_kN_k.$$
(13)
Apart from natural constants, $`C_k`$ is given as the product of $`g_k^2/I_k`$ (with $`g_k`$ denoting the $`g`$-factor)  and a factor accounting for the sensitivity of the resonance circuit. The latter would be difficult to determine absolutely. However, since only the ratio of $`P_kN_k`$ for the two spin species occurs in eq.(12), replacing this by the ratio of NMR signals via eq.(13), the common instrumental factor cancels out. This requires using the same, linear resonance circuit at the same frequency in the measurements of $`_\mathrm{d}`$ and $`_\mathrm{p}`$, which therefore have to be performed at different main magnetic fields to account for the different gyromagnetic ratios.

Combining the various measurements into ratios, we thus finally obtain

$$b_{\mathrm{i},\mathrm{d}}=\frac{1}{\sqrt{6}}\frac{g_\mathrm{d}^2}{g_\mathrm{p}^2}\frac{\varphi _2\varphi _3}{\varphi _1\varphi _2}\frac{_\mathrm{p}}{_\mathrm{d}}b_{\mathrm{i},\mathrm{p}}.$$
(14)
It is the gist of this method that we donot suppose any exact knowledge of neither $`N_\mathrm{d}`$, $`N_\mathrm{p}`$, $`P_\mathrm{d}`$, $`P_\mathrm{p}`$, $`d`$ or $`\lambda `$, nor of any absolute calibration factors.

## 4 Some practical comments

The choice of the sample is governed by several factors. First, we consider the isotopic composition. Best sensitivity is attained for $`\phi _\mathrm{d}^{}\phi _\mathrm{p}^{}`$. DNP keeps the spin temperatures of protons and deuterons equal . Using the Brillouin functions in the high-temperature limit,

$$\frac{P_\mathrm{d}}{P_\mathrm{p}}\frac{4\gamma _\mathrm{d}}{3\gamma _\mathrm{p}}0.2.$$
(15)
From eq.(9),

$$\frac{N_\mathrm{d}}{N_\mathrm{p}}5\frac{P_\mathrm{p}}{P_\mathrm{d}}\frac{\phi _\mathrm{d}^{}}{\phi _\mathrm{p}^{}}.$$
(16)
Thus, $`\phi _\mathrm{d}^{}\phi _\mathrm{p}^{}`$ for $`N_\mathrm{d}/N_\mathrm{p}4`$ %. On the other hand, to keep the systematic uncertainties induced by non-linearities of the NMR resonance circuit small, requires $`_\mathrm{d}_\mathrm{p}`$, which happens to be true for about the same ratio. Considering a typical material with number density $`6.7\times 10^{22}`$ cm-3 of hydrogen sites and the measurements being done at $`\lambda =0.4`$ nm, the corresponding pseudomagnetic precession angles are

$$\phi _\mathrm{d}^{}=90\text{rad }P_\mathrm{d},\phi _\mathrm{p}^{}=440\text{rad }P_\mathrm{p}.$$
(17)
Since the Ramsey technique is sensitive to precession angles of at least one degree, a nuclear polarisation of a few percent is sufficient already for a target only $`3`$ mm long.

Relaxation times are known to be strongly temperature-dependent. Using a magnetic field of $`2.5`$ T and a cryostat providing $`T100`$ mK, the nuclear spin systems are frozen. Typical spin-lattice relaxation times of $`500`$ hours and cross relaxation time of $`100`$ hours, obtained for (CH2OH)2 doped with Cr(V) paramagnetic centres (radicals), were published in ref. . These are sufficiently long, since, for the given neutron intensity, a single run of the experiment will take much less than one hour. Operating with an amount of radicals below the optimum value to reach maximum polarisation (which is not required here), spin relaxation may be suppressed even further.

Apart from stability requirements, keeping low the amount of radicals is also of interest to suppress systematic errors which may be due to their local magnetic fields. Nuclear spins close to the paramagnetic electrons contribute to pseudomagnetic precession but usually stay undetected by NMR, since their resonance frequencies are strongly shifted with respect to the NMR line of the nuclei in the bulk. Two measures may be taken. First, the radical concentration should be kept as low as possible. Second, deuterated radicals should be used, in order to keep the few protons of the sample away from the radicals. Also for practical reasons, a good candidate of a target material is a plastic deuterated to the required abundance and doped with deuterated nitroxyl radicals . Note in addition that, using a novel stroboscopic technique of simultaneous small angle neutron scattering and NMR measurements, detailed information can be obtained about the polarisation state of protons close to radicals dissolved in a deuterated matrix . That way one can determine a correction, if necessary. Note also, that the paramagnetic electrons are always polarised very close to $`100`$ %, and their fields cause a small magnetic neutron precession, which is included in the instrumental phase $`\phi _0`$.

The NMR detection scheme is a crucial part of the setup and deserves special consideration. An analysis of several possibilities is described in a separate publication . Apart from its special role in the measurement, NMR can also be used to monitor the polarisation state of the spin systems at any stage of the experiment. Many systematic variations of parameters can be performed to check the independence of the measurement with respect to the conditions which have dropped out in the derivation of eq.(14). Different samples with different chemical composition, varied degree of deuteration and polarisation will be used. A systematic variation of target and beam parameters should be able to demonstrate the reliability of the measurement (this strategy was also adopted in a recent determination of the spin-dependent n3He scattering length ).