text:
# Cooperativity in sandpiles: statistics of bridge geometries

## 1 Introduction

The surface of a sandpile cloaks within it a vast array of complex structures – networks of grains whose stability is interconnected, surrounded by pores and necks of void space. Bridges – arch-like structures, where mutual stabilisation is a principal ingredient – are prime among these, spanning all manner of shapes and sizes through a granular medium. They can be stable for arbitrarily long times, since the Brownian motion that would dissolve them away in a liquid is absent in sandpiles – grains are simply too large for the ambient temperature to have any effect. As a result, they can affect the ensuing dynamics of the sandpile; a major mechanism of compaction is the gradual collapse of long-lived bridges in weakly vibrated granular media, resulting in the disappearance of the voids that were earlier enclosed . Bridges are also responsible for the ‘jamming’  that occurs, for example, as grains flow out of a hopper.

The difficulty of even identifying, leave alone analysing, structures as complex as bridges in a three-dimensional assembly should not be underestimated; such an algorithm now exists, and its use has resulted in the identification and classification of a panoply of bridge configurations generated via numerical simulations . However, even given the presence of such data, it is a far from obvious task to cast it in theoretical terms. In the following, we present a theory for the formation of bridges, which, among other things, is able to explain at least some of the salient features of our numerical data.

We first define a bridge. Consider a stable packing of hard spheres under gravity, in three dimensions. Each particle typically rests on three others, which stabilise it, in the sense that downward motion is impeded. A bridge is a configuration of particles in which the three-point stability conditions of two or more particles are linked; that is, two or more particles are mutually stabilised. They are thus structures which cannot be formed by the sequential placement of individual particles; they are, however, frequently formed by natural processes such as shaking and pouring, where cooperative effects arise naturally. In a typical packing, upto 70 percent of particles are involved in bridge configurations. The detailed history of a stabilisation process identifies a unique set of bridges, since the stabilising triplet for each sphere is well defined. While it is impossible to determine bridge distributions uniquely from an array of coordinates representing particle positions, we are able via our algorithm to obtain ones that are the most likely to be the result of a given stabilisation process .

We now distinguish between linear and complex bridges via a comparison of Figures 1 and 2. Figure 1 illustrates a complex bridge, i.e., a mutually stabilised cluster of five particles (shown in green), where the stability is provided by six stable base particles (shown in blue). Of course the whole is embedded in a stable network of grains within the sandpile. Also shown is the network of contacts for the particles in the bridge: we see clearly that three of the particles each have two mutual stabilisations. Figure 2 illustrates a seven particle linear bridge with nine base particles. This is an example of a linear bridge. The contact network shows that this bridge has a simpler topology than that in Figure 1. Here, all of the mutually stabilised particles are in sequence, as in a string. A linear bridge made of $`n`$ particles therefore always rests on $`n_\mathrm{b}=n+2`$ base particles, whereas the number $`n_\mathrm{b}`$ of base particles of a complex bridge of size $`n`$ varies from one bridge to another, and always obeys $`n_\mathrm{b}<n+2`$, because of the presence of loops in the contact network of complex bridges.

An important point to note is that bridges can only be formed sustainably in the presence of friction; the mutual stabilisations needed would be unstable otherwise! Although we use Monte Carlo simulations (described below) which do not contain friction explicitly, the configurations we generate correspond to those in nature that include frictional effects. In particular, they generate coordination numbers in a range that is consistent with the presence of friction . Additionally, we have analysed the configurations of molecular dynamics simulations in the limit of high friction222We thank Leo Silbert for letting us use configurations generated from his molecular dynamics simulations, where we have found and analysed the distribution of bridges.; these generate distributions of bridges very similar to our own. Our simulations suggest strong analogies between bridge structure and force distributions in granular media, which have been measured both in experiments  and in simulations . This indicates that bridges might be complementary (and experimentally easier) probes of inhomogeneities in these systems.

## 2 Simulation details

We have examined bridge structures in hard assemblies that are generated by an established, non-sequential, restructuring algorithm , whose main modelling ingredients involve stochastic grain displacements and collective relaxation from them.

This algorithm restructures a stable hard sphere deposit in three distinct stages.

* The assembly is dilated in a vertical direction (with free volume being introduced homogeneously throughout the system), and each particle is given a random horizontal displacement; this models the dilation phase of a vibrated granular medium.
* The packing is compressed in a uniaxial external field representing gravity, using a low-temperature Monte Carlo process.
* The spheres are stabilised using a steepest descent ‘drop and roll’ dynamics to find a local minimum of the potential energy.

Steps (2) and (3) model the quench phase of the vibration, where particles relax to locally stable positions in the presence of gravity. Crucially, during the third phase, the spheres are able to roll in contact with others; mutual stabilisations are thus allowed to arise, mimicking collective effects. The final configuration has a well-defined network of contacts and each sphere is supported, in its locally stable position, through point contacts, by a set of three other spheres uniquely defined. In practice, the final configuration may include a few non-stabilised particles.

The simulation method recalled above builds a sequence of static packings. Each new packing is built from its predecessor by a random process and the sequence achieves a steady state. In the steady state thus obtained, structural descriptors such as the mean packing fraction and the mean coordination number fluctuate about well-defined mean values. The steady-state mean volume fraction $`\mathrm{\Phi }`$ is found to lie typically in the range $`\mathrm{\Phi }0.55\text{ – }0.61`$. This volume fraction depends on all the parameters of the simulation, and can therefore be adjusted to any desired value in the above range. The mean coordination number is always $`Z4.6\pm 0.1`$. Since for frictional packings the minimal coordination number is $`Z=d+1`$ , this confirms that our 3D configurations correspond to those generated in the presence of friction. In fact, the mean coordination number of molecular dynamics configurations of frictional sphere packings is slightly above 4.5, in the limit of a large friction coefficient . We recall that, in the absence of friction, the coordination number is that of an isostatic system, $`Z=2d`$ , hence $`Z=6`$ in the present three-dimensional situation.

To be more specific, simulations were performed in a rectangular cell with lateral periodic boundaries, and a hard disordered base. The simulations, performed serially on a desktop workstation, have a very time-consuming Monte Carlo compression phase with each simulation taking several days’ worth of CPU time. This is necessary in order to have simulation results which are reproducible, without appreciable dependence on Monte Carlo parameters or system size. Each of our configurations includes about $`N_{\mathrm{tot}}=2200`$ particles. We have examined approximately $`100`$ configurations from the steady states of the reorganisation process for the following two values of the packing fraction: $`\mathrm{\Phi }=0.56`$ and $`\mathrm{\Phi }=0.58`$. Segregation is avoided by choosing monodisperse particles; a rough base prevents ordering. A large number of restructuring cycles is needed to reach the steady state for a given shaking amplitude; about 100 stable configurations (picked every 100 cycles in order to avoid correlation effects) are saved for future analysis. From these configurations, and following specific prescriptions, our algorithm identifies bridges as clusters of mutually stabilised particles .

Figure 3 illustrates two characteristic descriptors of bridges used in this work. Along the lines of , the main axis of a bridge is defined using a triangulation of its base particles. Triangles are constructed by choosing all possible connected triplets of base particles: the vector sum of their normals is defined to be the direction of the main axis of the bridge. The orientation angle $`\mathrm{\Theta }`$ is defined as the angle between the main axis and the $`z`$-axis. The base extension $`b`$ is defined as the radius of gyration of the base particles about the $`z`$-axis (which is thus distinct from the radius of gyration about the main axis of the bridge).

## 3 Geometrical characteristics. Size and diameter distribution

In the following, we present statistics for both linear and complex bridges; however it is the former that we analyse in greater detail, both because they are conceptually simpler and because they are in fact more numerous in our simulations. The latter clearly show that as a linear bridge gets longer, it is increasingly likely to develop branches, whereupon it becomes complex .

