text:
# Small and Large Scale Granular Statics

## 1 Introduction

The modeling of granular materials has been a subject of ongoing research in the engineering community (see e.g., Levy01 ). In recent years, this subject has found renewed interest among physicists Jaeger96b ; Jaeger96a ; DeGennes99 ; Kadanoff99  (having been studied in the distant past by great physicists such as Coulomb, Faraday, Reynolds and others).

The behavior of “granular gases”, which are obtained by e.g., sufficiently strong shaking or shearing (so that the material behavior is dominated by interparticle collisions), has been quite successfully modeled using approaches based on extensions of the kinetic theory of gases Goldhirsch03 . However, the behavior of dense granular matter, which is dominated by prolonged interparticle contact, has proven more difficult for modeling. For the description of the quasi-static behavior, elasto-plastic models are commonly used by engineers Nedderman92 ; Savage98b .

This paper is concerned with the static behavior of granular systems. In elasto-plastic models, one often uses (linear) elasticity below yield (although parts of a static system are sometime assumed to be at incipient yield Savage98b ). However, in recent years a very different class of models has been proposed for describing the statics of granular materials, based on the notion of “force propagation”, suggested by the observation of force chains in experiments on granular materials Drescher72 , as well as simulations Cundall79 ; Radjai99 . These models (see e.g., Wittmer96 ; Bouchaud95 ; Wittmer97 ; Tkachenko99 ) typically yield hyperbolic partial differential equations for the stress field, in contrast with the elliptic, non-propagating nature of the classical equations of static elasticity. It has been claimed that the hyperbolic description tends to an elastic-like one at large scales Bouchaud01 ; Socolar02 ; Otto03  (however, the physical interpretation of the macroscopic fields in this case is not clear).

Recently, the response of granular slabs resting on a horizontal floor to a ‘point force’ applied at the center of the top of the system has been studied experimentally Silva00 ; Reydellet01 ; Serero01 ; Geng01 ; Geng03 ; Mueggenburg02 ; Moukarzel03 . In Silva00 ; Geng01 ; Geng03 , the intergrain force distribution has been measured in two-dimensional (2D) systems as a function of vertical and horizontal distance from the point of application of the force. In Moukarzel03 , the particle displacements for similar 2D systems have been measured. In Reydellet01 ; Serero01 ; Mueggenburg02 , the vertical force acting on the floor has been measured in three-dimensional (3D) systems. Prominent force chains have been observed in ordered 2D systems; these force chains fade out with increasing disorder. For pentagonal particles in 2D arrangements the measured force distribution is single peaked and the width of the peak is linearly related to the vertical distance, in conformity with elasticity. The results for cuboidal particles obtained in Silva00  appear to suggest a parabolic behavior, consistent with a diffusive model, although the systems studied were quite small. In Moukarzel03 , the width of the measured distribution of displacements, as a function of the vertical distance from the particle which is directly displaced, follows a square root dependence (as expected from a diffusive model) for small distances of a few particle diameters, crossing over to a linear dependence at larger distances (consistent with an elastic description). Ordered 3D packings exhibit multiple force peaks for shallow systems Mueggenburg02  and less structure for deeper ones. Somewhat larger (in terms of number of particles), disordered 3D systems Reydellet01 ; Serero01  exhibit a single peak in the force distribution measured at the floor, whose width is proportional to the depth of the system.

The experimental evidence appears to be contradictory: different experiments seem to support fundamentally different descriptions of the response of granular materials (in the case of 2D systems, it has been suggested that there may be a crossover from a hyperbolic to an elliptic behavior with increasing disorder Geng01 ). The thesis presented in this paper is that these seemingly contradictory experimental results (and theoretical explanations) are not necessarily at odds with each other. This thesis is based on the observation that most of the studies (perhaps all) rejecting the elliptic description have been devoted to small systems, of the size of a few dozen of particle diameters at most, whereas many engineering studies consider rather large granular systems. Since elasticity and other macroscopic descriptions are not valid on small scales, at which local anisotropies and randomness play a major role, one should not be surprised that such descriptions fail on small scales. Indeed, simulations Goldenberg02 ; Goldhirsch02  reveal the existence of a crossover from microscopic to macroscopic behavior of granular assemblies (as well as other systems Tanguy02 ) as a function of system size or resolution. We argue that such a crossover is observed in some of the experiments mentioned above. Strictly isostatic systems Moukarzel01  have been shown to be described by hyperbolic stress equations Tkachenko99 ; Head01 , and numerical simulations suggest that systems of frictionless spherical particles approach isostaticity in the limit of infinite rigidity Silbert01 ). However, we argue that since real granular systems have finite rigidity and usually experience frictional interactions, they cannot be generically isostatic (the same presumably holds even for frictionless non-spherical grains). The isostatic limit is a singular case, whose physical consequences for real systems are at best unclear. Therefore the controversy surrounding the correct description of granular statics is mostly a question concerning the behavior of small granular systems. The latter require a grain-scale (“microscopic”) description, rather than a macroscopic one.

A second point stressed below is the distinction between force and stress. Whereas interparticle forces can exhibit force chains which look like they contradict elasticity, the latter does not describe the nature of the forces but rather that of the stress field. The stress field involves an averaging over the forces (whose result is resolution dependent) and leads to less pronounced structure than the underlying force field. The small scale structure of the interparticle forces cannot be taken to consist an argument against an elliptic description or in favor of it, since it relates to small scales and it does not deal with the objects with which elasticity or plasticity are concerned. The large scale response of granular packing is shown to be consistent with a (possibly anisotropic) elastic description. The fact that in non-cohesive granular materials there is no significant attraction among the particles may lead to modifications of the contact network, which may strongly affect the anisotropy of the system. This effect is expected to be small for small applied forces (for non-isostatic systems) and for pre-stressed systems, in particular for disordered systems. Otherwise, it may be accounted for using a nonlinear, incrementally elastic model, with stress-history dependent elastic moduli.

The third point made in this paper is that while models employing frictionless particles can reproduce some properties of granular packings, friction can be of utmost importance for the description of granular matter (a rather intuitive fact). Results demonstrating the importance of accounting for frictional interactions are presented in Sec. 5.

## 2 The microscopic picture: forces

In attempting to describe granular materials in terms of continuum mechanics, by analogy to “regular”, atomic materials, one usually considers the “microscopic” scale to be that of the individual particles (whose internal dynamics should be well described by continuum mechanics).

One of the simplest granular systems is a collection of frictionless spherical particles. A typical microscopic (particle scale) description of such a system is given by the particle’s radii, $`\{R_i\}`$, their masses, $`\{m_i\}`$, center of mass positions, $`\{\stackrel{}{r}_i(t)\}`$, and velocities, $`\{\stackrel{}{v}_i(t)\}`$, at time $`t`$. It is typically assumed (e.g., in the context of simulations of granular materials Cundall79 ; Herrmann98 ; Wolf96 ) that the particles are quite rigid, so that the interaction between two particles (in the frictionless case) depends only on their respective distance, or, more conveniently, on their imaginary overlap $`\xi _{ij}(t)R_i+R_j|\stackrel{}{r}_{ij}(t)|`$, where $`\stackrel{}{r}_{ij}(t)\stackrel{}{r}_i(t)\stackrel{}{r}_j(t)`$. The contact interactions are usually modeled by treating the particles as macroscopic objects, described by the equations of continuum mechanics (see e.g., Gladwell80 ; Johnson85 ). For two frictionless elastic spheres, a classical result by Hertz (see e.g., Landau86 ) is that the force is proportional to $`\xi ^{3/2}`$, while for cylinders, it is linear in the overlap. For noncohesive particles, only repulsive forces are possible. Even for frictionless particles, internal dissipation as in e.g., viscoelastic particles, gives rise to a dependence of the force on the relative velocity $`\dot{\xi }`$ as well (for some examples of force schemes commonly used in simulations, see e.g., Schafer96 ; Sadd93 ; Walton95 ).

The interparticle forces for a given configuration of such particles subject to given boundary conditions (e.g., specified displacements of the particles on the boundary, or forces applied to them) and body forces such as gravity can be determined, for a static system, using the equations of equilibrium (Newton’s laws) and the force-displacement relation. We reiterate that when full force laws for particle interactions are known or modeled, the statics and dynamics of the system are fully determined (they may be history-dependent for history-dependent force laws, as commonly used for frictional interactions).