The formation of either type of bridge in shaken granular media is, of course, a dynamical process; we, however, adopt an ergodic viewpoint  here, inspired by the ensemble statistics of bridges generated by our simulations. As in polymer theory , we differentiate between linear and complex geometries; in our phenomenology, we regard a linear bridge as a random chain which grows in arc length $`s`$ by the sequential addition of links to its end, while complex bridges, being branched structures, cannot grow sequentially. This replacement of what is in reality a collective phenomenon by a random walk is somewhat analogous to the ‘tube model’ of linear polymers : both are simple but efficient effective pictures of very complex problems.

In this spirit, we first address the question of the length distribution of linear bridges. We define the length distribution $`f_n`$ as the probability that a linear bridge consists of exactly $`n`$ spheres. We make the simplest and the most natural assumption that a bridge of size $`n`$ remains linear with some probability $`p<1`$ if an $`(n+1)^{\mathrm{th}}`$ sphere is added to it, in the above sense. This simple effective picture therefore leads to the exponential distribution

$$f_n=(1p)p^n.$$
(3.1)
This exponential law is a common feature of all models for random ‘strings’ made of $`n`$ units, irrespective of the way they are formed. An illuminating example is provided by percolation clusters in one dimension, which are the prototype of random linear objects .

The above argument leading to an exponential distribution can alternatively be reformulated as follows, by means of a continuum approach. This will give us the opportunity to introduce a formalism that will be extensively used later. A linear bridge is now viewed as a continuous random curve or ‘string’, parametrised by the arc length $`s`$ from one of its endpoints. Within the above picture of a sequential process, we assume that a linear bridge disappears at a constant rate $`\alpha `$ per unit length, either by changing from linear to complex or by collapsing. The probability $`S(s)`$ that a given bridge survives at least up to length $`s`$ thus obeys the rate equation $`\dot{S}=\alpha S`$, and therefore falls off exponentially, according to $`S(s)=\mathrm{exp}(\alpha s)`$. The probability distribution of the length $`s`$ of linear bridges therefore reads $`f(s)=\dot{S}(s)=\alpha \mathrm{exp}(\alpha s)`$. This is the continuum analogue of (3.1).

Figure 4 shows a logarithmic plot of numerical data for the length distribution $`f_n`$ of linear bridges. The data exhibit an exponential fall off of the form (3.1), with $`p0.37`$, i.e., $`\alpha 0.99`$333Despite its proximity to unity, this value of $`\alpha `$ has no special significance, as can be seen from the corresponding value of $`p`$.. This exponential decay of the distribution of linear bridges is very accurately observed until $`n12`$.

Complex bridges begin to be dominant around $`n8`$. For a complex bridge, the number $`n`$ is now referred to as the ‘size’ of the bridge. As complex bridges are branched objects, it is natural to expect that their size distribution resembles that of generic randomly branched objects, such as e.g. lattice animals or critical percolation clusters . These prototypical examples have a power-law fall-off in their size distribution:

$$f_nn^\tau .$$
(3.2)
The numerical data shown in Figure 5 clearly show that the size distribution of bridges obeys a power-law behaviour of the form (3.2), and thus behave as generic randomly branched objects. The measured value of the exponent $`\tau `$ seems to coincide with the limiting value $`\tau =2`$, at and below which the mean size $`n=nf_n`$ is divergent.

We now turn to the diameter of linear and complex bridges, which is one of the most important of their geometrical characteristics. The typical diameter $`R_n`$ of a bridge of size $`n`$ is such that $`R_n^2`$ is the mean squared end-to-end distance over all the bridges of size $`n`$. We examine this quantity separately for linear bridges, as well as the ensemble of all bridges. In both cases we expect and find a power-law behaviour

$$R_nn^\nu .$$
(3.3)
For linear bridges, $`D_{\mathrm{lin}}=1/\nu _{\mathrm{lin}}`$ can be interpreted as their fractal dimension. The observed value (see Figure 6) $`\nu _{\mathrm{lin}}0.66`$, i.e., $`D_{\mathrm{lin}}1.51`$, lies between those for a self-avoiding walk in two ($`\nu _2=3/4`$) and three ($`\nu _30.59`$) dimensions . This seems entirely reasonable, since linear bridges are likely to behave as typical self-avoiding curves. Furthermore, in three dimensions, their effective dimensionality is expected to lie between 2 and 3: they may start off confined to a plane, and then collapse onto each other, due to the effects of vibration.

For complex bridges, we again find a power law of the form (3.3). The geometrical exponent $`\nu _{\mathrm{all}}0.74`$ (see Figure 7) is now naturally compared with the value $`\nu 0.875`$ of 3D percolation . (Although the figure shows statistics for ‘all’ bridges, these mostly reflect the behaviour of complex bridges, which dominate beyond about $`n8`$.)

Lastly, we plot in Figure 8 the probability distribution function of the base extension $`b`$, for linear bridges; this is clearly a measure of the spanning, and hence the jamming, potential of a bridge (see Figure 3). We plot also in Figure 8 distribution functions that are conditional on the bridge size $`n`$. This was done to compare our results with those of , bearing in mind of course that our results refer to sphere packings in three dimensions whereas theirs refer to disk packings in two. In our case, the conditional distributions are sharply peaked, and only extend over a finite range of $`b`$; in theirs the peaks are smaller and the distributions broader. Our results indicate that at least in three dimensions, and upto the lengths of linear bridges probed, bridges of a given length $`n`$ have a fairly characteristic horizontal extension; we can predict fairly reliably the dimension of the orifice that they might be expected to jam. The results of  indicate, on the contrary, that there is a wide range of orifices which would be jammed by disk packings of a given length, in two dimensions. Also, for our three-dimensional bridges, the cumulative distribution has a long tail at large extensions, reflecting chiefly the existence of larger bridges than in the data of . (We note that this long tail is characteristic of three-dimensional experiments on force chains  which we will discuss at greater length below.) All of this seems entirely reasonable given the more complex topology of three-dimensional space, as well as the greater diversity of possible sphere packings within it.

In the lower part of Figure 8 we have plotted the logarithm of the normalised probability distribution of base extensions, against the normalised variable $`b/b`$, where $`b`$ is the mean extension of bridge bases. This figure emphasises the exponential tail of the distribution function, and also shows that bridges with small base extensions are unfavoured. It is interesting to note that the form of the normalised distribution resembles that found in force distributions of emulsions , and that there are also similarities with normal force distributions in MD simulations of particle packings (see e.g. ). In our simulations, the sharp drop at the origin with a very tangible peak strongly resembles that obtained in the latter case, in the limit of strong deformations . Realising that the measured forces  propagate through chains of particles, we use this similarity to suggest that our bridges are really just long-lived force chains, which have survived despite strong deformations. We suggest also that with the current availability of 3D visualisation techniques such as NMR , bridge configurations might be an easily measurable, and effective, tool to probe inhomogeneous force networks in shaken sand.

## 4 Orientational distribution of linear bridges. Theory and simulation

Linear bridges predominate at small $`n`$, which is after all the regime of interest for finite-sized simulations and experiments. The theory we present below focuses on such bridges; the distribution of their orientation is predicted, and then compared with our simulation results.

As before, a linear bridge is modelled as a continuous curve, parametrised by the arc length $`s`$ from one of its endpoints. We propose the following effective description of the local and global orientation of this curve. For the sake of simplicity, we choose to focus on the most important degree of freedom, which is the tilt with respect to the horizontal; the azimuthal degree of freedom is therefore neglected. Accordingly, we define the local or ‘link’ angle $`\theta (s)`$ between the direction of the tangent to the bridge at point $`s`$ and the horizontal, and the mean angle made by the bridge from its origin up to point $`s`$,

$$\mathrm{\Theta }(s)=\frac{1}{s}_0^s\theta (u)du.$$
(4.1)
The local angle $`\theta (s)`$ so defined may be either positive and negative; it can even change sign along the random curve which represents a linear bridge. Of course, the orientation angle $`\mathrm{\Theta }`$ measured in our numerical simulations pertains to the three-dimensional realm of real sphere packings. It is there positive by construction, being defined at the end of Section 2 as the angle between the main bridge axis and the $`z`$-axis.