In the case of frictionless isostatic systems (in which the mean coordination number is exactly $`z=2d`$, where $`d`$ is the dimension of the system) the forces can be determined from the equations of equilibrium alone (and are therefore independent of the force-displacement law; however, the particle displacements certainly depend on this law). It has been suggested Moukarzel98c ; Roux00  that frictionless granular systems become isostatic in the limit of infinite rigidity (giving rise to a macroscopic behavior which is very different from elasticity), and this appears to be borne out by numerical simulations Makse00 ; Silbert01 . However, the relevance of this limit to real materials is questionable, since real materials cannot be infinitely rigid. Any additional contacts created if the rigidity is allowed to be finite will render the system hyperstatic (so that there is a “phase transition” to an isostatic behavior only at infinite rigidity Moukarzel98c ). The rigidity should of course be compared to the confining forces or body forces (in a system under gravity and confined by walls, the confining force is related to gravity). If the confining forces are very small, the system would indeed be expected to be close to marginal stability. As mentioned above, the static indeterminacy associated with hyperstatic systems simply means that the equations of equilibrium are insufficient for determining the forces, so that additional equations (e.g., force-displacement laws) are required. Static indeterminacy does not mean that there’s no unique solution for the forces in a real system. A similar situation occurs on the macroscopic, continuum level (see e.g., Savage98b ). The rigid limit can be approached in many different ways (e.g., the stiffness of each interparticle contact may be different), and, even if assuming that the same (isostatic) contact network is obtained for different distribution of the interparticle stiffness, yielding the same interparticle forces, the particle displacements will certainly be different, hence the rigid limit in not unique, at least in this sense.

In several experiments, photoelastic particles were used in order to measure the stress in granular systems Drescher72 ; Howell99 ; Silva00 ; Geng01 ; Geng03 . These measurements probe the intraparticle stress, i.e., the stress inside each particle. Following the above, these should be interpreted as measurements of microscopic fields (the macroscopic description of granular systems regards the particles as microscopic, and does not resolve any details below the particle scale). The microscopic fields corresponding to these measurements are the interparticle forces, which can be deduced from these internal stress measurements (as described in Geng03 ). As mentioned, these forces should be distinguished from the “macroscopic” stress field in the system.