Our simulations show that the angle $`\mathrm{\Theta }(s)`$ typically becomes smaller and smaller as the length $`s`$ of the bridge increases. Small linear bridges are almost never flat ; as they get longer, assuming that they still stay linear, they get ‘weighed down’, arching over as at the mouth of a hopper . Thus, in addition to our earlier claim that long linear bridges are rare, we claim further here that (if and) when they exist, they typically have flat bases, becoming ‘domes’.

We use these insights to write down equations to investigate the angular distribution of linear bridges. These couple the evolution of the local angle $`\theta (s)`$ with local density fluctuations $`\varphi (s)`$ at point $`s`$:

$`\dot{\theta }=a\theta b\varphi ^2+\mathrm{\Delta }_1\eta _1(s),`$ (4.2)
$`\dot{\varphi }=c\varphi +\mathrm{\Delta }_2\eta _2(s).`$ (4.3)
The effects of vibration on each of $`\theta `$ and $`\varphi `$ are represented by two independent white noises $`\eta _1(s)`$, $`\eta _2(s)`$, such that

$$\eta _i(s)\eta _j(s^{})=2\delta _{ij}\delta (ss^{}),$$
(4.4)
whereas the parameters $`a`$, …, $`\mathrm{\Delta }_2`$ are assumed to be constant.

The phenomenology behind the above equations is the following: the evolution of $`\theta (s)`$ is caused by the addition to the bridge of single particles to the bridge, along the lines of our effective sequential picture. The motion of particles within their cages gives rise to the fluctuations of local density at a point $`s`$; thus $`\varphi `$ can be regarded as an effective collective coordinate, with $`\theta `$ seen as an independent-particle coordinate .

The individual dynamics are simple: the first terms on the right-hand side of (4.2), (4.3) say that neither $`\theta `$ nor $`\varphi `$ is allowed to be arbitrarily large. Their coupling via the second term in (4.2) arises as follows: large link angles $`\theta `$ will be increasingly unstable in the presence of density fluctuations $`\varphi ^2`$ of large magnitude, which would to a first approximation ‘weigh the bridge down’, i.e., decrease the angle $`\theta `$ locally. Density fluctuations (otherwise known as ‘dilatancy’) also play a key role in the dynamics of the angle of repose of a sandpile, which are described in a concurrent paper . There, sandpile collapse is described in terms of an activated process, where an effective temperature competes against configurational barriers generated by dilatancy.

Reasoning as above, we therefore anticipate that for low-intensity vibrations and stable bridges, both density fluctuations $`\varphi (s)`$ and link angles $`\theta (s)`$ will be small. Accordingly, we linearise (4.2), obtaining thus an Ornstein-Uhlenbeck equation

$$\dot{\theta }=a\theta +\mathrm{\Delta }_1\eta _1(s).$$
(4.5)
Let us make the additional assumption that the initial angle $`\theta _0`$, i.e., that observed for very small bridges, is itself Gaussian with variance $`\sigma _0^2=\theta _0^2`$. The angle $`\theta (s)`$ is then a Gaussian process with zero mean for any value of the length $`s`$. Its correlation function can be easily evaluated to be

$$\theta (s)\theta (s^{})=\sigma _{\mathrm{eq}}^2\mathrm{e}^{a|ss^{}|}+(\sigma _0^2\sigma _{\mathrm{eq}}^2)\mathrm{e}^{a(s+s^{})}.$$
(4.6)
The characteristic length for the decay of orientation correlations is $`\xi =1/a`$, and the variance of the link angle reads

$$\theta ^2(s)=\sigma _{\mathrm{eq}}^2+(\sigma _0^2\sigma _{\mathrm{eq}}^2)\mathrm{e}^{2as},$$
(4.7)
where

$$\sigma _{\mathrm{eq}}^2=\frac{\mathrm{\Delta }_1^2}{a}$$
(4.8)
is the ‘equilibrium’ value of the variance of the link angle. Thus as the chain gets longer, the variance of the link angle relaxes from $`\sigma _0^2`$ (that for the initial link) to $`\sigma _{\mathrm{eq}}^2`$, in the limit of an infinite bridge.

We predict from the above that the mean angle $`\mathrm{\Theta }(s)`$ will have a Gaussian distribution, for any fixed length $`s`$. By inserting (4.6) into (4.1), we derive its variance:

$$\mathrm{\Theta }^2(s)=2\sigma _{\mathrm{eq}}^2\frac{as1+\mathrm{e}^{as}}{a^2s^2}+(\sigma _0^2\sigma _{\mathrm{eq}}^2)\frac{(1\mathrm{e}^{as})^2}{a^2s^2}.$$
(4.9)
The asymptotic result

$$\mathrm{\Theta }^2(s)\frac{2\sigma _{\mathrm{eq}}^2}{as}\frac{2\mathrm{\Delta }_1^2}{a^2s}$$
(4.10)
confirms our earlier statement (see below (4.1)) that a typical long bridge has a base that is almost flat. It can be viewed as consisting of a large number $`as=s/\xi 1`$ of independent ‘blobs’444This is qualitatively reminiscent of a similar problem in polymers, where the de Gennes picture of ‘blob’ dynamics in dilute solutions gives way to rigid rod-like behaviour – see e.g. ., each of length $`\xi `$.

The result (4.10) has another interpretation. As $`\mathrm{\Theta }(s)`$ is small with high probability for a very long bridge, its extension in the vertical direction reads approximately

$$Z=z(s)z(0)s\mathrm{\Theta }(s),$$
(4.11)
so that $`Z^2s^2\mathrm{\Theta }^2(s)2(\mathrm{\Delta }_1/a)^2s`$. Going back to the discrete formalism, we have therefore

$$Z_nn^{1/2}.$$
(4.12)
The vertical extension of a linear bridge is thus found to grow with the usual random-walk exponent $`1/2`$, whereas its horizontal extension exhibits the non-trivial exponent $`\nu _{\mathrm{lin}}0.66`$ of Figure 6. Thus, long linear bridges are domelike; they are vertically diffusive but horizontally superdiffusive. We recall that the two-dimensional arches found in  were diffusive in a vertical direction; our present results show that in three dimensions, the vertical diffusivities of bridge structures remain, and are enriched by a horizontal superdiffusivity. Evidently, jamming in a three-dimensional hopper would be caused by the planar projection of such a dome.

We now compare these predictions with data from our simulations. The numerical data shown in Figure 9 confirm that the mean angle has, to a good approximation, a Gaussian distribution in this typical situation of bridges of size $`n=4`$. We recall that the angle $`\mathrm{\Theta }`$ measured in simulations pertains to the three-dimensional world. Hence only the positive half of the Gaussian is used for comparison with results of simulations, and the $`\mathrm{sin}\mathrm{\Theta }`$ three-dimensional Jacobian is divided out. It has also been checked to high accuracy (less than 2 percent for a bin size $`\delta \phi =5`$ degrees) that the azimuthal angle $`\phi `$ of the main bridge axis has the expected flat distribution. Figure 10 shows the measured size dependence of the variance $`\mathrm{\Theta }^2(s)`$, for both volume fractions. The numerical data are found to agree well with a common fit to the first (stationary) term of (4.9), with a common value of the parameters $`\sigma _{\mathrm{eq}}^2=0.093`$ and $`a=0.55`$. The ‘transient’ effects of the second term of (4.9) are invisible with the present accuracy (see below). This observed agreement therefore provides a first confirmation of the validity of our theory, quite remarkable given that it is entirely independent of the simulations. We conclude that in spite of its simplicity, our theory captures the main features of linear bridges.