The distribution of force magnitudes in a static granular packings is a microscopic quantity which has been extensively studied in experiments Mueth98 ; Blair01  and simulations Radjai96 . An exponential behavior of the distribution at large forces appears to be quite universal in experiments on granular systems, independent of the degree of disorder Blair01 , the friction coefficient Blair01 , or the rigidity of the particles Erikson02 , and has also been observed in simulations of granular systems with different models for the interparticle forces (e.g., Nguyen00 ; Silbert02a ). The universality of the force distribution appears to extend to other systems such as foams, glasses, colloids etc. (see O'hern00  and references therein). The exponential tail of the distribution is reproduced in simple models such as the (parabolic) q-model Liu95 ; Coppersmith96 . The distribution for smaller forces appears to be less universal, and it has been suggested that the appearance of a peak in the force distribution near the mean force may signal the onset of jamming or a glass transition O'hern02 .

Interestingly, a qualitatively similar force distribution is obtained in purely harmonic networks: Fig. 1 shows the force distribution obtained for an ensemble of random networks of linear springs constructed as follows. Points are placed on a 2D triangular lattice with spacing $`d`$ (with square-shaped boundaries), and then their $`x`$ and $`y`$ coordinates are randomly displaced by $`\pm 0.04d`$. Points whose distance is less than $`1.02d`$ are connected by linear springs (whose equilibrium length is equal to this distance) with equal spring constants (this results in an average dilution of about $`12\%`$ of the springs compared to the perfect lattice). A uniform isotropic compression of $`1\%`$ is applied to the boundary particles, and the interparticle forces are calculated. The force distribution presented in Fig. 1 is obtained from an average over the force histograms of $`100`$ systems of $`1085`$ particles. The force was normalized by the mean force in the ensemble (a similar distribution is obtained for a normalization by the mean force for each system; the variation in mean force among different systems is relatively small, which may indicate that the system is far from ‘jamming’ O'hern02 ). The tail of the logarithm of the distribution is fit quite well with a line of slope $`3.8`$, similar to the slope obtained in experiments on highly compressed disordered packings of soft rubber spheres Erikson02  (similar distributions were obtained for a scalar harmonic network of unequal springs in Nguyen00 ). For the case of networks with no force dilution (the same connectivity as in the perfect lattice), the force distribution is Gaussian with a half-width of a few percent of the mean, i.e., a much narrower distribution). These results indicate that a random connectivity should be consequential for the force distribution, which may be the reason that even for highly compressed disordered spheres (whose contact network is still disordered), the distribution is qualitatively similar to that observed in less compressed systems Erikson02 . A similar effect has been observed in simulations of granular systems under different applied pressures Nguyen00 .

The forces in one of the realizations of the ensemble are shown in Fig. 2. Force chains are clearly observed (note that there are very few tensional forces, so that they do not significantly affect the force distribution in this case). Similar force chains have been observed in a polydisperse Lenard-Jones system Tanguy02  (incidentally, the concept of a force chain is not well-defined: in the case of a homogeneous strain applied to a uniform lattice, the forces are equal, so that it is reasonable to define the force chains to contain forces whose magnitude is larger than a uniform cutoff, e.g., the mean force, as used in Fig. 2; However, for a non-uniformly strained system, e.g., systems subject to gravity, in which the mean force increases with depth, such a global cutoff makes little sense). The results shown in Fig. 2 indicate that force chains are not specific to granular systems. Force chains are microscopic features of microscopically disordered systems (or even inhomogeneously strained ordered systems, as described below), and their presence does not necessarily indicate any macroscopic inhomogeneity, or inconsistency with a macroscopic elliptic, or elastic, description. It is quite certain that if one could observe the individual interparticle forces in atomic systems (which may not be quite well defined, since a quantum description is appropriate for such systems), one would also observe force chains.

It is important to note that a significant portion of the stress (even in a homogeneously strained system) is carried by forces which do not belong to the force chains. An example is provided by a system of frictionless polydisperse disks (with radii uniformly distributed within $`10\%`$ of the maximum radius) which is confined by side walls and a floor, with a uniform force applied to the particles of the “top” layer (without gravity). The interparticle forces are taken to be linear in the overlaps. Fig. 3 shows the forces in the system. Fig. 4 shows the fraction of the applied vertical force carried by the forces whose magnitude is greater than the mean (i.e., those belonging to force chains, using the definition mentioned above), compared to that carried by all the forces (which is of course equal to $`1`$), for forces in horizontal “slices” of the system, as a function of the vertical coordinate, $`z`$. As seen in Fig. 4, only about $`80\%`$ of the applied force is carried by the force chains. Furthermore, the force carried by the chains fluctuates with depth, so that the forces in the chains do not obey the conditions of force equilibrium. It is therefore questionable whether a model which describes the stress exclusively in terms of the force chains is justifiable.

The (near-)universality of the force distribution, in particular the fact that it is observed in simulations of random systems with harmonic interactions, does not make possible the differentiation between different models on the basis of the force distribution (in particular, the observation of such a distribution does not preclude an elliptic description). The same statement applies to the observation of force chains. A more sensitive and direct test should be rendered by the response of a granular system to inhomogeneous external forcing, such as that provided by localized forces. The latter seem to be consistent with elasticity, as described below.

## 3 Macroscopic fields and continuum equations in terms of microscopic quantities

Continuum descriptions of materials are often based on phenomenological arguments (usually motivated by experimental findings), rather than on derivations from the underlying microscopic dynamics. A unique feature of granular materials is that due to the typically large sizes of the constituents, it is relatively easy to access the “microscopic” scales experimentally. On the other hand, in most practical applications, the number of particles is such that a detailed particle-level description becomes intractable, and a continuum description is required. The fact that experiments on granular systems can yield both microscopic information and macroscopic information (possibly even in the same experiment) is useful to the elucidation of the connection between these two descriptions.

In order to obtain a macroscopic description of a system in terms of the microscopic fields, we employ a spatial coarse-graining approach Glasser01 ; Goldhirsch02 . We stress that the only averaging considered here is spatial (the approach can be extended to include temporal coarse-graining as well Glasser01 , but here we consider static configurations). Since static granular packings are typically found in metastable states, far from equilibrium, and thermal energy scales are negligible, such systems do not explore any phase space so that it is hard to justify the kind of ensemble average commonly used in statistical mechanics. An average over configurations (i.e., average over different disordered systems which are presumed to be prepared in the same way) is commonly performed when analyzing experimental data, due to the large fluctuations obtained in many experiments. However it is not clear a-priori if self-averaging occurs, i.e., that at least for large enough scales the macroscopic behavior of a single “typical” realization is the same as that of the average behavior over many realizations. Self-averaging may be valid for some quantities and not for others. Therefore we choose not to assume a-priory that any ensemble averaging is justified; instead we relate the macroscopic and microscopic fields in a way that is relevant for single realizations.

Following Glasser01 , define the coarse-grained (CG) mass density $`\rho (\stackrel{}{r},t)`$ and momentum density $`\stackrel{}{p}(\stackrel{}{r},t)`$ at position $`\stackrel{}{r}`$ and time $`t`$ as

$`\rho (\stackrel{}{r},t)`$ $``$ $`{\displaystyle \underset{i}{}}m_i\varphi [\stackrel{}{r}\stackrel{}{r}_i(t)],`$ (1)
$`\stackrel{}{p}(\stackrel{}{r},t)`$ $``$ $`{\displaystyle \underset{i}{}}m_i\stackrel{}{v}_i(t)\varphi [\stackrel{}{r}\stackrel{}{r}_i(t)],`$ (2)
where $`\varphi (\stackrel{}{R})`$ is a non-negative coarse-graining function (with a single maximum at $`\stackrel{}{R}=0`$) of width $`w`$, the coarse-graining scale, and $`\varphi (\stackrel{}{R})𝑑\stackrel{}{R}=1`$.

Upon taking the time derivative of the macroscopic fields $`\rho `$ and $`\stackrel{}{p}`$, performing straightforward algebraic manipulations Glasser01  and using Newton’s laws, one obtains the equation of continuity and the momentum conservation equation, respectively:

$`\dot{\rho }`$ $`=`$ $`\mathrm{div}(\rho \stackrel{}{V})`$ (3)
$`\dot{p}_\alpha `$ $`=`$ $`{\displaystyle \underset{\beta }{}}{\displaystyle \frac{}{r_\beta }}\left[\rho V_\alpha V_\beta \sigma _{\alpha \beta }\right],`$
where the velocity field is defined by $`\stackrel{}{V}\stackrel{}{p}/\rho `$, Greek indices denote Cartesian coordinates, and the explicit dependence of the CG fields on $`\stackrel{}{r}`$ and $`t`$ has been omitted for compactness. Since this paper focuses on the stress field, we have omitted the energy equation, which can be derived in a similar way Goldhirsch02 .

In addition to obtaining the standard equations of continuum mechanics from microscopic consideration, this coarse graining procedure provides an expression for the stress tensor $`\sigma _{\alpha \beta }`$ in terms of the microscopic entities:

$`\sigma _{\alpha \beta }(t)`$ $`=`$ $`{\displaystyle \underset{i}{}}m_iv_{i\alpha }^{}(\stackrel{}{r},t)v_{i\beta }^{}(\stackrel{}{r},t)\varphi (\stackrel{}{r}\stackrel{}{r}_i(t))`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij;ij}{}}f_{ij\alpha }(t)r_{ij\beta }(t){\displaystyle _0^1}𝑑s\varphi [\stackrel{}{r}\stackrel{}{r}_i(t)+s\stackrel{}{r}_{ij}(t)],`$
where $`\stackrel{}{v}_i^{}(\stackrel{}{r},t)\stackrel{}{v}_i(t)\stackrel{}{v}(\stackrel{}{r},t)`$ is the fluctuating velocity, $`\stackrel{}{f}_{ij}(t)`$ is the force exerted on particle $`i`$ by particle $`j`$, and $`\stackrel{}{r}_{ij}(t)\stackrel{}{r}_i(t)\stackrel{}{r}_j(t)`$.

The first term in Eq. (3) is the kinetic stress (which vanishes for static configurations), and the second term is known as the contact stress. Note that the standard Born-Huang expression Born88 : $`\sigma _{\alpha \beta }=\frac{1}{2V}_{ijV;ij}f_{ij\alpha }r_{ij\beta }`$ is equivalent to the expression for the contact stress in Eq. (3) if the coarse-graining function is taken constant inside a volume $`V`$ and zero outside it, provided that the interparticle separation is much smaller than the coarse-graining length scale (typically $`\sqrt[3]{V}`$).

The above expressions can be used to calculate the macroscopic fields from the microscopic ones (obtained e.g., from simulations or experiments), and compare them to the predictions of macroscopic models or direct experimental results. In order to close the set of continuum equations \[Eqs. (3)\] the stress and energy flux (the latter is not related to the considerations below) need to be expressed as functionals of the pertinent macroscopic fields. As mentioned, such constitutive relations are often obtained empirically or conjectured. In some cases they are derived from the microscopic dynamics. The above exact expression for the stress field provides a framework for a systematic derivation of constitutive relations (as suggested for elastic networks in Goldhirsch02 ).

Here, we are concerned with the interpretation of experimental data in terms of microscopic variables and macroscopic fields. The fact that the contact stress includes a sum over all contacts for each particle (i.e., even for very small CG scales the stress components correspond to specific “averages” over the forces on each particle) already suggests that a “picture” of the forces in the packing does not correspond directly to the macroscopic stress field (they are certainly related, i.e., one would usually expect large stress components in regions where large force magnitudes are observed). In particular, as shown in Sec. 2, and further discussed below, force chains do not necessarily indicate macroscopic anisotropy or inhomogeneity.

## 4 Numerical results for model frictionless systems

Consider a two-dimensional system of uniform disks (arranged on a triangular lattice) subject to a vertical external force at the center of the top layer Goldenberg02 . Experiments on such systems are described in Geng01 ; Geng03 . Consider first the case of nearest neighbor harmonic interactions, i.e., the disks are coupled by equal linear springs (whose rest length is the diameter of a disk). Clearly, real cohesionless particles do not experience any significant attractive interactions; however, there are a few insights to be obtained from the study of this system. Fig. 5 presents the forces in the system. Force chains are evident.

A contour plot of the “vertical stress component” $`\sigma _{zz}`$ \[computed using Eq. (3)\] for the same system is shown in Fig. 6 (with $`\varphi (𝐫)=\frac{1}{\pi w^2}e^{(|𝐫|/w)^2}`$, and $`w=d`$, the particle diameter, i.e., a fine resolution). The force chains are not evident any more. The model described above corresponds, in the continuum (long-wavelength) limit, to an isotropic 2D elastic medium GoldenbergUP . The observed force chains, which break isotropy, can be attributed to the fact that the local environment of a particle in contact with a finite number of other particles cannot be isotropic. Under homogeneous macroscopic deformation, all forces would be equal in a lattice configuration. However, the concentrated applied force yields an inhomogeneous deformation, which leads to the local anisotropy being reflected in the distribution of the forces. The elastic continuum description of the stress (to linear order in the strain) is isotropic, and cannot be expected to reflect this microscopic anisotropy. For small system sizes (in which the strain gradients on a particle scale are relatively large), this anisotropy can be observed in the stress field (a very clear example is shown in Goldenberg02  for a macroscopically isotropic 3D system, whose microscopic symmetry is cubic). These results, as well as those presented in Sec. 1 for disordered elastic systems, show that force chains do not necessarily indicate anisotropy or inhomogeneity of the material on sufficiently large scales; more importantly their existence does not require a non-elastic (microscopic) interaction.

Note that only the forces between the particles and the floor (a single such force per particle) are used in the calculation of the stress at the bottom of the packing. Hence on the bottom (but not in the bulk of the system), the spatial distribution of $`\sigma _{zz}`$ is equivalent (up to coarse-graining) to that of the microscopic forces. For sufficiently large systems, the distribution of forces on the bottom corresponds closely to the stress calculated using linear elasticity Goldenberg02 , even “almost without coarse-graining”, i.e., for a microscopic CG scale.

A more realistic force model consists of ‘one-sided’ springs, i.e., springs that snap when in tension. Fig. 7 presents the forces obtained for the same system presented in Fig. 7, but with ‘one-sided’ springs. Compared to the system of regular springs, the application of the concentrated force at the top of the packing leads to rearrangements in the contact network: some horizontal springs in the region under the point where the force is applied are disconnected (as also observed in Luding97  for a pile geometry) but the force chains in both systems are qualitatively similar. The force distribution vs. the horizontal coordinate at different depths is in good agreement Goldenberg02  with experiment Geng01 ; Geng03 . For slightly disordered systems Goldenberg02 , the force chains are qualitatively similar, though somewhat shorter.

The corresponding vertical stress field $`\sigma _{zz}`$ is shown in Fig. 8. The stress field in this case is clearly quite different from that obtained using ‘two-sided’, harmonic springs: the response for ‘one-sided’ springs is double-peaked. This is obviously related to the disconnected springs below the point of application of the external force. In Goldenberg02 , it has been shown that a model with harmonic springs in which the spring constant for the horizontal springs, $`K_1`$, is different from that of the oblique springs, $`K_2`$, corresponds (in the continuum limit) to an anisotropic elastic system. For sufficiently large $`K_2/K_1`$, the response of such an elastic system has two peaks Goldenberg02  (see also Otto03  for a more detailed analysis of the case of an infinite half-plane; the results presented in Goldenberg02  are for a finite slab on a rigid floor). The absence of horizontal springs corresponds to the limit $`K_2/K_1\mathrm{}`$, the extreme anisotropic limit, which corresponds to an isostatic system. Note that the stress field, but not the displacement, depends only on $`K_2/K_1`$. In the case considered here, $`K_1=0`$ and $`K_2`$ is finite; for $`K_2\mathrm{}`$, the rigid limit, the displacement is zero. The double peaked stress distributions are similar to those obtained from hyperbolic models. It follows that hyperbolic-like behavior can be obtained using an anisotropic (yet, still elliptic) elastic model (which becomes formally ‘hyperbolic’ in the limit of very large anisotropy; see also Cates99 ; Otto03 ).

The ‘stress-induced anisotropy’ Goldenberg02  observed in the case of ‘one-sided’ springs can be thought of as a nonlinear extension of the linear elastic continuum behavior obtained in a network of harmonic springs. While the (possibly position-dependent) elastic moduli in linear elasticity are time-independent material properties, a possible extension would be to introduce a stress history dependence of the elastic moduli (i.e., the anisotropy induced by the breaking of contacts in certain regions may be considered a result of a tensile stress in those regions). A similar type of stress-induced anisotropy has been suggested in the context of plastic models for soil mechanics Oda93 . If the particle positions do not change significantly, so that only the contact network is modified in response to the applied stress, the behavior can possibly modeled as ‘incrementally elastic’. Under certain condition (corresponding to plastic yield), the system would no longer be able to support the applied stress without a major rearrangement of the particles. Incipient plastic yield may possibly be related to a local extreme anisotropy typical of a marginally stable isostatic configuration.

The force chains obtained both for the harmonic case and the ‘one-sided’ case are quite similar Goldenberg02  to those observed experimentally Geng01 ; Geng03 , as an average over different realizations. This averaging is required due to experimental variations: although the particles are arranged on lattice, there is still some disorder present due to some variability in particle diameter, and possibly also in the contact properties GengPC . Indeed, while perfect atomic lattices may be obtained at low enough temperatures, since all atoms of the same isotope are exactly identical, macroscopic particles are never truly identical, so that perfect periodicity can never be obtained. It also appears that the force chains obtained using the two models are quite similar (similar chains are also obtained in slightly disordered systems Goldenberg02 ). The stress field appears to be more sensitive to the anisotropy induced by the applied force.

In the experiments reported in Reydellet01 ; Serero01 ; Mueggenburg02 , the forces on the floor were measured. In Reydellet01 ; Serero01 , the width of the pressure probe (which would correspond to the CG scale of the measured stress) was $`1030`$ particle diameters. The bottom stress profiles measured are quite consistent with continuum elasticity (note that the depths of the systems studied were $`20300`$ particle diameters). Experimental deviations from the predictions of isotropic elasticity Serero01  can be reproduced by anisotropic elasticity Goldenberg02 . Narrower or wider peaks than those obtained for isotropic systems can be obtained for small anisotropy, while for very large anisotropy, two peaks are expected (see also Otto03 ). An additional possible cause for deviations from the isotropic elastic calculations presented in Serero01  is finite rigidity of the floor GoldenbergUP . The more shallow systems used in the experiments may even be small enough for the finite size effects Goldenberg02  to be significant. Any anisotropy in these experiments is obviously much weaker than the strong anisotropy observed in the model ordered system of ‘one-sided’ springs. Several effects may explain this: first, the systems used in the experiment are highly disordered, so that inhomogeneous, random anisotropy may be expected on intermediate (already macroscopic) scales, presumably averaging out to an isotropic, or nearly isotropic, behavior at sufficiently large scales. In this case, the large-scale effect of contacts breaking due to applied forces would be significantly less pronounced than in the ordered system described above. A second possibility is the effect of frictional forces, which may either prevent contacts from breaking, or reduce the anisotropy of the response. Third, the model systems discussed above were unstressed before the application of the force, while the experimental ones are pre-stressed by gravity, which may compress some of the contacts such that the tension due to the applied force is insufficient to break them.

In Mueggenburg02 , individual forces on the floor were measured, and the results were averaged over realizations (which, as mentioned, is not necessarily equivalent to spatial coarse-graining). The regular packings used in Mueggenburg02  (FCC and HCP) are macroscopically anisotropic. The fact that some of the horizontal contacts (contacts among particles in the same layer) may be absent increases the anisotropy further (possibly in an inhomogeneous way; as mentioned above, a granular packing cannot be perfectly periodic). The extreme limit in which there are no such horizontal contacts corresponds to an isostatic system. Such anisotropy (possibly further enhanced by the applied force) may explain the discrete peaks observed for relatively shallow systems composed of $`9`$ layers of particles (and the fact that they appear to be consistent with a picture of “force propagation” appropriate for isostatic systems). However, for deeper system (about $`20`$ particle diameters), there appears to be a crossover to a smoother behavior, which should correspond to the crossover to the continuum limit (note that the depth of the systems used in Mueggenburg02  was smaller than the depth required in our calculations on 3D systems Goldenberg02  for reaching the continuum limit, so deeper systems may still show dependence on the depth).

## 5 Effects of Friction

As shown, some features of granular response may be reproduced using models employing frictionless and even harmonically interacting particles. However, it is clear that friction is consequential for granular materials.

For frictional spheres, the microscopic description, as described in Sec. 2, must be extended to include (at least) the orientations of the particles, and interparticle torques in addition to the forces. The description of static and kinetic friction requires the use of more complicated force models, which depend on the particle orientations and their relative tangential velocities, and possibly on the history of contact deformation (see e.g., Schafer96 ; Sadd93 ; Walton95 ; Vu-Quoc99 ).

In experiments performed on regular 2D packings of photoelastic disks Geng03 , the directions and “strengths” of the force chains observed upon application of a localized force to the top of the packing appear to depend quite strongly on the angle of the applied force with respect to the horizontal (in the following, all angles are given with respect to the horizontal). A particularly intriguing effect is that for some angles, force chains appear not only in the lattice directions ($`0,\pm 60^{},\pm 120^{},180^{}`$ for a triangular lattice), but also, apparently, in new directions which can be identified as $`\pm 30^{}`$ (in fact, in individual realizations, rather than their average as reported in Geng03 , it appears that force chains appear also at $`\pm 90^{}`$, i.e., the vertical direction GengPC ). These directions correspond to next-nearest neighbor directions in the triangular lattice. Since interactions among the particles only exist for particles in contact, there is no direct next-nearest neighbor interaction. The fact that the forces themselves, and not just the contact points, appear to be aligned with these $`\pm 30^{}`$ directions, suggests that frictional forces among the particles (tangential to the contact normals, which result in interparticle torques) are necessary for obtaining forces (and chains) at angles different from the lattice directions. For an applied force at $`\pm 90^{}`$, it appears that the frictional forces are small enough such that the results obtained in this case Geng01 ; Geng03  are described quite well by a model of frictionless particles with linear force-displacement laws Goldenberg02 .

In order to elucidate the role of frictional forces and torques in the quasi-static response of granular materials in general, and in particular in order to gain an understanding of the experimental results mentioned above Geng03 , we performed discrete element simulations with normal and tangential linear spring-dashpot forces among the particles (see e.g., Cundall79 ; Herrmann98 ), possibly the simplest model for frictional disks. The simulation parameters were chosen to correspond to those of the experimental system GengPC . Experimentally, the force-displacement law for the photoelastic disks was found to be fit quite well by $`f\xi ^{3/2}`$ GengPC , as predicted by the standard Hertz theory for elastic ellipsoids in contact (see e.g. Landau86 ), rather than the linear relation (with logarithmic correction) expected for cylinders in contact Schwartz71 , which appears to imply that the contact region between the “disks” is elliptic rather than rectangular. The simulation model described above employs a linear force-displacement law, so that an effective mean spring constant was estimated on the basis of the range of forces used in the experiments. The tangential spring constant was taken to be one-half the normal spring constant (a rough estimate consistent with the Hertz-Mindlin model Mindlin53  for oblique contact forces). The normal and tangential spring constants used are $`k_n=3000\overline{m}g/\overline{R}`$ and $`k_t=1500\overline{m}g/\overline{R}`$, where $`\overline{R}`$ and $`\overline{m}`$ are the mean particle radius and mass, respectively, and $`g`$ is the gravitational acceleration. The friction coefficient used is $`\mu =0.94`$ for particle-particle contacts and $`\mu ^{\mathrm{wall}}=0.35`$ for particle-wall contacts. The systems studied here are composed of polydisperse disks, with radii distributed uniformly in the interval $`[R\delta R,R]`$, where $`\delta R/R=810^3`$ (i.e., a small polydispersity).

The system is first relaxed to a static state under gravity (until the total kinetic energy per particle is less than $`10^9\overline{m}g\overline{R}`$), and then relaxed again with an external force applied at the center of the top layer (in some cases the force was increased linearly with time from zero to prevent the “buckling” of the top layer which leads to major rearrangements; these are beyond the nearly elastic behavior considered here).

For comparison with the experiments presented in Geng03 , an external force of magnitude $`F=150\overline{m}g`$ was applied to the center top particle at angles of $`15^{},30^{},45^{},60^{},75^{}`$, and $`90^{}`$. The simulated systems consisted of $`29`$ rows of $`80`$ particles, which is similar to the size of the systems used in the experiments Geng01 ; Geng03 ; GengPC . Fig. 9 presents the forces obtained for different applied force angles. The same particle configuration was used in all cases. No significant particle rotation occurred except for the particles adjacent to the one on which the force is applied. For a force at an angle of $`15^{}`$, buckling occurred in the top row, causing major rearrangements. Such buckling was also observed in the experiments, where it was apparently stabilized, limiting the rearrangements to a small region near the point of application of the force, but we have not been able to prevent major rearrangements in the simulation. As mentioned, tangential forces such as friction give rise to interparticle torques. Simulations with an applied torque (in addition to the applied force) show that this torque does influence the observed force chains GoldenbergIP .

The results are quite similar, qualitatively, to those observed in the experiment Geng03 . Note that the results shown in Fig. 9 are for a single configuration, while the results presented in Geng03  are for an average over configurations. The results obtained in simulations for different realizations of the disorder are qualitatively similar GoldenbergIP . The agreement of the results obtained using a relatively simple force model with the experiments is encouraging. A more detailed study of the effects of friction on the forces and the stress field will be presented elsewhere GoldenbergIP .

## 6 Concluding remarks

We have shown that the seemingly inconsistent results of different kinds of experiments studying the static response of granular packings to a concentrated force can all be understood within the same framework of an essentially elastic (elliptic) picture once the distinction between forces and stress is made and the possible consequences of small system size, as well as anisotropy, are taken into account. The effect of applied stresses on the contact network may be modeled as a nonlinear, incrementally elastic model (which may be further extended to describe yielding).

Somewhat surprisingly, many aspects of the response of such systems can be understood using models of frictionless particles. However, some effects do require the introduction of friction, as in the example of the force chains obtained for oblique applied forces described in this paper. We note, however, that the model for the friction used in the simulations described in Sec. 5 consists of tangential springs (with the additional Coulomb condition). This indicates that even for static frictional systems (below yield) an elastic continuum model, which probably includes rotational degrees of freedom (e.g., a Cosserat continuum model Eringen68 ), may be appropriate.

Another important issue which requires further study is the effect of disorder, and the relation of spatial averaging to averaging over the disorder.


format: text

text:
# Spin current shot noise as a probe of interactions in mesoscopic systems

## Abstract

It is shown that the spin resolved current shot noise can probe attractive or repulsive interactions in mesoscopic systems. This is illustrated in two physical situations : i) a normal-superconducting junction where the spin current noise is found to be zero, and ii) a single electron transistor (SET), where the spin current noise is found to be Poissonian. Repulsive interactions may also lead to weak attractive correlations (bunching of opposite spins) in conditions far from equilibrium. Spin current shot noise can also be used to measure the spin relaxation time $`T_1`$, and a set-up is proposed in a quantum dot geometry.

Non-equilibrium noise plays a key role in mesoscopic physics BB . Low-temperature correlations of the time fluctuations of the electronic current indeed give unique information about the charge and the statistics of quasiparticles. For non-interacting electrons, the scattering approach is very powerful for a variety of systems BB ; scattering1 ; scattering2 . The reduction of shot noise from the Schottky value originates from the Pauli exclusion principle which forbids two wavepackets with the same quantum numbers to be superimposed MartinLandauer . On the other hand, Coulomb interactions also act in correlating wavepackets, and noise is indeed more sensitive to interactions than the conductance. Coulomb interactions may decrease or increase noise correlations, depending on the physical regimes martin\_buttiker ; superpoisson ; loss\_sukho . Full counting statistics are also promising as a probe of interactions fullcounting . Yet, in a given mesoscopic structure, the effects on the shot noise of Fermi statistics and of interactions are intimately mixed. In contrast, we propose in this Letter that spin-resolved shot noise can unambiguously probe the effects of electronic interactions. The basic idea is that the Pauli principle acts only on electrons with the same spin. Therefore currents wavepackets carried by quasiparticles with opposite spins are only correlated by the interactions.

Spin-resolved shot noise has received very little attention up to now, contrarily to the total current shot noise. For instance, spin shot noise was recently considered in absence of charge current wang . On the other hand, the effect of a spin-polarized current on charge and spin noise was investigated, with complex behaviours due to spin accumulation Bulka . Noise is also an efficient probe for testing quantum correlations in two-electron spin-entangled states loss\_entangle ; entanglermetal ; fazio ; balents  or electron spin teleportation TP . In contrast, we consider here simple and general mesoscopic structures in which the average current is not spin-polarized, but where the currents carried by quasiparticles with different spins can be separately measured. A possible set-up for this purpose will be described at the end of this Letter. To clarify our statement, let us first consider a general mesoscopic device made of a normal metal with non-interacting electrons, non magnetic terminals $`i,j`$,.., and one channel for simplicity (generalization is obvious). In absence of magnetic fields and spin scattering of any kind, the scattering matrix is diagonal in the spin variable and spin-independent, $`s_{ij}^{\sigma \sigma ^{}}=\delta _{\sigma \sigma ^{}}s_{ij}`$. This trivially leads to spin-independent averaged currents $`I_i^\sigma =I_i^\sigma `$. In a similar way, spin-resolved noise, defined as $`S_{ij}^{\sigma \sigma ^{}}(tt^{})=\frac{1}{2}\mathrm{\Delta }I_i^\sigma (t)\mathrm{\Delta }I_j^\sigma ^{}(t^{})+\mathrm{\Delta }I_j^\sigma ^{}(t^{})\mathrm{\Delta }I_i^\sigma (t)`$, where $`\mathrm{\Delta }I_i^\sigma (t)=I_i^\sigma (t)I_i^\sigma `$, can be evaluated. One easily finds that at any frequency the noise power between terminals $`i`$ and $`j`$ is diagonal in the spin variables, $`S_{ij}^{\sigma \sigma ^{}}(\omega )=\delta _{\sigma \sigma ^{}}S_{ij}(\omega )`$. Thus, choosing an arbitrary spin axis $`𝐳`$, the total (charge) current noise $`S_{ij}^{ch}=S_{ij}^{}+S_{ij}^{}+S_{ij}^{}+S_{ij}^{}`$ and the spin current noise $`S_{ij}^{sp}=S_{ij}^{}+S_{ij}^{}S_{ij}^{}S_{ij}^{}`$, defined as the correlation of the spin currents $`I_i^{sp}(t)=I_i^{}(t)I_i^{}(t)`$, are strictly equal. On the contrary, in presence of Coulomb interactions, one expects that $`S_{ij}^{}=S_{ij}^{}0`$, or equivalently $`S_{ij}^{sp}S_{ij}^{ch}`$. This can happen for instance if the scattering matrix couples carriers with opposite spins, as Andreev scattering at a normal- superconductor (NS) interface, or in quantum dots in presence of strong Coulomb repulsion. The sign of the correlation $`S_{ij}^{}`$ (bunching or antibunching) is of special interest.

Let us first consider a NS junction, where S is a singlet superconductor and N a normal metal. The scattering matrix coupling electron (e) and holes (h) quasiparticles in the metal is composed of spin-conserving normal elements $`s_{ee}^{\sigma \sigma }`$, $`s_{hh}^{\sigma \sigma }`$, and Andreev elements $`s_{eh}^{\sigma \sigma }`$, $`s_{he}^{\sigma \sigma }`$ coupling opposite spins. The calculation of the total zero- frequency noise $`S^{ch}=_{\sigma \sigma ^{}}S^{\sigma \sigma ^{}}`$, using the unitarity of the scattering matrix, reduces at zero temperature to the well-known result $`S^{ch}=\frac{4e^3V}{\pi \mathrm{}}Tr[s_{he}^{}s_{he}(1s_{he}^{}s_{he})]`$, where the trace is made on the channel indexes beenakkerNS ; martinNS . We remark here that it is easy to calculate the spin-resolved correlations $`S^{\sigma \sigma }`$ and $`S^{\sigma \sigma }`$, and to check that they are exactly equal. The result of this observation is that for a NS junction, at $`T=0`$, the spin current shot noise is strictly zero, $`S^{sp}=0`$. The current correlation between electrons with opposite spins is $`S^{}=S^{}`$, therefore positive. This ”bunching” of opposite spins carriers is an obvious consequence of the Andreev process, since each spin-up quasiparticle crossing the junction is accompanied by a spin-down quasiparticle. This nearly instantaneous correlation is due to the conversion of Cooper pair wavepackets in S, into pairs of normal wavepackets which carry no spin, therefore the spin current noise is zero. It has been discussed in a three-terminal geometry in Ref. FazioSFF, .

Let us now turn to a very different situation, that of a quantum dot in the Coulomb blockade regime. Here, instead of the attractive correlations manifested by the NS junction, repulsive correlations are expected. Let us consider a small island connected by tunnel barriers to normal leads $`L`$ and $`R`$ with electrochemical potentials $`\mu _{L,R}`$, such as $`eV=\mu _L\mu _R`$ (Fig. 1). The spectrum of this quantum dot is supposed to be discrete, e.g. the couplings $`\mathrm{\Gamma }_{L,R}2\pi |t_{L,R}|^2N_{L,R}(0)`$ to the leads verify $`\mathrm{\Gamma }_{L,R}<<\delta \epsilon `$, the level spacing. We also assume that max$`(eV,k_BT)>>\mathrm{}\mathrm{\Gamma }_{L,R}`$ and that only one level of energy $`E_0`$ sits between $`\mu _R`$ and $`\mu _L`$. The dot can be in three possible charge states, depending on whether the level is occupied by zero, one or two electrons (Fig. 1). These states will be indexed as $`N=0`$, $`N=1`$ (with spins $``$, $``$) and $`N=2`$. Let us denote as $`U(N)`$ the Coulomb energy for the state $`N`$, $`\mathrm{\Delta }E_{L,R}^+(N)=E_0\mu _{L,R}+U(N+1)U(N)`$ the energy to add an electron to state $`N`$ from leads $`L,R`$, and $`\mathrm{\Delta }E_{L,R}^{}(N)=E_0+\mu _{L,R}+U(N1)U(N)`$ the energy to remove an electron from state $`N`$ towards $`L,R`$. Let us further assume that $`\mathrm{\Delta }E_L^+(0)`$, $`\mathrm{\Delta }E_R^{}(1)<<k_BT`$. This implies that the transitions from $`N=0`$ to $`1`$ involve electrons coming only from $`L`$, and the transitions from $`N=1`$ to $`0`$ involve electrons going only into $`R`$. Let us allow the Coulomb energy to vary and consider the possibility of transitions from $`N=1`$ to $`2`$, only from $`L`$, e.g. $`\mathrm{\Delta }E_R^{}(2)<<k_BT`$. Yet, $`\mathrm{\Delta }E_L^+(1)`$ can take any value. This describes the following situation : if $`\mathrm{\Delta }E_L^+(1)>>k_BT`$, the transition to state $`N=2`$ is forbidden and one has the simple SET case with only two charge states $`N=0`$, $`1`$, in the resonant regime at low temperature (Fig. 1a). If on the contrary $`\mathrm{\Delta }E_L^+(1)<<k_BT`$, then the three charge states $`0`$, $`1`$, $`2`$ are involved in the charge transport (Fig. 1b). The physical situation under consideration corresponds for instance to fixing the gate voltage such as $`U(1)=U(0)`$, and varying the ratio between $`k_BT`$ and the Coulomb excess energy $`U(2)U(1)`$.

Let us write the master equation describing this system mastereq1 ; mastereq2 ; mastereq3 . Assuming a constant density of states in the reservoirs and defining $`x`$ as the Fermi function $`x=[1+exp(\beta \mathrm{\Delta }E_L^+(1))]^1`$, the non-zero transition rates are $`\mathrm{\Gamma }_L^+(0)=\mathrm{\Gamma }_L`$, $`\mathrm{\Gamma }_R^{}(1)=\mathrm{\Gamma }_R`$, $`\mathrm{\Gamma }_L^+(1)=x\mathrm{\Gamma }_L`$, $`\mathrm{\Gamma }_L^{}(2)=(1x)\mathrm{\Gamma }_L`$ and $`\mathrm{\Gamma }_R^{}(2)=\mathrm{\Gamma }_R`$ (Fig. 1b). Then the populations $`p_0`$, $`p_{}`$, $`p_{}`$ and $`p_2`$ verify

$`\dot{p}_0=2\mathrm{\Gamma }_Lp_0+\mathrm{\Gamma }_R(p_{}+p_{})`$
$`\dot{p}_{}=(\mathrm{\Gamma }_R+x\mathrm{\Gamma }_L)p_{}+\mathrm{\Gamma }_Lp_0+((1x)\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)p_2`$
$`\dot{p}_{}=(\mathrm{\Gamma }_R+x\mathrm{\Gamma }_L)p_{}+\mathrm{\Gamma }_Lp_0+((1x)\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)p_2`$ (1)
$`\dot{p}_2=2((1x)\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)p_2+x\mathrm{\Gamma }_L(p_{}+p_{})`$
Let us first consider the limit $`x=1`$. Then the transition rates from $`L`$ or into $`R`$ do not depend on the charge state, which means that this limit is equivalent to that of a resonant state without Coulomb charging energy. The solution of Eqs. (Spin current shot noise as a probe of interactions in mesoscopic systems) factorizes in this case, e.g. $`p(n_{},n_{})=p(n_{})p(n_{})`$, so that for each spin component, the probabilities $`p(0)`$ and $`p(1)`$ of empty and occupied states verify $`\dot{p}(0)=\dot{p}(1)=\mathrm{\Gamma }_Lp(0)+\mathrm{\Gamma }_Rp(1)`$. From this one derives the average current $`I=2e\frac{\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{\mathrm{\Gamma }_L+\mathrm{\Gamma }_R}`$ and the zero-frequency shot noise $`S_{ij}(\omega =0)=2eI(1\frac{2\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{(\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)^2})`$, independently of the couple of junctions $`i,j=L,R`$ ChengTing . Moreover, it is simple to check that spin $``$ and $``$ currents are uncorrelated, thus $`S_{ij}^{}=S_{ij}^{}=0`$, or equivalently $`S^{sp}=S^{ch}`$, as can also been derived by the scattering method in the quantum coherent regime. We thus have another example of the general behavior for uncorrelated transport.

Let us now consider the SET case $`x=0`$, where charge transport is maximally correlated. The charge noise is given by the expression $`S_{ij}(\omega =0)=2eI(1\frac{4\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{(2\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)^2})`$ nazarovstruben . Apart from an effective doubling of the rate $`\mathrm{\Gamma }_L`$, this result is qualitatively similar to that obtained without interactions. Therefore the charge noise is not the best possible probe of interactions. We now show that, on the contrary, the behaviour of the spin noise is completely different. Indeed, using the method by Korotkov korotkov , we find that

$`S_{ij}^{\sigma \sigma }=eI(1{\displaystyle \frac{2\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{(2\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)^2}})S_{ij}^{\sigma \sigma }=eI{\displaystyle \frac{2\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{(2\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)^2}}`$
$`S_{ij}^{sp}=2eI`$ (2)
The striking result (Eq. (2)) for $`S^{sp}`$ resembles a Poisson result and corresponds to maximal fluctuations. The correlations between currents of opposite spins are negative, like a partition noise. Yet spin-up and spin-down channels are separated energetically rather than spatially, and wavepackets with up or down spins exclude each other because of interactions rather than statistics.

The above result, obtained at zero temperature with perfect spin coherence inside the island, can be interpreted in the following way : electrons come from reservoir $`L`$ with random spins. Even though the average spin current is zero, each junction is sequentially crossed – due to Coulomb repulsion – by elementary wavepackets with well-defined but uncorrelated spins. This implies very short time correlations (on the scale of tunneling through one of the barriers) therefore the spin current exhibits Poisson statistics. On the contrary, charge current wavepackets are correlated on times $`\mathrm{}/\mathrm{\Gamma }_i`$, leading to the reduction as compared to the Poisson value. Notice that Eqs. (2) is a consequence of the restriction to two charge states : as can be easily checked, the analysis of the SET involving only $`N=1`$ and $`2`$ states (instead of $`0`$, $`1`$) yields exactly the same result.

The general solution of Eqs. (Spin current shot noise as a probe of interactions in mesoscopic systems) offers an interpolation between the uncorrelated and the maximally correlated regimes. We find that the average current is given by $`I=e\frac{2\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{\mathrm{\Gamma }_R+(2x)\mathrm{\Gamma }_L}`$. The spin current noise components $`S_{ij}^{\sigma \sigma ^{}}`$ (i,j=L,R) can also be calculated and do not depend on the couple $`(i,j)`$ of junctions chosen. The expression for the spin noise is

$$S_{ij}^{sp}=2eI(1\frac{2x\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{(\mathrm{\Gamma }_R+\mathrm{\Gamma }_L)(\mathrm{\Gamma }_R+x\mathrm{\Gamma }_L)})$$
(3)
The expression for the total (charge) noise $`S^{ch}`$ is too lengthy to be written here. Figs. 2, 3 shows the variation with $`x`$ of the charge and spin current noise. The spin noise is maximum for $`x=0`$, decreases monotonously and merges the charge noise at $`x=1`$. The role of the asymmetry of the junctions is very striking. First, if $`\mathrm{\Gamma }_R>\mathrm{\Gamma }_L`$, we find that $`S^{sp}`$ is always larger than $`S^{ch}`$ (Fig. 2), like in the ideal SET ($`x=0`$). On the other hand, if $`\mathrm{\Gamma }_R<\mathrm{\Gamma }_L`$, $`S^{sp}`$ happens to be smaller than $`S^{ch}`$ for $`x>x_c\mathrm{\Gamma }_R/\mathrm{\Gamma }_L`$ (Fig. 3). This implies that $`S^{}>0`$, contrarily to the naive expectation for repulsive interactions. This unexpected behavior can be explained as follows : if $`\mathrm{\Gamma }_R<\mathrm{\Gamma }_L`$, the low charge states are unfavored and the high ones favored, despite of Coulomb repulsion : for $`x>x_c`$ , $`p(2)`$ becomes larger than $`p(0)`$. When the SET occasionnally reaches the state $`N=0`$, a first transition leads to state $`1`$, but then the most probable transition is to state $`2`$ since $`\mathrm{\Gamma }_L^+(1)=x\mathrm{\Gamma }_L>\mathrm{\Gamma }_R^{}(1)=\mathrm{\Gamma }_R`$ : two electrons enter the dot successively, with opposite spins, leading to a certain degree of bunching. Here the anomaly is due to a kind of ”population inversion” (the most energetical state is favoured), manifesting a strong departure from equilibrium (Fig. 3). Yet, the effect is rather weak, less than $`10\%`$, contrarily to the NS junctions where attractive correlations are $`100\%`$.

Let us now consider how the above results are modified by spin relaxation, due for instance to spin-orbit scattering. Let us simply focus on the SET with two charge states $`0`$, $`1`$ and introduce a spin-flip rate $`T_1^1=\gamma _{sf}`$. The master equations will be written in this case

$`\dot{p}_0=2\mathrm{\Gamma }_Lp_0+\mathrm{\Gamma }_R(p_{}+p_{})`$
$`\dot{p}_{}=\mathrm{\Gamma }_Rp_{}+\mathrm{\Gamma }_Lp_0+{\displaystyle \frac{1}{2}}\gamma _{sf}(p_{}p_{})`$ (4)
$`\dot{p}_{}=\mathrm{\Gamma }_Rp_{}+\mathrm{\Gamma }_Lp_0+{\displaystyle \frac{1}{2}}\gamma _{sf}(p_{}p_{})`$
The introduction of spin relaxation obviously changes neither the average curent nor the charge shot noise. On the contrary, the spin shot noise is altered. For instance, opposite spin noise correlations between junctions $`L`$ and $`R`$ become $`S_{LR}^{\sigma \sigma }=eI(\frac{2\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{(2\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)^2}\frac{\gamma _{sf}}{2(\mathrm{\Gamma }_R+\gamma _{sf})})`$, and can even become positive. This results in a spin noise

$$S_{LL}^{sp}=2eI;S_{LR}^{sp}=2eI\frac{\mathrm{\Gamma }_R}{\mathrm{\Gamma }_R+\gamma _{sf}}$$
(5)
The reduction of the spin noise $`S_{LR}^{sp}`$ from the ”Poisson” value is a fingerprint of spin relaxation. Remarkably enough, the spin noise on junction $`L`$ is not affected, since the spins of entering wavepackets are uncorrelated whatever happens in the island. While transient current measurements have allowed to measure $`T_1`$ in presence of Zeeman splitting T1 , our result suggests an alternative method which does not require a magnetic field. Notice that noise was recently proposed to test spin flip in absence of Coulomb repulsion spinflip .

Let us now propose a set-up for the measurement of spin current correlations in a single-electron transistor. One may consider a four-terminal configuration spinflip , where the two left terminals $`L1`$, $`L2`$ are ferromagnetic metals with opposite spin polarizations, having the same chemical potential $`\mu _L`$ (Fig. 4). Similarly, terminals $`R1`$ and $`R2`$ have opposite polarizations, respectively parallel to those of $`L1`$, $`L2`$, and the same chemical potential $`\mu _R`$. If the junction parameters are the same for $`L1`$, $`L2`$ on one hand, and for $`R1`$ and $`R2`$ on the other hand, then the net current flowing through the SET is not spin polarized. Yet, it is possible to measure separately the spin current components in each of the four terminals, e. g. measure the noise correlations $`S_{L1L1}`$, $`S_{L1L2}`$, $`S_{L1R1}`$, $`S_{L1R2}`$, etc… If each terminal generates a fully spin-polarized current, the analysis of this set-up can be mapped onto the above model, and the previous results hold. In the more realistic case where polarization is not perfect, the above measurement would yield a mixing of the spin noise with the charge noise. If those are sufficiently different (strong repulsive correlations), they could still be distinguished, which allows to probe the Coulomb correlations by the method of spin current noise.

In summary, we have proposed to probe the attractive or repulsive correlations induced by interactions by measuring the noise correlations of the spin components of the current. This requires not to break the spin symmetry in the device, e.g. the total current is not spin-polarized. We have illustrated this trend on two simple and classical mesoscopic devices. First, a NS junction shows opposite spin bunching due to attractive correlations. Second, a SET in the sequential regime shows in general repulsive correlations (antibunching), but those can be weakly attractive far from equilibrium. Extensions to other regimes or multiple dot systems is quite promising.

###### Acknowledgements.

The authors are grateful to Th. Martin for fruitful discussions concerning the ”partition noise” analogy. LEPES is under convention with Université Joseph Fourier.


text:
# Magnetohydrodynamic Models for the Structure of Pulsar-Wind Nebulae

## 1 Abstract

The Crab Nebula is well-described at optical wavelengths and above by a steady-state magnetohydrodynamic model due to Kennel and Coroniti (1984). Can this class of model describe other pulsar-wind nebulae? I exhibit simple generalizations of KC models, for various values of $`\sigma ,`$ the ratio of magnetic to particle flux input at the wind shock. I calculate the evolution of the electron spectrum and synchrotron emissivity in the nebula assuming spherical symmetry, a steady state, and a purely toroidal magnetic field. Emission profiles depend on the initial magnetic field $`B_0`$, the electron spectral index $`s`$, and the angle of the toroidal axis with the line of sight, $`\varphi `$. I show integrated spectra and radial profiles for various cases, along with predicted variations of photon index $`\mathrm{\Gamma }`$ with radius in X-rays. Most models predict much smaller sizes in X-rays, and curves of $`\mathrm{\Gamma }(r)`$ which do not resemble observations. Further elaborations of the dynamics of PWNe seem necessary.

## 2 Introduction

Pulsar-wind nebulae (PWNe) are bubbles of relativistic particles and magnetic field inflated by a pulsar, forming center-brightened radio and X-ray sources of synchrotron radiation Weiler80 , RC84 , KC84a , KC84b . Most early work on PWNe assumed they were surrounded by shell supernova remnants (SNRs), seen or unseen, and interacted with the interior of those remnants. Before the passage of a reverse shock, the PWN would expand into unshocked ejecta in the SNR interior; afterwards, it would be compressed by the reverse-shock passage and expand much more slowly into the thermalized, shocked ejecta RC84 .

The most elaborate early model for the high-energy (optical and above) emission from PWNe was that of Kennel & Coroniti KC84a , KC84b  who solved the steady-state relativistic MHD equations in spherical symmetry with an azimuthal magnetic field to obtain velocity and density profiles of the outflowing relativistic material. They showed that a strong (high compression) MHD shock required a fluid dominated by particles, not magnetic field, as parameterized by $`\sigma (B/4\pi )/n\gamma ^2\beta mc^2`$, the ratio of magnetic (plus electric) energy flux to particle flux just upstream of the wind shock. Here $`n`$ is the comoving density and $`\gamma `$ the flow Lorentz factor. Applying the Rankine-Hugoniot jump conditions at the wind shock (radius $`r_s`$) showed that unless $`\sigma 1`$, the postshock fluid would remain at high speed, within a factor of 3 of the upstream (relativistic) speed – clearly inconsistent with the observed expansion of the Crab Nebula at $`v2000`$ km s<sup>-1</sup>. Behind the shock, pure flux-freezing meant that the toroidal magnetic field strength evolved as $`B\rho r`$, with $`\rho `$ the matter density. Kennel & Coroniti KC84b  were able then to follow the evolution of electron energies and predict both the integrated spectrum and the brightness profile of remnants. They specialized to the Crab Nebula and showed calculations primarily for that case. They were able to describe the optical through hard X-ray spectrum well, as well as the gross variation of nebular size with X-ray energy (“size” being typically the FWHM of emission measured with relatively crude imaging instruments). However, the high value of preshock Lorentz factor of the wind they found ($`\gamma 10^6`$) made it difficult to explain the radio emission which requires much lower energies, and more electrons, than would be produced by thermalizing this extremely relativistic wind.

Since the advent of the newer generation of X-ray telescopes, beginning with ASCA and continuing with the Chandra X-ray Observatory and XMM-Newton, we have seen an explosion in the number and variety of PWNe in X-rays, and have observed a much wider range of morphological type. It is of interest to ask if the class of steady-state MHD models proposed by KC can be generalized appropriately to describe other PWNe, and if not, to get an idea of what kind of extensions might be necessary. This is the project described in this paper.

One major problem in modeling several larger, fainter PWNe such as 3C 58 Torii00  appears in comparing radio and X-ray images. In Kennel & Coroniti’s model for the Crab, synchrotron losses on outflowing particles cause the nebular size to decrease with frequency at optical and higher frequencies. In fact the Crab is considerably smaller at optical than radio frequencies, and continues to diminish in size roughly as $`R\nu ^{0.15}`$ up to 40 keV Ku76 . However, 3C 58 shows X-ray emission extending to near the edges of the radio emission, though it is faint there (Slane, private communication). These results call into question the MHD picture of transport of particle energy and magnetic flux throughout the volume of a PWN.

It is straightforward to show from the MHD equations KC84a  that postshock flow solutions have accelerating and decelerating branches, with the latter beginning at $`v_0=c/3`$ and initially decelerating as $`r^2`$ (pure hydrodynamic, isobaric flow; magnetic energy is unimportant since $`\sigma `$ must be small to fit the outer nebular boundary conditions). However, in this case $`B\rho r1/rvr`$, so the magnetic energy rises until it becomes dynamically important, slowing the deceleration and causing the velocity to approach a constant value, with $`B1/r`$. The full hydrodynamic solution can be well approximated by a simple analytic form: defining $`uv/v_0,`$

$$u(r)=\frac{1}{r^2}(1u_{\mathrm{}})+u_{\mathrm{}},$$
(1)
where the asymptotic velocity is $`v_{\mathrm{}}/c=\sigma /(1+\sigma )`$ (Fig. 1).

## 3 Particle evolution and luminosity calculation

Electron energies $`E`$ evolve due both to radiative and adiabatic losses: $`\dot{E}=(dE/dV)(dV/dt)aB^2E^2`$, where for synchrotron losses $`a=1.57\times 10^3`$ cgs. This equation can be integrated to give KC84b , R98

$$E(t)=\frac{E_0\alpha ^{1/3}}{1+aE_0B^2\alpha ^{1/3}𝑑t}=\frac{E_0\alpha ^{1/3}}{1+E_0\left(\frac{aB_0^2r_0}{v_0}\right)z^{8/3}u^{10/3}𝑑z}$$
(2)
where $`\alpha \rho /\rho _0,`$ $`zr/r_0,`$ and $`uv/v_0`$. The quantity $`E_f(aB_0^2r_0/v_0)^1`$ is a “fiducial energy”, that an initially infinitely energetic electron would have after radiating in a field $`B_0`$ for a time $`r_0/v_0`$ (basically a parameterization of $`B_0`$). From this result, the evolution of an arbitrary distribution can be calculated:

$$N(E)=N(E_0)\frac{dE_0}{dE}\frac{dV_0}{dV}=N(E_0(E))\frac{dE_0}{dE}\alpha =N(E_0(E))\left(\frac{E_0^2}{E^2}\right)\left(\frac{\rho }{\rho _0}\right)^{4/3}.$$
(3)
Given the approximate dynamics above, the electron distribution and synchrotron emissivity can be calculated numerically at each radius. If magnetic field lines are circles in an equatorial plane, an elementary rotation can give brightness profiles as a function of aspect angle $`\varphi `$ between the polar axis and the line of sight. The synchrotron emissivity was integrated numerically to generate the brightness profiles shown above, for an injected flat power-law $`N(E)=KE^s`$ with $`s<2`$ as deduced from observed PWN radio spectra. The flat spectra cause a “bump” to appear just below the high-energy cutoff that appears as a result of radiative losses.

## 4 Conclusions

Steady-state MHD models of pulsar-wind nebulae can be well approximated by a simple velocity law, which then dictates the evolution of the injected electron distribution and the integrated spectrum and profile of emission, depending on $`\sigma `$. Spectral properties also depend on B and spectral index $`s`$. Varying the aspect angle changes the flux normalization but has no effect on the integrated spectrum and little on the profile shape at photon energies for which losses are important. For very flat injected spectra ($`s1`$), the effects of the energy-loss “bump” are significant, causing concave-up curvature over a broad frequency range.

The models are in significant conflict with the observed radial dependence of photon index $`\mathrm{\Gamma }`$ in X-rays (e.g., G21.5-0.9 Slane00 ), which appears roughly linear with $`r`$. All models drop in size by at least a factor 2 between radio and X-rays, but some observed PWNe do not show this. Unless magnetic-field strengths are unrealistically low, most X-ray profiles have HWHM radii only 3–6 times the injection radius – also perhaps at odds with observations.

We conclude that there are significant discrepancies between the predictions of the simple MHD models and X-ray observations of PWNe. Nonsteady or nonspherical flows, and/or electron transport by diffusion as well as convection (e.g., RJ91 ), may be necessary to explain the observations.

## Index