We end this section with the following remarks. First, more subtle effects, including the effects of transients via the second term of (4.9), and the dependence of the parameters $`\sigma _{\mathrm{eq}}^2`$ and $`a`$ on the packing fraction $`\mathrm{\Phi }`$, are deserving of further investigation. As a preface to our second remark, we recall that in our phenomenological picture linear bridge formation proceeds sequentially, while complex bridges with their branched structures form collectively, in space.555We reiterate once again that this effective picture replaces the complexity of the full dynamical picture, where all bridges form cooperatively. We might expect that with increasing density $`\mathrm{\Phi }`$, branched structures would become more and more common; linear bridge formation, with its ‘sequential’ progressive attachment of independent blobs would then become more and more rare. Our theory should therefore cease to hold at a limit packing fraction $`\mathrm{\Phi }_{\mathrm{lim}}`$, which is reminiscent of the single-particle relaxation threshold density , at which the dynamics of granular compaction crosses over from being single-particle (sequential) to collective. Finally, while we have allowed for the local variability of density fluctuations $`\varphi `$ via (4.3), we have so far not explicitly coupled these to bridge orientations. All of these issues form part of ongoing work, since much more numerical data is needed before a detailed comparison with theory is possible.

## 5 Discussion

In the above, we have classified bridge structures in granular media as either linear or complex, analysed their structures as obtained in our simulations, and obtained values for some of their more important geometrical exponents. We have also presented a theory for their orientational correlations, which is in reasonably good agreement with our simulations, and more importantly, makes predictions that can be observed experimentally. Do linear bridges really predominate at short lengths? Do they really cross over to being domelike at large lengths, if they survive? The answers to these and other questions, if obtained experimentally, will not just satisfy what is arguably mere scientific curiosity. Our investigations above suggest that long-lived bridges are natural indicators of sustained inhomogeneities in granular systems. Consequently, experimental and theoretical explorations of geometry and dynamics in stable bridge networks will have immediate implications for issues such as hysteresis and cooperative stability, which are some of the most visible manifestations of complexity in granular media.

AM warmly thanks SMC-INFM (Research and Development Center for Statistical Mechanics and Complexity, Rome, Italy), and the Service de Physique Théorique, CEA Saclay, where parts of this work were done. GCB acknowledges support from the Biotechnology and Biological Sciences Research Council UK and thanks Luis Pugnaloni for help with the figures. Kirone Mallick is gratefully acknowledged for fruitful discussions.

## References


format: text

text:
# Short-Bright GRBs: spectral properties

## 1 Introduction

Evidences supporting the possible different nature of short and long GRBs are: (1) their bimodal duration distribution (Kouveliotou1993 ) with mean duration of $`20`$ sec and $`0.3`$ sec for long and short events, respectively; (2) their different temporal properties (e.g. number and width of pulses in the light curve, Norris2000 , Nakar2002 , McBreen2001 ); (3) the larger hardness ratio of short bursts (Kouveliotou1993 , Tavani1998 , Paciesas2001 ). Moreover, theoretical models for the progenitors of GRBs associate short bursts to the merger of compact objects in a binary system (Goodman1986 , Meszaros1992 ) while long GRBs seem to be connected to the core–collapse of massive stars (Woosley1993 , Vietri1998 ).

Despite the increasing understanding of the nature of long duration $`\gamma `$–ray bursts, the population ($`1/3`$) of sub-second short GRBs is still largely not understood. This is mainly due to (i) the low signal-to-noise ratio which strongly limits the analysis of the prompt (temporal and spectral) emission of short events, and (ii) the lack of any firm afterglow measurement for short GRBs, which indeed represented a major advance in unveiling the mystery of long bursts.

## 2 Short vs. Long: spectral properties

The emission properties of long GRBs have been studied in details (Band1993 , Preece2000 ) by fitting their time average (Band1993 ) and time resolved (Preece2000 , Ghirlanda2002 ) broad-band high-resolution spectra. Nonetheless, the comparison with the spectral properties of short bursts has been based mainly on the analysis of their fluence hardness ratio (e.g. Cline1999 , Yi-ping Qin2001 ) which is marginally representative of the effective spectral shape. For this reason we analyzed (Ghirlanda2003a ) a sample of bright–short BATSE bursts (selected with peak flux $`10`$ phot/cm2 sec between 50-300 keV) by fitting their $`\gamma `$–ray spectra ($`301800`$ keV) with the standard spectral functions (e.g. Preece2000 ). The set of spectral parameters (namely the low energy photon spectral index $`\alpha `$ and the EFE peak energy $`E_{peak}`$) were compared with those of long bright bursts (Preece2000 ) both considering time integrated and time resolved spectra.

As shown in fig.1 (Left) short and long GRBs present different $`\alpha `$ distributions and only marginally different $`E_{peak}`$ distributions (with a small Kolmogorov-Smirnov test probability - i.e. 4% and 10%, respectively - for the two population to be similar). The average values are $`\alpha =0.84\pm 0.15`$, $`E_{peak}=305\pm 22`$ keV and $`\alpha =0.58\pm 0.10`$, $`E_{peak}=355\pm 30`$ keV for long and short GRBs, respectively. This suggests that short bursts have (average) harder spectra than long events, as also indicated by the analysis of their hardness ratio, *but*, intriguingly, this difference is due to a considerably *harder low energy spectral component* present in short bursts rather than to a different peak energy.

The comparison of short GRBs with the time resolved spectral properties of long bursts shows that short events have spectra similar to the first 1 sec of long GRBs, as shown in fig.1 (Right), with a high KS probability of 23% and 80% for $`\alpha `$ and $`E_{peak}`$ to be similar. This might suggest that a similar mechanism operates for the complete duration of short bursts and in the first 1 sec of long GRBs.

## 3 Short vs. Long: energy and luminosity

The analysis of the intrinsic properties of (still few) long bursts with known redshift (Amati2002 ) highlighted a possible correlation ($`E_{iso}E_{peak}^{1.93}`$) between is the equivalent isotropic burst energy $`E_{iso}`$ and the spectral peak energy $`E_{peak}`$. This correlation has been also confirmed by Hete–II (Lamb2003 ) and a similar relation has been found between $`E_{peak}`$ and the burst isotropic luminosity $`L_{iso}`$ (Yonetoku2003 ).

If short bursts were similar to long events they might satisfy these correlations between $`E_{iso}`$ ($`L_{iso}`$) and $`E_{peak}`$. However, no redshift of short GRBs has been measured so far. Nonetheless, we can still verify (Ghirlanda2003a ) the above hypothesis assuming a variable redshift (between 0.1 and 10) and scaling the observed spectral properties of the sample of bright short bursts in the source rest frame.

Fig.2-a shows that for any assumed $`z`$ short bursts (solid curves) populate a region below the $`E_{peak}`$ \- $`E_{iso}`$ correlation of long bursts (Amati2002  \- solid line). On the other hand (fig.2-b) the luminosity of short events is consistent with the proposed relation for long bursts (Yonetoku2003  \- solid line). In conclusion short and long bursts seem to emit a *similar equivalent isotropic luminosity* but different (lower in short bursts) energy due to their different duration. Nonetheless, short and long bursts might still have similar emitted energies if short bursts are less collimated than long events although, in this case, short bursts would have a much larger luminosity. Under the first hypothesis we can further extract from the $`L_{iso}`$-$`E_{peak}`$ relation the possible redshift distribution of short bursts (fig.2-right, solid line). This is consistent with the observed $`z`$ distribution of long events (fig.2-right, dotted line) and again supports a possible similarity of short and long GRBs.

## 4 Discussion

The comparison of the spectral properties of short and long bursts pointed out that short bursts are harder than long events due to a (average) harder low energy spectral component (rather than to a different peak energy). The spectra of short bursts are similar to the first 1 sec of the emission of long bursts. Intriguingly, short bursts present a similar intrinsic luminosity but a lower (isotropic) energy than long bursts if the recently found correlations between these quantities were true also for short GRBs. Under such hypothesis the implied redshift distribution of short bursts results similar to that observed in long GRBs, and this prediction will be tested in the forthcoming Swift era. Nonetheless, short GRBs might still have a similar energy to that of long bursts if their collimation angle is much larger than that of long events, but in this case their intrinsic luminosity should be much larger than that of long bursts. These results suggest the presence of a common mechanism operating at the beginning of short and long bursts which could explain their similar spectral properties and luminosity. If this is the case a possible difference in the burst dynamical evolution (e.g. fallback of the pre-GRB ejected material) might play a crucial role in distinguishing these two classes.

G. Ghirlanda would like to thank the organizing committee for the full grant support for the participation to this conference.


text:
# Novel pinning phenomena in a superconducting film with a square lattice of artificial pinning centers

## Abstract

We study the transport properties of a superconducting Nb film with a square lattice of artificial pinning centers (APCs) as a function of dc current, at a temperature close to the superconducting transition temperature of the film. We find that, at low dc currents, the differential resistance of the film shows the standard matching field anomaly, that is, the differential resistance has a local minimum at magnetic fields corresponding to an integer number of flux lines per APC. However, at higher dc currents, the differential resistance at each matching field turns to a local maximum, which is exactly opposite to the low current behavior. This novel effect might indicate that the flux lines in the APC system change their flow mode as the dc current is increased.

Superconducting films with a lattice of artificial pinning centers (APCs) have been investigated for several years fiory ; victor1 ; vitali ; pannetier ; reichhardt1 . It is found that the vortex lattice is highly ordered at certain matching fields, $`H_n`$, where the number of flux lines equals an integer multiple $`n`$ of the number of APCs. This ordering effect provides additional pinning interactions that result in a local minimum resistance and maximum critical current at the matching fields. Interest has been growing in such systems in the presence of driving forces reichhardt2 . Although much work has been done in systems with random disorder vinokur ; giamachi , the behavior of APC systems is only beginning to be explored. In this letter, we investigate the transport properties of a superconducting Nb film with a square lattice of APCs as a function of dc current through the film. Our results show novel behavior in a specific range of the dc current, which has not been observed in previous experiments.

Figure 1(a) shows a schematic of the device discussed in this letter. This device was patterned onto an oxidized Si substrate using laser interferometric lithography. The entire device is made of a single layer of Nb film (100 nm thick) with a square lattice of APCs with period $`d`$$``$1 $`\mu `$m and diameter $`D`$$``$0.3 $`\mu `$m. A scanning electron microscopy (SEM) image of the Nb film is shown in Fig. 1(b) with a 2 $`\mu `$m size bar. Technical details of the patterning procedure can be found in Ref. .

Figure 2 shows the differential resistance of the Nb film as a function of magnetic field normal to the film with different dc bias currents ($`I_{dc}`$) through the sample. The resistance was measured by conventional four-terminal techniques using an ac resistance bridge. The temperature of the system was fixed at $`T`$$`=`$7.07 K, close to the superconducting transition temperature, $`T_c`$$`=`$7.13 K, of the film. All the resistance curves show oscillations as a function of magnetic field, with a fundamental period of $``$ 18 G, corresponding to one superconducting flux quantum $`h/2e`$ per unit cell of the square lattice, or alternately, one flux line per APC. At low dc bias currents ($`I_{dc}`$$``$100 $`\mu `$A), the differential resistance curve has a local minimum at magnetic fields $`H_n`$ corresponding to an integer number $`n`$ of flux lines per APC. This behavior is the one that is conventionally observed. It is associated with strong pinning of flux lines at the matching fields $`H_n`$, leading to lower dissipation associated with flux line motion, and hence a lower resistance. At higher dc bias currents ($`I_{dc}`$$``$110 $`\mu `$A), however, the local minimum becomes a maximum at each matching field, exactly opposite to the results at the low dc current. The current at which this change occurs is between 100 $`\mu `$A and 110 $`\mu `$A, where the behavior of the differential magnetoresistance shows a dramatic change in a small current range of 10 $`\mu `$A. At even higher dc bias currents ($`I_{dc}`$$`>`$160 $`\mu `$A), not plotted out in Fig. 2, the magnetoresistance oscillations are washed out.

The novel behavior observed in Fig. 2 can also be seen in Fig. 3(a), the differential resistance measurements as a function of dc current at two different magnetic fields of 18 G and 9 G, corresponding to the first matching field $`H_1`$, and $`H_1/2`$, respectively. This measurement has been done by the same four-terminal technique used for Fig. 2, and at the same temperature $`T`$$`=`$7.07 K. Figure 3 illustrates that the $`dV/dI\text{vs}.I_{dc}`$ behavior can be divided into three regimes: I, II and III yaron . In regime I, the measurements show zero differential resistance, which indicates that the vortex lattice is static due to the strong pinning force. In regime II, the differential resistance increases gradually when we increase the dc current. This regime is called the plastic flow regime, where the vortex lattice has been partly disordered and shows plastic flow behavior. In regime III, the measurements show a sharp resistance jump to the normal state of the film, which indicates that the vortex lattice has been completely disordered and shows avalanche motion. In our measurements, there exist in addition three sub-regimes in the plastic flow regime, which are labeled II(1), II(2) and II(3) in Fig. 3. In II(1) the measurements show the standard matching field anomaly, i.e., the 18 G curve, where each flux line is expected to be tied to one APC, has a smaller resistance compared with the 9 G curve. In II(2), the differential resistance measured at 18 G is larger than the one measured at 9 G, which is opposite to the results in regime II(1). In II(3), the 18 G curve and the 9 G curve merge together. No oscillations appear in this sub-regime. It should also be pointed out that the crossing point between regime II(1) and II(2) is quite sensitive to the system temperature. As shown in Fig. 3(b), this crossing point moves from 103 $`\mu `$A to 152 $`\mu `$A when we cool our system down just 0.01 K, from 7.07 K to 7.06 K. In other words, the novel behavior in the regime II(2) only happens in a very small temperature range, which might be the reason why it has not been observed in previous work.

To our knowledge, no current theory can explain this unusual behavior. We propose here a phenomenological picture to explain our experimental results. The total force f acting on each flux line in a superconducting film with APCs can be expressed as $`\text{f}=\text{f}^{vp}+\text{f}^{vv}+\text{f}_T+\text{f}_d`$, where $`\text{f}^{vp}`$ is the pinning force due to the trapping of the flux line by a APC, $`\text{f}^{vv}`$ is the repulsive interaction with other flux lines, $`\text{f}_T`$ the force associated with thermal fluctuations, and $`\text{f}_d`$ the Lorentz force, present when transport current is sent though the film reichhardt1 ; reichhardt2 . The first two terms in f would result in a highly ordered vortex lattice configuration, especially at each matching field; the last two terms drive the system from the highly ordered vortex lattice state to a disordered state. The motion of flux lines is determined by the balance between $`\text{f}^{vp}+\text{f}^{vv}`$ and $`\text{f}_T+\text{f}_d`$. In our experiments, we fix the system temperature at a value slightly below the superconducting transition temperature of the Nb film, that is, we set $`\text{f}_T`$ to be constant, and vary $`\text{f}_d`$ by tuning the dc current through the device.

The behavior demonstrated in Fig. 3 might be associated with a change of the flow mode of the flux lines when the dc current is increased. At very low dc currents (regime I), the pinning force dominates the system, $`(\text{f}^{vp}+\text{f}^{vv})`$$``$$`(\text{f}_T+\text{f}_d)`$. As a consequence, the flux lines are static and the measurements show zero differential resistance. This corresponds to the situation shown in Fig. 4(a), which shows a schematic of the potential landscape seen by the vortex lines at $`H_1`$. There is a minimum of potential at each APC that pins exactly one vortex line. Any vortex line that is displaced from its pinning center by thermal fluctuations is immediately forced back due to the repulsive interactions between vortices, leading to a minimum in the resistance due to vortex motion. For half-filling (corresponding to $`H_1/2`$, the situation shown in Fig. 4(b)), the distance between the vortices is larger, leading to a smaller vortex-vortex interaction. Each vortex is therefore not so rigidly pinned as in Fig. 4(a), giving rise to a larger resistance due to vortex motion.

The effect of applying a dc current is to add a finite slope to the potential profile, as shown in Figs. 4(c)-(f). The larger the current, the steeper the slope. At low $`I_{dc}`$, (corresponding to region II(1)), the pinning potential is still strong enough that the situation described above for Figs. 4(a) and (b) still holds, and the resistance at $`H_1`$ (corresponding to Fig. 4(c)) is less than the resistance at $`H_1/2`$ (Fig. 4(d)). At still higher $`I_{dc}`$ (corresponding to region II(2)), however, the situation changes. Consider the case of full filling, at $`H_1`$, shown schematically in Fig. 4(e). Here, the effective pinning potential at one end of a pinning center is effectively reduced by the slope of the potential induced by the dc field, so that thermal fluctuations can easily kick a flux line out of the potential well. If this happens in the center of the superconducting film, the repulsive interactions between vortices are still likely to push the vortex line back into its pinning site. If the vortex pops out of a pinning site at the edge of the film, however, it will leave the film altogether. At this point, the flux line at the pinning site immediately adjacent to the now empty site will no longer feel the repulsive interaction due to the absent vortex line. If the repulsive interaction from its remaining vortex neighbors is sufficient to overcome the reduced pinning potential, it will jump to the empty pinning center. This jump mechanism will be continued down the linear chain of vortices, leading to a mass movement of vortices in the lattice, and hence a higher resistance. For half-filling ($`H_1/2`$, shown in Fig. 4(f)), the smaller repulsive interaction between vortices, arising from the greater inter-vortex distance, is not sufficient to overcome the effective pinning potential, so that the flux lines do not move. Hence the resistance of the sample at $`H_1/2`$ is lower than its resistance at $`H_1`$. This mechanism clearly depends on a fine balance between the dynamic pinning potential, which depends on $`I_{dc}`$, and the repulsive vortex-vortex interactions, and hence can only be seen in a narrow range of $`I_{dc}`$. At even higher currents, corresponding to region II(3), the potential profile is tilted so much by the applied dc current that the vortex lattice flows at both $`H_1/2`$ and $`H_1`$, so that no clear oscillations are seen in the magnetoresistance. A similar picture also holds at higher magnetic field, where there is more than one vortex per unit cell. One might also point out that vortex-antivortex pairs are present even at zero field, as evidenced by the presence of a finite resistance.

An alternative explanation is based on the remarkable similarity between the data shown in Fig. 2, and the oscillations in conductance observed in an Aharonov-Bohm interferometer formed from a quantum dot yacoby . In that experiment, a similar change in the phase of the oscillations is observed as the phase difference between the two arms of the interferometer is varied by means of an external gate potential. Indeed, in that experiment, the oscillations at the fundamental period $`\mathrm{\Phi }_0`$ disappear at the cross-over point, and only oscillations at $`\mathrm{\Phi }_0/2`$ remain. Hints of this behavior are also seen in our data (see curve (d) in Fig. 2). This suggests that the transport current in our system is modulating the phase of the interference around each antidot in the lattice. The exact nature of this modulation is not clear to us at the moment, and further investigation is required to understand the origin of this unusual behavior.

This work is supported by the NSF under grant number DMR-0201530.

Fig. 1, Z. Jiang *et al.*

Fig. 2, Z. Jiang *et al.*

Fig. 3, Z. Jiang *et al.*

Fig. 4, Z. Jiang *et al.*


text:
# Eclipsing binaries in open clusters. I. V615 Per and V618 Per in h Persei

## 1. Eclipsing binaries in open clusters

Theoretical stellar models are increasingly sophisticated and accurate. Different input physics and different sets of models often produce very similar predictions for the physical parameters of stars. This means that observational tests of stellar model physics require high quality data to determine the success of a particular model compared to other models. This is exacerbated by the number of physical properties which are not observable, allowing several quantities to be adjusted freely to fit the observed data, for example interior helium and metal abundances and the effect of different theories and treatments of convection.

Detached eclipsing binaries (dEBs) can can be analysed to provide accurate and absolute values for the masses and radii of two stars of the same age and chemical composition (Andersen 1991), but the age and composition are not constrained by observations so can be varied freely when attempting to fit models to data. Photometric studies of open clusters are also commonly used to test theoretical predictions. Analyses of colour-magnitude diagrams can test quite subtle physics, but this method is usually limited in its usefulness due to degeneracies between the distance to the cluster, its age, metallicity and reddening effects, and the presence of unresolved multiple stars and field stars.

Eclipsing binaries in open clusters provide a way of combining these two methods to increase the number of physical properties which are directly observed. Accurate masses and radii of two stars are augmented by some knowledge of their age and chemical composition due to their membership of the cluster. A further advantage of having accurate radii is that the distance to each component of the dEB can be found via empirical calibrations of surface brightness (e.g., Barnes & Evans 1976). This allows the determination of the distance to a cluster entirely independently of main sequence fitting or theoretical models.

The dEBs V615 Persei and V618 Persei were discovered by Krzesiński, Pigulski & Kołaczkowski (1999) in the well-studied open cluster h Persei. Their spectral types are B7 V and A2 V and their periods are 13.7 and 6.4 days respectively. They have proper motions (Uribe et al. 2002), systemic velocities and photometric characteristics consistent with membership of the cluster.

The Perseus Double Cluster is composed of the open clusters h Per and $`\chi `$ Per, and contains a large number of unevolved early-type stars. The question of whether the clusters are physically related has been studied frequently throughout the last century. The results of recent photometric studies are converging to identical values of age ($`\mathrm{log}\tau =7.0\pm 0.1`$ years), and distance modulus ($`11.70\pm 0.05`$ mag), but final agreement has not yet been reached on whether the clusters are co-evolutionary (see Southworth, Maxted & Smalley 2004).

## 2. Observations and analysis

Grating spectra with a resolution of 0.2 Å were obtained with the Isaac Newton Telescope on La Palma. The spectral window included Mg ii 4481 Å, known to give good radial velocities for early-type stars, and He i 4471 Å, useful for the determination of effective temperatures. Exposure times of 1800 s gave an approximate signal to noise per pixel of 50 for V615 Per and 15 for V618 Per.

The observed and disentangled spectra were analysed with synthetic spectra calculated using uclsyn (Smalley et al. 2001) and Kurucz atlas9 model atmospheres. Effective temperatures, rotational velocities, and light ratios were derived by fitting observations using least squares. All stars are slow rotators, with only the most massive star here rotating faster than synchronous.

Radial velocities were derived using the two-dimensional cross-correlation technique todcor (Zucker & Mazeh 1994) and synthetic templates generated using uclsyn. Circular orbits were fitted to the derived radial velocities using the sbop code (Fig. 1). Uncertainties due to synthetic template mismatch were found using templates with a large range of effective temperature, rotational velocity and microturbulence velocity, and added in quadrature to the formal errors of the parameters of the spectroscopic orbit calculated by sbop.

The discovery light curves (Krzesiński et al. 1999) are the result of approximately 200 hours of observations in the broad-band $`UBVI`$ filters. Additional photometric data were obtained with the Jakobus Kapteyn Telescope, on La Palma, in $`BVI`$ for V615 Per and in the Strömgren $`b`$ and $`y`$ filters for V618 Per.

The light curves were fitted using the ebop code (Fig. 2). A spectroscopic light ratio was used to break the degeneracy between the light ratio and the ratio of the radii for V615 Per. The quoted uncertainties were derived from solutions of different-filter light curves and increased to include uncertainties in other quantities, e.g., orbital period and (theoretical) limb darkening coefficients.

## 3. Discussion and conclusions

We have derived absolute dimensions (Table 1) for two early-type detached eclipsing binaries in the young open cluster h Persei (Southworth et al. 2004). All four stars are slow rotators and both orbits have negligible eccentricity, despite the timescales for rotational synchronisation and orbital circularisation being much greater than the age of the cluster. This strongly supports the hypothesis of binary star formation by ‘delayed break-up’ (Tohline 2002).

The stars have been plotted in the mass–radius plane in Fig. 3 with the Granada theoretical models (Claret & Giménez 1998 and references therein) and a good fit is found for a chemical composition of ($`X`$,$`Z`$) = (0.63,0.01). Whilst several studies have derived atmospheric metal abundances of B stars in h and $`\chi `$ Per using high-resolution spectroscopy, this appears to be the first determination of the interior metallicity of stars in h Per using theoretical models. This metallicity is significantly below the solar metallicity assumed in recent photometric studies, suggesting that such studies could be systematically biased.

Definitive light curves are needed for V615 Per and V618 Per to derive accurate (1%) photometric parameters. Such data would allow a precise determination of the age, empirical distance, and metal and helium abundances, of h Per and a discriminate test of different stellar evolutionary models.

## References

Andersen, J. 1991, A&ARv, 3, 91

Barnes, T. G. & Evans, D. S. 1976, MNRAS, 174, 489

Claret, A. & Giménez, A. 1998, A&AS, 133, 123

Krzesiński, J., Pigulski, A. & Kołaczkowski, Z. 1999, A&A, 345, 505

Smalley, B., Smith, K. C. & Dworetsky, M. M. 2001, uclsyn Userguide

Southworth, J., Maxted, P. F. L. & Smalley, B. 2004, MNRAS, accepted

Tohline, J. E. 2002, ARA&A, 40, 349

Uribe, A., et al. 2002, PASP, 114, 233

Zucker, S. & Mazeh, T. 1994, ApJ, 420, 806


text:
# Prospects for Detecting a Neutrino Magnetic Moment with a Tritium Source and Beta-beams

## 1 Introduction

Physicists have long speculated on the existence of a neutrino magnetic moment. Since recent evidence for neutrino oscillations provides indirect evidence for a neutrino mass, this suggests that neutrinos have small magnetic moment, although the size will depend on the nature of the neutrino mass, Dirac or Majorana. In the case of a Dirac mass, standard model interactions give the neutrino a magnetic moment of $`3\times 10^{19}(m_\nu /\mathrm{eV})`$ in units of Bohr magnetons, $`\mu _B`$, although beyond the Standard Model interactions could give the neutrino a larger magnetic moment. A positive measurement of the neutrino magnetic moment would therefore provide valuable information for understanding the neutrino mass mechanism.

At present, limits on the neutrino magnetic moment come both from direct measurements and indirect considerations. The most recent experiments using reactor neutrinos set the present upper limit of $`\mu _\nu <1.0\times 10^{10}\mu _B`$ at $`90\%`$ C.L. by the MUNU collaboration at the Bugey reactor  and $`\mu _\nu <1.3\times 10^{10}\mu _B`$ at $`90\%`$ C.L. by the TEXONO collaboration at the Kuo-Sheng reactor . Previous measurements set similar upper limits. From the pioneering experiment at Savannah River  and the following experiments at the Kurtchatov and Rovno reactors, the limits $`\mu _\nu <24\times 10^{10}\mu _B`$ , $`\mu _\nu <2.4\times 10^{10}\mu _B`$ at $`90\%`$ C.L. , $`\mu _\nu <1.9\times 10^{10}\mu _B`$ at $`95\%`$ C.L.  have been derived respectively. From solar neutrino-electron scattering at SuperKamiokande a limit of $`<1.5\times 10^{10}\mu _B`$ at $`90\%`$ C.L. has been obtained . From astrophysical and cosmological considerations, namely Big Bang nucleosynthesis, star lifetime and cooling and supernova explosions , one infers upper limits in the range $`10^{11}10^{12}\mu _B`$. The exact values for these limits are model dependent. Since neutrinos mix, and different flavors will contribute in different cases, the reactor, solar and astrophysical limits cannot be compared directly. However, for experiments which make use of terrestrial neutrinos and have a short baseline relative to the energy of the neutrinos, those originally produced as $`\nu _e`$ ($`\overline{\nu }_e`$) will have had little opportunity to change flavor. Whether the baseline (L) is long or short depends on the energy of the neutrinos. For mixing associated with the solar neutrino solution, there will be little transformation as long as $`(L/\mathrm{m})25(E/\mathrm{keV})`$, and for transformation governed by the third mixing angle and mass squared difference, the baseline is short as long as $`(L/\mathrm{m})0.4(E/\mathrm{keV})`$. In the latter situation the transformation will not be large even if the baseline is long, because the mixing angle is small and, for this reason, as yet unknown.

Since the neutrino magnetic moment interaction operates in the mass eigenstate basis, an incoming electron neutrino will scatter via a photon into at least three possible final states. For Dirac neutrinos, we can write, $`\mu _e^2=_i|U_{ie}|^2\mu _i^2`$, where $`\mu _i`$ are the magnetic moments in the mass eigenstate basis, and $`|U_{ie}|^2`$ are entries in the Maki-Nakagawa-Sakata (MNS) matrix. However, if the magnetic moment were not diagonal the above equation should be expanded to include transition moments. Similarly in the case of a Majorana neutrino, one needs to take into account a sum over all transition moments. In all cases we consider here however, the measurement is of the scattering of an electron neutrino or antineutrino into any type of neutrino by way of the magnetic moment, so we can compare the case of reactor neutrinos, a tritium source and beta-beams directly without the need to decompose a neutrino magnetic moment into its possible components in the mass basis.

The scenario with a tritium source that we consider is a 200 MCurie, 20 kg source which is used in conjunction with a large Time Projection Chamber (TPC) detector. This has just been proposed as a way to explore many neutrino properties  one of which is the neutrino magnetic moment. This idea in principle is very attractive because of the low threshold of the detector and the large number of neutrinos from the source.

“Beta-beams” is a novel method to produce neutrino beams which consists in boosting radioactive nuclei, that decay through beta-decay, to obtain an intense, collimated and pure neutrino beam . The use of this new concept to have a beta-beam facility for low-energy neutrinos has been proposed by Volpe  and has the important benefit that the neutrino spectra will be well understood. Interest in this method is fairly recent and studies of the feasibility of such beams are planned .

In this letter we investigate the potential for setting new limits on the neutrino magnetic moment using a tritium source or a low-energy beta-beam. We use the specific example of a TPC detector to observe neutrino-electron scattering. We compare the sensitivities which could be obtained by these methods with what can be achieved with reactor neutrinos.

## 2 Calculations and Results

Each of the sources we consider emits a neutrino flux which peaks at different energies. For the tritium source and the low-energy beta-beams, the neutrino flux is just that of a standard beta decay spectrum from a nucleus at rest. To calculate the spectrum, we use a Fermi function to approximate the final state Coulomb interaction. Tritium is a low energy source, since the maximum energy of the antineutrinos it produces is 18.6 keV. The ions typically considered for the beta-beams, such as $`{}_{}{}^{6}\mathrm{He}`$, a $`\beta _{}`$ emitter and $`{}_{}{}^{18}\mathrm{Ne}`$, a $`\beta _+`$ emitter, have Q-values of a few MeV, although other ions with smaller Q-values may be considered in the future. In table 1 we give the fluxes for these two beta-beam emitters. Reactor antineutrinos have a different spectrum since they are produced by many pathways in the course of the fissioning of the nuclei. The understanding of the exact nature of the reactor flux at low energies is therefore less certain than in the case of the beta decay of a single isotope, but the spectrum peaks at around 1-2 MeV. We use the antineutrino fluxes as given by .

In all cases considered the neutrinos are detected through recoil electrons from neutrino-electron scattering. The weak and electromagnetic cross section as a function of recoil energy $`T`$ and neutrino energy $`E_\nu `$ is given by

$`{\displaystyle \frac{d\sigma }{dT}}`$ $`=`$ $`{\displaystyle \frac{G_F^2m_e}{2\pi }}\left[(g_V+g_A)^2+(g_Vg_A)^2\left(1{\displaystyle \frac{T}{E_\nu }}\right)^2+(g_A^2g_V^2){\displaystyle \frac{m_eT}{E_\nu ^2}}\right]`$
$`+`$ $`{\displaystyle \frac{\pi \alpha ^2\mu _\nu ^2}{m_e^2}}{\displaystyle \frac{1T/E_\nu }{T}}`$
where $`g_V=2\mathrm{sin}^2\theta _W+1/2`$, $`\theta _W`$ is the Weinberg angle, $`g_A=1/2`$ (-1/2) for $`\nu _e`$ ($`\overline{\nu }_e`$), $`m_e`$ is the electron mass and $`G_F`$ is the Fermi coupling constant. Following Eq. (1) we consider scattering of neutrino on free electrons. The effect of atomic binding is discussed in . These cross sections have to be averaged by the neutrino flux from the relevant source :

$$\frac{d\sigma }{dT}=\frac{\frac{dN_\nu }{dE_\nu }\frac{d\sigma (E_\nu )}{dT}𝑑E_\nu }{\frac{dN_\nu }{dE_\nu }𝑑E_\nu }$$
(2)
where $`dN_\nu /dE_\nu `$ is the number of neutrinos per unit energy emitted by the neutrino source. Note that the electron recoil energy is restricted to the values :

$$T\frac{2E_\nu ^2}{2E_\nu +m_e}$$
(3)
At lower electron recoil energies the flux-averaged weak cross section becomes constant while the part proportional to the magnetic moment diverges. Experimentally one searches for the contribution from the magnetic moment by looking for an excess of counts at low electron recoil energy. The lack of such an observation has produced the current limit on the magnetic moment e.g. . In tables 2 and 3 we give flux-averaged electron-neutrino cross sections (Eq.2) as a function of recoil energy.

In order to consider the number of counts as a function of electron recoil energy we take the specific example of a TPC detector as in Ref.. For both the beta-beams and the tritium source, the geometry is the same $`4\pi `$ configuration where the source is located at the center. The radius of the detector is 10 meters and contains 20 tons of Xenon. In the comparison with the reactor experiments we consider a detector located 18 m from the source and having 1 m3 volume as used in the MUNU experiment , although we take Xenon as the detection medium.

## 3 Toward New limits on the Magnetic Moment: Tritium and Beta-Beams

Here we compare the sensitivities of the recoil electron spectrum to a non-zero magnetic moment for the three different neutrino sources. In Figs. 2-3 we show the electron scattering rates as a function of electron recoil energy for cases with and without a neutrino magnetic moment. The event rate is obtained by combining the flux-averaged cross sections with both the number of electrons in the detector and the number of neutrinos per unit time and area in the detector. A 100 % efficiency has been assumed in all cases.

Figure 2 shows the electron rate for the tritium source. The detector is sensitive to an electron energy above 100 eV and the maximum electron energy from this source is 1.26 keV. We find about 900 events/year in the absence of a neutrino magnetic moment. For a magnetic moment of $`10^{11}\mu _B`$, this number is increased to ten thousand events per year. Since the cross section is diverging, most of these events occur at the lower energies, where the points are so large that they are not shown on the figure. As can be seen from the figure the limit $`\mu _\nu =10^{11}\mu _B`$ is clearly distinguishable while $`\mu _\nu =10^{12}\mu _B`$, which would produce about one hundred events in this energy range, seems more challenging. From this, one expects a sensitivity of a few $`\times 10^{12}\mu _B`$ to be achievable, with the exact value to be determined by detailed detector simulations. Fig 2 clearly shows that the use of a 200 MCurie tritium source, which has a well understood spectrum coupled with a low threshold detector, e.g., can be extremely attractive from the point of view of setting limits on the neutrino magnetic moment.

Scattering rates for neutrinos produced by reactors are shown in Figure 2. For the purposes of this figure, we assume a 2750 MW reactor, like the Bugey reactor, and a fuel composition of $`55.6\%`$ $`{}_{}{}^{235}\mathrm{U}`$, $`32.6\%`$ $`{}_{}{}^{239}\mathrm{Pu}`$, $`7.1\%`$ $`{}_{}{}^{238}\mathrm{U}`$ and $`4.7\%`$ $`{}_{}{}^{241}\mathrm{Pu}`$. For these values the number of neutrinos produced is about $`4\times 10^{20}\nu /\mathrm{s}`$ and as stated before, we assume that the detector is located about 18 meters away from the source. Since the reactor neutrino flux peaks at 1-2 MeV there are significant electron recoil events in the tenths of MeV range. Between 0.1 MeV and 1.2 MeV, there would be about 400 events in three months with this detector and this number would increase to 500 with a neutrino magnetic moment of $`5\times 10^{11}\mu _B`$. Therefore, if the electron recoil is measured in this range, it is in principle sensitive to a magnetic moment of $`5\times 10^{11}\mu _B`$ or even smaller. In order to reach this it is crucial to understand the low energy part of the neutrino spectrum . In fact, there are several sources of uncertainty, like for example $`\beta `$-decays, which are not well known, either of the fission products or induced by neutrons absorbed by the materials. A detailed discussion of the influence of such uncertainty on the $`\mu _\nu `$ sensitivity, for different energy ranges of the electron recoil energy, is presented in . The knowledge of the low energy spectra is a limiting factor in the MUNU experiment  which quotes a limit based on data taken above 900 keV.

Let us now consider the case of a beta-beam source. Similarly to the case of a static tritium source, an advantage of the beta-beams is that the neutrino fluxes can be very accurately calculated. Figure 3 shows the electron-neutrino scattering events in the range of 0.1 MeV to 1 MeV and 1 keV to 10 keV respectively. (In Figure 3b we have rounded to the nearest integer number of counts). The shape of the flux-averaged cross sections is very similar to the reactor case as reflected in the event rates shown in the figures. As can be seen, by measuring electron recoils in the keV range with a beta-beam source one could, with a sufficiently strong source, have a very clear signature for a neutrino magnetic moment of $`5\times 10^{11}\mu _B`$. These figures are for Helium-6 ions, however, similar results can be obtained using neutrinos from $`{}_{}{}^{18}\mathrm{Ne}`$. The results shown are obtained for an intensity of $`10^{15}\nu /s`$ (i.e. $`10^{15}`$ ions/s). If there is no magnetic moment, this intensity will produce about 170 events in the 0.1 MeV to 1 MeV range per year and 3 events in the 1 keV to 10 keV range per year. These numbers increase to 210 and 55 respectively in the case of a magnetic moment of $`5\times 10^{11}\mu _B`$.

Although the first feasibility study shows that an intensity of a few $`\times 10^{13}\mathrm{ions}/\mathrm{s}`$ can be reached for $`{}_{}{}^{6}\mathrm{He}`$ , higher intensities such as $`10^{15}\mathrm{ions}/\mathrm{s}`$ may be obtained for specific ions and/or designs . An intensity at least that high is needed to improve the present limit on the neutrino magnetic moment.

## 4 Conclusions

We have discussed different routes for obtaining improved limits on the neutrino magnetic moment. The best scenarios for improving the limits, or indeed for detecting a neutrino magnetic moment, will involve three main aspects. One needs to understand well the spectrum of the neutrino source, to be able to measure electron recoil at as low energy as possible, and to have an intense neutrino flux.

The main advantage of a large tritium source when combined with a low threshold detector such as a large TPC detector, is that the neutrino flux is well known and many counts can be obtained at low recoil energy. Such a configuration is very attractive since one can obtain sensitivities to the neutrino magnetic moment of a few $`\times 10^{12}\mu _B`$.

Reactor neutrinos have provided the best terrestrial limits on the neutrino magnetic moment to date $`\mu _\nu 10^{10}\mu _B`$. Reactor experiments have the advantage that there is an intense source that can be utilized. In order to improve the current limits using reactor neutrinos, one needs to be able to measure and understand well the electron recoil spectrum below an MeV, for which it is necessary to understand the low energy component of the neutrino flux.

Beta-beams provide a new and as yet unexplored possibility for measuring a neutrino magnetic moment. As in the case of the tritium source, the knowledge of the spectrum of neutrinos is very good. The limiting feature of this scenario is the intensity of the neutrino fluxes. One needs at least $`10^{15}\mathrm{ions}/\mathrm{s}`$ in order to improve the current limit and reach a few $`\times 10^{11}\mu _B`$. For higher fluxes, one obtains better limits. For a low energy application such as this, it would be extremely helpful to evaluate ways to attain high production rates for beta-beam ions.

We thank J. Bouchez, I. Giomataris, M. Lindroos, A. Villari, H. Weick for useful discussion. G.C.M. thanks IN2P3 for their support while this work was being completed and also acknowledges support from DOE Grant DE-FG02-02ER41216.